Transcriptomic Analysis of Female Panicles Reveals Gene Expression Responses to Drought Stress in Maize ( Zea mays L.)

: Female panicles (FPs) play an important role in the formation of yields in maize. From 40 days after sowing to the tasseling stage for summer maize, FPs are developing and sensitive to drought. However, it remains unclear how FPs respond to drought stress during FP development. In this study, FP di ﬀ erentiation was observed at 20 and 30 days after drought (DAD) and agronomic trait changes of maize ears were determined across three treatments, including well-watered (CK), light drought (LD), and moderate drought (MD) treatments at 20, 25, and 30 DAD. RNA-sequencing was then used to identify di ﬀ erentially expressed genes (DEGs) in FPs at 30 DAD. Spikelets and ﬂorets were suppressed in LD and MD treatments, suggesting that drought slows FP development and thus decreases yields. Transcriptome analysis indicated that 40, 876, and 887 DEGs were detected in LD / CK, MD / CK, and MD / LD comparisons. KEGG pathway analysis showed that ‘biosynthesis of other secondary metabolites’ and ‘carbohydrate metabolism’ were involved in the LD response, whereas ‘starch and sucrose metabolism’ and ‘plant hormone signal transduction’ played important roles in the MD response. In addition, a series of molecular cues related to development and growth were screened for their drought stress responses.


Introduction
Under the influence of global warming, changes in climatic conditions are creating unusual weather phenomena worldwide, often imposing drought stress on crops [1,2]. From the agricultural perspective, drought often results in decreased crop productivity and growth [3][4][5], especially for cereal crops. Maize (Zea mays L.) is one of the most important cereal crops and has the most extensive planting area globally [6,7]. One of the most important factors limiting maize growth and development around the world is a lack of water [8][9][10][11]. Accordingly, improving tolerance of maize to drought stress is essential for achieving high and stable yields in cereal crops.
As a multidimensional stress, water limitation triggers a wide variety of plant responses; these range from responses at the physiological and biochemical levels to the molecular level [12][13][14][15][16]. When external drought stimuli are perceived and captured by sensors on cell membranes, the signals are transmitted through multiple signal transduction pathways. Then, plant can regulate the expression of drought-responsive genes to protect themselves from the harmful effects of external stimuli [17,18]. The expressed products of drought-responsive genes are mainly proteins involved in signaling cascades and transcriptional regulation (such as protein kinase, protein phosphatase, and transcription factors) and functional proteins [19]. With the rapid development of high-throughput sequencing technologies, transcriptome analyses have been conducted to identify stress-mediated differences at the level of gene expression. Previous research has shown that many significantly differentially regulated genes that were associated with drought tolerance are induced in different organs of maize plants [20][21][22][23][24][25][26]. For example, 249 and 3000 differentially expressed genes (DEGs) were involved in root tissues after 6 h of light and severe drought stress, respectively [23]. In leaves, a total of 619 DEGs and 126 transcripts had their expression levels altered by drought stress at flowering time [24]. In tassels, 1902 DEGs were found after 5-7 days of drought stress [25]. In young ears, a total of 1825 DEGs were identified on the 5th day of drought stress at the V9 stage [26].
The panicle stage (from jointing to flowering) is the key stage for panicle differentiation and development in maize, the number of rows per ear and the number of grains per row are dependent on spikelet and floret differentiation at this time [27]. Therefore, the growth and development of female panicles (FPs) play an important role in the formation of maize yields. Although great advances in understanding differentiation of FPs and how drought stress affects genes transcription in FPs have been achieved in the past few decades [26,[28][29][30][31], so far, progress in understanding the general molecular basis of FP development in response to long-term drought stress across the panicle stage has not been reported.
Accordingly, in this study, maize inbred line PH6WC (6WC) was used as drought-sensitive experimental material [32], soil water was controlled by means of drip irrigation for 30 days, and the gene expression dynamics of developing FPs at 30 days after drought (DAD) were investigated using transcriptomic analysis. The DEGs were identified and assigned to functional categories to reveal the various metabolic pathways in FPs that are involved in responses to long-term drought stress at different levels. Furthermore, differences in transcription factors between treatments were also analyzed. Overall, the exploration and function prediction of drought-response genes in maize FPs represent an efficient approach to improving the molecular breeding of drought-resistance maize cultivars.

Plant Material and Growth Conditions
Field experiments were conducted at field experiment stations (34 • 31 N, 115 • 35 E, 50.7 m above sea level) during the maize growing season (June-October, 2018) in Shangqiu (Henan, China). The maize inbred line 6WC was grown in 9 experimental plots (each plot was 2 m wide, 3.3 m long, 1.8 m deep) which were under a movable awning and filled with luvo-aquic soils, with a 20 cm sand filter layer at the bottom. Maize were planted into four rows, with 40 cm between rows in each plot. Two to three seeds were sown at each acupoint, with subsequent thinning to one seedling conducted at the trifoliate stage (V3). The final stand density was 8 plants m −2 . During the jointing and tasseling stages, topdressing fertilizer was applied, and weeds, insects, and diseases were controlled throughout the experiment. The top soil (0-40 cm layer) had a pH (water) of 7.3, mean mineral P content of 3.24 g/kg, and inorganic N at sowing of 3.60 g/kg. The average daily maximum and minimum temperatures of the field experiment during the trial were 32.98 • C and 20.71 • C, respectively.

Drought Stress Treatments
During FP development, soil moisture was controlled by means of drip irrigation at 80 ± 5% of the field water capacity (FWC) (well-watered, CK), 60 ± 5% of FWC (light drought, LD), and 45 ± 5% of FWC (moderate drought, MD). The meter was checked every morning and evening throughout the growth period to guide adjustments of the soil moisture. When the treated soil moisture dropped towards its lower limit, moderate drip irrigation was carried out until the upper limit level was reached, and the irrigation volume was measured by a water meter. After 30 DAD, the drought treatment plots were rehydrated to the CK level. Other field management measures reflected standard field management practices.

Measurement of Morphology and Microscopic Observation of Female Panicles
Plant height was measured from the ground to the top of the leaves in their natural growth state at 20, 25, and 30 DAD. The length and width of all the green leaves were measured by ruler in order to calculate leaf area, and the leaf area index (LAI) was determined according to this method [33]. Dry matter accumulation in stalks, leaves, tassels and ears of maize plants were measured at 20, 25, and 30 DAD. After 30 min of defoliation at 105 • C, dry weight was determined after being dried at 75 • C until a constant weight was reached. The percentage of drought limitation was calculated as (T2−T1)/T1. Here, T2 was shoot dry matter under the MD or LD treatment, while T1 was shoot dry matter under the control or LD treatment [34]. FPs were dissected with a dissecting needle at 20 and 27 DAD and analyzed under a stereomicroscope (Guanpujia, SMZ-B2, Beijing, China) to observe the differentiation of developing female inflorescences. There were three biological replicates in each group.

RNA Isolation and Illumina Sequencing
Three plants with consistent growth were selected from each treatment at 30 DAD, FP bracts were sampled, and the upper, middle, and lower parts of the ears were mixed evenly and then frozen at −80 • C prior to RNA-sequencing (RNA-seq) analysis. Total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion, Inc., Austin, TX, USA). RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Shanghai OE Biotech (Shanghai, China) conducted RNA-Seq library construction and high-throughput sequencing based on total RNA from the female inflorescence. The libraries were sequenced on the Illumina sequencing platform (HiSeq 2500l, Illumina, San Diego, CA, USA) with 125 bp paired-end reads. The raw RNA-seq data have been uploaded to NCBI SRA (BioProject ID: PRJNA604094).

Read Mapping and Differential Expression Analysis
Base calling was conducted using the raw image data generated by sequencing to obtain sequence data, and the called raw data (raw reads) were stored in fastq format. Raw data (raw reads) were processed using Trimmomatic [35]. The reads containing poly-N runs and low-quality reads were removed to obtain the clean reads. Then, the clean reads were mapped to the reference (NCBI_B73_v4) genome [36] using HISAT2 (version 2.2.1.0) [37]. The Fragments Per kb Per Million Reads (FPKM) values were calculated using cufflinks (version 2.2.1) [38,39], followed by differential expression analysis using DESeq (version1.18.0) [40]. Genes with |fold change| > 2 and p < 0.05 were identified as differentially expressed genes with p presented as raw p-values rather than FDR adjusted p-values.

GO and KEGG Enrichment Analysis
The significantly expressed Gene Ontology (GO) terms were selected by GO enrichment analysis according to the GO database (http://geneontology.org/). The differences in the frequency of assignment of GO terms in the DEG set were compared with the expressed genes in the CK, LD, and MD samples (p < 0.05). Functional groups encompassing DEGs were identified based on GO analysis, and pathway analysis was conducted according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/), with manual reannotation based on several databases and a literature search.

Differential Expression Verification by Quantitative Real-Time PCR (qRT-PCR)
Transcriptome sequencing data were validated by qRT-PCR. Total RNA was reverse transcribed using EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TRANS, Beijing, China). The qRT-PCR experimental methods used HiScript Q RT SuperMix for qPCR (Vazyme, Nanjing, China). The primer sequences were designed using Primer 5 and are listed in Supplementary Materials Table S1. The relative quantification 2 −∆∆Ct method was used to calculate the expression level of target genes in different treatments.

Statistical Analysis
Data collation and graphic rendering were conducted with SIGMAPLOT 14.0 (Systat Software Inc., San Jose, CA, USA), Microsoft Excel, and Microsoft PowerPoint 2016 software (Microsoft Corp., Redmond, WA, USA). All data are expressed as the mean ± SD value from three independent experiments unless otherwise stated. Data were analyzed by one-way ANOVA using Duncan's multiple range test at a p < 0.05 significance threshold in SPSS (IBM Corp., Armonk, NY, USA).  Table S2).

Overview of RNA Sequencing and Mapping
A total of 49.42 million raw reads were obtained from PH6WC transcriptome libraries (Supplementary Materials Table S3). More than 96.72% of them (47.80 million clean reads) remained after discarding low-quality reads and reads containing adaptor sequences, which were then used for downstream analyses. The clean reads were mapped to the B73 reference genome (ZmB73_RefGen_v4). Overall, 94.34-94.52% of clean reads from nine samples were mapped onto the reference genome (Supplementary Materials Table S3). On average, approximately 43.93 (91.29%), 43.73 (91.13%), and 43.70 (91.17%) million reads from the CK, LD, and MD treatments, respectively, were uniquely mapped onto the reference genome.
Compared with the CK treatment, only 40 genes (log 2 foldchange > 1 and p < 0.05), including nine up-regulated and 31 down-regulated genes, showed significantly differential expression in the LD treatment, and a total of 212 up-regulated and 664 down-regulated genes were identified in the MD treatment. A total of 887 DEGs, including 208 up-regulated and 679 down-regulated genes, were identified in the MD versus LD comparison (Figure 3a). A Venn diagram of the DEGs illustrated that there were 10 common genes that appeared in the LD vs. CK and MD vs. CK comparisons, five genes shared between the LD vs. CK and MD vs. LD comparisons, and 565 genes shared between the MD vs. LD and MD vs. CK comparisons. However, there were no DEGs commonly expressed in all three comparisons (Figure 3b).

Gene Expression Validation by qRT-PCR
To investigate the changes in gene expression at the mRNA level, eight randomly selected genes and three specific genes cuc2 (LOC103629107), TE1 (LOC541683), DLF1 (LOC100037791) were analyzed using quantitative real-time RT-PCR for validation of RNA-seq. The level of expression of the genes amplified is shown in Figure 4. Raw data were compared to transcriptomics data (Supplementary Materials Table S4), which closely resembled each other, validating the differential expression of the genes identified as being under drought stress.

GO Annotation and Enrichment
A total of 26, 629, and 630 DEGs were assigned by GO analysis conducted based on the genes from the LD vs. CK, MD vs. LD, and MD vs. LD comparisons, respectively. The most significantly regulated 20 terms among biological processes from the MD vs. CK and MD vs. LD comparison genes are shown in Figure 5a,b, but there were no significant terms (i.e., terms with gene number > 2 and p < 0.05) resulting from the LD vs. CK comparison. The up-regulated terms from the MD vs. CK comparison are involved in "regulation of timing of plant organ formation," "developmental process," and "regulation of cell proliferation." The down-regulated terms "response to water deprivation" and "post-embryonic plant morphogenesis" were also enriched. The up-regulated terms from the MD vs. LD comparison are "regulation of timing of plant organ formation," "regulation of cell proliferation," and "gibberellin biosynthetic process." The down-regulated terms "reductive pentose-phosphate cycle," "phosphate ion homeostasis," and "ethylene-activated signaling pathway" were enriched among genes from the MD vs. LD comparison. The genes associated with GO terms related to development, growth, and responses to stimulus were also significantly different among the three comparisons. These genes that changed in their levels of transcriptional expression were basically the same in the MD vs. CK and MD vs. LD comparison groups, but lower in expression differences in the LD vs. CK comparison group (Figure 5c,d, Supplementary Materials Table S5).

Metabolic Pathways Related to Soil Drought Stress
To further characterize genes affected by drought stress, we performed a KEGG pathway classification analysis to identify functional enrichment of DEGs. Thus, 8, 72, and 74 terms were significantly enriched in the transcriptome profile comparisons of LD vs. CK, MD vs. CK, and MD vs. LD groups (Supplementary Materials Table S6). The significant differences in the top 20 enriched KEGG pathways in the MD vs. CK and MD vs. LD comparisons are shown in Figure 6. In the MD vs. CK comparison, genes associated with the pathway "Starch and sucrose metabolism" were most enriched followed by those associated with "plant hormone signal transduction" (Figure 6a, Supplementary Materials Table S7). In the MD vs. LD group, "Plant hormone signal transduction," "Starch and sucrose metabolism," and "Glycosphingolipid biosynthesis" were the three most enriched terms (Figure 6a, Supplementary Materials Table S7). Pathway entries with the corresponding number of genes (among those pathways with more than two genes) are shown, and the corresponding -log 10 p-value of each entry is sorted in descending order. The number of DEGs in each pathway is positively related to the size of plot, and the p-values shown in red are more significant.

Responses of Plant Growth and Female Panicle Differentiation to Soil Drought Stress
The panicle stage is the most important productive stage in corn development, and soil drought stress in this stage can affect the plant growth rate, prolong the growth processes of the panicle stage, hinder the normal differentiation and development of ears, and ultimately lead to decreased crop seed setting rates and yields [42][43][44][45][46]. Further, the drought response depends on the time and intensity of water loss as well as the developmental stage [47,48]. In this study, LD and MD treatments compared with the CK treatment decreased green leaf area and significantly suppressed shoot dry matter accumulation over the prolonged drought treatment, and relative to the LD treatment, the MD treatment affected plant growth much more (Figure 1), which is consistent with previous research by Boonjung et al. [49].
FPs are the precursor to maize ears, and the differentiation and development of FPs mainly occurs from the V9 to VT stages, which include growth cone extension (V9), spikelet differentiation (V11), floret differentiation (V12), and formation of the sexual organs (VT). Developing organs are sensitive to drought, especially during their early phases [50]. When drought occurs between the V9 and VT stages, how does the degree of drought affect the formation of ears? In our study, spikelet and floral differentiation, as observed under stereomicroscope, were significantly inhibited and thus delayed by soil drought (Figure 2). Some studies have shown that the number of kernel rows is determined at the spikelet differentiation stage, and the floret differentiation period is the key period that affects grain number [51][52][53]. Here, mature ears in the MD and LD treatments were much shorter and thinner than those in the CK treatment; in addition, the bald tip length and number of unfilled grains were both increased under drought treatments. As indicated above, drought affected grain yields as well (Figure 2c).

Genes Involved in Development and Growth in Response to Soil Drought Stress
Drought treatments affected the expression of genes associated with development and growth of the inflorescence (Figure 4). Terminal ear 1 (te1) maize mutants have shortened internodes, abnormal phyllotaxy, leaf pattern defects, and partial feminization of tassels [54]. Similarly, cup-shaped cotyledon 2 (cuc2) mutants have been reported to have abnormalities in the regulation of the shoot meristem boundary and formation and subsequent development [55,56]. Here, cuc2 were up-regulated under moderate drought stress (Figure 4, Supplementary Table S5), combined with developmental change (Figure 2a), implying that the differential expression of the gene under MD treatments may be related to the mature delay of the FPs tissue. DLF1 was also up-regulated under MD stress (Figure 4, Supplementary Materials Table S5), suggesting that the trans-activator protein encoded by, this gene plays an important role in the signal transduction pathway and the regulation of plant growth at the FP development stage [57].

Genes Involved in Auxin Signaling in Response to Soil Drought Stress
Auxin is an important phytohormone that is closely related to plant resistance to adverse environmental conditions, and it can induce rapid and transient expression of some genes, including auxin response factor genes (ARF) and primary auxin response genes (Aux/IAA, GH3, SAUR and LBD); the protein products of these genes can specifically bind to ARFs to activate or inhibit downstream gene expression under drought [58][59][60][61]. In the current study, auxin signaling genes were involved in the response to drought, as the expression of IAA-conjugating genes (GH3) was up-regulated, and the expression levels of auxin biosynthesis genes were down-regulated after MD stress (Figure 7), leading to the reduction of auxin levels (Supplementary Materials Figure S1a). This implies that drought improves GH3 transcription to help maintain endogenous auxin at an appropriate level in maize [62,63]. However, when the concentration of auxin increases, auxin combines with transport inhibitor response 1 (TIR1), causing Aux/IAA ubiquitination and degradation; then, ARF is released, which further activates the expression of small auxin-up RNA (SAUR) genes [64]. SAUR genes are early auxin-responsive genes involved in plant growth, and the SAUR family regulates a series of cellular, physiological, and developmental processes in response to environmental signals [65][66][67]. In our study, three SAUR genes were down-regulated under the MD treatment (Figure 7), which may explain MD-induced inhibition of maize growth.

Reactive Oxygen Scavenging System and Ion Channel in Response to Soil Drought Stress
Limited water supply enhances the production of reactive oxygen species (ROS) [68,69], and plants are protected by glutathione S-transferase (GST) and other antioxidant enzymes scavenging excessive ROS from the damage caused by ROS [70]. This is because GST comprises a large superfamily of multifunctional protein and participates in ascorbic acid (AsA)/glutathione (GSH) cycling pathways [68]. Here, probable glutathione S-transferase GST12 was significantly down-regulated under the MD treatment compared with the CK treatment. However, GSTU6 (LOC103637303) and GST activity was up-regulated under MD stress (Supplementary Materials Table S8 and Figure S1b), which may scavenge ROS and protect both plant cell membrane structure and protein activity [71,72], implying that GST is involved in responding to drought stress [73,74].
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/2/313/s1, Table S1: The primer sequences for qRT-PCR, Table S2: Ear and grain yield traits at harvest time under CK, LD and MD treatments, Table S3: Number of reads based on RNA-Seq data of CK, LD and MD treatments, Table  S4: RNA-Seq expression levels of the eight genes for qRT-PCR, Table S5: Expression pattern of the differently expressed genes about development progress (A) and growth (B) under CK, LD and MD treatments, Table S6 (FIRI2017-19).

Conflicts of Interest:
The authors declare no conflict of interest.