Identification and Analysis of NaHCO3 Stress Responsive Genes in Wild Soybean (Glycine soja) Roots by RNA-seq

Abstract

Soil alkalinity is a major abiotic constraint to crop productivity and quality. Wild soybean (Glycine soja) is considered to be more stress-tolerant than cultivated soybean (G. max), and has considerable genetic variation for increasing alkalinity tolerance of soybean. In this study, we analyzed the transcriptome profile in the roots of an alkalinity tolerant wild soybean variety N24852 at 12 and 24 h after 90 mM NaHCO₃ stress by RNA-sequencing. Compared with the controls, a total of 449 differentially expressed genes (DEGs) were identified, including 95 and 140 up-regulated genes, and 108 and 135 down-regulated genes at 12 and 24 h after NaHCO₃ treatment, respectively. Quantitative RT-PCR analysis of 14 DEGs showed a high consistency with their expression profiles by RNA-sequencing. Gene Ontology (GO) terms related to transcription factors and transporters were significantly enriched in the up-regulated genes at 12 and 24 h after NaHCO₃ stress, respectively. Nuclear factor Y subunit A transcription factors were enriched at 12 h after NaHCO₃ stress, and high percentages of basic helix-loop-helix, ethylene-responsive factor, Trihelix, and zinc finger (C2H2, C3H) transcription factors were found at both 12 and 24 h after NaHCO₃ stress. Genes related to ion transporters such as ABC transporter, aluminum activated malate transporter, glutamate receptor, nitrate transporter/proton dependent oligopeptide family, and S-type anion channel were enriched in up-regulated DEGs at 24 h after NaHCO₃ treatment, implying their roles in maintaining ion homeostasis in soybean roots under alkalinity. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis showed "phenylpropanoid biosynthesis" and "phenylalanine metabolism" pathways might participate in soybean response to alkalinity. This study provides a foundation to further investigate the functions of NaHCO₃ stress-responsive genes and the molecular basis of soybean tolerance to alkalinity.

Keywords: Gene Ontology (GO) enrichment analysis; RNA-seq; alkalinity; differentially expressed gene (DEG); ion transporter; wild soybean.

FIGURE 1

Ion concentrations and ratios in N24852 and Lee 68 under control (CK) and alkaline stress for 12 and 24 h. (A) Na⁺ concentration in roots, (B) Na⁺/K⁺ ratio in roots, (C) Na⁺/Ca²⁺ ratio in roots. (D) Na⁺ concentration in leaves, (E) Na⁺/K⁺ ratio in leaves, (F) Na⁺/Ca²⁺ ratio in leaves. Alkaline stress was 90 mM NaHCO₃ (pH = 8.5) treatment, and control was 0 mM NaHCO₃ (pH = 6.5). ^∗ and ^∗∗ indicate significant difference at 0.05 and 0.01 level by student’s t-tests between two soybean genotypes, respectively. DW, dry weight.

FIGURE 2

Venn diagram of the differentially expressed genes (DEGs) between alkaline stress (90 mM NaHCO₃, pH = 8.5) and control. A12h, 12 h of alkaline stress; A24h, 24 h of alkaline stress.

FIGURE 3

Heatmap of the overlap DEGs between two time points (12 and 24 h). Heatmap was plotted using heatmap.2 function of the R/Bioconductor package gplots. Hierarchical clustering of the DEGs was done by complete method with Euclidean distance. The gene expression levels were transformed by log₂ (FPKM+1) and the values were centered and scaled in row direction. X-axis, samples; Y-axis, differentially expressed gene names.

FIGURE 4

Enriched Gene ontology (GO) terms for up-regulated genes at 12 h (A) and 24 h (B) after 90 mM NaHCO₃ (pH = 8.5) stress. The Singular Enrichment Analysis (SEA) was carried out by agriGO. The GO terms with their IDs are written in boxes. The significant (P < 0.01, FDR < 0.05) GO terms are in colored boxes (the degree of color saturation is positively correlated to the enrichment level of the GO term), and non-significant terms are in white boxes.

FIGURE 5

Comparison of RNA-seq results and qRT-PCR analysis of gene expression levels. (A) Log₂ fold change of 14 genes at 12 h after 90 mM NaHCO₃ (pH = 8.5) stress. (B) Log₂ fold change of 14 genes at 24 h after 90 mM NaHCO₃ (pH = 8.5) stress. (C) The comparison of log₂ fold change obtained by RNA-seq (x-axis) and qRT-PCR (y-axis). ^∗∗Differentially expressed genes at FDR ≤ 0.01 and | log₂FoldChange|≥ 1 from RNA-seq.

FIGURE 6

Transcription factor families of the DEGs at 12 h (A) and 24 h (B) of NaHCO₃ treatment. A total of 130 and 173 transcription factors (TFs) were differentially expressed at 12 and 24 h after 90 mM NaHCO₃ (pH 8.5) stress, respectively, and the percentages of TF families in DEGs and all genes in soybean genome were shown as blue bars and red bars, respectively. ^∗Significantly enriched TF family identified by hypergeometric Fisher exact test (P < 0.01) and Benjamini and Hochberg method (FDR < 0.05).