A total of 44 diverse Miscanthus genotypes comprising 41 M. sinensis (NO. 1 to NO. 41), 2 Miscanthus floridulus (NO. 48 to NO. 49) and 1 M. × giganteus (NO. 50) species used in this study were mainly collected from the Southwest of China (Supplementary Table S1). Each genotype was propagated to six individuals (three for control and three for treatment) by dividing rhizome and grown in plastic pots (20 cm × 25 cm) with mixture soil (50% loam with 50% fine sandy). All materials were grown in greenhouse with a completely randomized design from April to July in 2014. Before subjected to water stress, all the plants grew with sufficient conditions. The experiment was designed to evaluate the genotypic difference in physiological and photosynthetic aspect responding to water stress. Prior to treat, all the replicates were well watered with the average of 90% SWC measured by Soil Moisture Equipment TDR300 (Santa Barbara, CA, USA). Naturally water stress was applied by stopping water for 30 days. The photosynthetic traits were measured at the end of the drought stress in control and treatment first, and then leaf samples were collected for determination physiological indexes from three replications.

The Pn, Gs, Ci, and Tn were determined by photosynthetic system (LI-6400, Li-Cor, Lincoln, NE, USA) with the CO₂ concentration was 400 μmol mol^-1 and photosynthetic active radiation (PAR) was 1,200 μmol⋅m^-2s^-1. Gs and Ci were measured under the condition that with a saturated light intensity of 1,000 μmol m^-2s^-1. Ten randomly selected leaves were measured to determine related traits for each replicate (Zhao et al., 2014).

The RWC of leaves was measured by drying method (Barrs and Weatherley, 1962). Approximate 0.2 g (FW) fresh leaves randomly selected from each pot were wrapped in absorbent paper and put into the 50mL tube with filled water. Tubes were placed in the dark place for 24 h, and then take out the leaves, dry the surface moisture, and weight SW. Followed the leaf samples were placed in the oven with 75°C for 3 days to constant weight, weight DW. The RWC was calculated as follows: RWC (%) = (FW - DW/SW - DW) × 100. For the extraction of MDA, a total of 0.1 g leaf samples were quickly frozen in liquid nitrogen (2 min), sufficiently grind in 2 volume solutions containing 50 mM phosphate buffer (pH = 7.8) and 1% (w/v) polyvinyl polypyrrolidone. The supernatant was then collected after the homogenate centrifuged at 12,000 rpm for 30 min and the content of MDA was examined (Dhindsa and Matowe, 1981).

Analysis of variance of the physiological and photosynthetic index was conducted by using SPSS 17.0 software (IBM, Armonk, NY, USA) and the Pearson correlation analysis of all traits was performed. Effects of both treatment and genotypes were determined by using the model of Least Significant Difference test (LSD).

A total of six RNA samples from M. sinensis under drought stress were extracted based on the manufacturer’s instructions of the RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA). RNA purity, concentration, and integrity were checked by using the RNA Nano 6000 Kit for the Agilent 2100 Bioanalyzer 2100 System (Agilent Technologies, USA). The qualified samples were primary stored at –80°C after isolation and quality assessment. Each cDNA library was constructed using 5 μg of total RNA and six cDNA sequencing libraries of M. sinensis were constructed using the NEB Next^® Ultra TM RNA Library Prep Kit for Illumina (New England Bio-labs, USA). Initially, the total RNA was treated with RNase-free DNase I (NEB) for 30 min at 37°C and Poly (A) mRNA was isolated using poly-T oligo-linked magnetic beads. After purification, the poly (A)-containing mRNA was fragmented into 200–250 bp pieces using fragment buffer (Ambion), and the first-strand cDNA was synthesized using random hexamer primers and the short fragments as templates. Then the RNase H and DNA polymerase I were applied to the products for second-strand cDNA synthesizing (16°C for 2 h). Finally, before hybridization, NEB Next Adaptor with hairpin loop structure was ligated to the cDNA and the 3′ ends of the DNA fragments were adenylated. Subsequently, the purity and quality of cDNA fragments library were evaluated using the Agilent Bioanalyzer 2100 system. After generating the clusters on a cBot System using the TruSeq PE Cluster kit v3-cBot-HS (Illumina), Illumina sequencing (paired-end technology in Illumina HiSeq 2000 platform) of the six libraries was performed for RNA-seq analysis. The de novo assembly of RNA-seq was conducted using Trinity¹ with employing the sequence splicing and redundancy removing. In total, 114,747 assembled unigenes were directly used to identify SSRs based on the software MISA² (MIcroSAtellite). SSR search criteria were conducted based on perfect mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide motifs minimum number of six, five, five, four, and four repeats, respectively. Primer3 v2.23 was used to design primers in the flanking regions of the SSRs.

A total of 93 SSR markers (Table 1) were randomly selected to validation in 54 plant materials including 44 M. sinensis, 5 M. floridulus, 1 M. sacchariflorus, 1 M. × giganteus, and 3 Hemarthria cultivars (Supplementary Table S1). The genomic DNA was extracted from the young leaves using plant DNA Genomic Extraction Kit (Tiangen Biotech, China). PCR amplifications were performed in 20 μl total volume reactions, including 3 μl genomic DNA (15 ng/μl), 0.4 μl Golden DNA Polymerase, 10 μl 2X Reaction Mix, 0.5 μl of each primer, and 5.6 μl ddH₂O. The PCR program was as follows: an initial denaturation at 94°C for 1 min, followed by 5 cycles of 1 min annealing at 60°C, and 1 min extension at 72°C, then followed by 35 cycles of 1 min denaturation at 94°C, 1 min annealing at 55–58°C, and 1 min extension at 72°C. A final 10 min extension at 72°C was used before the PCR reactions were completed. The obtained PCR products were separated on 6% PAGE under 350 volts for 2 h and gel was stained by AgNO₃ solutions.

Population structure analysis was performed using STRUCTRE v2.3.4 software (Pritchard et al., 2000) with the “admixture model”. Briefly, the software default model parameters were set with 10,000 iterations and 10,000 replications of MCMC. 20 independent runs were applied in this study with varying the assumed number of genetic groups (K) from one to six for each run. Maximum likelihood and delta K (ΔK) values were used to determine the optimum number of subgroups (Pritchard et al., 2000; Evanno et al., 2005). In clustering analysis, the UPGMA dendogram were construct by using similarity coefficients based on SAHN module in the NTSYS-pc version 2.10 software.

The association analysis was conducted by using TASSEL 2.1 software to reveal the association between drought related traits and marker alleles (Bradbury et al., 2007). Four models (simple linear, Q, K, and Q+K) were tested by using Quantile–quantile (QQ) plots to identify the best model fitting drought tolerance related traits for association mapping in M. sinensis populations. The significant threshold for selecting the associations between alleles and traits was set at P < 0.01. Phenotypic variation explained (R²) indicated the fixed marker effects.

Relative water content determines plant water status when exposed to drought stress. As a physiological indicators that directly response to the extent of plant water deficiency in drought environments, all plants demonstrated high level of leaf RWC in control conditions with an average RWC of 88.27% while evident effect in drought stress with a decrease to 43.21% (Table 2). The leaf RWC significantly difference among genotypes based on ANOVA analysis under drought stress which ranged from 12.44 to 66.11%, indicating that leaf RWC could well separate the drought resistant and drought sensitive genotype of Miscanthus. The MDA content directly reflected the stability of cellular membrane in the leaves. After a period of drought stress treatment, the MDA content largely increased than control from 11.71 to 54.63 nmol/g FW and showed an extensive variation under drought conditions (19.62–143.78 nmol/g). For the photosynthetic productivity, Pn, Gs, Ci, and Tn were all decreased under drought stress, and significant differences among genotypes were observed under control and drought stress (Table 2).

As an undomesticated new bioenergy crop, several characteristics related to drought tolerance successfully distinguish Miscanthus from many other crops. Since high biomass yield is the primary goal by meeting Miscanthus when all above ground biomass is harvested annually, the drought tolerance assessment criterion seems clearly establishment. However, for optimum and stable productivity, at least a 3 years establishment is required after transplanting of Miscanthus, which is a typical perennial grass species (Clifton-Brown and Lewandowski, 2000; Anzoua et al., 2015; Nie et al., 2016). Directly evaluating the biomass productive of Miscanthus in the field is time and labor consuming under drought conditions. Greenhouse experiment is favorable for detecting the mechanism responding to drought stress which can be identified some reliably measured trait to predict future performance for M. sinensis resistant breeding programs.

Leaf RWC reflecting the plant water status, MDA content representing the extent of cellular membrane peroxidation and Pn, Gs, Ci, and Tn comprehensively showing the plant biosynthesis and catabolism process, which provide rapid and easy measurements for screening drought-tolerant plant materials and characterize drought tolerance response mechanism (O’Neill et al., 2006; Jansen et al., 2009; Jiang et al., 2009). Significant correlations were observed between Pn, Gs, Ci, and Tn under well-watered and drought conditions, respectively (Table 3). As expect, Pn significant positively correlated with Gs and Tn, and negatively correlated with Ci under drought stress. Leaf RWC significantly correlated with Pn, Gs, Ci, and Tn under well-watered conditions, but such correlations were not observed under drought stress. MDA had no significant correlation with any other measured traits under control, but significantly correlated with RWC under drought stress. The correlation analysis results showed that there is no significant correlation between the physiological and photosynthetic responds to drought stress, indicating that varied and different regulating mechanism was existed in the panel of Miscanthus accessions.

In order to avoid damage caused by drought stress, one hand the plant reduced water transpiration and improved water use efficiency to protect the photosynthetic mechanism through stomatal closure, on the other hand, adjusted osmotic pressure to increase water uptake and maintained the integrity of the cell structure and function (Yu et al., 2013). In this study, a significant correlation between RWC and MDA content were detected under drought condition and the decreasing of RWC and increasing of MDA content significantly different from that of the control, potentially indicating the occurrence of drought stress, suggesting that they are sufficient parameters of drought tolerance in greenhouse experiment. Generally, drought induced ABA synthesis, which causing stomatal closure to prevents water loss through transpiration and play an important role in maintaining leaf water status (Yu et al., 2013). Previously studies had been showed that M. × giganteus exhibited little stomatal regulation ability under mild drought stress compared to M. sinensis (Clifton-Brown and Lewandowski, 2000). However, the significant correlation between RWC and Pn, Gs, Ci, or Tn were just detected in well watered condition in this study, indicating that Miscanthus holding moisture could not through the stomatal closure pathway but antioxidant pathways to maintain the cell structure and function under drought stress.

In addition, different phenotypic responses of Miscanthus to drought stress are important for marker-trait association. The extensive variation observed among the evaluated genotypes regarding physiological and photosynthetic respond indicating that the large range of comprehensive drought tolerance among the evaluated genotypes. The genotypes tested may become important candidates for investigating of mechanisms of plant respond to drought stress and could possibly be used in breeding programs. However, the potential variability in short-term versus long-term drought response should be taken into account in order to further develop new varieties of drought tolerant.

Simple sequence repeats are considered as one of the most useful molecular markers due to their broad range of applications in genotype identification, genetic map construction, and marker-trait association. Using the Illumina HiSeq 2000 platform, a total of 349,393,396 raw reads were generated (accession numbers: SRP095822 in the NCBI SRA database) and a total of 316,200,846 high-quality clean reads were assembled into 114,747 unigenes after filtering and trimming the raw reads. Using MISA script, the microsatellites were searched across all unigenes and a total of 16,566 SSRs were identified with motif lengths ranging from one to six bp (Table 4). The tri-nucleotide repeats, di-nucleotide repeats and mono-nucleotide repeats were the three most abundant SSRs types in this study and were responsible for 54.63%, 23.84% and 12.81% of the total amount of SSRs, respectively (Table 4). However, the mono-nucleotides were excluded in the subsequent analysis owing to the presence of homopolymer affecting the sequencing quality (Gilles et al., 2011). The number of the given repeat unit of SSRs ranged from 4 to >10 and with the number of repeat units increased, the frequency of the given SSR structure progressively decreased (Table 5). As for the two most abundant repeat motif types (di- and tri-nucleotides), the frequency of AG/TC motif type account for 52.00% in di-nucleotide repeat motifs, and that the frequency of CCG/GGC was the most abundant motif type in the tri-nucleotide, accounting for 41.80% (Table 6).

The results were very consistent with previous studies that the tri-nucleotide was the most abundance repeat among the types of SSRs (Kaur et al., 2012; Yates et al., 2014; Huang et al., 2016; Tian et al., 2017). In comparison to other nucleotide repeats, the tri-nucleotide expansions always lead to a homopolymeric amino acid run in an unaltered protein due to the retention of the original reading frame in translated regions (Kaur et al., 2012). However, the amplification of repeats other than tri-nucleotide is prevented in coding DNA, resulting in frame-shift mutations (Metzgar et al., 2000), which could be a possible explanation for the abundance of tri-nucleotide repeats. In the other hand, the AG/TC of di-nucleotide SSRs and CCG/GGC of tri-nucleotide SSRs were reported to be the largest number of motif in Hemarthria and grain crops (Temnykh et al., 2000; Kantety et al., 2002; Huang et al., 2016). The inverse GA/CT motif can describe multiple codons when transcribed into mRNA. Furthermore, a high proportion (8 and 10%) of amino acids Ala and Leu were translated from inverse GA/CT motif contributing it widely distribution across the genome (Kantety et al., 2002).

To validate the SSRs markers identified by RNA-seq, a total of 93 SSR primer pairs were randomly selected for validation including different motifs. Validation of these SSRs was based on the different size among 51 different Miscanthus genotypes and 3 Hemarthria species (Supplementary Table S1). In total, 70 (75.3%) SSRs were successfully amplified whereas 30 have polymorphic, including 9 for di-nucleotide, 16 for tri-nucleotide and 5 SSRs with motif lengths greater than 3 (Table 7). A total of 110 polymorphic loci generating from 30 SSRs were determined to be used for cluster analysis. The UPGMA dendrogram constructed by genetic similarity (GS) data obviously showed that three clusters were formed when the GS coefficient was 0.52. The group 1 accessions were mainly Miscanthus sinensis collected from the south of Sichuan and Chongqing which latitude was lower than N 31°. The M. sinensis genotypes of group 2 were primarily collected from the North of Sichuan and other genotypes of M. floridulus, M. sacchariflorus, M. × giganteus were also assigned into this groups. All three Hemarthria cultivars were classed to Group 3 as a foreign species (Figure 1). These results from the cluster analysis suggested that the Miscanthus SSR markers identified in this study have a high level of polymorphism clearly differentiating the Miscanthus species. More importantly, they had been proved having high transferability in the Andropogoneae tribe.

Previous studies on Miscanthus marker developing are limited to the number of available for use, and always markers from Tribe Andropogoneae (Poaceae) C4 crops such as sorghum (Sorghum bicolor L. Moench), maize (Zea mays L.), and sugarcane (Saccharum officinarum L.) were used for genetic studies (Ho et al., 2011; Kim et al., 2012; Chae et al., 2014; Yook et al., 2014; Nie et al., 2016). Although these large number of SSRs have been successfully applied to M. sinensis studies and proven with high transferability (Hernandez et al., 2001; Zhao et al., 2011; Xu et al., 2013; Chae et al., 2014; Yook et al., 2014; Nie et al., 2016), the specificity and function of markers seems less than satisfactory. The SSRs developed from Miscanthus transcriptome under drought conditions were aimed to the function regions among expressed unigenes, which could be facilitate better discovering of the marker-trait associations and higher efficiency.

Population structure of 44 genotypes was analyzed by using STRUCTURE V2.3.3 software. A total number of 110 filtered marker loci (by dropping the loci with minor allele frequency less than 5%) retained for structure and association analysis. Two optimum subgroups were determined by maximum likelihood and delta K (ΔK) values (Figure 2). Accordingly, 24 genotypes were assigned to group 1, 14 genotypes were to group 2, and 6 genotypes were referred not clearly assigned to the specific groups based on a significant membership threshold (Q-value) of 0.60 (Figure 3).

The association analysis was performed by taking Q (population structure) and K (relative kinship) into consideration using GLM and MLM model. Drought resistant related traits measured in this study under drought and control condition were used to test the model of Q, K, Q+K and simple linear model in QQ plots (Figure 4). Generally, in most conditions, the behavior of Q+K model was the best approximation to the excepted cumulative distribution of P-values. By implementing the Q+K models, a significant improvement was observed when compared with the other models. Based on Q+K model, a total of seven loci were detected significantly associate with RWC, Pn, Tn, Gs, and Ci (P < 0.01) under drought stress. Among them, 2 SSRs developing from CL12059.Contig2_All_5665_1 and CL10784.Contig3 _All_5453_1 were significantly associated with RWC, 1 SSRs from CL10592.Contig1_All_5404_1 associate with Pn and Gs simultaneously, 1 SSRs from CL11502.Contig2_All_5573_1 associate with Tn, 1 SSRs from CL10112.Contig3_All_5301_4 associate with Gs and 1 SSRs from CL10164.Contig1_All_5323_2 associate with Ci (Table 8). The contribution of single significant associated markers to the phenotypic variation was varied ranging from 11.89 to 25.30% (Table 8).

Molecular markers have been widely used to analysis the genetic relationship of Miscanthus all around its distributed areas (Xu et al., 2013; Clark et al., 2014; Yook et al., 2014; Anzoua et al., 2015; Nie et al., 2016). Although some genetic maps (Kim et al., 2012; Ma et al., 2012; Swaminathan et al., 2012; Liu et al., 2015) and association mapping (Zhao et al., 2011; Slavov et al., 2014; Nie et al., 2016) with high density and resolution have been constructed, mapping of QTLs using the genetic map are limited owing to its absent of genomic information. Zhao et al. (2013) used 23 SSR markers transferred from Brachypodium distachyon to conduct marker-trait association analysis and 115 loci were generated. Among them, nine markers were detected to be significantly (P < 0.01) associated with heading date and biomass yield. Slavov et al. (2014) conducted genome-wide association studies in a 138 M. sinensis population by using single-nucleotide variants markers and a total of 17 significant associations were detected with phenology, morphology, and cell wall composition traits. Nie et al. (2016) used 104 pair of markers to conducted the marker-trait association analysis and among 1059 generated loci, 12 significant associations of biomass yield related traits were identified. In this study, among 110 loci generated from 30 randomly selected SSRs, 7 significant associations were identified with RWC and photosynthetic trait under drought stress indicating that the developed SSR markers have a relative stable linked efficiency to the drought tolerance related traits.

By searching for the function annotation of 7 unigenes accounting for the associated markers in the protein database (Supplementary Table S2), an definite answer was appeared. Especially, unigene CL10592.Contig1_All_5404_1 and CL10112.Contig3_All_5301_4 associated with Pn and Gs were annotated function related to vegetative cell wall protein, while unigene CL12059.Contig2_ All_5665_1 and CL10784.Contig3_All_5453_1 associated with RWC were annotated related to arginine/serine-rich protein and ubiquitin-like-specific protease. With the abundant genetic basis information, the developed markers could be potential candidates for improving Miscanthus selection for enhanced drought tolerance in a breeding program.