RNA-seq Analysis of Cold and Drought Responsive Transcriptomes of Zea mays ssp. mexicana L.

The annual Zea mays ssp. mexicana L. is a member of teosinte, a wild relative of the Zea mays spp. mays L. This subspecies has strong growth and regeneration ability, high tiller numbers, high protein and lysine content as well as resistance to many fungal diseases, and it can be effectively used in maize improvement. In this study, we reported a Zea mays ssp. mexicana L. transcriptome by merging data from untreated control (CK), cold (4°C) and drought (PEG2000, 20%) treated plant samples. A total of 251,145 transcripts (N50 = 1,269 bp) and 184,280 unigenes (N50 = 923 bp) were predicted, which code for homologs of near 47% of the published maize proteome. Under cold conditions, 2,232 and 817 genes were up-regulated and down-regulated, respectively, while fewer genes were up-regulated (532) and down-regulated (82) under drought stress, indicating that Zea mays ssp. mexicana L. is more sensitive to the applied cold rather than to the applied drought stresses. Functional enrichment analyses identified many common or specific biological processes and gene sets in response to drought and cold stresses. The ABA dependent pathway, trehalose synthetic pathway and the ICE1-CBF pathway were up-regulated by both stresses. GA associated genes have been shown to differentially regulate the responses to cold in close subspecies in Zea mays. These findings and the identified functional genes can provide useful clues for improving abiotic stress tolerance of maize.


INTRODUCTION
In order to survive in nature through evolution, plants developed multiple strategies to adjust to various abiotic stresses by promoting a series of physiological and metabolic processes such as stomatal closure, repression of cell growth and photosynthesis, and activation of respiration (Mishra et al., 2014). Stress induced damages are reduced or removed by improving sugar content, increasing antioxidants, taking mechanical action and inducing molecular chaperone chemicals in plants. The gene regulation networks involved are reviewed elsewhere (Zhu, 2016).
Low temperature and drought stresses are major environmental factors that impact on the geographical distribution and composition of plant species. They can also lead to a decrease in crop quality and productivity (Viswanathan and Zhu, 2002). A number of genes and alternative spliced isoforms responding to these stresses at the transcriptional level have been reported (Leyva et al., 1995;Shinozaki et al., 2003;Thatcher et al., 2016). Transcription factors make up of a large group which extensively involve in these processes (Zhu, 2016). Among them a class of DREB/CBF transcription factors which bind to the DRE/CRT element (A/GCCGAC) play roles in the regulation of the expression of target genes in response to cold and drought in Arabidopsis (Nakashima et al., 2009). DREB/CBFlike genes were significantly up-regulated by both stresses of cold and drought (Agarwal et al., 2006). AaDREB1 from Adonis amurensis, as an example, was capable to enhance tolerance to salt, drought, and low temperature in transgenic Arabidopsis and rice (Zong et al., 2016). Transcription factors including NAC, bZIP, MYB, and MYC also play important roles in cold and drought stresses. Overexpression of MlNAC5 (Miscanthus lutarioriparius), SlNAC1 (Suaeda liaotungensis K.), and VaCBF4 (Vitis amurensis) enhanced drought and cold stress tolerance of Arabidopsis, respectively Yang et al., 2015;Zong et al., 2016). MAPK related genes that are involved in signal transduction are significantly induced by both cold and drought (Pitzschke et al., 2009). A MAPK gene (GRMZM2G174170_T01) from maize, for example, was activated under cold and drought stresses and was suggested to contribute to stress tolerance in maize (Shan et al., 2013). The expression levels of most ZmVQ genes (41 out of 61 members) encoding VQ motif-containing proteins were changed by the drought stress, and half of ZmVQ genes were co-expressed with ZmWRKY genes in maize (Song et al., 2015).
Phytohormones such as abscisic acid (ABA), ethylene, cytokinin (CK), auxin (IAA), gibberellin (GA), and jasmonate (JA) play important roles in regulating plant growth and development and also in the responses to various biotic and abiotic stresses (Peleg and Blumwald, 2011). ABA synthesis is one of the fastest responses of plants to abiotic stress, triggering ABAinducible gene expression and causing stomatal closure, thereby reducing water loss via transpiration and eventually restricting cellular growth. 9-cis-epoxycarotenoid dioxygenase (NCED) is a key enzyme in the ABA biosynthetic pathway in several plants which cleaves carotenoids to form the phytohormone ABA (Riahi et al., 2013). Cold or drought stress can induce ABA biosynthesis and exogenous ABA can improve cold or drought tolerance in tomato, Arabidopsis and tobacco (Thompson et al., 2000;Wan and Li, 2006;Zhang et al., 2008). In maize, thousands of genes are thought to be involved in abiotic stress. In particular, the gibberellin (GA) metabolic genes. Thus, understanding their expression changes could help to detect molecular mechanisms of the GA pathway under stress conditions (Colebrook et al., 2014). The response of GA metabolic genes to abiotic stress has been investigated and these genes participate in the CBF1-mediated stress-response pathway (Niu et al., 2014) and the GA20ox gene that is responsible for GA biosynthesis is down-regulated by cold treatment in Zea maize (Shan et al., 2013).
Zea mays ssp. mexicana L. is a member of teosintes and a close wild relative of cultivated maize (Almeida et al., 2011). Although teosintes have not yet been widely used in maize breeding, the high genetic diversity shows that Zea mays ssp. mexicana L. is a genetic reservoir for the improvement of agronomic characteristics of cultivated maize and teosintes (Wang et al., 2008). Since its introduction to China from Japan in 1979, the planting areas of teosinte have increased extensively in southern China, including Guangdong, Fujian, Guangxi, and Sichuan Provinces (Song et al., 2005). As an important forage or silage source, teosinte has a larger plant biomass, higher tiller number and better resistances to various stresses than cultivated maize (Niazi et al., 2015). Similar to maize, most lines of Zea mays ssp. mexicana L. originate from high altitude in northern and central Mexico and can adapt to acid soil with the best growing temperatures ranging from 25 • to 35 • C (Fukunaga et al., 2005). It is considered to be a very cold-sensitive crop, especially during the germination and early autotrophic growth stages, despite it having originated at altitudes of about 1,000-2,000 m (Hincha and Zuther, 2014).
RNA-seq (RNA sequencing) has become the most convenient and cost effective tool for understanding gene structure and quantitative transcriptome profiling. It does not depend on the existence of known genomic sequences, and shows a great ability for detecting differentially expressed genes with a broader dynamic range of expression levels (Wang et al., 2009). To discover genetic bases and molecular mechanisms in response to cold or drought in Zea mays ssp. mexicana L., the RNA-seq technique was employed in this study and the transcriptome of Zea mays ssp. mexicana L. was presented for the first time. In total, 414,232,462 high quality clean reads were obtained by RNAseq and were used for de novo assembly and annotation of genes from Zea mays ssp. mexicana L. Differentially expressed genes (DEGs) under cold or drought were also identified. These data will be valuable for the exploration of genetic and molecular mechanisms in response to stresses in Zea mays ssp. mexicana L. and provide gene resources for breeding programs.

Plant Materials and Growth Conditions
Plump seeds of Zea mays ssp. mexicana L. variety "8493" were used in this study. Seeds were washed 3 times with distilled water, then soaked in 75% ethanol and 2% sodium hypochlorite for 10 min and 3 min, respectively, and washed 3 times with distilled water before planting in plastic boxes (54 × 28 × 7 cm) containing soil substrates (Jiffy, Netherlands, http://www. jiffygroup.com/jiffy_product_category/substrates/). The boxes were transferred into a climate control box (RXZ 500-C, JIANGNAN Instrument), allowing the seeds to germinate and grow under a photoperiod of 10-h light/14-h dark and a humidity of 60% at 25 • C. Thirteen-day-old seedlings were assigned into 3 groups for different treatments. For cold treatment (cold), seedlings were treated for 12 h at 4 • C under 10-h light/14h dark while for drought treatment (drought), seedlings were treated with Hoagland solution containing PEG2000 (20%) for 3 h under the same light condition. Seedlings of the control group (CK) grew under the same photoperiod and were assigned neither of the above treatments. After treatments, fresh tissues (roots, stems and leaves) were sampled and pooled, followed by quick-freezing with liquid nitrogen and storing at −80 • C for further analysis. Two replicates for each treatment were applied and plant materials from at least three seedlings were pooled for each sample.

RNA Extraction
Total RNA was extracted using TRIzol R reagent (Invitrogen, Carlsbad, USA) according to the manufacturer's instructions. RNA was treated with RNase-free DNase I (Takara, Japan) to remove any possible DNA. The integrity was checked by gel electrophoresis and the Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA). Concentration of total RNA was determined using NanoDrop 8000 spectrophotometer (NanoDrop, Wilmington, DE). Total RNA with RIN values ≥7.3 and 28S:18S ratios ≥1.2 were accepted for the following RNA-seq analysis.

RNA-seq Analysis
The Illumina HiSeq2500 platform was applied for RNA-seq analysis. To obtain a comprehensive overview of the Zea mays spp. mexicana L. transcriptome and the transcript profiles under cold and drought, six libraries were constructed and pairedend sequencing was carried out according to the manufacturer's instructions (Illumina, San Diego, CA). In brief, poly(A)-tailed mRNA was enriched using biotin-Oligo (dT) magnetic beads and fragmented into short fragments of 200-700 bp, followed by converting into double-stranded cDNA. The cDNA fragments were then purified with a QiaQuick PCR extraction kit and adapters were added to both ends of the short fragments. cDNA fragment pools were loaded to Illumina HiSeq2500 platform for sequencing. The generated raw data were processed for further analysis. Library construction and RNA-seq were carried out by Novogene Bioinformatics Technology Co. Ltd, China.

De novo Assembly and Sequence Clustering
The quality of raw data was controlled by using the FastQC tool (http://www.bioinformatics.babraham.ac.uk). Reads containing substantial unknown nucleotides (the proportion of N > 5%) and low-quality reads (>50% of the bases with a quality score ≤5) were discarded, and the first 10 bases of processed data were removed. The generated clean data from each sample were merged to perform de novo assembly using Trinity with default parameters and an optimized k-mer length of 25. A set of transcripts were generated and served as a reference transcriptome. Based on the shared sequence, transcripts from a gene were clustered into a subclass and the longest transcript was regarded as the unigene. The above processes were performed by Trinity software with default parameters (http://trinityrnaseq.sf.net, updated on 2014/07/17; Haas et al., 2013). To identify possible plant transcription factors (TFs), transcriptional regulators (TRs) and protein kinases (PKs), Plant Transcription factor and Protein Kinase Identifier and Classifier (iTAK) V1.5 was employed to analyse unigenes with the best match result.
Raw Illumina sequences and assembled sequences are available in the Gene Expression Omnibus (GEO) database of the National Center for Biotechnology Information (NCBI) (accession number: GSE76939).

Identification of Differentially Expressed Unigenes (DEGs)
The paired-reads from each sample were mapped to an assembled reference transcriptome by Bowtie software v1.1.1 with default parameters, and the number of mapped reads were calculated by RSEM (v 2.15) (Langmead et al., 2009;Li and Dewey, 2011). FPKM values were assigned to each unigene, representing normalized expression level by eliminating the influences from library construction and the length of genes (Trapnell et al., 2010). Unigenes with FPKM ≥ 1 were used for further analysis. Fold changes for each unigene under cold or drought condition were determined by comparing the FPKM value to that in control sample, and those unigenes with more than 2-fold change and adjusted P < 0.005 were identified as differentially expressed genes (DEGs) by DESeq 2 (V 1.2.10) package of the R statistical program (Love et al., 2014).

Annotation and Function Classification of Transcripts and Unigenes
Annotations were assigned to each unigene based on the top hit in BLASTX search against the protein databases, with the non-redundant (Nr) protein database at GeneBank (http://www.ncbi.nlm.nih.gov) as the highest priority, followed by Swiss-Prot (http://www.expasy.ch/sprot), KEGG (http://www.genome.jp/kegg) and eukaryotic KOG (http://www.ncbi.nlm.nih.gov/KOG) in that priority order. The significant thresholds of E-value were set at ≤10 −6 . CDS (5 ′ ->3 ′ ) predictions of unigenes were confirmed with the priority result of BLAST searched protein.
The accession number and the GI code of the top hit at the BLASTX search against the Nr database were retrieved and ID mapping was performed to obtain GO annotations for the queried unigenes. To further annotate the unigenes, the Blast2GO (v2.5.0) program was employed to get GO annotations based on molecular function (MF), biological process (BP), and cellular component (CC) features (Conesa et al., 2005). GO enrichment analysis were carried out by AgriGO software with FDR < 0.05.

Quantitative Real-Time PCR Analysis
To further verify the expression profiles of genes in our Illumina RNA-seq data, total RNAs were extracted from three independent groups, which were different from those used for RNA-Seq. First-strand cDNA synthesis and qRT-PCR were carried out with the PrimeScript TM RT reagent Kit with gDNA Eraser (TaKaRa) and TransStart R Tip Green qPCR SuperMix (Transgen), respectively. qRT-PCR was carried out in the BIO-RAD CFX96 sequence detection system according to the manufacturer's instructions. A total of 12 genes were randomly selected and gene-specific primers were designed with online Primer-blast of NCBI. The unigene annotated maize actin-2 was used as the endogenous control (GeneBank accession number NP_001146931.1). Each PCR reaction (20 µL) contained 10 µL 2 × Green qPCR SuperMix, 0.2 µM of each primer and Frontiers in Plant Science | www.frontiersin.org appropriately diluted cDNA. The thermal cycling conditions were 94 • C for 30 s, followed by 40 cycles of 5 s at 94 • C and 30 s at 60 • C then 72 • C for 30 s. At the second dissociation stage, 95 • C for 10 s, 65 • -95 • C with increment of 0.5 • C for 0.05 s were used. All reactions were performed in triplicate, including the nontemplate controls. The relative expression level was calculated with the 2 − CT method (Livak and Schmittgen, 2001).

De novo Assembly and Functional
Annotation of the Zea mays ssp. mexicana L. Transcriptome To comprehensively understand the transcript profiles of Zea mays spp. mexicana L. under cold and drought, a total of six samples, including two biological repeats for each condition, were collected and two replicated runs of Illumine sequencing were performed for each sample. Raw reads of 153,219,582, 156,576,828 and 150,941,988 were generated from the control, cold and drought group, respectively. After quality check, the adapters, low-quality sequences and ambiguous reads were removed. A total of 136,233,336 (control), 141,732,126 (cold), and 136,267,000 (drought) clean paired-end reads, corresponding to 17.04 Gb (control), 17.72 Gb (cold), and 17.02 Gb (drought) clean bases, respectively, were used for future assembly. De novo assembly was carried out using the Trinity software and 251,145 transcripts and 184,280 unigenes with N50 values of 1,269 and 923 bp were obtained, respectively ( Table 1). Approximately 47% of assembled transcripts and 35% of unigenes had a length of more than 500 bp (Supplemental File 1).
Homologs were also identified by searching against the protein database of the most closed relative Zea mays L. and the model species, Arabidopsis thaliana. A total of 60,583 and A total of 2,157 and 512 unigenes were identified as transcription factors (TF) or transcription regulators (TR) based on the Blast search against the Nr protein database at NCBI, covering 62 TF families and 25 TR families, respectively ( Table 2). In addition, predicted protein kinases (PK) that play crucial roles during signal transduction and other biological processes were also obtained, including 1,598 unigenes covering 67 PK families ( Table 3).

Expression Quantification of Unigenes and Identification of Differentially Expressed Genes (DEGs)
The expression levels of unigenes were evaluated using FPKM values. A total of 43,432 (23.57%), 42,849 (23.25%), 55,928 (30.35%), 48,945 (26.56%), 47,327 (25.68%), and 43,879 (23.81%) unigenes had FPKM ≥ 1 from the sample of control-1, control-2, cold-1, cold-2, drought-1, and drought-2, respectively. However, the percent of unigenes with FPKM ≥ 100 from each sample was less than 1%, ranging from 0.78 to 0.84% (Supplemental File 1). Compared with the control, the cold treatment at 4 • C resulted in a larger number of genes with high FPKM values than did the drought treatment with PEG, indicating that the responses induced by cold was more dramatic than that by drought stress in Zea mays spp. mexicana L. at least with the relative levels of stress used in this trial.
To investigate differential responses of Zea mays spp. mexicana L. to cold or drought stress, data from replicated treatment groups were merged to perform pair-wise comparisons, and unigenes that matched the criteria of |log 2 Change|> = 2 and adjusted p < 0.005 were identified as differentially expressed genes (DEGs). Totally, 5,338 DEGs were obtained, with 3,049 and 614 genes differentially responding to cold and drought stresses, respectively. Among those coldresponsive DEGs, 2,232 were up-regulated and 817 were down-regulated ( Figure 2C, Supplemental File 4). Notably, annotations for 41.1% of the up-regulated and 25.7% of the down-regulated DEGs were not retrieved from the public maize database (Supplemental File 4), suggesting some special FIGURE 1 | Functional classification of assembled unigenes. Functional classification of the assembled unigenes based on Gene Ontology (GO) categorization. GO annotations were assigned to unigenes based on the best BLAST hits, and a total of 98,947 GO terms were retrieved for 28,382 unigenes which were grouped into 49 categories. The left Y-axis represents the percentages of unigenes in each main category. The right Y-axis indicates the numbers of unigenes in each GO category.
regulation paths in Zea mays spp. mexicana L. For those DEGs which responded to drought, however, the number of DEGs was less, with only 532 and 82 genes being activated and inhibited by drought, respectively, of which 14.47 (77) and 4.2% (17) were unannotated (Supplemental File 4). The expression level of the whole set of DEGs in all six samples were also evaluated by plotting the Log 10 -transformed FPKM values in descending order (Figure 2A), and the distributions of DEGs' expression suggested activation was present in both cold-and droughttreated samples, with more dramatic transcriptional changes in the cold-treated samples. Taken together, these results indicated that Zea mays spp. mexicana L. had distinct gene networks in response to cold and drought.
Venn analysis indicated that 194 unigenes were regulated by both factors and 2,855 and 420 were specifically induced by either cold or drought ( Figure 2B, Supplemental File 4).Of the 194 shared DEGs, 179 genes had similar response patterns to both cold and drought, including 152 up-regulated genes and 27 down-regulated genes. Two unigenes (c52482_g3 and c52482_g6) similar to protein phosphatase 2C3 and an unannotated unigene (c50928_g2) gave the greatest upregulation, while the c29390_g1 (thioredoxin H4), c63024_g4 (hypothetical protein ZEAMMB73_373831) and an unannotated unigene (c55190_g2) comprised were the three most downregulated. Fifteen genes showed the opposite response patterns under cold and drought stresses with 13 genes being downregulated by cold but up-regulated by drought.

GO Enrichment of DEGs
GO enrichment analysis was performed to investigate overrepresented biological processes or molecular functions involved in responses to cold or drought stress. A total of 45 and 89 GO terms were identified based on FDR < 0.05 as overrepresentations under cold and drought, respectively (Tables 4, 5 and Supplemental File 7).
Apparently, terms related to the response to abiotic stimuli, especially those to temperature and light, were significantly enriched in cold-treated Zea mays spp. mexicana L. seedlings. Genes for heat shock factor proteins (HSF) such as hsf1 (c59237_g1), hsp3 (c60986_g3), hsf4 (c46749_g1), hsf7 (c54684_g2), hsf26 (c46143_g1), hsp81 (c56395_g1), and hsp101 (c58173_g1) were down-regulated under cold stress. With respect to light responsive genes those for phosphatase 2C50 (c57868_g1) and chloroplastic-like stress response protein (c62415_g1) were up-regulated under cold stress. Additionally, terms related to the regulation of gene expression, mostly at the level of transcription, were over-represented. Several genes for encoding transcription factors, AP2/EREBP (c48755_g1 and c45827_g1), NAC (c60421_g1 and c47289_g2), bZIP (c61581_g3 and c59299_g1), WRKY (c59902_g1, c38115_g1 and c63357_g2 as examples) and ethylene-responsive transcription factor (c56755_g1 and c63899_g2), were up-regulated under cold stress. Together, these observations emphasized the roles of transcriptional regulation of cold-related responses in this species.
For the drought associated responses, however, there were a few differences. Under drought stress, DEGs were enriched for a total of 89 processes. Responses to stimuli, especially those to water deprivation, and regulation at the level of transcription were the most dramatic. Genes for the response to water deprivation were up-regulated under drought. These genes include dehydration-responsive element-binding protein 1D (c49699_g2 and c55715_g7), galactinol synthase gene (c60364_g1) and aquaporin PIP1-5 (c56750_g3 and c35748_g1). Genes involved in abscisic acid responses, such as EID1-like F-box protein (c57232_g1, c57232_g2 and c57232_g3), NAC domain transcription factor (c52045_g1), protein phosphatase 2C (c57868_g1, c50507_g1 and c60116_g4), and stachyose synthase (c55229_g2), were up-regulated under drought.
To further identify the differences between responses induced by cold and drought in Zea mays spp. mexicana L., a direct comparison was carried out between cold and drought-treated samples. A total of 2,514 assembled unigenes were identified by |log 2 change|> = 2 and adjusted p < 0.005 ( Figure 2C). GO enrichment analysis for these DEGs indicated that the differences between cold and drought induced responses were in relation to abiotic stresses of temperature and water, transcriptional regulation and protein phosphorylation processes (Supplemental File 5). The FDR value in relation to the abscisic acid related signaling pathway was 0.011 which was very close to the criteria of 0.01 (Table 5), suggesting that this pathway was also a major difference between the two datasets. Together, the enrichment patterns identified by directly comparing responses to cold and drought were in agreement with the above results obtained from comparing each in turn with the untreated control.

qRT-PCR Verification
In order to validate the RNA-seq results, 12 DEGs with different expression patterns were selected for qRT-PCR analysis using gene-specific primers (Supplemental File 6). Among these DEGs, nine expressed under both cold and drought stresses, one regulated by cold only and two regulated by drought only were included. Eight of 10 cold-responsive DEGs and night of 11 drought-responsive DEGs were confirmed by qPCR (Figure 3).

DISCUSSION
RNA-Seq technology has been successfully used at the level of transcriptomes or whole genomes analysis in many plants and crops (Lu et al., 2010;Xia et al., 2011;Xu et al., 2012Xu et al., , 2014Zhao et al., 2012;Garcia-Seco et al., 2015). Zea mays ssp. mexicana L. is a subspecies of Zea mays, and it has a genome size of 2.58 Gb (2n = 20) which is more than 20 times greater than the Arabidopsis thaliana (0.115 Gb) genome, 6 times greater than the rice (0.43 Gb) genome and equivalent to the Zea mays ssp. mays genome (2.73 Gb) (Hufford et al., 2012). The variety 8493, which serves as a representative of Zea mays spp. mexicana L., is reported to grow well under 25 • C and can adapt to semi-arid and semi-humid environments (Huang et al., 1992). In this study, RNA-seq analyses were performed for de novo assembly of the transcriptome of Zea mays spp. mexicana L. and identification of the differentially expressed genes (DEGs) in response to cold or drought stress.
An Overview of the Zea mays ssp. mexicana L. Transcriptome By pooling the clean reads from all samples, we obtained a limited transcriptome which provided a profile the transcript dynamics under cold and drought stresses. The transcriptome identified 184,280 unigenes, and 251,145 transcripts which Blast analysis divided in protein coding and other genes. When growing under cold (4 • C) condition, the change of gene expression was dramatic at the global level of transcriptome (Figures 2A,C), indicating that cold induces an extensive activation of transcription. Moreover, GO analysis (Table 4) suggested there was an activation of gene networks involving in cold/heat-induced responses. The drought stress used here, however, induced fewer transcriptional changes, (only 15% as many) than cold ( Figure 2C). These observations conformed the geographic distribution and the ideal growth conditions of Zea mays ssp. mexicana L. (Huang et al., 1992;Hufford et al., 2012). Similar observations have also been found in the elite maize inbred line Zheng 58 of Zea mays in which it was reported that the number of differentially expressed genes under drought stress was approximate 42% of that under cold (Shan et al., 2013), suggesting that the more sensitive response to cold rather than drought would be a conserved mechanism in Zea. It is well known that gene expression in response to cold stress is somewhat different from that to drought stress (Shinozaki et al., 2003). In this study, only 194 DEGs are shared in stressed Zea mays ssp. mexicana L., (6.4 and 31.1% of coldand drought-related DEGs, respectively). Among these, nearly 90% are regulated in a similar manner by both stresses, while only 15 DEGs are oppositely regulated (Supplemental File 4), indicating that there is a shared network to regulate the cold and drought induced responses in Zea mays ssp. mexicana L. On the other hand, specific regulations in response to cold or drought stress are also clearly visible in Zea mays ssp. mexicana L. The DEGs in response to cold are mostly enriched in processes of "Response to abiotic stimulus (GO:0009628), " especially the "Response to heat (GO:0009408), " and "Response to high light intensity (GO:0009644)" (Table 4), while the processes of "Response to water (GO:0009415), " "Response to water deprivation (GO:0009414), " and "protein amino acid dephosphorylation (GO:0004722)" are dominant in the collection of drought-induced DEGs ( Table 5).

Transcriptional Regulatory Networks and Signaling Pathways involved in Cold and
Drought Induced Responses in Zea mays ssp. mexicana L.
Transcription factors play crucial roles in the regulation of target gene expression by specifically binding to cis-acting elements in the promoter regions (Agarwal et al., 2006; Van Buskirk and Thomashow, 2006). In this study, WRKY, CBL, MYB, NAC, bHLH, bZIP, and AP2/EREBP families that are known to be involved in stress tolerance in plants were identified. Members of WRKY (c56800_g2 and c56800_g1), CBL (c51372_g1 and c49176_g1), MYB (c54045_g1), bHLH (c48928_g1), bZIP (c59299_g1), AP2/EREBP (c51992_g7 and c48755_g1) family are up-regulated by both stresses of cold and drought ( Table 2), suggesting shared upstream pathways for signal transduction and regulation under these stimuli. From those down-regulated TF proteins, only one NAC (c57823_g1) protein is consistently down-regulated by both stresses, which is similar to the findings in Arabidopsis and rice (Lu et al., 2007;Takasaki et al., 2010) and suggests a negative regulation mechanism in Zea mays ssp. mexicana L. It is notable that some previously uncharacterized TFs are significantly up-or downregulated in response to cold or drought only. For example, a member of TUB family (c61146_g4), was found in coldinduced DEGs set but not in the drought-induced set. This is similar to what has been observed in V. amurensis (Xu et al., 2014). In addition, some members in the ULT (c53036_g1), BIM2 (c64573_g2), and DBP (c51832_g3) families are found to be up-regulated by cold only ( Table 2). The AP2/ERF transcription factor family comprises four major subfamilies, namely AP2, ERF, RAV, and dehydration-responsive element binding proteins (DREB/CBF) (Shinozaki et al., 2003). The roles of DREBs/CBFs in stresses induced response have been well established in numerous plants including A. thaliana , rice (Challam et al., 2015), maize (Liu et al., 2013), and soybean (Kidokoro et al., 2015). The DREB1s/CBFs can activate the downstream cold-responsive genes via specific binding to the DRE/CRT cis-acting element in their promoters, and expressions of these genes are positively correlated with cold tolerance (Mizoi et al., 2012). In Zea mays ssp. mexicana L., three DREB genes (c51992_g7, c49699_g1 and c49699_g2) were found to be significantly induced by both cold and drought stresses, suggesting that these proteins might functionally act in Zea mays ssp. mexicana L. during cold and drought stresses.
Signaling components such as kinases participate in numerous processes, including cell division, developmental programs, hormonal responses, and drought, salinity and reactive oxygen species (ROS) signaling responses . In particular, the MAPK cascade, serving as a signaling joint for various abiotic stress, positively regulates cold acclimation in plants (Chinnusamy et al., 2010). In Arabidopsis, the AtMEKK1/ANP1 (MAPKKK)-AtMKK2 (MAPKK)-AtMPK4/6 FIGURE 3 | Verification for expressions of 12 selected DEGs. Relative expressions of the 12 genes in cold-or drought-treated Zea mays spp. mexicana L. seedlings were tested by qRT-PCR. Transcription levels were normalized to that of actin and compared to the untreated control. Values of fold change are showed in mean ± SEM. Three independent replicates were performed for each gene. Fold changes of these genes detected in RNA-Seq are also showed for reference. Relative transcription levels were calculated by log2 − CT method. The names of the closest putative orthologs are indicated for each gene, and the annotations are as following: c52482_g3, uncharacterized protein LOC_100192073; c61291_g1, calcium-binding protein CML16; c57152_g1, transcription factor FER-LIKE IRON DEFICIENCY-INDUCED; c54045_g1, MYB DNA-binding domain superfamily protein isoform X1; c59856_g7, heat shock factor protein HSF30; c51372_g1, CBL-interacting serine/threonine-protein kinase 11 isoform X1; c57309_g1, E3 ubiquitin-protein ligase EL5-like; c54089_g2, heat shock protein 82; c56552_g1, 17.0 kDa class II heat shock protein; c58118_g1, photosystem I reaction center subunit N, chloroplastic-like; c58084_g2, protein phosphatase 2C37; c60754_g2, homeobox-leucine zipper protein ATHB-6.
(MAPK) mitogen-activated protein kinase cascade is necessary for cold acclimation, and AtMKK2 (MAPKK) was identified as a key signal transducer. In addition, overexpression of AtMKK2 induces the cold responsive genes CBF2 and CBF3 and improves cold tolerance in Arabidopsis (Teige et al., 2004). In this study, one unigene coding for a MAP kinase (c14147_g1) and another one for a MAPKK protein (c62016_g1) were up-regulated by cold stress only, suggesting that the MAPK cascade may be involved in response to cold in Zea mays ssp. mexicana L.

ICE1-CBF pathway Genes Involved in Cold and Drought Stresses
Cold induced responses come down to different metabolic pathways, gene regulations and cell compartments. The ICE1-CBF-COR transcriptional cascade is a highly conserved module which is activated in response to cold stress (Zhu et al., 2007;Lissarre et al., 2010), and has been demonstrated to have essential roles in response to cold, drought, and other abiotic   (Zhou et al., 2011;Peng et al., 2015). ICE1 is a MYCtype transcription factor that binds to the MYC recognition cis-elements (CANNTG) in the promoter of CBF3/DREB1A protein to activate the expression of CBF3/DREB1A, which consequently regulates 40% of COR genes and 46% of coldregulated transcription factor genes (Lee et al., 2005). In Arabidopsis, three CBF proteins, namely CBF1, 2, and 3, are induced by low temperature but not by abscisic acid (ABA) or dehydration (Medina et al., 1999). CBF1 and CBF3 can respond to cold at the same time, but the response of the CBF2 lags behind. CBF2 negatively regulates the expression of CBF1 and CBF3 during cold acclimation, while CBF1 and CBF3 are not involved in regulation of other CBF genes. CBF1/CBF3 have lower cold acclimation ability and play a subsidiary function in the CBF/DREB1 pathway and cold acclimation (Gilmour et al., 2004;Novillo et al., 2007;Kim et al., 2015). Differences in the functions of CBFs in mangrove Avicennia marina, European bilberry (Vaccinium myrtillus) and Brassica rapa have been observed. For example, CBF4 is regulated by drought and ABA but not by cold (Haake et al., 2002;Lee et al., 2012;Oakenfull et al., 2013;Peng et al., 2013). In this study, two putative ICE1 homologs (c60792_g1, c13919_g1) were identified in this subspecies. However, expressions of these genes were not affected by cold. A possible explanation is that the ICE1 homolog may be induced rapidly and transiently, and is attenuated in a short time, making it undetectable. In addition, CBF6 has been found to play a role in freezing tolerance in Festuca pratensis (Ergon et al., 2016). In our study, a CBF6 (c57018_g4) homologous gene was found to be up-regulated in response to cold and drought stresses, suggesting that this CBF6 homologous gene might play roles in stress resistance in Zea mays ssp. mexicana L.
Trehalose May be a Common Regulator of the Response to Cold and Drought Stresses in Zea mays ssp. mexicana L.
Trehalose is widely found in various organisms and is a nonreducing disaccharide formed by linking two glucose units with an a,a-1,1-glycosidic linkage (Wang et al., 2005). The synthetization of trehalose involves two-steps. The first step is the formation of trehalose-6-phosphate (T6P) from UDP-glucose and glucose-6-phosphate, which is catalyzed by trehalose-6phosphate synthase (TPS). The second step is the generation of trehalose through dephosphorylation which is catalyzed by trehalose-6-phosphate phosphatase (TPP) (Zang et al., 2011). The genes coding for the TPS and TPP enzymes are thus important for the regulation of trehalose production (Ilhan et al., 2015). Transgenic plants harboring TPS or TPP genes have enhanced tolerance to abiotic stresses by stabilizing dehydrated enzymes, as well as by protecting protein and cellular membranes structures (Nuccio et al., 2015). In our study, a key trehalose-6phosphate phosphatase (TPP) gene (c45624_g2) was also found to be up-regulated under both cold and drought stresses in Zea mays ssp. mexicana L. This is in agreement with the positive roles of trehalose in tolerating low temperature and drought stresses.
Phytohormone Dependent Pathways in Zea mays ssp. mexicana L.
Plant hormones play vital roles in adaption to stressful environments by regulating growth, development, nutrient allocation and source/sink transitions. Previous studies showed that different types of plant hormones, including abscisic acid (ABA), cytokinins (CKs), jasmonic acid (JA), gibberellins (GAs), brassinosteroids (BRs), ethylene (ETH), and auxin, have functions in regulating freezing tolerance by CBF-dependent or -independent pathways (Shi et al., 2015). ABA is the key phytohormone used in plants for regulation of the responses to abiotic stresses. As an example, the secondary hyperosmotic stress derived from drought or salt stress causes the accumulation of ABA which in turn rapidly elicits downstream responses in plants (Zhu, 2016). ABA accumulation can cause stomatal closure, thereby reducing water loss and eventually restrict cellular growth. At the molecular level, ABA regulates gene transcription, protein synthesis, signaling pathways, ion transport (and the transport of other organic molecules) and the production of important protectants against dehydration and photoinhibition (Riahi et al., 2013). Cold stress is also reported to induce ABA biosynthesis, and the exogenous application of ABA improves the cold tolerance of plants (De Zelicourt et al., 2016). ABA biosynthesis is one of the important steps used by plants in responding to stresses via the ABA-dependent pathway. The 9-cis epoxycarotenoid dioxygenase 1 (NCED1) is regarded as a key rate-limiting enzyme in ABA biosynthesis. It was originally identified from maize viviparous 14 mutants, and is responsible for cleavage of the ABA precursor C40-cis-epoxycarotenoids, to either 9-cis-neoxanthin, 9-cis-violaxanthin or both to produce xanthoxin, the direct C15 precursor of ABA (Raghavendra et al., 2010;Boursiac et al., 2013). In Zea mays ssp. mexicana L., a homolog of NCED1 (c63716_g5) gene was identified. It was upregulated by drought (Supplemental File 4), which is consistent with results reported previously in Arabidopsis (Behnam et al., 2013) and Phaseolus vulgaris L. (Qin and Jan, 1999), suggesting conserved regulation of this drought induced response in Zea mays ssp. mexicana L. It is notable that this NCED1 gene was also activated by cold stress (Supplemental File 4), which is similar to recent observations for OsNCED1 and OsNCED3 (Maruyama et al., 2014), suggesting that the ABA-dependent pathway plays roles in response to cold and drought in Zea mays ssp. mexicana L. PP2C proteins are also reported to be involved in ABAdependent regulations of stress induced responses, such as the 76 phosphatase 2C (PP2C-type) proteins which have been identified from Arabidopsis. Most of these (PP2C-type) proteins fall into 10 groups (A-J) (Schweighofer et al., 2004). ABA signaling can be mediated by type 2C protein phosphatases (PP2Cs) such as HAB1 and ABI2, which inhibit stress-activated SnRK2 kinases (Nakashima et al., 2009). When directly comparing the expression levels under cold with those under drought conditions in Zea mays ssp. mexicana L., four candidate PP2C genes, including c52482_g3, c52482_g6, c62830_g2, and c57868_g1, showed up-regulated patterns. A putative protein kinase SnRK1 (c57228_g1), however, was down-regulated only under cold stress. These results taken together suggest that PP2C protein is involved in ABA-dependent signaling and affects the regulation of cold induced responses in Zea mays ssp. mexicana L.
Previous studies on Zea mays ssp. maize L. transcript profile under cold, salt, and drought stresses indicate that the abundance of bioactive GAs was reduced and that the negative regulatory factor DELLAs were accumulated (Shan et al., 2013), suggesting a negative effects of gibberellin in stress induced responses. This is in disagreement with what we observed in cold-stressed Zea mays ssp. mexicana L. Genes for the GAs receptor GID1L2 precursor (c47232_g1), the GA 2-beta-dioxygenase (c57017_g1, c33506_g1 and c57017_g2) and the GA20 oxidase (c63822_g4, c49254_g1 and c56760_g2) that are involved in GAs synthesis or signaling were up-regulated by cold, suggesting a positive role of gibberellin in this subspecies in response to cold stress. Thus, GAs could differentially regulate the responses to cold in these two close Zea mays subspecies. Future comparative study will provide better understanding of the molecular basis of the GA associated response to cold in Zea mays.
In summary, this study provided an overview of the transcriptome of Zea mays ssp. mexicana L., and the changes of transcript abundance in response to cold or drought stress. Less than a half of the reported Zea mays proteome have homologs in Zea mays ssp. mexicana L., indicating that specific regulation networks or mechanisms exist in these close relatives. DEGs responsive to cold or drought conditions serve as candidates for further study on tolerance related gene networks in Zea mays ssp. mexicana L. The effects of drought and cold reported here have arisen from a limited range of potential types and severities of stress. Thus, it should be noted that a greater range of treatments for (e.g., timing, severity, frequency, source) need to be examined in future studies to provide more clues for understanding the adaptation and tolerance mechanisms in this species.

AUTHOR CONTRIBUTIONS
XL conducted the bulk of experiments, data analysis, and wrote the paper draft. XZ and YC participated in the treatments of experimental materials. MZ and DM contributed to data interpretation and writing. CY and SL were responsible for experimental design, data interpretation and writing.

FUNDING
This research was supported by the Guangdong Science and Technology Department (2015A020209155) and the National Natural Science Foundation of China (U1201212).

ACKNOWLEDGMENTS
We thank the technical support of Qi Zhao and Yue-yuan Zheng from Sun Yat-Sen University, Hao-ran Shu from Guangzhou Magigen Biotechnology Corporation.