Transcriptome analysis reveals the impact of arbuscular mycorrhizal symbiosis on Sesbania cannabina expose to high salinity

Arbuscular mycorrhiza can improve the salt-tolerance of host plant. A systematic study of mycorrhizal plant responses to salt stress may provide insights into the acquired salt tolerance. Here, the transcriptional profiles of mycorrhizal Sesbania cannabina shoot and root under saline stress were obtained by RNA-Seq. Using weighted gene coexpression network analysis and pairwise comparisons, we identified coexpressed modules, networks and hub genes in mycorrhizal S. cannabina in response to salt stress. In total, 10,371 DEGs were parsed into five coexpression gene modules. One module was positively correlated with both salt treatment and arbuscular mycorrhizal (AM) inoculation, and associated with photosynthesis and ROS scavenging in both enzymatic and nonenzymatic pathways. The hub genes in the module were mostly transcription factors including WRKY, MYB, ETHYLENE RESPONSE FACTOR, and TCP members involved in the circadian clock and might represent central regulatory components of acquired salinity tolerance in AM S. cannabina. The expression patterns of 12 genes involved in photosynthesis, oxidation-reduction processes, and several transcription factors revealed by qRT-PCR confirmed the RNA-Seq data. This large-scale assessment of Sesbania genomic resources will help in exploring the molecular mechanisms underlying plant–AM fungi interaction in salt stress responses.

Sesbania cannabina, a soil-improving legume is used as green manure to increase the many crops yield in saline soils. Therefore, exploiting S. cannabina -AM fungi symbiosis represents a good strategy for reclamation of saline soils. We recently observed that AM fungi could improve salt-tolerance in S. cannabina and increase biomass 14,15 . In order to comprehend the underlying mechanisms of mycorrhizal plant response to saline soils, in this study, we investigated the effects of Funneliformis mosseae inoculation on the growth of S. cannabina exposed to short and long-term salt stress and identified plant genes with key roles in the response to AM fungi by transcriptome profiling in S. cannabina shoot and root using Illumina sequencing. We focused on characterization the transcriptome of mycorrhizal plant under salt stress to provide insights into the molecular mechanism that allows AM fungi to induce salt tolerance in plants.

Results
transcriptome sequencing and assembly. In order to confirm the alleviation effect of AM fungi on S. cannabina under salinity, plant biomass was measured in the presence or absence of 200 mM NaCl with F. mosseae inoculation. In 10 days, the growth of S. cannabina inoculated with F. mosseae was faster than un-inoculated plants under saline treatment (Fig. S1). Compared with non-inoculated (NM) plants, AM-inoculated (AM) plants exhibited 61.5% and 64.7% increases in shoot and root dry weights in saline conditions, respectively, but the difference was not significant in water-treated groups (Fig. S1b,d). Similar patterns were observed in fresh biomass accumulation (Fig. S1c,e). These results suggested that F. mosseae promoted plant growth and biomass under saline conditions.
To understand the mechanisms underlying AM fungi-induced plant salt tolerance, 12 S. cannabina RNA-Seq libraries containing two biology replicates for each treatment were constructed from tissues (shoot and root) of AM and NM plants exposed to 200 mM NaCl for 0, 3 and 27 h. The sequencing data were deposited in the NCBI GEO database under accession number GSE99532. Altogether, 1,254,209,094 Illumina paired-end raw reads were generated (Table 1). After adaptor sequences, ambiguous nucleotides and low-quality sequences were removed, 1,199,463,482 clean reads were retained. As shown in Table S1, the clean reads assembly resulted in 202,621 unigenes in the range from 222 to 15,720 bp, with a N50 (2204 bp). The fragments per kilobase of exon per million fragments mapped (FPKM) data were tested by correlation analysis to evaluate sampling between each pair of bioreplicates. The minimum value of all correlation coefficient between each biological pairs was 0.74. (Fig. S2). sequence annotation. As shown in Table 2 www.nature.com/scientificreports www.nature.com/scientificreports/ (44.01%) in Pfam. Altogether, 215,546 unigenes (74.07%) were successfully annotated in at least one of the seven databases and 1505 unigenes (0.51%) annotated in all seven databases.
Comparative transcript expression analysis of AM and NM plants in response to salt stress. Differentially expressed genes (DEGs) were defined as unigenes that had significantly enriched or depleted transcripts after the treatments. Based on pairwise comparison, differential expression analysis between AM and NM plants exposed to saline treatment for different times (0, 3 and 27 h) was performed (  Table 2. BLAST analysis of non-redundant unigenes against public databases. www.nature.com/scientificreports www.nature.com/scientificreports/ differentially expressed in the shoots and roots of AM plants, respectively, compared with NM plants. Moreover, 1188 (0.66%) and 684 (0.38%) transcripts were differentially expressed in the shoots and roots of AM plants, respectively, after 3 h of salt treatment, whilst 1049 (0.58%) and 1200 (0.67%) transcripts were differentially expressed in the shoots and roots of AM plants, respectively, after 27 h of salt treatment. Commonly up-and downregulated unigenes were identified between AM and NM plants in response to short-and longterm salt stress to evaluate the overlap degrees. There were 238 upregulated and 280 downregulated transcripts in AM plant shoots comparing 3 h with CK ( Fig. 2a,b), while 459 transcripts were upregulated and 438 transcripts were downregulated in AM plants shoots comparing 27 h with CK ( Fig. 2c,d).
In order to reveal co-regulation patterns among the DEGs in mycorrhizal S. cannabina in response to short and long periods of saline treatment, the expression profiles of the unigenes were analysis by a hierarchical clustering algorithm. It showed distinct expression patterns for DEGs in AM and NM plants in response to saline treatments (Fig. S3).

Weighted Gene Coexpression Network Analysis (WGCNA).
To explore the acquired salinity tolerance in AM S. cannabina, 10,371 genes were selected to construct a scale-free coexpression network using WGCNA. WGCNA is a systemic approach especially for understanding biology networks instead of individual genes. In order to confirm that the network was biologically relevant, the scale-free topology model fit and the mean connectivity of the network was evaluated over a range of soft threshold powers (β), before selecting β = 9 (Fig. S4). A dynamic hierarchical tree algorithm was used to divide the clustering tree constructed from the DEGs, resulting in five coexpression modules, which were labeled turquoise (3,999 genes), blue (2,832 genes), brown (2,328 genes), yellow (676 genes), and green (263 genes) (Figs 3 and S5). The gene cluster modules were enriched in specific GO functional terms, and statistically significantly enriched terms (p < 0.05) were selected for further analysis. The genes in the turquoise module were mainly associated with organic substance biosynthetic processes (p = 5.1087E-16), organonitrogen compound metabolic processes (p = 3.2913E-20), cell parts (p = 5.0073E-18), intracellular regions (p = 9.10E-13), and structural molecule activity (p = 2.7081E-60); genes www.nature.com/scientificreports www.nature.com/scientificreports/ in the blue module were mainly involved in oxidation-reduction processes (p = 8.1242E-28), single-organism metabolic processes (p = 8.40E-19), membrane protein complexes (p = 0.028578), oxidoreductase activity, acting on other nitrogenous compounds as donors, with NAD or NADP as acceptor (p = 0.032091), metal ion binding (p = 0.000449) and catalytic activity (p = 4.03E-11); genes in the brown module were involved in single-organism metabolic processes (p = 0.000863) and oxidoreductase activity (p = 0.012702); the genes of the yellow module were involved in organic substance catabolic processes (p = 9.97E-05), lipid metabolic processes (p = 0.001364) and hydrolase activity, acting on glycosyl bonds (p = 0.000244); and genes of the green module were mainly involved in trehalose biosynthetic processes (p = 0.036706), disaccharide biosynthetic processes (p = 0.036706) and oligosaccharide biosynthetic processes (p = 0.64995).
At last, the results of eigengene-trait correlation analysis indicated that the blue module was positively correlated with both salt stress and AM inoculation after multiple testing correction (p < 0.05) (Fig. S6). Identifying the genes in this module and their biological roles in response to salinity and AM inoculation was of particular interest, so functional annotation was performed. Detailed information on the GO enrichment in the modules is provided in Table S3 with P-value < 0.05. According to it, the blue module was significantly enriched with unigenes functioning in oxidation-reduction processes, oxidoreductase activity, photosynthetic membranes, and iron ion binding. These biological activities may play central roles in mycorrhizal plants under saline stress. KEGG pathways (http://www.genome.jp/kegg/) enrichment analysis (Table S4) showed that photosynthesis, carotenoid biosynthesis and glycine, serine, alanine, aspartate and glutamate metabolism were the most enriched physiological processes.
Hub genes are genes which have the most connections in a network. As shown in the blue module, 34 of the 85 hub genes (the top 3%) encoded transcription factors (Tables S5 and S6). The majority of transcription factors  Table S5).

Validation of RNA-seq data by qRt-pCR.
To verify our bioinformatics analyses, the expression levels of 12 selected genes which encode WRKY transcription factor, MYB-related transcription factor LHY, ethylene-responsive transcription factor, transcription factor TCP, photosystem I subunit VI, ribose 5-phosphate isomerase A, light-harvesting complex II chlorophyll a/b binding protein, copper/zinc superoxide dismutase (SODC), nickel-containing superoxide dismutase, catalase, glutathione peroxidase, peroxidase were quantified by qRT-PCR. The results showed comparable expression patterns of all the selected genes in real-time PCR analysis and in the RNA-Seq data (Fig. S7a,b). The correlation between RNA-Seq and qRT-PCR platforms was also calculated using SPSS V17.0 (SPSS, Inc., Chicago, IL, USA), the pearson correlation coefficients were 0.903 and 0.939 (Fig. S7c), confirming that our experimental results were valid.

Discussion
It has been proposed that mycorrhizal colonization enhances plant salt tolerance by improving photosynthetic ability, water and nutrient uptake, ion balance and osmolyte concentrations, among other characteristics [16][17][18] . However, information about the molecular bases of such effects is limited. The present study demonstrated that F. mosseae alleviated salinity stress in S. cannabina as shown in previous studies. Illumina sequencing and WGCNA analysis were used to reveal the gene network stimulated by AM fungi in S. cannabina plantlets under salt treatment.
Improvements in photosynthetic activity have been reported in mycorrhizal plants growing under saline conditions. It has been found that AM inoculation improved the photosynthetic efficiency of maize, primarily though increasing gas exchange capacity and the PSII efficiency and regulating the energy distribution between photochemical and non-photochemical processes 19 . Similar conclusions were drawn by Lin et al. 20 in experiments with Leymus chinensis seedlings. Another study showed that AM inoculation enhanced photosynthetic efficiency by improving stomatal conductance and protecting photochemical events of PSII against salt stress 10 . In this study, photosynthesis was one of the most enriched pathways and included genes related to photosystem II, photosystem I and the light-harvesting chlorophyll-protein complex (Fig. S8a,b; Table S4). Enzymes involved in carbon assimilation such as 5-bisphosphate carboxylase/oxygenase (Rubisco) also play a key role in photosynthesis process 21 . Mo et al. 22 discovered that AM inoculated watermelon seedlings had higher water-use efficiency and Rubisco activity than NM seedlings under drought treatment. In contrast, data about Rubisco enzyme in AM plants under salt treatment is limited. The present study showed genes related to carbon fixation in photosynthetic organisms were also significantly enriched ( Fig. S8c; Table S4). All of the above, plus we found that AM inoculation increasing photosystem II efficiency (ΦPSII) values, carbohydrate content, and the activities of ADP-glucose pyrophosphorylase (AGPase) and starch synthase, and a decrease in non-photochemical quenching www.nature.com/scientificreports www.nature.com/scientificreports/ of chlorophyll fluorescence (NPQ) (Fig. S9), led to the conclusion that AM fungi enhance salinity tolerance of S. cannabina by maintaining photosynthetic ability.
Salinity stress can induce the formation of ROS, its excessive accumulation cause oxidative stress 23 . Several reports have proved that AM inoculation enable host plants to tolerate saline stress by elevating the activities of antioxidant enzymes including superoxide dismutase (SOD), catalase (CAT), glutathione reductase (GR), and peroxidase (POX) 16,[24][25][26] . In the current study, the most enriched GO functional term (p < 0.05) in the module that was positively correlated with both saline treatment and AM inoculation was 'oxidation-reduction process' . Some of the unigenes were verified by qRT-PCR, and the increased expression levels of genes related to SOD, CAT, GR and POX in AM plants compared with NM plants indicated that AM fungi enhanced the ROS-scavenging capacity in S. cannabina plants. To date, the studies about effects of AM inoculations on nonenzymatic antioxidants accumulation in the host plant are limited. Ascorbate is the most important antioxidant in plants and, in association with other components of the antioxidant system, protects plants against oxidative damage resulting from aerobic metabolism and photosynthesis. Carotenoids function as photoprotectants by quenching ROS before oxidative damage can occur, or by active non-photochemical quenching (NPQ), or heat dissipation, of excess light energy 27,28 . Our data showed that pathways associated with ascorbate and aldarate metabolism and carotenoid biosynthesis were significantly enriched ( Fig. S10; Table S4), which could reflect key mechanisms of acquired systemic salt tolerance in AM S. cannabina.
To sum up, this study showed that tolerance of S. cannabina to salinity stress was improved by AM inoculation. In order to identify genes essential in the acquired salinity tolerance of AM plants, we investigated transcriptome data from S. cannabina tissues under salt stress using Illumina sequencing and WGCNA analysis. This analysis revealed one coexpression module that responded to salt stress and AM inoculation in shoots. The results showed that genes related to photosynthesis, ROS scavenging and specific transcription factors played central roles in the plant salt tolerance acquired through AM symbiosis. Our work may facilitate plant functional genomic studies, such as by suggesting candidate genes as potential markers of tolerance to salt stress.

Materials and Methods
plant material and treatments. S. cannabina (Retz.) Pers. seeds from a single population (Shandong Academy of Agricultural Sciences, China) were surface sterilized (10% hydrogen peroxide, 10 min) before use. The seeds were pregerminated and sown in autoclaved zonolite. Glomus mosseae (BGC NM03D) was propagated with Trifolium repens plants for 60 days, infected roots, hyphae, spores and substrate was used as inoculum. After 2 weeks, each seedlings were transferred to a PVC pots (1 L) with a partition in the middle, and the roots of each seedling were divided and buried on both sides of the partition with autoclaved zonolite; on one side the roots were inoculated with 10 g inoculums containing approximately 121 spores. S. cannabina plants were cultured in 12 h photoperiod, 25/17 °C (day/night), and 600 μmol m −2 s −1 light intensity. Three-week-old plantlets were used for all treatments with 200 mM NaCl solution. The 27-h group was treated first and then the 3-h group was treated 24 h later. After the treatments were performed, all samples from the three treatments (0, 3, and 27 h) were harvested at the same time point, frozen immediately in liquid N 2 and stored at −80 °C before use.

Quantification of Arbuscular mycorrhizal fungal colonization.
Percentage of mycorrhizal roots colonization was measured according the gridline intersection method after samples were stained with trypan blue as previously described 14 . RNA isolation and sequencing. Total RNA extraction and purification, and cDNA library construction were performed as previously described 29 . Samples from three seedlings were pooled for a single biological replicate; total RNA from two biological replicates was isolated. RNA-Seq library was constructed using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ispawich, USA). Then, Illumina HiSeq. 4000 platform were used for RNA-Seq libraries sequencing, 150-bp paired ended reads were generated. It was carried out by Novogene Bioinformatics Technology Co. Ltd (Tianjing, China).
Assembly and annotation. The clean reads were obtained by removing low-quality regions and adapter by NGS QC Toolkit (v2.3) 30 , for all following analyses. Trinity was used for transcriptome assembly with min_ kmer_cov set to 2 by default 31 . Functional annotation of unigene were based on seven public databases, Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Swiss-Prot (a manually annotated and reviewed protein sequence database), KOG/COG (Clusters of Orthologous Groups of proteins), KO (KEGG Ortholog database), Pfam (Protein family) and GO (Gene Ontology), using the Blast2GO software with a cutoff E-value of 10 −6 .
Quantification and differential analysis of gene expression levels. First, gene expression levels were estimated by RSEM 32 , the clean data were mapped back onto the transcriptome assembly using bowtie2 by default (mismatch 0). The read count for unigenes was obtained and normalized to fragments per kilo base of transcript sequence per millions base pairs sequenced (FPKM). Second, differential expression between two treatments was analyzed using the DESeqR package (1.10.1) 33 . The P-values were adjusted by Benjamini's approach 34 for controlling the false discovery rate (FDR); FDR adjusted P < 0.05 and absolute value of log2 ratio > 1 were set as the threshold for significant differential expression. GO enrichment analysis of DEGs was carried out using the GOseq R package based on Wallenius' non-central hyper-geometric distribution 35 , which can adjust for gene length bias in DEGs. The statistical enrichment of DEGs in KEGG pathways 36,37 was performed by the KOBAS software.
www.nature.com/scientificreports www.nature.com/scientificreports/ Construction and Visualization of Co-expression network. Gene co-expression networks were constructed by WGCNA package (v1.29) in R 38 . Genes with averaged NRPKM values from two replicates > 2 were used for the WGCNA unsigned network construction, and the average NRPKM values were imported into WGCNA. Modules were obtained using the automatic network construction function blockwiseModules with default settings, except that the power of nine was chosen to make the networks showed an approximate scale-free topology (model fitting index R 2 = 0.9), TOM-Type was signed, minModuleSize was 100, and merge Cut Height was 0.25. The edge weight (from 0 to 1) of any two genes was determined by the topology overlap measure provided in WGCNA. Connectivity was defined as the sum of the weights across all the edges of a node, and the top 3% of genes with the highest connectivity in the network were defined as hub genes 39 .
To identify modules that were significantly associated with the acquired salt tolerance in AM plants, the module eigengene was calculated in each module, and then correlated with each sample (saline treatments, AM fungi inoculation and tissues). Only modules with P < 0.05 were considered related modules. Coexpression patterns and interactions of hub genes were visualized using Cytoscape 40 . qRt-pCR analysis. Total RNA was extracted from shoots using TRIzol (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions, with three biological replicates each sample. The reverse transcription reactions were performed using the GoScript ™ Reverse Transcription System (Promega, Madison, WI, USA).
The quantitative real-time (qRT)-PCR analysis was conducted with an iCycler iQ real-time PCR detection system (BIO-RAD). Each reaction was completed in a total volume of 20 μl including 0.2 μM primers, 2 μl diluted cDNA, and 10 μl 2 × SYBR Green PCR Master Mix (TaKaRa Bio Inc., Dalian, China). Gene-specific primers were showed in Table S5. Actin and tubulin genes were used for normalization (Table S7). Relative gene expression levels of each gene were calculated using the 2 −△△ Ct method. photosynthetic parameters. Gas exchange and modulated chlorophyll fluorescence parameters were simultaneously detected using a LI-6400XTR Portable Photosynthesis System equipped with a Fluorometer (Li-Cor, Lincoln, NE, USA). Efficiency of Photosystem II (ΦPSII) and non-photochemical quenching of chlorophyll fluorescence (NPQ) were calculated according to a previously described method 15 .
Measurement of total carbohydrate content. Root tissue was used to measure the total carbohydrate content according to previously described method 29 . In brief, non-structural carbohydrates were extracted from 100 mg non-AM root tissues. The extract was hydrolyzed by adding HCl (3 M) boiling for three hours. Then it was neutralized by adding Na 2 CO 3 . Total carbohydrate content was measured using anthrone-sulfuric acid colorimetric method, the absorbance was monitored at 630 nm.
Measurement of enzyme activities. The extraction and investigation of the activities of the AGPase and starch synthase were performed as previously described 29 , with three replicates for each treatment.

Data Availability
Raw sequencing data, as well as the assembly was deposited in NCBI GEO database under accession number GSE99532. All data generated or analysed during this study are included in this published article [and its supplementary information files].