Changes in ASTTI in two switchgrass genotypes
The alkaline salt tolerance trait index (ASTTI) was used to evaluate the physiological responses of Alamo (alkali-tolerant genotype) and AM-314/MS-155 (alkali-sensitive genotype) under alkaline salt stress for 0, 3, 6, 12 and 24 h. Ten physiological traits, including the relative water content (RWC), relative electrical conductivity (REC), the contents of malondialdehyde (MDA), free proline, soluble protein, soluble sugar, and reduced glutathione (GSH), and the activities of superoxide dismutase (SOD), peroxidase (POD) and catalase (CAT), were divided into three groups for subsequent analysis. The RWC, REC and MDA contents were used as indicators of the membrane system, and the ASTTI value of the root RWC (Fig. 1A) showed an upward trend, while the ASTTI values of the REC and MDA content showed a downward trend (Fig. 1B and C). The root REC and MDA contents in AM-314/MS-155 were significantly higher than those in Alamo at 4 h post alkaline salt stress. The contents of free proline, soluble protein and soluble sugar were used as indicators of the osmotic adjustment system, and the osmotic regulation level of the two cultivars first increased with increasing stress duration and then decreased (Fig. 1 D-F). The osmotic regulation ability of Alamo was significantly higher than that of AM-314/MS-155 and was highest at 0, 6 and 24 h. For the indicators of the antioxidant system (GSH, SOD, POD, CAT), the activities of SOD, POD and CAT of Alamo were significantly higher than those of AM-314/MS-155 (Fig. 1H-J), while the content of GSH, a nonenzymatic substance, was significantly lower than that of AM-314/MS-155 at the early stage (Fig. 1G). In addition, we observed that there were significant physiological differences between the two varieties at different time points under alkali stress, especially at 6 and 24 h. Therefore, we chose these two time points for further transcriptome analysis.
De novo assembly and annotation of unigenes
Six cDNA libraries were constructed from total RNA extracted from E5 stage roots of Alamo and AM-314/MS-155 treated with alkaline salt stress for 0, 6 and 24 h. A total of 114.03 Gb of clean data were obtained, with 6.05 Gb of clean data for each sample and a Q30 base percentage of 93.21% or more (Additional file 1: Table S1). Since the whole-genome sequence of switchgrass is not currently publicly available, valid readings from the six libraries were combined for reassembly (Additional file 2: Table S2). A total of 108,319 unigenes were obtained after reassembly, including 73,636 unigenes in AM-314 and 65,492 unigenes in Alamo. The total length, N50 length and mean length were 106,429,710 nt, 1,751 nt, and 982.56 nt, respectively. Functional annotation of assembled sequences was primarily based on BLAST homology searches with the NCBI nonredundant protein (Nr), Evolutionary Genealogy of Genes: Nonsupervised Orthologous Groups (eggNOG), Gene Ontology (GO), Pfam, Swiss-Prot, Eukaryotic Orthologue Groups (KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Clusters of Orthologous Groups (COG) databases, and a total of 66,253 genes were annotated (Additional file 3: Table S3). Among these annotated sequences, 63,663 (96.09%) sequences had significant matches in the Nr database. On the basis of the homology among sequences of different species, 25,805 (40.56%) sequences were found against Setaria italica, and 5,124 (8.05%) sequences had clandestine hits with Sorghum bicolor, followed by Zea mays (4,732, 7.44%), Oryza sativa (1,590, 2.50%), Phaeosphaeria nodorum (1,564, 2.46%), Verticillium dahliae (950, 1.49%), V. alfalfae (839, 1.32%), Phytophthora sojae (680, 1.07%), Togninia minima (627, 0.99%) and Hordeum vulgare (584, 0.92%). Only 21,129 (33.21%) annotated sequences were similar to those of other plant species (Additional file 4: Figure S1A). The gene expression levels in response to alkaline salt stress were evaluated by converting the mapped read count for each gene into the expected number of fragments per kilobase of reproduction per million mapped reads (FPKM) to eliminate the effects of transcript length and sequencing differences on computational expression (Additional file 5: Table S4). The gene expression levels were not evenly distributed in the different stress environments in the boxplot diagram of the FPKM values (Additional file 4: Figure S1B). Then, the correlations of each biological sample were evaluated by Pearson correlation coefficients, and a value of r2 close to 1 indicated a strong correlation between two replicate samples (Additional file 4: Figure S1C). Finally, all the unigenes were used for further identification of DEGs after the exclusion of abnormal samples.
Identification of DEGs associated with alkaline salt stress
To identify the DEGs of the two different genotypes under alkali stress, the expression patterns of DEGs were analysed by comparing 6-h and 24-h libraries with the control library for Alamo and AM-314/MS-155. A total of 10,219 DEGs whose expression was up- or downregulated between samples (fold change ≥ 2 and false discovery rate (FDR) < 0.01) at any pair of alkaline salt-treated points were identified (Additional file 6: Table S5). In the early stage of alkaline salt stress (0-6 h), 4483 genes (58.02% of the DEGs) were upregulated in AM-314/MS-155, and 867 genes (61.36% of the DEGs) were upregulated in Alamo (Table 1). In the late stage of alkaline salt stress (24 h), 2942 genes (56.17% of the DEGs) were upregulated in AM-314/MS-155, while Alamo had 2732 upregulated genes (61.61% of the DEGs). The proportion of upregulated genes in the tolerant genotype (Alamo) was higher than that in the sensitive genotype (AM-314/MS-155) in both the early and late stages of stress. In addition, the DEGs of AM-314/MS-155 were mainly concentrated in the early stage, while Alamo showed greater advantages in the late stage.
To further analyse the effects of alkaline salt treatment in the two genotypes over time, the expression trends of DEGs were clustered into 16 modules (Fig. 2A). Then, the mainstream gene expression trends (6 groups) were screened over the duration of the stress (Fig. 2B). The first group contained 1,831 enriched genes, which mainly reflected the functional classification of genes expressed abundantly in AM-314/MS-155 but with little or no expression in Alamo (Fig. 2C). A large number of ribosomal proteins and a small number of cytochrome and energy-related genes were enriched in the ribosome and oxidative phosphorylation pathways and gradually downregulated with increasing stress time (Additional file 7: Table S6). The expression level of genes in group 2 (1,492 genes) decreased sharply in the early stage of alkaline salt stress in both switchgrass genotypes. KEGG analysis of the genes in group 2 revealed that most were involved in the ribosome, phenylalanine metabolism, plant hormone signal transduction, and ribosome biogenesis in eukaryotes pathways (Additional file 8: Table S7). The expression of genes in group 3 (1,261 genes) was completely opposite to that in group 1, mainly reflecting the functional classification of genes expressed abundantly in Alamo but with little or no expression in AM-314/MS-155. Genes in this group functioned mostly in the glycolysis/gluconeogenesis, plant-pathogen interaction, amino sugar and nucleotide sugar metabolism, and ubiquitin-mediated proteolysis categories (Additional file 9: Table S8). The expression of genes in group 4 (1,072 genes) increased strongly in the early stage of alkaline salt stress, and most of these genes functioned in starch and sucrose metabolism, valine, leucine and isoleucine degradation, plant hormone signal transduction, galactose metabolism, and fatty acid degradation (Additional file 10: Table S9). Although there was a fast decrease in gene expression levels in both group 5 (718 genes) and group 6 (736 genes), the decrease in the expression of genes in group 5 was mainly in AM-314/MS-155, while the decrease in the expression of genes in group 5 was mainly in Alamo. Genes involved in oxidative phosphorylation, protein processing in the endoplasmic reticulum, glycolysis/gluconeogenesis, and glycerophospholipid metabolism were enriched in group 5 (Additional file 11: Table S10), and genes in group 6 were enriched in ribosome, oxidative phosphorylation, glycolysis/gluconeogenesis, and DNA replication (Additional file 12: Table S11).
Identification of TFs and PKs
In this study, 1,480 TF genes were predicted using ITAK software (http://itak.feilab.net/cgi-bin/itak/index.cgi) and were classified into 64 TF families (Additional file 13: Table S12). The most abundant TF families were the C2H2 (153 genes), bZIP (108 genes), bHLH (93 genes), C3H (91 genes), MYB-related (82 genes), NAC (81 genes), WRKY (74 genes), GRAS (57 genes) and AP2/ERF-ERF (52 genes) families. The genes in the bHLH, C2H2, NAC and WRKY families were mainly upregulated under alkali stress, while the genes in the AP2/ERF-ERF, bZIP and MYB-related families showed a relatively balanced number of upregulated and downregulated members (Additional file 14: Figure S2). There were 32 TF genes that showed differential expression in both Alamo and AM-314/MS-155 under alkaline salt stress, including 7 TF genes that showed differential expression only at 6 h, 10 TF genes that showed differential expression only at 24 h, and 15 TF genes that showed differential expression at both 6 and 24 h (Additional file 15: Table S13). Interestingly, the TF genes only differentially expressed at 6 h were all downregulated, and most of them belonged to the MYB-related family. Nine of ten TF genes differentially expressed only at 24 h were unregulated, and most of them belonged to the NAC family. There were 27 upregulated TF genes and 16 downregulated TF genes specifically expressed in Alamo, and many of them were differentially expressed in the late stage of alkali stress (24 h) and belonged to the bZIP, MYB-related and bHLH families (Additional file 16: Table S14).
In addition, a total of 1,718 PKs, such as calcium/calmodulin-dependent protein kinase (CaMK), CDPK, MAPK and receptor-like kinase (RLK), were predicted in AM-314/MS-155 and Alamo under alkaline salt stress for 6 and 24 h (Additional file 17: Table S15). The CBL-interacting protein kinases (CIPKs), which play a role as Ca2+ sensors, were upregulated and highly expressed in Alamo under alkaline salt stress for 24 h. RLKs accounted for the largest proportion among many protein kinase families, and the genes belonging to the RLK family maintained a high expression level at 6 h and were then downregulated at 24 h.
Weighted gene coexpression network analysis (WGCNA) of DEGs in response to alkaline salt stress
To further study the patterns of association between the DEGs responding to alkaline salt stress and physiological indicators across the two genotypes, WGCNA was performed to explore the gene modules for synergistic expression. After the filtering of low-expression genes (FPKM <5), the remaining 7,485 genes were classified into 20 different modules (combined modules) (Fig. 3A).
The blue module showed the highest positive correlation with the REC (R=0.82, P=3e-5) (Fig. 3B). This module was characterized by a significant upregulation of DEGs in the alkali-sensitive genotype AM-314/MS-155 at 6 and 24 h (Fig. 3C). Genes in this module were significantly enriched in the peroxisome, amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism and plant-pathogen interaction pathways by KEGG (Additional file 18: Figure S3A). The brown module showed a positive correlation with soluble protein (R=0.64) and the RWC (R=0.49). This module represented the downregulated genes in Alamo under alkaline salt stress. The functions of these genes were mainly concentrated in the plant hormone signal transduction, phenylpropanoid biosynthesis and biosynthesis of amino acid pathways (Additional file 18: Figure S3B). The dark magenta module showed a high positive correlation with MDA, GSH and SOD, with R values of 0.70, 0.68 and 0.63, respectively. This module showed the features of significantly upregulated genes in both Alamo and AM-314/MS-155 under alkaline salt stress for 6 h. Genes in this module were mainly enriched in the amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, and other functional pathways (Additional file 18: Figure S3C). The light steel blue 1 module was positively correlated with proline, soluble sugar and CAT, with R values of 0.80, 0.71 and 0.64, respectively (Fig. 4-B). The module mainly showed genes in the alkali-tolerant Alamo that were significantly upregulated under alkaline salt stress for 6 h. These genes were mainly enriched in the functions of phenylalanine metabolism and phenylpropanoid biosynthesis (Additional file 18: Figure S3D).
Verification of RNA-Seq sequencing data by qRT-PCR analysis
The DEGs associated with alkaline salt stress were selected for qRT-PCR assays to verify the accuracy of the RNA-Seq sequencing data. Ten genes were selected randomly from the DEGs that were coexpressed in all sequencing samples. The expression patterns of these ten genes were the same as those of the RNA-Seq assay (Fig. 4A). The expression levels indicated by the qRT-PCR results were strongly correlated with the differential gene expression levels identified by mRNA-seq according to the Pearson correlation coefficients, which were 0.8934 and 0.9144 under alkaline salt stress for 6 and 24 h, respectively (Fig. 4B), and demonstrated the reliability of the RNA-seq data.