De novo assembly of the Japanese lawngrass (Zoysia japonica Steud.) root transcriptome and identification of candidate unigenes related to early responses under salt stress

Japanese lawngrass (Zoysia japonica Steud.) is an important warm-season turfgrass that is able to survive in a range of soils, from infertile sands to clays, and to grow well under saline conditions. However, little is known about the molecular mechanisms involved in its resistance to salt stress. Here, we used high-throughput RNA sequencing (RNA-seq) to investigate the changes in gene expression of Zoysia grass at high NaCl concentrations. We first constructed two sequencing libraries, including control and NaCl-treated samples, and sequenced them using the Illumina HiSeq™ 2000 platform. Approximately 157.20 million paired-end reads with a total length of 68.68 Mb were obtained. Subsequently, 32,849 unigenes with an N50 length of 1781 bp were assembled using Trinity. Furthermore, three public databases, the Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-prot, and Clusters of Orthologous Groups (COGs), were used for gene function analysis and enrichment. The annotated genes included 57 Gene Ontology (GO) terms, 120 KEGG pathways, and 24 COGs. Compared with the control, 1455 genes were significantly different (false discovery rate ≤0.01, |log2Ratio |≥1) in the NaCl-treated samples. These genes were enriched in 10 KEGG pathways and 73 GO terms, and subjected to 25 COG categories. Using high-throughput next-generation sequencing, we built a database as a global transcript resource for Z. japonica Steud. roots. The results of this study will advance our understanding of the early salt response in Japanese lawngrass roots.


Background
Plant growth, development, and production depend largely on soil; however, more than 800 million hectares of land worldwide are subjected to salt stress. High-salinity soils prevent the reabsorption of water from plant roots, the first organ that responds to salt stress (Zhu, 2002;Bartels and Sunkar, 2005;Munns, 2005;Petricka et al., 2012); therefore plants can be poisoned by the excess uptake of salts. The level of sensitivity of the plant salt stress response plays a crucial role in determining resistance to high salt stress (Tracy et al., 2008;Schmidt et al., 2013;Wei et al., 2013). Given the current scarcity of water resources, recycled water, and water use have been the focus of many studies (Marinho et al., 2013). Improving the abiotic stress resistance of plants is an effective method for overcoming water shortages and for the cultivation of large areas of saline land.
A focus on early detection would assist the determination of plant stress responses (Schmidt et al., 2013). Model plants such as rice have evolved a regulatory network of early responses to chilling stress (Yun et al., 2010). As a vital signaling component in many biological processes (Mittler et al., 2004;Schippers et al., 2012), reactive oxygen species (ROS) can rapidly induce abiotic stress responses-H 2 O 2 levels induced by salt stress rise within several minutes in rice (Hong et al., 2009). Understanding the physiological and molecular mechanisms underlying salt stress tolerance in halophytes is therefore of great importance. Several halophytes, including Thellungiella halophila, Mesembryanthemum crystallinum, Suaeda, and Populus, have been studied (Dyachenko et al., 2006;Sun et al., 2010;Fukao et al., 2011;Yoon et al., 2014); however, our understanding of the unique tolerance mechanisms of halophytes is limited. Thus, further study of various halophytes may provide new information regarding the ability of plants to resist salt damage.
Zoysia Willd (family Poaceae, subfamily Chloridoideae, tribe Zoysieae) are perennial grasses, consisting globally of about 10 recognized species (Tsuruta et al., 2011). Zoysia grasses are widely used as a warm-season turfgrass for home lawns, golf courses, athletic fields, and parks (Ge et al., 2006). The species belonging to this genus are the most salt and cold tolerant of the C 4 grass species in the family Poaceae. They are indigenous to China, Japan, and Korea. As one of the three most important commercial species, Zoysia japonica Steud. also exhibits marked tolerance to abiotic stress. It is distributed naturally in mountainous areas, along riversides, and in coastal areas, and it shows moderate tolerance to weak shade. Further, it can survive in various soils, ranging from infertile sands to clays (Tsuruta et al., 2011).
Few genomic sequence resources are available for Zoysia grasses. Only 11 mRNA sequences (complete cDNAs) are available in the database of the National Center for Biotechnology Information (NCBI) for Z. japonica (as of April 2015). This grass has an advanced transformation system (Asano, 1989;Inokuma et al., 1998;Ge et al., 2006), but the lack of sequence resources has limited the exploitation of Z. japonica's genetic resources.
To determine the molecular mechanism underlying the salinity tolerance of plants (Uddin et al., 2012), a large-scale analysis of differentially expressed genes (DEGs) is necessary.
Sequencing technology, which is used to analyze the transcriptome of a species without complete genome information, has undergone rapid progress in recent years with respect to the number of base calls and cost. Therefore, the technology is used in studies of many species and by researchers in many fields; for example, it has been used to study tissue regeneration in newts (Looso et al., 2013), glucosinolate metabolism in radish , and the leaf transcriptome of Camelina sativa .
In this study, we used the HiSeq ™ 2000 platform to perform RNA sequencing (RNA-seq) of Zoysia grass roots. We compared the transcriptomes of plants grown under saline and normal conditions to identify differences in gene expression, and to identify the function of key transcripts and genes in the initiation of ROS-related signal transduction. Significant expression differences were found among genes involved in many metabolic pathways. Many novel genes were also identified and inferred to be expressed, specifically in the salt-treated plants.
To the best of our knowledge, this is the first report of a transcriptome of Zoysia grass. We used H 2 O 2 as a marker to identify DEGs between the control and salt-treated plants. This will improve our understanding of the mechanisms involved in short-term stress responses of Zoysia. These sequence data may also enhance our understanding of the molecular mechanisms in plants under salt stress, and provide a public dataset for use in future studies of Zoysia.

Plant Materials and Treatment
Z. japonica Steud. cv. Zenith was grown from seed in soil in a propagation tray. Three weeks after germination, individual seedlings were transplanted to pots (diameter: 15 cm, depth: 14.5 cm) filled with a mixture of topsoil and coarse river sand (1:1) in a greenhouse (25 • C during the day/20 • C at night, 16 h of light/8 h of dark, 800 µmol m −2 s −1 photosynthetically active radiation, and 75% relative humidity). The plants were irrigated with water during the first 3 weeks to keep the soil moist. After that, they were transferred to pots and maintained with 1/2 Murashige and Skoog cultivation solution once a week and water irrigation twice a week. NaCl treatment (150 mM) was initiated 3 months after germination.

H 2 O 2 Assay
To visualize H 2 O 2 accumulation, samples (from control and salt-treated plants) were immediately placed in a 1 mg/ml 3,3 ′diaminobenzidine (DAB)-HCl solution in 10 mM phosphate buffer (pH 3.8) at 25 • C for 4 h in the dark. After staining, the root of each plant was boiled in 95% (v/v) ethanol for 30 min and rehydrated in 70% (v/v) ethanol for 48 h at 25 • C. H 2 O 2 was visualized as a reddish-brown coloration. Each experiment was repeated using 10 plants.

RNA Preparation and Library Preparation for Analysis
Total RNA was extracted using TRIzol reagent according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, USA).
The extracted RNA was treated with RNase-free DNase I (Takara Inc., Kyoto, Japan) for 45 min at 37 • C to remove residual DNA. The quality of the RNA was evaluated using a NanoDrop 2000. The cDNA was prepared by pooling 10 µg of RNA each from the control and salt-treated samples.
Total RNA was extracted from normal and salt-treated plant roots. Poly A+ mRNA was obtained using a NEBNext Poly(A)mRNA Magnetic Isolation Module. Then, according to the instructions of the NEBNext mRNA Library Prep Master Mix Set for Illumina and NEBNext Multiplex Oligos for Illumina, a mixed cDNA library of salt-stressed (Case) seedlings and control (CK) plants was prepared.

Sequencing, De novo Assembly, and Annotation
The cDNA library was sequenced using the Illumina HiSeq ™ 2000 platform. After cleaning the raw reads and discarding lowquality reads, we ran Trinity (Yoon et al., 2014) to assemble the clean reads into transcripts. The clean reads were used to construct a K-mer dictionary and each K-mer was used as an initial contig. The most frequent K-mer, which was in the dictionary with K-1 overlaps with the current contig end, was extended until neither direction could be extended further. The contigs that shared at least one K-1-mer were then read across the junction sites to build a pool. Each contig was used to construct a de Brujin graph with trimmed spurious edges and compacted linear paths. It then reconciled the graph with reads and pairs, and output one linear sequence for each splice form and/or paralogous transcript represented in the graph. The final selection of the most important section from the partial transcript was the unigene. To annotate sequences obtained by de novo assembly, we used BLASTX with a significance threshold of E ≤ 10 −5 . The assembly unigenes were aligned against the plant protein datasets of the non-redundant protein (NR), Swiss-prot protein, TrEMBL, Gene Ontology (GO), Clusters of Orthologous Groups (COGs), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.

Identification and Functional Annotation of DEGs
We calculated the expression levels as reads per kilobase exon model per million mapped reads (RPKM) for each root sample. DEGseq software was then used to determine significant DEGs defined as a fold change ≥2 and false discovery rate (FDR) <0.01. Using BLAST, DEGs were aligned against the NR, GO, COGs, and KEGG databases.

Gene Comparison and Sequence Alignment within Relative Species
We previously analyzed published results on the salt response of model plants that are most homologous to Zoysia, including Setaria, barley, rice, and maize (Ueda et al., 2004(Ueda et al., , 2006Qing et al., 2009;Puranik et al., 2011). We used Blastn with a significance threshold of E ≤ 10 −10 , sequence identified ≥30% and the length ≥30 aa. The DEGs were blasted against the salt-response genes that were identified in relative plants.

Quantitative Reverse Transcription PCR (RT-qPCR) Verification
We subjected nine ERF salt stress-related unigenes (Comp132815_c0, Comp211710_c0, Comp196695_c0, Comp 233109_c5, Comp223421_c0, Comp217707_c0, Comp206 878_c0, Comp206878_c1, and Comp195676_c0) to RT-qPCR analysis. Root samples were collected from additional plants grown and tested under identical conditions. As per the Invitrogen protocol, total RNA was extracted from the root samples and then subjected to reverse transcription. The primers were designed using an online tool (http://www.idtdna.com/site). ZjActin (GenBank: GU290545.1) was used as a housekeeping gene. RT-qPCR was performed on a Bio-Rad RCR platform using a SYBR R Green Real-Time PCR Mix (TransStart Green qPCR SuperMix) to detect transcript abundance. Amplification was performed as follows: denaturation at 95 • C for 30 s, followed by 40 cycles of denaturation at 95 • C for 5 s, annealing at 60 • C for 15 s, and extension at 72 • C for 10 s. All amplifications were performed with three replicates. We calculated the relative expression levels of the selected unigenes using the 2 − Ct method. The data from the reactions were analyzed using Bio-Rad software.

Phylogenetic Tree Construction
Each of nine ERFs had more than one homologous genes in the PlnTFDB database. We selected 21 homologous genes (Supplementary File 1) from that database to perform phylogenetic analysis. Multiple sequence alignments were performed using the Clustal X program with default parameters. Phylogenetic analysis was performed by aligning a maximumlikelihood phylogenetic tree using MEGA 6 software, which was supported by a bootstrap test with 1000 iterations.

Putative Molecular Markers
To investigate the distribution of simple sequence repeats (SSRs), the MIcroSAtellite identification tool (MISA version 1.0; http:// pgrc.ipk-gatersleben.de/misa) was used to search the assembled transcripts. We scanned and counted the SSRs in all selected unigenes. The parameters were adjusted for the identification of perfect mono-, di-, tri-, tetra-, and penta-nucleotide motifs with a minimum of 12, 6, 5, 5, and 4 repeats, respectively.

Physiological Responses in the Roots of Salt-treated Plants
Most previous studies focused on long-term adaptations to salt stress in grasses (Wang and Jiang, 2007;Bian and Jiang, 2009;Du et al., 2009;Hu et al., 2012); few have investigated the initiation of the stress response that plays an important role in this type of stress signaling. Changes in the environment of plants increase the level of ROS (Hong et al., 2009), with the production of ROS being an early response to abiotic stress (Zhu, 2002;Mittler et al., 2004). H 2 O 2 is the predominant salt-induced ROS (Yang et al., 2007); thus, it was used as a marker of the degree of salt stress in root samples.
As shown in Figure 1, the H 2 O 2 concentration in roots increased shortly after salt treatment. Compared with the control plants, a yellow substance was visible on the roots. After 30 min of salt treatment, oxidative stress caused by H 2 O 2 was visible throughout the roots. Thus, we selected representative samples at 0 and 30 min for transcriptomic sequencing.
Illumina Sequencing and De novo Assembly of Zoysia Root Transcripts RNA-seq technology is an indispensable tool for the wholegenome analysis of complex stress treatments. Compared with traditional large-scale sequencing, de novo whole-genome analysis is less costly and more efficient. It is suitable for those plants whose genomic sequences are unknown. We used RNA sequencing to analyze the transcriptome of Z. japonica Steud. roots.
To obtain a comprehensive transcriptome, two cDNA libraries denoted "CK" and "Case" prepared from three repeat RNA samples from normal and salt-treated roots were subjected to paired-end read sequencing using the Illumina platform. Pairedend read technology increases the depth and improves de novo assembly efficiency. Sequencing using the Illumina HiSeq ™ 2000 platform resulted in the generation of 21.30 G raw reads. More than 80% had Phred-like quality scores at the Q30 level (error < 0.1%). After removing reads with adaptors, reads with unknown nucleotides, and low-quality reads, we obtained 105.47 million clean reads with an average GC content of 55.83% (Table 1A). The Trinity method can generate full-length transcripts without reference genomes (Grabherr et al., 2011;Haas et al., 2013). Using this method for de novo assembly, all highquality clean reads were assembled into 719,182 contigs (Supplementary File 2) with an average length of 201.56 bp. Contigs of 100-500 bp were most frequent, accounting for 94% of the total. Subsequently, the contigs were clustered into 32,849 unigenes of which the mean length was 1107 bp and the N50 value was 1781 bp (Table 1B). There were 20,499 unigenes of ≥500 bp, and 5224 unigenes of ≥2000 bp. Most unigenes were in the range 200-500 bp (37.59%). Among these genes, the longest and shortest were 15,772 and 201 bp, respectively. The unigene lengths facilitated annotation and classification. The random distribution of the unigenes is presented in Figure 2.

Functional Annotation and Classification of Assembled Unigenes
The assembled sequences were first searched against the NR protein database, and analysis indicated that 55% of the sequences ranged from 1.0E −5 to 1.0E −50 , while 18,006 of sequences with an E < 10E −50 displayed strong homology (Figure 3). The distribution of identity is shown in Figure 3B. The majority pattern was 60-80% similarity (12,111) and 40-60% similarity (9812). The closest species was Setaria itallica, with 14,950 genes (45%) matched. The next closest species was Sorghum bicolor, which showed 17% homology with Zoysia. This implies that the transcripts were assembled and annotated correctly (Dang et al., 2013;Wang et al., 2013). Figure 3C shows homologies with plants in the family Poaceae. Recent research (Ahn et al., 2015) supports our data showing that the most homologous model species to Zoysia is Setaria itallica, followed by Sorghum bicolor and Zea mays. Due to the differences in tissue between cultivar and treatment, our results were not completely consistent with that research.

Identification of DEGs
We used the RPKM method (RPKM > 0.1) to calculate the expression levels of unigenes in the control and salt-treated samples. Using DESeq software (FDR < 0.01, FC > 2) to analyze the DEGs between salt-and control-treated root samples, a total of 1648 differentially expressed unigenes (Figure 4) were obtained. Among those genes, 948 were upregulated and 700 were downregulated (Supplementary File 4).

GO Classification
To categorize unigenes functionally, we assigned GO terms (P < 0.01). In total, 1,455 unigenes were enriched in 73 GO terms (Figure 5, Supplementary File 5). Among the unigenes, 798 were upregulated, and 647 were downregulated. After analyzing the statistically enriched GO functions related to DEGs, the majority of GO terms were assigned to biological processes distributed in 67 subcategories (91.78%), followed by cellular components and metabolic processes, with three subcategories each (4.11%). The most over-represented GO terms in response to salt stress were "response to stimulus" (294 unigenes), "response to stress" (229 unigenes), "response to chemical stimulus" (224 unigenes), "response to abiotic stimulus" (187 unigenes), and "response to an organic substance" (187 unigenes). Many genes in these GO terms responded to short-term salt stress, suggesting that they play important roles in the early salt response. We also found nine significantly enriched hormone-related GO terms, including ABA-related ("response to abscisic acid stimulus"), JA-related ("response to jasmonic acid stimulus" and "cellular response to jasmonic acid stimulus"), and SA-related ("response to salicylic acid stimulus, " "cellular response to salicylic acid stimulus, " "salicylic acid-mediated signaling pathway, " "salicylic acid metabolic process, " "salicylic acid biosynthetic process, " "systemic acquired resistance, " and "salicylic acid-mediated signaling pathway") terms. These hormones were upregulated by salinity, and they induced genes involved in alleviating salt stress (Wang et al., 2001). A GO terms analysis may help in understanding the salt tolerance mechanism involved when plants suffer a sudden increase in soil salinity.

Calcium Ion (Ca 2+ ) Signaling
Ca 2+ acts as a secondary messenger in salt stress responses (Xuan et al., 2013). The cytosolic free Ca 2+ concentration ([Ca 2+ ] cyt ) changes within seconds in plants subjected to salt stress (Knight et al., 1997), leading to downstream signaling. RSA1, a nuclear Ca 2+ -binding protein, can sense changes in the free Ca 2+ level elicited by salt stress in the nucleus and transduce this signal by activating the SOS1 promoter (Guan et al., 2013). P-type Ca 2+ channels can sense changes in the [Ca 2+ ] cyt and transduce the signal downstream by activating specific targets (Huda et al., 2013). Among the DEGs, we identified 12 unigenes belonging to a major Ca 2+ sensor family; of these, 10 were upregulated and two were downregulated (Table 2A).

ROS Scavenging
Salt stress can cause the rapid accumulation of ROS, including superoxide, H 2 O 2 , and hydroxyl radicals. ROS can perturb cellular redox homeostasis leading to DNA damage and changes in protein or membrane function. Plants possess ROS scavenging systems to control ROS levels and cope with oxidative stress. Z. japonica has a similar ROS scavenging system to that of other plant species. Cold treatment significantly increases antioxidant enzyme levels (Xuan et al., 2013), which may provide protection against oxidative damage in this species (Xu et al., 2013a). Among our DEGs, we identified 14 unigenes belonging to the ascorbate-glutathione cycle, GPX pathway, and Prx/Trx pathway.

Comparison of Salt-responsive Genes in Zoysia and Relative Species
To further characterize the identified genes, we used BLAST to compare 255 DEGs ( Table 3, Supplementary File 6) with previously identified salt-responsive genes in four relative of the model plants. Genes in barley and rice had similar DEGs, and more than 60% of these DEGs were homologous sequences with Zoysia. Only 27 of 296 genes in maize were matched with those in Zoysia. These results suggest that barley and Zoysia use a similar mechanism to respond to salt stress, whereas maize has a different system for this response.

Major Regulators of the Salt-stress Response Transcriptome
Over the past decade, genes related to abiotic stress have been discovered and functionally characterized. Among these genes, TFs are master regulators that control gene clusters . A single TF can regulate the expression of downstream genes by binding specifically to a cis-acting element in the promoter of a target gene. Members of the AP2/ERF, MYB, WRKY, and NAC families were shown to regulate salt tolerance. Enhanced expression of DREB2A can improve saltstress tolerance in rice. Three rice NAC TFs (SNAC1, SNAC2, and NAC5) act as positive regulators of the salt stress response (Schmidt et al., 2013). A total of 3751 unigenes ( Table 4) were identified as potential TFs with an average length of 1224 bp. The length distribution of the TFs is shown in Table 3. More than 43.05% of the TFs were longer than 1000 bp; the shortest TF was 201 bp, while the longest was 11,427 bp. In total, 80 TF families were identified (Supplementary File 7). The largest group was the AP2/EREBP family (187 unigenes), followed by the bHLH family (158 unigenes) and the WRKY family (137 unigenes). There were nine TF families containing more than 100 unigenes.
We selected the TF family that was most closely related to abiotic stress from the DEGs (Supplementary File 8).

The AP2/EREBP Family
Increasing evidence indicates that ERF proteins are involved in the salt response. They adjust the expression of downstream genes and fine-tune crosstalk between signaling pathways (Agarwal et al., 2006;Nakano et al., 2006;Fukao et al., 2011). Many studies have shown that ERFs enhance salt tolerance (Gilmour et al., 2000;Park et al., 2001;Jung et al., 2007). Among the unigenes we identified were 187 ERF genes, including 34 DEGs. About 33 DEGs were upregulated and only 1 DEG (comp223372_c0) was downregulated. In both libraries, comp232359_c0 had the highest expression level, but the expression of this gene was slightly increased in the case libraries. The largest change was found for comp234321_c1, which was upregulated more than four-fold (log 2 fold change).

The bZip Family
Members of the bZIP family of TFs respond to abiotic stress. They have important roles in the ABA response and the regulation of oxidative and pathogen defense responses (Yun et al., 2010). In total, 115 genes were identified in our database, with seven DEGs. Comp219239_c0 had a higher expression level, and its expression differed slightly between the two sample libraries.

The NAC Family
The NAC domain, which was identified based on consensus sequences from the Petunia NAM, Arabidopsis ATAF1/2, and CUC2 proteins, was plant-specific and contained the highly conserved NAC DNA-binding domain and variable C-terminal domains. The NAC family is considered to be important in plant development, senescence, auxin responses, and biotic stress responses (Aida et al., 1997;Kim et al., 2004;Olsen et al., 2005;Lu et al., 2007;Hao et al., 2011;Nakashima et al., 2012;Sun et al., 2012). After salt treatment, the expression of 131 NACs differed, and only 11 unigenes were considered to be DEGs. The transcript of comp231035_c0 was present at a high level in FIGURE 5 | GO classification of the unigenes. The results are summarized for three main categories: biological process, cellular component, and molecular function. In total, 1445 DEGs with BLAST matches to known proteins were assigned to GO categories. The WRKY Family As one of the best-characterized TF families, the WRKY TF family has been suggested to play an important role in plant stress responses. The WRKY protein family contains a conserved amino acid sequence motif, WRKYGOK, at the N-terminus and a novel zinc-finger-like motif at the C-terminus. Numerous studies have demonstrated that WRKY family regulates various plant processes from development to biotic and abiotic stress responses (e.g., wounding, drought, and salinity). Eight WRKYs in wheat respond to low temperature, salt stress, and heat treatment (Wu et al., 2008). Enhancing the expression of WRKY25 and WRKY33 can increase the salt tolerance of Arabidopsis thaliana (Jiang and Deyholos, 2009;Li et al., 2011). The overexpression GmWRKY54 increases salt and drought tolerance in soybean. H 2 O 2 treatment significantly induces the expression of AtWRKY30, AtWRKY75, AtWRKY48, AtWRKY39, AtWRKY6, AtWRKY53, and AtWRKY22 (Vanderauwera et al., 2005;Miao and Zentgraf, 2007;Zhou et al., 2011). In our study, we identified 137 WRKYs, including 20 DEGs. Of these, seven WRKYs (comp231105_c0, comp224918_c0, comp232021_c2, comp230256_c0, comp216034_c0, comp218505_c0, and comp190635_c0) had high transcription levels in both libraries.
The MYB Family MYB proteins are the largest class of TFs in plants. Nearly 9% of the TFs in Arabidopsis are MYBs, and more than 1600 TFs have been identified (Riechmann et al., 2000). MYB proteins  (Cominelli et al., 2008). The overexpression of MdoMYB121 in tomato and apple confers improved tolerance to drought (Cao et al., 2013). TaMYB73 has been shown to enhance salt resistance by regulating the expression of related genes (He et al., 2012). A total of 129 genes were identified; among them, only 14 showed significant differences in expression.

The bHLH Family
TFs belonging to the bHLH family are important in plant development, circadian regulation, and stress responses. Here, we found 158 bHLH genes; 19 genes from this family were DEGs.

Putative Molecular Markers
Molecular markers are important for molecular breeding (Liu et al., 2012;Chen et al., 2013Chen et al., , 2014Dai et al., 2013;Xu et al., 2013b;Li et al., 2014). SSRs are 1-6-bp iterations of DNA sequences that occur only in noncoding regions. The occurrence of SSRs in transcribed sequences has been established. The roles of SSRs in plants such as rice (Cho Yg et al., 2000), bread wheat (Gupta et al., 2003), and sugarcane (Cordeiro et al., 2001) have been reported. In our study, 32,849 unigenes were used to detect SSRs, and a total of 4842 SSRs were identified in 3951 unigenes using MISA (Supplementary File 9). Among them, 702 unigenes contained at least 2 SSRs. The largest fraction was mononucleotide SSRs (2305), followed by tri-nucleotide SSRs (1694), and di-nucleotide SSRs (778). This phenomenon corresponds to natural selection. Molecular markers are a useful resource for determining functional genetic variation. Our results will facilitate the prediction of molecular markers for Zoysia.

RT-qPCR
We identified nine ERFs with significant changes; their lengths varied from 260 to 873 bp (Figure 6). PCR amplification showed that all RT-qPCR primers (Supplementary File 10) used produced only single fragments of the expected lengths (100-250 bp; 100% success rate). All positive clones from the validation studies were subjected to Sanger sequencing, and the results confirmed those obtained using the Illumina method. Our qPCR results for the nine unigenes are in agreement with those from the DESeq analysis of our RNA-seq data.

Phylogenetic Analysis
To further predict and distinguish the function of the nine unigenes, phylogenetic analysis was performed using 21 known plant AP2/EREBP protein sequences that were homologous to these unigenes. The analysis revealed the presence of five distinct clusters (Supplementary File 11), suggesting that these five groups might differ functionally in some respect. Comp206878_c1, Orza_sativa_subsp._indica_OsIBC009072, and Zea_mays_GRMZM2G129674_P01 belong to class I, whereas comp196695_c0, Zea_mays_GRMZM2G175525_P01, and Sorghum_bicolor_5041828 belong to the class II. In the phylogenetic tree, comp132815_c0, comp223421_c0, and comp211710_c0 were not closely related to the known AP2/EREBP genes, and comp233109_c5, comp217707_c0, comp195676_c0, and comp206878_c0 had a low bootstrap value. If these seven unigenes are novel ERFs or novel genes, the full lengths of these unigenes can be used further to inform our phylogenetic tree.

Conclusions
Here, we presented the first comprehensive transcriptome data of Z. japonica Steud. roots. In total, 32,849 unigenes were identified. The large number of transcripts identified will serve as a global resource for future studies. TFs play a crucial role in the early response to salt stress; thus, we identified candidate TFs related to salt stress, which will be the subject of future studies. In addition, a total of 4842 SSRs were identified. This information will facilitate future studies of plant biology and molecular breeding.

Author Contributions
HL conceived and designed the experiments and contributed the reagents. XQ and NJ wrote the manuscript. XQ, ZY, and FB performed the experiments. XX conducted the TF analysis. YS conducted the DEG analysis. ZL and SX conducted the H 2 O 2 assay. XL and LX provided the materials and analytic tools.

Availability of Supporting Data
The RNA sequence dataset supporting the results in this article is available from the NCBI Sequence Read Archive. The accession number is SRX838917.