Transcriptomic characterization of candidate genes responsive to salt tolerance of Miscanthus energy crops

Given the growing need for biofuel production but the lack of suitable land for producing biomass feedstock, development of stress‐tolerant energy crops will be increasingly important. We used comparative transcriptomics to reveal differential responses to long‐term salt stress among five populations of Miscanthus lutarioriparius grown in the natural habitats and salinity experimental site. A total of 59 genes were found to be potentially responsive to the high‐salinity conditions shared by the five populations, including those involved in detoxification, plant defense, photosynthesis, and signal transduction. Of these genes, about 70% were related to abiotic stress response. Among five populations, the most contrasting performance between relatively high survival rates and the relatively weak growing traits was in accordance with the down‐regulation of genes involved in growth and up‐regulation of genes related to plant stress tolerance in one of the populations. These results might reveal a potential tolerance‐productivity trade‐off, where resources were allocated from growth to stress resistance. The comparative transcriptomics of different populations among different environments will provide a basis for breeding and domestication of energy crops.


Introduction
With the increasing demand for fuel production from renewable resources, the development of second-generation energy crops capable of growing on marginal land becomes increasingly urgent (Sang & Zhu, 2011;Allwright & Taylor, 2016). Given that salinity is a major adverse environmental factor affecting plant growth and productivity (Boyer, 1982), it is important to improve salt tolerance of energy crops. Considerable efforts have been undertaken to elucidate salt-responsive mechanisms of plants (Flowers & Yeo, 1995;Zhu, 2001;Zhang et al., 2004;Brinker et al., 2010;Cherel et al., 2014;Bushman et al., 2016). However, understanding of plant responding to long-term salt stress in the field conditions that could facilitate energy crop development remained limited.
Over the past decades, intense studies were devoted to understanding the mechanisms of plant responding to salt stress using model plants under short-term (several hours) salt stress in controlled laboratory or greenhouse conditions (Munns & Termaat, 1986;Zhu et al., 1998;Taji et al., 2004;Fujii & Zhu, 2009;Sun et al., 2010;Wang et al., 2015). Although this has been a powerful approach for revealing detailed molecular mechanisms, mechanistic study of salt tolerance of plants growing in salinity soil under long-term field conditions is needed to close the gap for the development of salt-tolerant crops (Zhu et al., 1998;Brosch e et al., 2005;Vicente et al., 2016). In such efforts, plants grown in the field conditions do represent a valuable resource for elucidating the tolerance mechanisms (Brosch e et al., 2005). The importance of acclimation process to long-term salt stress in the field is even more crucial for perennial grasses, as they face repeated episodes of abiotic stresses during their life cycles (Brosch e et al., 2005;Moinuddin et al., 2014). Thus, a better understanding of long-term salt tolerance mechanisms under field conditions should benefit breeding salt-tolerant perennial grasses.
Generally, plants responding to various abiotic stresses often show a notable alteration in gene expression at the transcriptional level (Kreps et al., 2002). Highthroughput RNA sequencing has been increasingly adopted to monitor global gene expression changes and identify candidate genes by comparative analyses of transcriptional profiles between normal and stressed conditions (Walia et al., 2005;Wang et al., 2009;Beritognolo et al., 2011;Alvarez et al., 2015). It is particularly useful in studying nonmodel plants whose reference genome sequences are not available (Trick et al., 2009;Libault et al., 2010a,b). Moreover, according to recent researches exploring the performance of native plants in adverse environments, transcriptomic analysis has been proven to be an effective approach to discover candidate genes with drastically altered expressional levels (Morozova & Marra, 2008;Wang et al., 2009;Champigny et al., 2013).
Several studies aiming to determine which physiological parameters best reflecting the response of Miscanthus to salt stress showed that salinity resulted in a reduction in plant growth and photosynthetic rates (Plazek et al., 2014). Nevertheless, Miscanthus could still maintain a relatively high growth rate and biomass accumulation compared to other plants (Plazek et al., 2014). Miscanthus lutarioriparius, an endemic species in central China, is considered to be a promising candidate of second-generation energy crops due to its capability of producing high biomass on the semiarid marginal land (Sang & Zhu, 2011;Yan et al., 2012;Liu & Sang, 2013;Mi et al., 2014;Fan et al., 2015;Xu et al., 2015). It is self-incompatible and reproduces through both cross pollination and clonal growth by rhizomes (Chen & Renvoize, 2006). The adaptation of M. lutarioriparius to the semiarid environment could be related to the change of expressional patterns of candidate genes responsible for abiotic stress resistance and the improvement of water use efficiency (Fan et al., 2015;Xu et al., 2015).
In this study, we sequenced transcriptomes of 50 oneyear-old M. lutarioriparius individuals grown in the salinity experimental field in Dongying, Shandong Province of China (DS). These individuals were collected from five natural populations (NP), and seeds were planted in the experimental field. The gene expression of 50 individuals was compared with those of 40 individuals whose leaves were directly sampled from these five natural populations for transcriptomic studies. The comparative transcriptomic study allowed for the identification of candidate genes responsive to long-term salt tolerance and one underlying mechanism of toleranceproductivity trade-off.

Sample collection
A total of 40 individuals of M. lutarioriparius from five populations across the natural habitats of the species (NP) and 50 individuals from same populations grown at the salinity experimental site (DS) were sampled in June and July of 2013, respectively. Our previous study had predicted a high degree of adaptability in M. lutarioriparius based on the genetic analysis of 644 individuals from 25 populations along the Yangtze River, which span the whole geographic range of the species (Yan et al., 2016). To investigate whether these populations have salt tolerance and reveal the underlying mechanism of salt tolerance at the transcriptomic level, we focused on five of these populations of M. lutarioriparius which represent different habitats in the middle reaches of Yangtze River in this study. They are LU5, LU7, LU10, LU14, and LU9, which are from Hekou, Jianli, Honghu, and Jiayu of Hubei Province and Junshan of Hunan Province, respectively. The population identities used here are corresponding to codes in Yan et al. (2016). Seeds of the five populations of M. lutarioriparius were collected in NP and planted in the experimental field in Dongying, Shandong Province (DS), with high salt content in April 2012 (Fig. 1). The five populations of M. lutarioriparius were planted with randomized block design to reduce the variety of salinity in the saline field. We investigated the survival rate, plant height, and net photosynthetic rate of each sampled individual in DS. To reveal the intrinsic expression patterns of underlying adaptation in natural populations, the five populations in NP had been sampled for transcriptome sequencing around noon of 18-22 June, 2013, with eight individuals for each of the populations (J. Yan, Z.H. Song, J. Greimler, J.Q. Li, T. Sang, unpublished). To obtain samples of DS harvested at a similar growing phase with that of NP, the collecting time was accordingly adopted. Due to the seasonal climate conditions, the same temperature was found about one month later in DS than in NP. Ten individuals from each of the populations grown in DS were sampled for transcriptome sequencing around noon on 24 July, 2013. For all individuals, the fourth leaf from the top of each individual was cut and immediately placed in liquid nitrogen and stored at À80°C for further analysis.

RNA-Seq, preprocessing, and annotating RNA-Seq data
Total RNA of leaves was extracted using Qiagen Plant Mini Kit (Qiagen, Stanford, CA, USA) and purified using the RNase-free DNaseI (TaKaRa, Otsu, Shiga, Japan) following the manufacturer's protocol. The concentration of purified RNA was quantified by NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and stored at À80°C until next step. The mRNA was isolated from 5 lg purified total RNA using one round of purification with oligo d (T) beads Dynabeads â mRNA Purification Kit (Invitrogen, Carlsbad, CA, USA). The cDNA libraries were prepared with the NEB-Next Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). The first strand cDNA was synthesized using random hexamer-primed reverse transcription. Following the second-strand cDNA synthesis and adaptor ligation, approximately 450-bp cDNA fragments were isolated by the selection of Ampure XP beads (Beckman Coulter, Brea, CA, USA). The isolated cDNA fragments were amplified by 10 cycles of PCR. Library integrity and quality were estimated with Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and Qubit 2.0 fluorometer (Life Technologies, Grand Island, NY, USA), respectively. The libraries were subsequently sequenced on the Illumina HiSeq 2500 system as 2 9 100-bp paired-end reads.
Raw data were sorted into the corresponding individual according to their indexed nucleotides using bcl2fastq-1. 8.4 (http://support.illumina.com/downloads/bcl2fastq_conver sion_software_184.html). The unstable sequences content in the first nine bases of the 100-bp reads were trimmed in all samples.
To control the quality, raw reads were filtered based on quality scores (Q = 20) and trimmed the indexed using FASTQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/) and FASTX (http://hannonlab.cshl.edu/fastx_toolkit). Although there was no reference genome sequence of M. luparioriparius, Xu et al. (2015) assembled a high-quality reference transcriptome (TSA accession no. GEDE00000000). To measure the expression level of each gene in each sample, trimmed reads of each individual were mapped to the Bowtie-build indexed reference transcriptome using TopHat v2.0.0 with default settings (Langmead et al., 2009;Trapnell et al., 2009;Xu et al., 2015). Expression level of each sample was calculated using FPKM standing for fragments per kilobase of exon per million fragments mapped with Cufflinks v2.0.2 (Trapnell et al., 2010). The trimmed sequence data are available at Sequence Read Archive (SRA) at NCBI under Project ID (SRP068901). By mapping reads from each sample onto the reference transcriptome, we obtained the number of genes expressed in each sample and the expression level for each gene.
SAMtools was used to identify single-nucleotide polymorphisms (SNPs) with default settings (Li et al., 2009). We eliminated SNPs with quality score ≤10 and minor allele frequency (MAF) ≤0.05 to ensure the accuracy of SNPs. After the lowquality reads were discarded, the remaining SNPs were retained for analysis. Genetic differentiation (F ST ) was estimated based on those SNPs (Wright, 1978).

Expression pattern analysis
We used Xu et al.'s method (Xu et al., 2015) to describe the level and diversity of gene expression at the population level. Population expression level (E p ) was calculated as the mean of FPKM values of the given individuals to descript the mean level of gene expression in each population. Considering the unequal number of samples collected from two sites, the extent of variability in relation to the mean of the population of gene expression in each population (expression diversity, E d ) was calculated using the formula (E d ¼ P n i¼1 jEiÀEpj nÀ1 ð Þ Ep ). The above n represents the number of individuals sampled and E i represents the FPKM of a given gene of the ith individual. The change of E p was further examined by the ratios of E p s (the E p in DS divided by the E p in NP).

Identification of population-specific and shared responsive genes
To more accurately detect population gene expression either shared or responding differently to long-term salt stress, we used both t-test and a fold change ranking (Cui & Churchill, 2003;Lee et al., 2012;Parker et al., 2016). The differentially expressed genes between the two sites were identified using the parametric t-test (normal bimodal distribution) or the nonparametric Wilcoxon test (nonnormal unimodal distribution). Gene expression change with P < 0.05 and larger than twofold change was considered to be statistically significant. Those genes significant differentially expressed in all of the five populations were considered as shared responsive genes. To further estimate the pattern of gene expression of shared responsive genes (down-regulated and up-regulated), we compared the E p s and E d s of these genes with those of all genes that expressed in both sites, respectively.

Functional annotation of population shared and differently responsive genes
Hierarchical cluster analysis of shared responsive genes in all individuals was carried out using MULTI EXPERIMENT VIEWER (MEV) v4.9 software using Spearman's rank correlation (Saeed et al., 2003). These normalized genes' FPKM values of individuals of two sites were put together as input data for MEV v4.9. The sequences of shared responsive genes were searched against the Arabidopsis database, using the defaulted TAIR BLASTn parameters to analyze the detailed functional categorization. The expected value (e-value) threshold was set at 10 À10 when we did the BLAST search. Functional classifications of genes expressed significantly different in each population were achieved using AgriGO (Du et al., 2010). KOBAS 2.0 was used to analyze the pathway annotation and enrichment of the differentially expressed genes in each population and chose false discovery rate (FDR) <0.05 as the threshold (Xie et al., 2011).

Sequencing quality, gene expression of five M. lutarioriparius populations
Based on RNA-seq, 50 M. lutarioriparius individuals collected from five populations in the experimental field DS generated approximately 1 369 882 768 raw reads and 128.9 Gb of 100-bp paired-end reads after quality control (Table S1). The detail data size information for each individual in DS is shown in Table S1. After excluding genes with average FPKM value of zero in each population, the number of expressed genes for five populations in DS ranged from 17 315 to 17 560.
Then, we compared the data with that of 40 individuals of the five populations in NP. The further population transcriptome analyses were based on 16 678 genes expressed in five populations in both environments. We compared the E p s for all genes from NP to DS, and we found that E p s in DS shifted significantly toward to higher levels, suggesting an overall increase in expression from NP to DS (Wilcoxon's test, P < 0.01; Fig. S1a). We also compared the change of E p s for each of the five different populations from NP to DS separately. From NP to DS, the proportion of genes that up-regulated was highest in LU5 and lowest in LU9 (67.4%, 62.8%, 62.1%, 61.4%, and 48.2% genes were up-regulated in LU5, LU10, LU7, LU14, and LU9, respectively; Fig. S1bf), suggesting the majority decrease in expression from NP to DS in LU9.

Responsive genes shared in five populations of M. lutarioriparius
The salinity of NP ranges from 0.55& to 1.65& (Harmonized World Soil Database v1.2, FAO/IIASA/ISRIC/ ISSCAS/JRC 2012), while the salinity of DS is around 7&, about seven times higher than that of NP on average. To accurately identify salt stress-responsive genes shared by different populations of M. lutarioriparius, we conducted t-test/Wilcoxon's test and a fold change method to screen genes significantly up or down-regulated in five populations (Fig. S2). As a result, 9945 genes were significantly differentially expressed in only one of the five populations (P < 0.05 and twofold change; Fig. 2a), while 375, 172, and 89 genes were significantly differentially expressed together in two, three, and four of the five populations, respectively (P < 0.05 and twofold change; Fig. 2a). Moreover, there were 59 differentially expressed genes shared by five populations of M. lutarioriparius, in which 46 genes were upregulated and 13 genes were down-regulated. The E p s of the 59 responsive genes shared in five populations of M. lutarioriparius tended to be significantly higher than that of all genes in both sites (Wilcoxon's test, P < 0.01; Fig. 2b, c). E d s of 59 responsive genes shared in five populations of M. lutarioriparius tended to be lower compared with that of the all genes in both environments (Fig. 2d, e). The ratio of E d s between DS and NP for the 59 responsive genes shared in five populations was significantly lower than that of all genes (P < 0.01; Fig. S3). Together these results demonstrated the conserved expression of these shared genes in response to salt stress.
Three distinct clusters were obtained based on hierarchical cluster analysis of those genes using Spearman's rank correlation ( Fig. 3a; Table 1). In addition, all 59 responsive genes shared in five populations were placed into several main functional categories based on the results of BLASTn against the Arabidopsis database. The major functional category was involved in plant defense (23.73%), followed by photosynthesis (15.25%), cellular metabolism (13.56%), signal transduction (8.47%), and detoxification (8.47%). Some of the responsive genes shared in five populations (5.08%) were unable to assign functions (Fig. 3b). A detailed functional categorization for 59 responsive genes shared by five populations is shown in Table 1. Of these genes, about 61%, 69%, and 90% genes were reported abiotic stress-responsive genes in clusters I, II, and III, respectively, with an average of about 70% (Table 1). Interestingly, of these genes, about 33%, 11%, and 56% of reported photosynthesis genes were within clusters I, II, and III, respectively (Table 1).

Population-different phenotype and candidate genes responding to salt stress
Pair-wise F ST estimated with SNPs from the transcriptomes between the five populations ranged from 0.028 to 0.049 (Table S2). Then, we monitored the survival rate and trait performances of each population respectively in response to salt stress to screen the variation of salt stress tolerance. The average survival rate for the majority of populations was smaller than 50% (Fig. 4). And there was significant difference for their survival rates among the five populations, with the survival rate 45.9, 38.3, 31.6, 29.4, and 26% in LU14, LU9, LU5, LU7, and LU10, respectively (Fig. 4). Moreover, plant height and photosynthetic rate of five populations were significantly different from each other (P < 0.05). The average plant height of LU10, LU5, LU14, LU7, and LU9 was 255, 243.3, 218.3, 209, and 201.7 cm, respectively (Fig. 4). The average photosynthetic rates of LU7, LU10, LU14, LU5, and LU9 were 35.0, 33.7, 33.0, 31.7, and 29.6 lm m À1 s À1 , respectively (Fig. 4). Therefore, among the five populations, LU9 showed the most contrasting performance between relatively high survival rates and relatively weak growing traits including plant height and photosynthetic rates in salinity soil. LU10 showed the contrary performance with relatively lower survival rates and relatively high growing traits for survival individuals in salinity soil.
We further monitored the expression performance of each population respectively in response to salt stress to screen the candidate genes. Among the 16 678 shared expressed genes in all five populations in two environments, 1774 (10.6%), 239 (1.4%), 981 (5.9%), 2015 (12.1%), and 912 (5.5%) population-specific responsive genes (P < 0.05 and twofold change) were detected in LU5, LU7, LU9, LU10, and LU14, respectively. The E p s and E d s of these genes were quite different among populations (P < 0.05; Fig. S4a-d). The E p s of LU9 exhibited the widest range in both two environments. What's more, for LU9 the E pS were largest in NP and the smallest in DS (P < 0.05; Fig. S4a, b), and the E dS were largest in NP and smallest in DS (P < 0.05; Fig. S4c, d). LU9, with most of the genes down-regulated, was found to be with the significantly widest range of E p ratio from NP to DS (P < 0.05; Fig. 5a). LU9 was also found to be with the significantly lowest E d ratio from NP to DS (P < 0.05; Fig. 5b).
The enriched GO terms in population-specific responsive genes were listed in Table S3. The significantly enriched GO terms of population-specific responsive genes in LU9 were involved in the up-regulation of defense response, the down-regulation of photosystem, the down-regulation of isopentenyl diphosphate biosynthetic process, mevalonate-independent pathway, and the pentose-phosphate shunt in the category of biological process ( Fig. 6; Table S3). Enriched GO terms of genes between the up-and down-regulated genes were similar among other four populations (Table S3). Thus, LU9 showed the most typical expression coordination, in which the decreased expression of genes involved in normal growth was accompanied by the increased expression of genes related to stress tolerance. To improve our understanding of the important biochemical pathways of population-specific responsive genes, we further analyzed the KEGG pathway of those genes in each population. With the KO annotation of expressed genes in each population as background, four pathways that were photosynthesis (ko00195), protein processing in endoplasmic reticulum (ko04141), phenylpropanoid biosynthesis (ko00940), and circadian rhythm plant (ko04712) had significant enrichment between two environments in LU9 (Table 2). Two pathways that were ribosome (ko03010) and aminoacyl-tRNA biosynthesis (ko00970) had significant enrichment between two sites in LU5 and LU10 (Table 2). One pathway that  Table 1. Table 1 The accession, annotation, and potential functional groups of responsive genes shared in five populations was performed using blastn of TAIR. The clusters were grouped using hierarchical clustering method. The direction of expression regulation is shown by plus sign (up) and minus sign (down), respectively. Different stress types and corresponding reference of the reported abiotic stress-responsive genes in those genes were also listed Gene               was phenylpropanoid biosynthesis (ko00940) had significant enrichment between two environments in LU7 (Table 2). One pathway that was ribosome (ko03010) had significant enrichment between two sites in LU14 (Table 2).

The important role of shared responsive genes in stress responses
To avoid getting false-positive results of responsive genes, only those showing significantly differential expression in all the five populations (P < 0.05 and twofold change) were identified as candidate genes. This is necessary for studying plants responding to stress in the field conditions where variables are relatively difficult to control. The finding that responsive genes shared by five populations tended to have lower E d s and higher E p s compared with all genes (Fig. 2b-e) may indicate that those shared responsive genes tend to be more conserved among individuals than others. They are expected to have experienced strong positive or purifying selective constraint (Hannah et al., 2006;Des Marais et al., 2012;Rengel et al., 2012;Lasky et al., 2014). This result was also in agreement with previous reports that E d s of genes with higher E p s tended to be more conserved in changing environment , which had been explained that highly expressed genes experienced stronger purifying selection than those at lower level. The results that the decrease of E d s of shared responsive genes following M. lutarioriparius being transplanted from its native habitats NP to DS in high saline field reflected the dominant trend of narrowed range of gene expression in the new environment. Thus, these shared responsive genes may be considered as consistently expressed salt-responsive genes and likely represent adaptive responses to common saline conditions (Lasky et al., 2014). Indeed, these shared genes are also found to play an important role in stress response in the related species. As we found that the functions of those genes were mainly involved in detoxification enzymes, plant defense and photosynthesis, and up to 70% of these genes were reported to be stress-responsive genes. For example, detoxification enzymes can reduce excessive reactive oxygen species (ROS) accumulation as various environmental stresses such that salinity universally causes oxidative stress with an excessive accumulation of ROS in plants such as Soldanella alpina and Arabidopsis thaliana (Kiyosue et al., 1994;Sunkar et al., 2003;Engel et al., 2006). Therefore, the complex regulation of the mRNA of catalase (CAT) (Engel et al., 2006), aldehyde dehydrogenase (ALDH) (Sunkar et al., 2003), and epoxide hydrolase (EH) (Kiyosue et al., 1994) seemed to constitute a detoxification mechanism that limits excessive ROS accumulation in M. lutarioriparius leaves. Plant defense was found to frequently respond to various stress environments in other plants (Park et al., 2005;Lockwood et al., 2010), which is according to our result that genes encoding DJ-1/PfpI proteins family and heat shock proteins (HSPs) were up-regulated in DS. Photosynthetic response to salinity stress is extremely complex because salinity stress indirectly affects photosynthesis and includes expression of those genes connected to growth inhibition or leaf shedding. Through limiting water consumption, plant could help to maintain the carbon assimilation (Chaves et al., 2009). In this case, down-regulated genes in M. lutarioriparius encoding photosynthesis-related proteins such as lightharvesting chlorophyll a/b-binding protein 5 (LHCP5) may function in the same way. Consistently, the previous studies also showed that plants responding to salt stress displayed complex molecular responses including the production of ion transporters and detoxification enzymes to re-establish homeostasis and resume growth in Arabidopsis (Zhu, 2001).

Fig. 6
Gene Ontology (GO) term enrichment result for biological process based on GO terms for population-specific responsive genes in LU9. The directed graph shows the relationship of the GO terms over-represented. All colored boxes are enriched with adjusted P value < 0.05. The degree of enrichment was shown from red (adjusted P value < 10 À7 ) to white (10 À1 < adjusted P value). Arrows stand for relationships between different GO terms. Red arrow stands for the up-regulated expression, and green arrow stands for the down-regulated expression. The potential tolerance-productivity trade-off mechanism in response to salt stress The different performances existing among different populations or species may be used for dissecting different salt tolerance mechanisms (Walia et al., 2005;Platten et al., 2013). This study analyzed the expression pattern of population-specific responsive genes in five populations. The E d s' range of these genes narrowed down only in LU9 when M. lutarioriparius was transplanted from NP to DS (Fig. S4c, d), suggesting that the expression for these genes in LU9 experienced a purifying selection (Brawand et al., 2011) or responded to long-term salt stress. Another possible explanation for this might be that these genes would be under a more systematic and strict regulation when encountered with salt stress environment (Zhu, 2001;Chen et al., 2002).
In addition, the different performances of phenotype were also observed among the five populations (Fig. 4). The strongest growth traits and lowest survival rate of LU10 might be a result of selection under which a relatively small portion of individuals with outstanding phenotype was able to survive (Randerson & Hurst, 2001). Therefore, the surviving individuals with such phenotype could have a great potential for energy crop development in salinity soil. In further studies, salt tolerance and growth traits of M. lutarioriparius can be selected separately for pyramiding breeding (Handa et al., 2014). The detailed response mechanisms could be revealed more clearly by associating the expression pattern with phenotype (Atwell et al., 2010). First of all, the responsive genes were analyzed for GO enrichment in each population separately, and the results showed that the up-regulation of stress response, the down-regulation of photosystem II assembly, isopentenyl diphosphate biosynthetic process, mevalonate-independent pathway, and the pentose-phosphate shunt were enriched in LU9 ( Fig. 6; Table S3). Moreover, pathway enrichment analysis showed that photosynthesis, protein processing in endoplasmic reticulum, and phenylpropanoid biosynthesis were enriched in LU9 (Table 2). Phenylpropanoid biosynthesis (ko00940) (Vogt, 2010) and the up-regulation of defense response might have contributed to the salt tolerance of LU9. The pentosephosphate shunt pathway that was well known as the fundamental metabolic pathway was down-regulated (Wood, 1986). And the down-regulation of photosystem II assembly, isopentenyl diphosphate biosynthetic process, mevalonate-independent pathway, and photosynthesis (Lichtenthaler et al., 1997) might have contributed to the lowest photosynthetic rate in LU9 (Fig. 4). The down-regulation of function in photosynthesis and cellular metabolism that consequently limited carbon dioxide uptake and compromised plant growth might be one cause of the weak phenotype in LU9 under salt stress (Chaves et al., 2009).
The contrasting performance of LU9 between relatively high survival rate and relatively weak growing traits including plant height and photosynthetic rate was in accordance with gene expression patterns in such a way that the decreased expression of genes involved in normal growth was accompanied by the increased expression of genes related to plant stress tolerance. These results revealed the underlying mechanism of tolerance-productivity trade-off, that is, the molecular basis for redirecting resources from growth to stress resistance (Dupont-Prinet et al., 2010;Muller-Landau, 2010). The coordination of decreasing the expression of genes related to normal growth and increasing the expression of genes involved in plant stress-responsive mechanisms could have contributed to allocation of resources from rapid growth to stress protection (Sabreen & Sugiyama, 2008;Lopez-Maury et al., 2009;Zakrzewska et al., 2011) and resulted in higher survival rates at a cost of other physiological attributes (Dupont-Prinet et al., 2010).

The findings of response of M. lutarioriparius to long-term salt stress in natural environment
We compared the findings of our study with those studies of poplar and Arabidopsis under short-term controlled environment salt stress. The CAT gene family involved in antioxidant defense in poplar (Ding et al., 2010) and SOS pathway that salt stress induced calcium-signaling pathway in Arabidopsis (Chinnusamy et al., 2004) were the well-known salt stress-associated genes which were revealed under short-term (0-48 h) controlled environment salt stress. We also found genes such as CAT involved in antioxidant defense and CBL that was identified as one of the calcium sensors in calcium-signaling pathway, which is in coincidence with the findings in short-term controlled environment (Engel et al., 2006;Zhang et al., 2008). And the expression of catalase was only up-regulated under long-term salt stress treatments, which was consistent with the previous study (Hernandez et al., 2010). Compared with findings of studies under short-term salt stress, the long-term study could also have a better observation of survival rates and growing traits.
We also compared the findings of our study with those studies of Miscanthus, pea, and wheat under controlled long-term salt stress. We found in M. lutarioriparius leaves, genes involved in the function of detoxification, plant defense, photosynthesis, and signal transduction were potentially responsive to the highsalinity conditions. One of the M. lutarioriparius populations had the potential tolerance-productivity trade-off, in which down-regulation of cellular metabolisms, photosynthesis genes, and the up-regulation of stressrelated genes were coordinated with the contrasting performance between relatively high survival rate in salinity soil and relatively weak growing traits including plant height and photosynthetic rates ( Fig. 4; Fig. 6). The previous studies found that Miscanthus x giganteus and Miscanthus sinensis tended to have a reduction in plant growth and photosynthetic rate in response to salt stress under controlled hydroponic conditions (Plazek et al., 2014), which were in accordance with the relatively weak growing traits. In the expression profile, the induction of antioxidant defense (Hernandez et al., 2000) and the reduction in photosynthesis-related genes (Kiani-Pouya, 2015) were also found in pea and wheat responding to long-term (15-60 days) salt stress under controlled environment. Those findings focused on investigating the different responses between salt-tolerant and salt-sensitive genotypes facing long-term salt stress. Salt-tolerant genotypes had stronger growing traits or higher expression of antioxidant defenserelated genes than salt-sensitive genotypes. However, only our study did find the link between gene expression and growth traits.
In this study, we aimed to understand how plants respond to long-term salt stress by conducting comparative population transcriptomic analyses between native habitat and salinity field. Through identifying candidate stress-responsive genes shared by populations, we demonstrated that consistently expressed responsive genes were likely to represent adaptive responses to common stress conditions and genes shared for stress responses appeared to have been subject to purifying selection. Based on the population which had the most contrasting performance between survival rates and growing traits, we found this contrasting performance was in accordance with the up-regulation of genes involved in stress tolerance and down-regulation of certain growth-related genes. The finding suggested that there might be a tolerance-productivity trade-off when plants were exposed to long-term salt stress in the field conditions, where resource allocation from growth to stress resistance might be the mechanism at the expressional level for stress tolerance. Fine tuning of the tolerance-productivity trade-off could potentially contribute to the development of second-generation energy crops that have satisfactory establishment rates and productivity in salinity land.

Supporting Information
Additional Supporting Information may be found online in the supporting information tab for this article: Figure S1. Comparison of E p s of M. lutarioriparius between field sites. (a) Distributions of E p s in all individuals sampled from native habitat (NP) and transplanted site (DS). The E p was calculated as the mean of FPKM values of all individuals in one site.
(b-f) Distributions of E p s in five populations between two sites, respectively. The E p was calculated as the mean of FPKM values of the individuals in each population. The black line indicates genes expressed in DS, while the black dotted line indicates genes expressed in NP. Figure S2. The volcano plot displays the relationship between fold-change and significance between the two sites in five populations respectively (a-e). The y-axis is the -log 10 P values (a higher value indicates greater significance) and the x-axis is the log 2 E p ratio between two experimental fields. Horizontal red line indicated the significance threshold (P < 0.05). The vertical red line identified as the two-fold change threshold. Figure S3. Comparison of E d s ratio of responsive genes shared in five populations and all genes expressed in both sites. Figure S4. Comparison of E p s and E d s of genes expressed significantly different in populations and all genes in the populations of M. lutarioriparius between two sites. (a, b) Distributions of log 2 E p of population-specific responsive genes and all genes of LU5, LU7, LU9, LU10, LU14 in DS and NP, respectively. (c, d) Distributions of E d s of population-specific responsive genes and all genes of LU5, LU7, LU9, LU10, LU14 in DS and NP, respectively. Each group was further divided based on location in DS (white box) or NP (gray box). Table S1. Data size and coverage of reference transcriptome of the M. lutarioriparius individuals used in our study. Table S2. The population genetic differentiation (F ST ) based on SNP data. Table S3. Significantly enriched GO terms of down-or up-regulated population-specific responsive genes in each population. Functional classifications of differentially expressed genes of each population were achieved using AgriGO. The P value cut-off of 0.05 was chosen in order to control the false discovery rate.