Developmental and transcriptional responses of maize to drought stress under field conditions

Abstract Drought is a common abiotic stress which significantly limits global crop productivity. Maize is an important staple crop and its yield is determined by successful development of the female inflorescence, the ear. We investigated drought stress responses across several developmental stages of the maize B73 inbred line under field conditions. Drought suppressed plant growth, but had little impact on progression through developmental stages. While ear growth was suppressed by drought, the process of spikelet initiation was not significantly affected. Tassel growth was reduced to a lesser extent compared to the observed reduction in ear growth under stress. Parallel RNA‐seq profiling of leaves, ears, and tassels at several developmental stages revealed tissue‐specific differences in response to drought stress. High temperature fluctuation was an additional environmental factor that also likely influenced gene expression patterns in the field. Drought induced significant transcriptional changes in leaves and ears but only minor changes in the tassel. Additionally, more genes were drought responsive in ears compared to leaves over the course of drought treatment. Genes that control DNA replication, cell cycle, and cell division were significantly down‐regulated in stressed ears, which was consistent with inhibition of ear growth under drought. Inflorescence meristem genes were affected by drought to a lesser degree which was consistent with the minimal impact of drought on spikelet initiation. In contrast, genes that are involved in floret and ovule development were sensitive to stress, which is consistent with the detrimental effect of drought on gynoecium development and kernel set.

Understanding the underlying physiological and genetic processes controlling maize reproductive development under water-limited conditions is crucial to this goal.
Among the major cereal crops, maize is the only monoecious plant bearing unisexual flowers. The male inflorescence, or tassel, develops from the shoot apical meristem at the top of the plant, whereas the female inflorescence, or ear, develops from lateral meristems in the axil of leaves. The initial steps of development in both inflorescences are very similar; however, sex related differences appear at later developmental stages (Cheng, Greyson, & Walden, 1983). The basic unit of the ear inflorescence is the pistillate spikelet, composed of a floret subtended by a pair of glumes.
Development of the ear is defined by several meristem transitions starting with spikelet pair meristems (SPMs) arising from the flanks of each inflorescence meristem (IM). Each SPM gives rise to two spikelet meristems (SMs) which terminate as a floral meristem (FM) (Cheng et al., 1983;Tanaka, Pautler, Jackson, & Hirano, 2013). The FM produces different organ primordia, including the gynoecium (female reproductive structure) which terminates in formation of an ovary composed of the embryo sac and the silk, a stigmatic structure. The ear is protected by several layers of husks, which are modified leaves. When silks exert from the husk, the female flowers are fertilized by pollen shed from the tassel and kernel development commences.
In the US Corn Belt the most severe consequences of drought occur at flowering time (Shaw, 1977) and result in multiple detrimental effects on the female inflorescence and to a lesser extent negative effects on the male inflorescence (Herrero & Johnson, 1980a). Silk elongation is dramatically reduced under limited water conditions due to a reduction in cell division and cell expansion (Fuad-Hassan, Tardieu, & Turc, 2008). The delay of silk extrusion results in a longer ASI (anthesis to silking interval) and can reduce pollination efficiency (Araus, Serret, & Edmeades, 2012). Even short periods of water deficits can lead to abnormalities in embryo sac development (Moss & Downey, 1971), as well as zygotic and early kernel abortion (Schussler & Westgate, 1990). For these reasons drought conditions occurring just before and shortly after pollination have the most profoundly negative effect on kernel set and final grain yield (Araus et al., 2012). In contrast, drought episodes during the grain filling period typically have a less severe impact on yield (Barker et al., 2005).
Despite the importance of maize ear growth and development under drought conditions, molecular and genomic data for this organ are surprisingly limited. Several RNA profiling expression studies under drought conditions have been conducted on maize seedlings or leaves (Avramova et al., 2015;Badicean, Scholten, & Jacota, 2011;Fernandes, Morrow, Casati, & Walbot, 2008;Shan et al., 2013;Yue, Zhuang, Li, Sun, & Zhang, 2008;Zheng et al., 2010). There are reports on the drought stress response in maize reproductive tissues including pre-pollinated ears (Zinselmeier et al., 2002), ears, and silks Zhuang et al., 2007), as well as ovaries and developing kernels (Kakumanu et al., 2012;Luo, Liu, Lee, Scully, & Guo, 2010;Marino et al., 2009). Most of these studies were conducted in a controlled environment and led to important discoveries of stress response in maize. However controlled environments lack fluctuating meteorological conditions such as solar radiation, heat, wind, and vapor pressure deficit, which play key roles in plant growth and development. The complexity of natural environments emphasizes the importance of conducting drought experiments in the field (Blum, 2014). To address this issue, we designed and conducted two drought stress experiments under field conditions. The first experiment was located at Woodland, California (USA), which has been developed as a managed-stress environment for phenotyping the drought response of maize (Campos, Cooper, Habben, Edmeades, & Schussler, 2004). Drought stress was imposed over 1 month with the maximum intensity of stress occurring around flowering time. A second experiment was conducted in Johnston, Iowa (USA), within the Corn Belt region. Drought conditions were created by growing maize in pots located in the field and precisely limiting water application. For both experiments, we collected vegetative and reproductive phenology data to understand the tissue specific response to drought stress.
To complement the observations of phenotypic responses, RNAseq profiling of leaves, ears, and tassels was conducted to determine relative differences in transcriptomic responses between fully irrigated and drought stressed plants. Overall, these field drought experiments revealed a tissue specific response to stress at both the phenotypic and transcriptional levels, as well as significant flexibility of gene expression in response to changing environmental conditions.

| MATERIAL S AND ME THODS
2.1 | Experimental design and drought stress management in the field at Woodland, CA | 3 DANILEVSKAYA Et AL.
B73 seeds were planted in 4-row plots in a randomized complete block design with four replications. The length of the plot was 4.4 m, 29 plants per plot. Plant density was 82 K plants per ha. Standard agronomic practices were used, as described in Gaffney et al., 2015. Well-watered (WW) plots were fully irrigated throughout the entire growth cycle of the crop. In contrast, plots that were designated to receive the drought treatment (DRT) were fully irrigated until plants reached the V7-V8 stage, when irrigation was terminated for the remainder of the experiment. Due to the high water holding capacity of the soil (Yolo silt loam) at this location, plants in this drought stress treatment continued to develop for several weeks after the irrigation was stopped. Daily high and low temperatures were obtained from weatherunderground.com. The nearest weather station (ID: KCAWOODL9) is located 3.6 miles away from the Woodland field.

| Tissue sampling for RNA-seq profiling in the field study
Tissues were sampled from plants of the middle two rows of each plot in both the WW and DRT treatments. The middle portions of the 10th-11th leaves from three plants were cut and these tissues were pooled to make a biological replicate. Leaves were sampled at the peak of stress in the afternoon between 14:00-15:00 p.m. from both well-watered and drought plots during the diurnal stress period. Within each plot, ears of uniform size were selected and pooled as the biological replicate. Prior to sampling ears, shoot bags were used to cover exposed prophylls prior to silk emergence. For the first sampling at V10-V12 (designated as V12), 32 ears <0.5 cm in length were pooled per plot. For the second sampling at V12-V14 (designated as V14), 16 ears 0.8-1.5 cm in length were pooled per plot.
For the third sampling at V16-V18(designated as V18), 3-4 ears that were 3-4 cm in length under WW and 2-3 cm in length under DRT were pooled per plot. For the final sampling at V18-R1(designated as R1), 3-4 ears that were 7-8 cm in length under WW and 3-4 cm in length under DRT were pooled. In the drought stressed plants, the R1 designation based on the timing of silk emergence in wellwatered plants, as the drought stressed ears did not exert silk beyond the husk at this stage. Visible silks were manually removed from ears at the time of sampling and whole ears were used. Tassels were also sampled from the same plants used for ear sampling. The 10 cm middle section of the main tassel spike was removed for the tassel tissue sample. Tassels were still hidden in the whorl at the first (V12) and second (V14) samplings. Tassel emergence occurred at the third (V18) sampling. Shedding was initiated just prior to the fourth (R1) sampling. For all samples, tissues were immediately frozen in liquid nitrogen in the field.

| Experimental design and drought stress management in the field pot study
The second field study was conducted in Johnston, IA (41°40′N, 93°42′W) during the summer of 2013. Plants were grown in pots that were arranged in double rows in wooden racks. The plant distance within a row was 20 cm. The row space between the two rows in the rack was 24 cm. The row space between racks was 110 cm. This arrangement produced 26 plants per 4.4 m which is close to the number of plants in a field plot (29 plants/4.4 m). To prevent root growth into soil, pots were supported above the ground on the wooden racks. Drainage holes in pots were covered by paper towels. Pots were covered with plastic bags to prevent rain water entry. The pots were laid out in a split plot (strip-plot) design with 24 reps per stage equally divided between WW and DRT treatments. Square tree pots (Stuewe&Sons TP818) of 10 L were filled with potting soil (Fafard ® 3B mix). All pots for the DRT treatment were filled by weight with the same amount of soil. Water was provided to each pot via irrigation drip tubes and Peters ® Excel fertilizer was used with a final N concentration of 100 ppm delivered to the plants. The fertilizer contained 15% N, 5% P 2 O 5 , 15% K 2 O, 5% Ca, 2% Mg, 0.0187% B, 0.0187% Cu, 0.075% Fe, 0.0375% Mn, 0.0075% Mo, and 0.0375% Zn. Plants were grown under full irrigation until they reached the V7-V8 stage. Irrigation was stopped from the DRT pots until 50% of the plants showed obvious leaf wilting. Pots were then re-watered to full soil capacity. This drying/rehydration cycle was repeated eight times sequentially from stage V8 to R3 over the course of the experiment.

| Phenotypic data collection in the field pot study
Growth stages were defined according to the book "Corn growth and development" (Abendroth, Elmore, Boyer, & Marlay, 2011).
Plant height was measured from the plant base to the collar of the top fully expanded leaf. Plant height and leaf number were recorded weekly. Pollen shed and silk emergence were recorded on a per plant basis. Ears were dissected from V10 to V14 in the lab and analyzed under a dissecting microscope. At the R1-R3 stages, measurements were done without a microscope. Ear length was measured from the base to the tip. Two rows of spikelets at the opposite sides of the ear were counted at every stage to determine total ovaries/ kernels per row.

| Growing degree day
Growing degree day was calculated using the standard formula GDD C = (T max + T min )/2) − T base where T max and T min are the daily high and low temperatures (Celsius). The base temperature (T base ) of maize was 10°C (Abendroth et al., 2011;Yang, Logan, & Coffey, 1995).

| Statistical analysis
A linear regression formula with R 2 value (Microsoft Excel) was used to determine the plant growth rate, the leaf appearance rate, the ear elongation rate, and spikelet initiation rate. T tests were calculated using the Microsoft Excel 2015 Analysis ToolPak add-in. All data were analyzed with the t test: Two-Sample Assuming Equal Variances using a two-tail approach and a p-value of 0.01.

| mRNA-seq library preparation and transcriptome data analysis
Total RNA isolation and deep cDNA sequencing was performed on the Illumina HiSequation 2500 system as described previously (Thatcher et al., 2014(Thatcher et al., , 2016. The generated RPKtM (Relative parts per kilobase per 10 million) data matrix was visualized and analyzed in Genedata Analyst ™ software (Genedata AG, Basel, Switzerland).
Genedata Expressionist ® for Genomic Profiling was used for RNAseq data analysis and visualization. Several statistical applications such as t test, ANOVA, linear models, and Principal Components Analysis (PCA) are built in this software. RPKtM values were normalized base-2 logarithm before all statistical analysis. The total number of expressed genes was detected by mapping reads to maize B73 reference genome sequence V2 as described in (Thatcher et al., 2016). One-way ANOVA analysis was performed to detect differentially expressed (DE) genes at a Q-value of 1E-6 with a combined effect of developmental stages and drought treatment. False discovery rate was corrected for by using a multiple hypothesis testing method (Benjamini, Drai, Elmer, Kafkafi, & Golani, 2001). A default of 1E-6 FDR was used for selecting of genes that are significantly affected by drought and stage. DE genes with a q-value of 1E-6 were designated as significant. PCA, a statistical procedure that models the variation in terms of its principal components, was used to reveal the impact of drought stress and developmental stage on DE genes by tissue. Student's t tests were applied to identify genes significantly affected by drought stress at each of the four developmental stages where a false discovery rate of 1E-2 was used as a cutoff.
The stage specific DE genes were subjected to K-means clustering with positive correlation distance to identify up-or down-regulated genes. Enriched function analysis was performed with GO Fisher's exact test (GOFET, p-value of 0.01) with biological processes as the default main ontology. GOFET is an integrated part of Transcriptome Data Analysis of the software package Genedata Analyst ™ .

| Quantitative RT-PCR
Quantitative RT-PCR amplifications were performed using TaqMan ® probe based detection system (Applied Biosystems). The quantity of target genes was determined by the standard curves of cDNAs pooled from the tissues where the genes are highly expressed. The relative expression level of genes was calculated by their quantification normalized to ubiquitin 5. The processes were performed following the User Bulletin #2 from Applied Biosystems at http:// www3.appliedbiosystems.com/cms/groups/mcb_support/documents/generaldocuments/cms_040980.pdf.

| Accession numbers
Sequence data from this article can be found in National Center for Biotechnology Information Gene Expression Omnibus under accession number GSE71723.

| Phenology of plants under well-watered and drought stressed conditions
The initial field experiment was conducted during the summer of 2012 in Woodland, CA in a randomized complete block design.  thereafter, in the drought (DRT) treatment plots irrigation was terminated. The well-watered (WW) plots were fully irrigated during the remainder of the experiment, and no rainfall occurred during the treatment portion of the study (Supporting Information Figure S1).
The first sampling was 11 days after the irrigation was terminated.
During these first 11 days, the temperature reached a high of 40°C, imposing heat stress on the plants. At the first sampling (June 19th), no visible signs of stress were observed in the WW plots ( Figure 1a).
However, a water deficit response was observed in the DRT plots in the form of plant leaf rolling and leaf wilting (Figure 1b) which is a typical manifestation of the low leaf water potential (Fernandez & Castrillo, 1999). The week before the second sampling (June 27th), the temperature was milder, dropping to approximately 26°C. The temperature slowly increased after the second sampling and continued to rise during the third and fourth samplings.
Vegetative and reproductive traits were collected for plants grown in both treatments. Plants in the DRT treatment had a final height reduction in 40% compared to the WW treatment, indicating the severity of the stress (Table 1). Plants in the DRT treatment produced 1-2 fewer leaves than WW plants (Table 1). Leaf appearance rate was also slightly delayed under DRT (Table 1 and Supporting Information Figure S2A). Plants flowered (pollen shedding and silk emergence) on the same day (July 10th) in the WW plots. In contrast, plants in the DRT treatment shed pollen 3 days after the WW plots and no silk exertion was observed. DRT treatment ears were dissected at the end of the experiment (July 10th) and showed minimal elongated silks compared to WW ears (Figure 1d,e), resulting in ears with no kernel set.
The length of the main tassel spike was reduced by ~20% in DRT plants compared to WW plants ( Table 1). The ear elongation rate was reduced approximately 1.6-fold under DRT relative to WW, which resulted in 30% smaller ears under stress (Table 1 and Supporting Information Figure S2B). Remarkably, the final spikelet number per row was not statistically different between treatments.
This observation suggests that ear elongation is more susceptible to drought stress than total spikelet number.
To further substantiate this observation, a precision phenotyping experiment was executed using a field pot approach (Supporting Information Figure S4) in Johnston, IA during the summer of 2013. In general, the weather was hot and dry with little precipitation during the experiment (Supporting Information Figure S3 Height and leaf number were collected from 28 well-watered (WW) and 32 drought (DRT) treated plants at GDD C 836. Staging notes were taken from V7 to R1. Tassel length was measured for 11 (WW) and 16 (DRT) plants at R1. Leaf appearance rate and ear elongation rate was calculated using linear regression based on data in Supporting Information Figure S2. Measurements represent mean ± SD. One day is approx. 11.4 GDD C , growing degree days. a DRT means are statistically different from WW at p < 0.01.
slightly reduced under DRT (Table 2 and Figure 2b). There was no significant difference in final leaf number between DRT and WW treatments (Table 2). However, there was a significant increase in the anthesis silking interval (ASI) in the DRT treatment compared to WW (Table 2).
To investigate ear growth and development, primary ears were dissected at five developmental stages starting at V10 (Supporting Information Figure S5). Ear length and spikelet number (measured as the number of spikelets along the length of the ear) were recorded for each sampled ear (Figure 2c,d). Ear elongation rate after V14, as estimated by linear regression models, was decreased by twofold in DRT relative to WW (Table 3 and Figure 2e) which resulted in 50% shorter ears under DRT (Table 3 and Supporting Information Figure   S6). This was similar to the reduction in ear length observed in the Woodland experiment (Table 1 and Supporting Information Figure   S2).
Spikelet initiation plateaued around V14 under both conditions ( Figure 2d). The estimated spikelet initiation rate before V14 was slightly lower under DRT, 2.13 spikelets per day, compared to 2.43 spikelets per day for WW (Table 3 and Figure 2f). The maximum spikelet number was reduced by 13% under DRT, relative to WW, while ear length was reduced by 50%. This is consistent with the 2012 study in Woodland. To investigate whether the spikelet number is less susceptible than ear elongation to stress in germplasm other than B73, we measured nine Pioneer proprietary inbred lines in the 2013 experiment. The results confirmed that, in general, ear elongation is more susceptible to DRT than ear spikelet number ( Figure 2g,h).
Overall, the B73 response to drought stress in the field pot study was comparable to the Woodland field experiment, which suggested that the level of stress was similar in both field experiments. The DRT treatment began at the same developmental stages in both locations and the weather was fairly consistent over the course of both experiments. Data from both studies demonstrated that overall plant growth was more sensitive to drought than progression through developmental stages. The same pattern was observed for the ear, where length was more affected by drought than spikelet initiation.

| Analysis of RNA-seq expression in leaf, ear, and tassel samples
To conduct RNA-seq transcription profiling, leaves, ears, and tassels were systematically sampled after initiation of the DRT treatment in the Woodland experiment (Supporting Information Figure S1). At the first sampling plants were at V12 (ears < 0.5 cm), followed by a second sampling at V14 (ears 0.8-1.5 cm). The third sampling point was at V18 (ears 3-4 cm under WW), followed by a final sampling at R1 (ears 7-8 cm under WW). No silk emerged under DRT conditions even though it was formed inside of the husk. Ninety-six RNA-seq libraries were generated and were comprised of three tissues, four sampling times, four biological replicates, and the two treatments (WW and DRT). Quality control and mapping reads to the B73 reference genome are described in a recent publication (Thatcher et al., 2016). The total number of expressed genes was in the range of 31,000-37,000 by tissue (Supporting Information Figure S7). The number of DE genes varied substantially by tissue, with leaves having the fewest (~3,500), followed by ear (~7,000), and finally tassel (~20,000) (Supporting Information Figure S7).
Principal component analysis was used to identify the major experimental factors (components) accounting for all DE genes across tissues, developmental stages, and treatments. Leaf WW samples were primarily separated from DRT samples based on water treatment rather than development stage, suggesting that the majority of DE genes were responding to the drought treatment in mature segments of the leaf ( Figure 3a). In contrast, WW and DRT ear samples were primarily clustered by development stages rather than by water treatment, reflecting developmental processes occurring in the ears (Figure 3b). At the later stages, a treatment specific clustering began to emerge. This suggests that most DE genes in ears are initially driven by development, but begin to respond to treatment over time ( Figure 3b). Virtually no treatment effect was detected in the tassel samples ( Figure 3c) suggesting that most DE genes were related to development rather than a result of watering treatment.
In order to identify DE genes responding to drought, K-Means clustering was performed at every sampling stage for WW and DRT samples of each tissue with a Q-value cutoff of 1E-2 from the set of DE genes (Supporting Information Figure S8). K-Means clustering was well supported by four biological replicates per sampling (Supporting Information Figures S9-S11). The clustering was then used to populate Gene Ontology (GO) categories.
It is important to consider the weather pattern during the Woodland field experiment, which imposed heat as an additional stress factor to drought ( Figure 4a). During the first week of drought stress, daily temperatures were 35⁰C to 40⁰C, exposing plants at the V12 sampling to not only drought stress, but also heat stress. The week before the V14 sampling, daily temperatures TA B L E 2 Vegetative and flowering traits of B73 in the field pot study, Johnston, IA 178.4 ± 9.9 a 3.57 20.0 ± 0.9 0.29 750.5 ± 13.5 a 798.4 ± 11.0 a 4.8 ± 1.9 a Vegetative traits were collected weekly. Final measurements were collected at R2 stage. Growth rate and leaf appearance rate were calculated using linear regression based on data in Figure 2. Measurements represent mean ± SD. One day is approx. 11 GDD C , growing degree days. a DRT means are statistically different from WW at p < 0.01.
were 25⁰C and plants experienced drought stress but less heat stress. The daily temperature then progressively increased from 30°C to 40°C, suggesting the highest level of combined heat and drought stress occurred after 32 days of drought stress at the R1sampling date.
Every tissue displays unique changes at the transcriptional level in response to abiotic stress. Leaf transcriptional changes increased, in general, over the course of the drought treatment, reaching the maximum number of DE genes at 32 days of drought (R1) (Figure 4b).
At V14 when the temperature was mild, the number of DE genes was F I G U R E 2 Plant growth and development under WW and DRT in the field pot study, Johnston IA, 2013. Effect of water treatment on B73 plant growth (a) and leaf appearance (b) by GDD C depicted as linear trend lines. Effect of water treatment on ear length (c) and spikelet number (d) by developmental stage. Linear trend lines of ear elongation after V14 (e) and spikelet initiation before V14 (f) by GDD C in B73. Effect of water treatment on average ear elongation (g) and spikelet initiation (h) in nine pioneer proprietary inbred lines by GDD C . Trend lines were calculated using linear regression models for each trait and water treatment. The x value in the linear regression formula denotes average change in each trait per 1 GDD C . R 2 describes how well the data fit the trend line where 1.0 is a perfect fit. Data points represent means ± SD. Traits separated by WW (solid blue) and DRT (dashed red) treatments The heat stress level was the lowest at this time. However, at the highest level of drought and heat stress (third and fourth samplings) few DE genes were found in tassel samples. This pattern suggests that the tassel transcriptome is less responsive to environmental stress which agrees with the phenotypic observation that tassel development is less sensitive to drought stress relative to ear development.

| Gene ontology enrichment of biological processes for de genes in leaf samples
In order to gain insight into the functional categories of drought responsive genes, gene ontology (GO) Fisher's exact test was used to analyze DE genes separated by tissue and stage. The tissue-specific top 20 up-and down-regulated GO functional terms, excluding highlevel and synonym terms, are presented in Figure 5. The complete dataset is available in Supporting Information Table S1 (leaf), Table   S2 (ear), and Table S3 (tassel).

| Leaf up-regulated genes
A consistent enrichment across developmental stages was observed in the "protein folding" category, reaching a maximum at the R1 sampling ( Figure 5a). Genes in the "response to heat" category were enriched in the V12 and R1 samplings which followed the temperature  Ear elongation rate was calculated using data collected after V14 (GDD C 600). Spikelet initiation rate was calculated using data collected prior to V14 (GDD C 600). The maximum number of spikelets occurred at R1 (GDD C 750). Spikelet initiation rate and ear elongation rate were calculated using linear regression based on data in Figure 2. Measurements represent mean ± SD. One day is approx. 11 GDD C , growing degree days. a DRT means are statistically different from WW at p < 0.01.
pattern. Gene enrichments in "cellular amino acid metabolic process" and "proline biosynthetic process" categories were also enriched, suggesting enhancement of osmolyte biosynthesis, such as proline and other amino acids. Some enrichment was also observed in genes related to translation machinery, suggesting stabilization of protein biosynthesis under drought conditions.
Surprisingly, there was enrichment of up-regulated genes in functional categories related to photosynthetic machinery such as "photosynthesis light reaction and light harvesting", "electron transport in photosystem II", "chlorophyll biosynthetic process", "porphyrincontaining compound metabolic process", "tetrapyrrole biosynthetic process", "phylloquinone biosynthetic process", and "thylakoid membrane organization" (Supporting Information Table S1). At later stages, the level of enrichment gradually declined (Figure 5a). This suggests that there may be compensatory responses in leaves at the first sign of drought stress to maintain photosynthesis. As drought stress continues and less water is available, leaves begin to show daytime wilting (Figure 1b) which likely impedes the photosynthetic rate. A similar up-regulation of photosynthetic genes under drought stress in B73 leaves was recently reported (Avramova et al., 2015).
However, in that study, photosynthesis was inhibited by stress as well as stomata conductance. To reconcile this paradox, the authors provided some evidence that investment in the photosynthetic machinery under stress facilitates photosynthesis during plant recovery when water becomes available (Avramova et al., 2015).
Abscisic acid (ABA) is a primary abiotic stress response hormone and is expected to have a strong response to drought (Hauser, Waadt, & Schroeder, 2011). However, few genes were  and GRMZM2G389301, EID1-like F-box protein 3), which are homologs of the Arabidopsis gene EDL3 involved in the regulation of ABA-signaling (Koops et al., 2011). The ABA biosynthesis gene VP14 (GRMZM2G014392 viviparous14, 9-cis-epoxycarotenoid, NCED1 dioxygenase) was also up-regulated at the V12 sampling.

| Leaf down-regulated genes
The greatest number of down-regulated genes was found in the "organic anion transport" and "amino acid transport" categories ( Figure 5b). This indicates that with low water movement through the plant under DRT conditions a reduction in transport of solutes may occur. Down-regulation of genes in biological categories such as "response to chitin" and "salicylic acid stimulus" suggests that plant immunity to pathogen invasion may weaken under drought conditions. The trend of down-regulation of genes involved in the "response to reactive oxygen species" was also observed, which could cause reactive oxygen species to accumulate in leaves under drought stress.

| Ear up-regulated genes
The In the "response to abscisic acid stimulus" category, duplicated genes EDL3 (GRMZM2G073324 and GRMZM2G389301 EID1-like F-box protein 3), bZIP transcription factor ABI-5 (GRMZM2G077124), and the ABA receptor PYL8 gene (GRMZM2G165567) were all up-regulated. This is similar to what was observed in leaves, except that in leaves activation of the ABA receptor genes was not detected.
No enrichment of the typical stress response categories such as "protein folding" or "response to heat" were found in ear samples at the V12 or V14 samplings. However, these categories were enriched later at the V18 and R1 samplings when plants had been exposed to abiotic stress for a longer period of time. At the later stages, some enrichment of up-regulated genes was observed in the "cell cycle" and "DNA replication" categories which may be explained by delayed ear growth under drought conditions (Figure 1d,e and Supporting Information Figure S2B).

| Ear down-regulated genes
The categories with the greatest number of down-regulated genes in ears were "translation", "DNA replication", "cell division", and many related processes (Figure 5d and Supporting Information Table S2). These down-regulated genes were predominant at the V12 sampling, but were less impacted at later stages. These results suggest that fundamental processes such as protein biosynthesis, DNA replication and cell division are suppressed under a combination of heat and drought stress, resulting in a delay in ear growth.
As soon as environmental conditions were less stressful (V14 sampling), transcriptional activity of these key genes was restored.
Down-regulation of genes in the categories of "glucose metabolic process", "glycolysis", and "tricarboxylic acid cycle" was observed at the V12 and R1 samplings, indicating the sensitivity of energy related processes to stress conditions. Down-regulation of genes involved in carbohydrate metabolic processes was enriched at the R1 stage. Interestingly, genes in the "wax metabolic processes" were down-regulated (R1 stage) signifying that wax biosynthesis may be impaired.

| GO enrichment of biological processes for DE genes in tassel samples
In contrast to the ear and leaf, drought response in the tassel was much more limited (Figure 4b,c,d). The majority of the DE genes were seen at the V14 sampling. The most prominent enrichment among the up-regulated genes was found in processes related to "protein degradation", "telomere maintenance", "chromosome and chromatin organization", and "meiosis" (Figure 5e). Among the down-regulated genes, functional enrichments were observed in "translation" and "ribosome biogenesis" categories (Figure 5f). At the later V18 and R1 stages there were very few differences in gene expression between WW and DRT samples. Transcriptional differences appeared to reflect a slight delay in tassel development under DRT and not a direct effect of stress on the tassel transcriptome.

| Comparison of top ranking stress-induced categories in leaf and ear samples
Two functional categories "protein folding" and "heat response" were highly enriched in leaves and ears under prolonged drought stress (Figure 5a,c). Neither of these functional categories was enriched in the tassel. We detected 66 genes whose expression level was increased >twofold in at least one tissue at the R1 stage ( Figure 6). The majority of these were annotated as

| 13
DANILEVSKAYA Et AL. F I G U R E 6 Expression of genes in categories "protein folding" and "heat response" in leaves and ears at the R1 stage. The darkest color corresponds to the higher levels of induction under stress. Cutoff effective size was above twofold threshold at least in one tissue. A default of 10 −6 FDR was used heat-shock proteins (HSP). There were 61 HSP genes that were up-regulated in leaves and 36 HSP genes in ears. The subcellular localization of the predicted proteins suggested that about 17 of them are localized to the chloroplast, signifying their putative function in protecting the photosynthetic machinery. In leaves, the top up-regulated genes encode the peptidyl-prolyl cis-trans isomerase (GRMZM2G15468) and HSP26 (GRMZM2G149647) were induced by stress >110-fold. However, these genes were not expressed in ears. The top DRT induced genes in ears were the heat shock gene HSP90 (GRMZM5G833699) and

| Expression of developmental genes in stressed ears
It is important to note that there were few differentially expressed inflorescence development genes identified in stressed ears, even though the entire meristem spectrum (from IM to SPM to SM and FM) was sampled over the course of the experiment. Some degree of enrichment was observed among down-regulated genes at the V14 stage in the "inflorescence morphogenesis" category ( Figure 5d) and at the R1 stage in the "sexual reproduction", "cell morphogenesis", and "developmental process involved in reproduction" categories (Supporting Information Table S2). The small number of developmental genes identified in the stressed ear samples led us to conduct qRT-PCR of several characterized maize genes involved in the patterning of the maize inflorescences (Table 4). We tested six genes essential for axillary meristem initiation: VT2 (vanishing tassel2) (Phillips et al., 2011), SPI1 (sparse inflorescence1) (Gallavotti, Barazesh, et al., 2008), ZmPIN1a (pinformed) (Gallavotti, Yang, Schmidt, & Jackson, 2008), BAF1 (barren stalk fastigiated) (Gallavotti et al., 2011), BIF2 (barren inflorescence2) (McSteen et al., 2007), and BA1 (barren stalk1) (Gallavotti et al., 2004). The expression patterns of these genes were similar in WW and DRT ear samples (Figure 7). None were repressed by stress conditions and several were slightly up-regulated with a p-value of 0.01. Absence of repression under stress conditions is consistent with the phenotypic data that show a limited effect of stress on spikelet initiation and spikelet number.
We also investigated expression of MADS box transcription factors which are key regulators of floral development. Out of the 43 maize MIKC-type MADS box genes (Zhao et al., 2010), at least 15 genes were down-regulated in ears under stress (Table 4). Expression patterns were confirmed by qRT-PCR for eight selected MADS box genes (Supporting Information Figure S13). Non-MADS floral developmental genes such as TS6 (tasselseed6), AP2 (floral homeotic), and DL (drooping leaf homolog) were also down-regulated by stress (Table 4). These results suggest that drought stress may delay floral development due to the repression of floral developmental genes.
Down-regulation of genes related to ovule development and pollen tube receptivity were also found in this study (Table 4). Among them were IG1 (indeterminate gametophyte1) (Evans, 2007) and a homolog of the Arabidopsis FERONIA gene. FERONIA is receptor-like protein kinase that is essential for interaction between the synergids and the pollen tube during embryo sac fertilization (Escobar-Restrepo et al., 2007). Another function of FERONIA is the regulation of cell elongation (Haruta, Sabat, Stecker, Minkoff, & Sussman, 2014). Down-regulation of genes involved in embryo sac development suggests a negative effect of abiotic stress on the formation of functional ovaries.

| Plant growth and development is affected by stress in a tissue specific manner
Drought stress is an area of extensive research in many crops and model plants because of its negative impact on agriculture productivity (Blum, 2014). It has been well documented that the failure of the maize female inflorescence (ear) to develop under stress conditions is a primary reason for grain yield loss in the field (Araus et al., 2012;Barker et al., 2005;Boyer & Westgate, 2004;Campos et al., 2006). Previous genomic studies of maize reproductive tissues grown under water-limiting conditions provide insight into their general drought response, but fail to comprehensively correlate the maize ear development with transcriptome response to drought under field conditions.
We conducted experiments covering the sequence of ear development starting from initiation up to silk emergence under WW and DRT conditions. This sampling was done in the field to capture ear responses to drought under natural weather conditions. In addition, we systematically collected tassel and leaf samples as well as phenotypic traits for the comparison of tissue specific growth and developmental processes under drought stress.
In both field experiments, plants were exposed to significant drought stress. There were differences in irrigation and high temperature fluctuations between the two locations which may explain some of the small observed phenotypic differences. For example,  (Figure 2g,h).
Tassel growth (measured by final tassel length) is less sensitive to abiotic stress than ear growth, being reduced by 20% under drought stress relative to well water as compared to a 50% reduction in ear size under the same conditions (Table 1). This finding is consistent with previously published studies that tassels are less susceptible to drought than ears (Herrero & Johnson, 1980a). Overall, our phenotypic data demonstrate that maize organs have a differential response to abiotic stress. Ear growth is the most sensitive to Moreover, developmental processes such as leaf and spikelet initiation (organogenesis) are less affected by drought than organ growth, a pattern that matches our transcriptional analysis.  Table S3). Few DE genes were detected in the V18 and R1 samples, which were the periods of the most intense abiotic stress.

| Transcriptome response is tissue specific and adjusts to environmental conditions
Moreover, no activation of heat shock genes, which are biomarkers for stress responses, was detected in tassels. Our findings are in agreement with previous studies indicating that due to the shoot apical dominance, tassel growth is less sensitive to stress than ear growth (Herrero & Johnson, 1980a). However, heat, not drought is detrimental for pollen viability (Herrero & Johnson, 1980b;Schoper, Lambert, & Vasilas, 1985).
Activation of heat shock proteins (HSP) is a universal stress response in virtually all organisms. These proteins function as molecular chaperons, preventing proteins from misfolding, denaturing, and degrading in order to support vital functions during stress episodes (Al-Whaibi, 2011). Genes annotated in GO categories "protein folding" and "response to heat" were among the top up-regulated genes in the leaf and there was a significant overlap between stress-induced genes found in our study and those found in sorghum seedlings (Johnson et al., 2014), suggesting conserved stress responses in leaves of C4 grasses. Many leaf-expressed HSP have predictive subcellular localization in chloroplasts which indicates that they may function to protect chloroplasts from stress (Hu et al., 2015). In the ear, heat shock genes were activated two stages later than in leaves, with smaller numbers and lower fold changes. In total, only 30 heat shock genes were induced in ears, and the genes with the largest change were HSP90 and HSP70 which are also induced in the leaf under heat stress. Heat shock proteins of this type are thought to be molecular chaperons that stabilize newly synthesized proteins to support their function under unfavorable conditions (Al-Whaibi, 2011).
One possible interpretation of the different dynamics of heat shock proteins in the leaf versus the ear is related to different tissue functions. The leaf senses and then rapidly responds to the changing environment by the activation of heat shock genes to protect chloroplasts. In contrast, in the ear, it takes a longer time to activate heat shock genes where their function may be related solely to stabilizing proteins under stressful conditions. The stronger response of heat shock proteins in leaves may be explained by the importance of chloroplast maintenance and the threat of water loss through transpiration. No heat shock genes were induced in the tassel, suggesting that it may be insulated from stress better than the other organs.

| Retarded ear growth correlates with the down-regulation of DNA replication and cell-cycle genes and the up-regulation of oxylipin biosynthetic genes under stress
Like previous studies (Boyer & Westgate, 2004;Westgate & Grant, 1989) we have shown that ear growth is very sensitive to abiotic stresses. However, limited data are available to explain the molecular mechanisms of the inhibition of ear growth under stress at developmental stages prior to silk emergence. Our RNA-seq data show that the first stress response in early V12 stage ears is a massive down-regulation of genes involved in DNA replication, cell cycle, and cell division (Figure 5d). About 140 genes within these categories had a >1.5-fold decrease. Down-regulation of DNA replication and cell-cycle genes ultimately leads to arrested cell division and cessation of ear growth. A similar response was observed in maize leaves whereby both cell division and cell expansion were down-regulated under DRT in concert with a decrease in cell cycle gene expression (Avramova et al., 2015). In addition, down-regulation of cell cycle genes was observed in drought stressed ovules, which is thought to be a cause of early kernel abortion (Kakumanu et al., 2012).
At the V12 stage, the greatest number of up-regulated genes was in the "oxylipin biosynthetic process", "lipid biosynthetic process", and "response to abscisic acid stimulus" categories. Activation of ABA-signaling pathways was expected, but the up-regulation of the oxylipin pathway appears to be novel. The most studied oxylipin is the plant defense hormone jasmonate, with a broad range of functions in growth, development, biotic (Wasternack & Hause, 2013) and abiotic stress responses (Savchenko, Zastrijnaja, et al., 2014).
Upstream genes in the biosynthetic oxylipin pathway AOS, LOX1, LOX3, LOX6 and two lipoxygenase genes with no specific names (GRMZM2G015419, GRMZM2G009479) were up-regulated. Similar drought induced expression of LOX genes has been observed in the leaf elongation zone (Avramova et al., 2015). Induction of these genes by drought may lead to an accumulation of 12-oxo-phytodienoic acid (12-OPDA), a precursor of jasmonic acid as was shown in Arabidopsis (Savchenko, Kolla, et al., 2014). The up-regulation of three maize homologs of the COI genes, which encode the jasmonate receptor that plays a key role in JA-signaling (Wasternack & Hause, 2013), were also found. Jasmonates are known repressors of cell-cycle genes especially in dividing tissues (Pauwels, Inze, & Goossens, 2009;Pauwels et al., 2008;Shyu & Brutnell, 2015). The activation of jasmonate biosynthesis and signaling pathways in developing ears may lead to repression of DNA replication, cell cycle, and cell division genes resulting in ear growth inhibition. We propose that jasmonates (or their precursors) may be the factors controlling ear growth under abiotic stress conditions which is supported by the emerging function of jasmonates in abiotic stress responses beyond their wellknown roles in biotic stress responses (Kazan,2015;Liu et al., 2015).

| Spikelet initiation is tolerant to stress which correlates with the stable expression of inflorescence meristem genes
Robust spikelet initiation was observed in ears under abiotic stress conditions. In Woodland, there was no difference between treatments in the final number of spikelets formed on the ear. While in Johnston, the spikelet number was modestly reduced under drought as compared to a dramatic reduction in ear length. Moreover, a study of nine Pioneer inbred lines showed a limited effect of drought on spikelet initiation and final number (Figure 2h). Spikelet number is determined by local auxin biosynthesis and auxin signaling in the inflorescence meristems that produces SPM and SM Gallavotti, 2013). There are at least six maize genes with well-established functions in spikelet initiation: BAF1 is required for ear formation (Gallavotti et al., 2011), VT2 and SPI1 function in local auxin biosynthesis (Gallavotti, Barazesh, et al., 2008b;Phillips et al., 2011), and ZmPIN2a, BIF2, and BA1 regulate auxin transport and signaling (Barazesh, Nowbakht, & McSteen, 2009;Gallavotti, Yang, et al., 2008;Gallavotti et al., 2004). Our data show that none of these genes are repressed in ears under drought stress (Figure 7). This finding is consistent with the limited effect of stress on the spikelet number (Table 4). Additionally, no significant expression differences were detected in genes controlling the inflorescence meristem size and maintenance (KN1, FEA2, 4, TD1), or meristem determinacy (RA1, 2, 3, BD1) (

| Down-regulation of floral development genes is consistent with the formation of defective ovaries under stress
Each spikelet (axillary) meristem (SM) gives rise to a floral meristem (FM) which develops into the gynoecium and terminates in formation of the ovary. The ovary is composed of the embryo sac and silk, a stigmatic structure which accepts pollen and facilitates fertilization. Floral development is mostly governed by the MIKC-MADS box transcription factors (Coen & Meyerowitz, 1991). There are 43 known MIKC MADS genes in the maize genome (Zhao et al., 2010) and 16 of them were found to be down-regulated in stressed ears, which is about 40% of all maize MADS box genes (Table 4). MADS proteins work in quadruple complexes and their combination defines specification of the floral organs (Coen & Meyerowitz, 1991). To some extent this model is pertinent for flower morphogenesis in grasses as well (Ciaffi, Paolacci, Tanzarella, & Porceddu, 2011). Stress sensitive MADS genes represent all functional A, B, C, D, and E classes (Table 4). MADS box genes that encode proteins that form the heterodimers ZAG1 and ZAG3 (BDE1) (Thompson et al., 2009) were both down-regulated under stress. Interestingly, ZMM16 which was downregulated by stress in our study (Table 4) was identified as a drought QTL in tropical maize suggesting a putative role in stress adaptation (Almeida et al., 2014). In addition, ZMM16 was also identified as a potential QTL candidate for accumulation of the ABA metabolite phaseic acid in ears (Setter et al., 2011). Non-MADS floral genes APETALA2like transcription factors (TS6, tasselseed6) and DL1 (drooping leaves), homologs of the rice YABBY transcription factor, were also downregulated in stressed ears.
The key gene for embryo sac development, IG1 (indeterminate gametophyte) (Evans, 2007), was also down-regulated by stress. This finding sheds light on the observation that even short periods of water limitation can lead to abnormalities in embryo sac development (Moss & Downey, 1971). The down-regulation of the FERONIA homolog is intriguing. In Arabidopsis a receptor-like kinase FERONIA is required for the pollen tube interaction with the synergids (Escobar-Restrepo et al., 2007). The function of the maize homolog has yet to be shown, but we can speculate that down-regulation of the FERONIA homolog might impair efficacy of fertilization under abiotic stress conditions. Collectively, many genes involved in gynoecium development are suppressed under drought stress, which ultimately leads to the formation of defective ovaries and a failure to produce kernels.
Due to acropetal ear development, this process is most prominent at the ear tip where the youngest florets are positioned.

| CON CLUS ION
The monoecious nature of maize allowed us to test and identify differential effects of abiotic stress on three key plant organs: leaf, ear, and tassel. Our results indicate that each organ perceives different levels of stress which in turn drive a differential growth and developmental response. The tassel displayed the lowest transcriptional response to stress relative to the other organs and thus seems to be relatively well insulated from abiotic stress. This may be related to the evolutionary strategy prioritizing pollen which can be dispersed over wide distance over female organ development, which are immobilized in the drought-stressed environment. In contrast the ear has distinct and independent responses to drought in that organogenesis (spikelet initiation) appeared to be relatively stable under drought stress while organ extension (ear length) was significantly altered. The reduction in ear growth under drought stress was associated with the down-regulation of gene expression from pathways involved in DNA replication, cell-cycle, and cell division.
This is consistent with published observations that inhibition of cell division is a common response to drought stress in plant organs such as silks (Fuad-Hassan et al., 2008), ovaries (Kakumanu et al., 2012), endosperm (Setter & Flannigan, 2001), and leaf meristem (Avramova et al., 2015). Our results also suggest that jasmonates may be mediators of cell cycle suppression (Kazan, 2015;Pauwels et al., 2008), but this hypothesis would need to be tested.
It is important to emphasize that the B73 line used in this study was released in 1972 and represents "older genetics" that is more susceptible to abiotic stress. In recent years advanced drought tolerant germplasm has been created; for example, Optimum ® AQUAmax ™ hybrids (Gaffney et al., 2015). It would be interesting to contrast the phenotypic and transcriptional changes in the parents of these modern drought tolerant hybrids with that of B73 to determine which attributes have been altered by breeding for abiotic stress tolerance.

ACK N OWLED G M ENTS
Authors are thankful to Jeff Habben, Jeff Schussler, Nate Coles, Carl Simmons,Antoni Rafalski, and Mike Muszynski for valuable suggestions and stimulating discussion; to Renee Lafitte and Salim Hakimi for their excellent field management; Lila Bell and Miriam Bell for tissue sampling; and Jesse Ourada for help with qRT-PCR.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest associated with the work described in this manuscript.