Transcriptome analysis of maize pollen grains under drought stress during flowering

Keywords: anthesis, differentially expressed genes, drought stress, pollen development, pollen vitality, reporductive tissue, Zea mays L.

Maize (Zea mays L.) is one of the most important grain crops underpinning global food security, due to its wide cultivation and high productivity (Liu et al. 2020). However, maize growth, development and grain yield are susceptible to drought stress (Lobell et al. 2014; Danilevskaya et al. 2019). Major developmental processes including seedling emergence, canopy establishment, pollen and silk development, and grain-filling are all affected by drought stress (Zheng et al. 2010; Song et al. 2017; Wang et al. 2019; Jia et al. 2020). Maize flowering, marked by tassel emergence and anthesis, is extremely sensitive to drought stress (Andrade et al. 1999; Song et al. 2017). The anther is a unique organ for male gametogenesis, and its abnormal development may cause male sterility under drought stress (Geng et al. 2018; Li et al. 2021).

In order to cope with abiotic stresses, plants have evolved many strategies ranging from changes in starch and sucrose metabolism to hormone signal transduction and gene expression regulation (Lenka et al. 2011; Min et al. 2016; AbdElgawad et al. 2020). Genes are up- or downregulated in response to drought stress, and considerable differences in the expression patterns of many genes are observed through gene expression profiling studies (Kakumanu et al. 2012; Song et al. 2017; Jia et al. 2020). Functions of stress-responsive gene products such as dehydration-responsive element binding proteins (DREB), heat stress transcription factors (HSF), myeloblastosis (MYB) and serine/threonine protein kinases (SnRKs) have been well characterised (Qin et al. 2007; Bechtold et al. 2013; Baldoni et al. 2015; Zenda et al. 2019). Drought stress in maize during flowering results in pollen mortality and ovule abortion through lower carbon availability (Andrade et al. 1999; Aslam et al. 2015; Shin et al. 2015). As such, the quality of pollen grains is recommended as a selective trait for drought tolerance and is used as an important index for the success of reproductive stages (Aslam et al. 2015; Li et al. 2021).

Whole transcriptome sequencing (RNA-Seq) is a method for gene expression analysis with high sensitivity and high throughput (Jain 2011; Portwood et al. 2019). It has been widely employed in maize research (Davidson et al. 2011; Kakumanu et al. 2012; Kim et al. 2021) with increasing availability of the maize reference genome and considerable gene expression information (Danilevskaya et al. 2019; Portwood et al. 2019; Liu et al. 2020). Studies relating gene expression to water deficiency in maize have been performed with roots (Poroyko et al. 2007), seedlings (Humbert et al. 2013), and developing ears and tassels across multiple stages (Miao et al. 2017; Danilevskaya et al. 2019). The differentially expressed genes (DEGs) identified from specific tissues were found useful in predicting differences in tolerance between maize cultivars (Wang et al. 2019; Zenda et al. 2019; Liu et al. 2020). Although transcriptome profiling has been widely used to identify drought-responsive genes and characterise molecular mechanisms in maize (Kakumanu et al. 2012; Song et al. 2017; Danilevskaya et al. 2019; Wang et al. 2019; Jia et al. 2020), the effects of drought stress on mature pollens during flowering require more study.

In order to enhance our understanding of the molecular responses to drought stress during maize flowering, we used RNA-Seq to investigate transcriptional expression in pollen grains of plants subjected to well-watered (WW) and drought-stressed (DS) conditions at the flowering stage. A maize hybrid sensitive to flowering-stage drought stress was chosen for the study. We sought to identify drought-responsive genes known or predicted to be involved in different metabolic and regulatory pathways. This study should provide valuable information on the functional genes and key regulatory pathways controlling the maize pollen grain response to drought stress during flowering, which could then benefit molecular breeding for flowering drought tolerance.

A pot trial was conducted during summer 2017 at the experiment field of Anhui Agricultural University campus, China (117°24′E, 31°83′N). The site is in a subtropical humid monsoon climate zone. Three seeds of maize hybrid ZhongDan909 (ZD909), which is widely grown in the local area, were sown into each of 30 plastic pots (25 cm diameter and 30 cm height) filled with sieved yellow brown loam soil. Soil composition was as follows: organic matter 10.6 g kg⁻¹, total nitrogen 1.13 g kg⁻¹, available nitrogen 81.5 mg kg⁻¹, available phosphorus 33.1 mg kg⁻¹, available potassium 76.2 mg kg⁻¹, pH 7.15, and field capacity 27.32%. A single seedling of the three seedlings at the third-leaf stage was retained for culturing. Pots were arranged randomly. Irrigation was manually provided to maintain soil water content at field capacity, and pots were moved to a shelter in the field when rain was predicted. Otherwise the plants were grown in the open. At tassel emergence, half of the plants were retained at field capacity as for previous stages (WW), and the rest were subjected to drought stress at 50% field capacity (DS). The drought treatment lasted for 7 days, from 3 days before to 4 days after anthesis. Then, all remaining plants were watered normally until they were fully mature and harvested, and the kernel number per ear were calculated.

Samples were collected whilst pollen shedding at 09:00–10:00 immediately after the 7 days of treatment of the plants in well-watered and drought stress treatments respectively. Each treatment was with three replicates at different pots. The samples were respectively obtained from three treatments. The collected pollen grains were taken to the laboratory at 25°C and observed under a microscope (BX53+DP73+cellSens; Olympus, Japan). Three microscopic fields were used to estimate average number of abnormal pollen. The collected pollen grains were placed into pollen germination solution in the dark at 25°C for 4 h, then observed under a microscope in order to count pollen tube germination rate from three microscopic fields. The pollen tube germination medium consisted of the following (L⁻¹): MES 3.99 g, PEG4000 240 g, Ca(NO₃)₂.2H₂O 0.7 g, MgSO₄.7H₂O 0.2 g, KNO₃ 0.1 g, and H₃BO₃ 0.1 g.

The samples of pollen grains for transcriptome analysis were collected, and then immediately frozen in liquid nitrogen and stored in a freezer at −80°C.

Total RNA was isolated using TRIzol reagent [Invitrogen (Thermo Fisher Scientific), Waltham, MA, USA] according to a protocol provided by the manufacturer. The total amount of RNA was analysed with a Bioanalyser 2100 and the purity was assessed by the RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, CA, USA) with RNA integrity number (RIN) >7.0. Approximately 10 μg RNA was used to isolate poly(A) mRNA with poly-T oligo-attached magnetic beads (Invitrogen). Following purification, mRNA was fragmented into smaller pieces using divalent cations under increasing temperature. The cleaved RNA fragment was then reverse-transcribed to construct the final cDNA library in accordance with the protocol for the mRNA-Seq sample Reverse Transcription Kit (Illumina, San Diego, CA, USA). The mean size for the paired-end libraries was 300 ± 50 bp. Paired-end sequencing was performed with an Illumina Hiseq 4000 at LC Sciences (Houston, TX, USA) by following the protocol recommended by the provider.

Six pollen-grain RNA samples from three different plants of both WW and DS treatments were used for transcriptome sequencing analysis. Before assembly, low-quality reads (reads containing sequencing adaptors or sequencing primer, nucleotides with quality score (Q) <20) were removed. In total, 299 769 018 clean reads were produced, and the reads from WW and DS plants were mapped to the UCSC maize reference genome (http://genome.ucsc.edu/) using the TopHat package. First, TopHat was used to remove a portion of reads based on quality information accompanying each read, and then the reads were mapped to the reference genome (ftp://ftp.ensemblgenomes.org/pub/plants/release-34/fasta/zea_mays/dna/Zea_mays.AGPv4.dna.toplevel.fa.gz Zea mays L.). TopHat allows multiple alignments (up to 20 by default) and a maximum of two mismatches when mapping the reads to the reference. It then builds a database of potential splice junctions, confirmed by comparing the previously unmapped reads against the database of putative junctions (Trapnell et al. 2010).

The mapped reads for each sample were assembled using StringTie, and then all the transcripts from samples were merged in reconstructing a comprehensive transcript list using Perl script. After the final transcript was assembled, StringTie and Ballgown were used to estimate the expression levels of all transcripts. StringTie was used to show the expression level for cDNAs (identified from Phytozome 12) by calculating fragments per kilobase transcript per million mapped reads (FPKM) values (Trapnell et al. 2010; Pertea et al. 2015). DEGs were selected with log₂(fold-change) ≥1 and with statistical significance (P-value) ≤0.05 by the Ballgown R package (Frazee et al. 2015). The DEG correlation network was visualised using the OmicStudio tools in R version 3.6.1 (https://www.omicstudio.cn/tool).

Gene Ontology (GO) (http://www.geneontology.org/) and functional annotations were performed for all identified DEGs using the Goatools 2.0 software (https://github.com/tang haibao/goatools) (P ≤ 0.05). Pathway enrichment analysis was performed for all identified DEGs in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/genes.html), using Blastx/Blastp 2.2.24+ and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) (Mao et al. 2005).

Total RNA was extracted using the TRIzol kit following the manufacturer’s protocols. Synthesis of cDNA by purified RNA samples was done with PrimeScript RT reagent kit with gDNA Eraser (Takara, Dalian, China). The primers of genes for qPCR were designed using Primer Premier 5.0 software (PREMIER Biosoft, Palo Alto, CA, USA). Sequences of gene primers are shown in Supplementary Table S1. After setting up the parameters, qPCR was performed in an ABI 7500 FAST Real-Time PCR System (Thermo Fisher Scientific), and three biological replicates were used for all qPCR samples. The WW samples served as the control, and gene expression levels were calculated by the 2^−(ΔΔCt) algorithm, then each value was log₂ transformed (Wang and Wang 2021).

Under drought stress, the time of anthesis was postponed by 3–4 days due to anther indehiscence (Fig. 1a, b), and there was a lower rate of grain formation (Fig. 1c, d). Final kernel number per ear averaged 439.3 for WW plants but only 236.3 for DS plants, a 46.2% reduction (Fig. 1i). In the DS treatment, 93.7% of pollen grains were abnormal, compared with only 22.0% in the WW treatment (Fig. 1e, f, j). Furthermore, only 7.3% of DS pollen grains germinated compared with 69.3% of WW pollen grains (Fig. 1g, h, k).

Fig. 1.  Comparison of reproductive traits under well-watered (WW) control and drought-stress (DS) conditions for maize hybrid ZD909: (a) WW and (b) DS tassel; (c) WW and (d) DS mature ear; (e) WW and (f) DS pollen grains; (g) WW and (h) DS pollen tube germination; (i) kernel number per ear; (j) rate of abnormality; (k) germination. Bars = 40 μm in e–h. Note: tube length twice long as the diameter recommended as active pollen grains.

In total, 301.89 million raw reads were obtained and then trimmed, and 150.46 million clean reads for the DS sample (ZD909G) and 149.31 million clean reads for the WW sample (ZD909Z) were achieved (Table 1). All of the clean reads were matched to the maize reference genome, allowing two base mismatches. As a result, 132 884 619 mapped reads for the DS sample and 133 330 329 mapped reads for the WW sample were obtained, with an average matching rate of 88.8% (Table 1).

Table 1.  Number of reads sequenced and mapped to Zea mays ZD909 genome.

Principal component analysis (PCA) based on the transcriptome analysis exhibited an obvious separation among DS and WW pollen grains (Fig. 2a). The transcriptional variation that occurred in response to drought stress was characterised by using DEG analysis, with the criteria of P ≤ 0.05 and log₂(fold-change) ≥1 used to identify genes that were differentially expressed between DS and WW pollen grains. In total, 308 DEGs were identified (Supplementary material File S1), among which, 130 DEGs were upregulated, and the other 178 were downregulated in DS pollen (Fig. 2b, File S2). Therefore, the results indicate that the expression of most pollen grain genes (57.8%) was significantly inhibited by drought stress. A heatmap was constructed based on the log₁₀(FPKM + 1) values for the 43 DEGs related to pollen grain development for high or low expression levels, with similar distribution patterns verified after clustering analysis (Fig. 3).

Fig. 2.  Analyses of differentially expressed genes (DEGs) between ZD909G (DS) and ZD909Z (WW). (a) PCA with principal components 1 and 2: red points represent DS samples, and blue points WW samples. (b) Volcano diagram illustrating the fold-changes and P-values for genes (log₂(fold-changes) ≥1 and P ≤ 0.05); red points represent upregulated genes for DS samples, blue points downregulated genes, and grey points no DEGs.

Fig. 3.  Heatmap of the overlap DEGs between ZD909G (DS) and ZD909Z (WW), plotted using the heatmap.2 function of the R/Bioconductor package. Hierarchical clustering of the DEGs was done by the complete method with Euclidean distance. Gene expression levels were transformed by log₁₀(FPKM + 1) and the values were centred and scaled in the row direction. The x-axis indicates samples; y-axis indicates DEG names.

Based on sequence homology, 308 DEGs were annotated into the top 50 functional subcategories, including 25 biological processes, 15 cellular components and 10 molecular functions (Fig. 4, File S3). Among the biological process subcategories, ‘biological process’ was the main functional group, followed by ‘protein phosphorylation’ and ‘regulation of transcription, DNA-templated’. Among the cellular component subcategories, ‘nucleus’ was the main functional group, followed by ‘integral component of membrane’ and ‘plasma membrane’. Among the molecular function subcategories, ‘protein binding’ was the main functional group, followed by ‘molecular function’ and ‘ATP binding’.

Fig. 4.  GO annotation analysis of DEGs between ZD909G (DS) and ZD909Z (WW). The x-axis indicates subcategories; y-axis indicates number of genes in a category.

Transcriptional variations that occur in response to drought stress were further characterised by using the DEGs of the six transcript libraries to discover the genes with significant differences in expression via KEGG enrichment analysis. Scatter plots of significant enrichment items of KEGG pathway analyses associated with DEGs between ZD909G (DS) and ZD909Z (WW) are shown in Fig. 5. The analysis showed that ‘starch and sucrose metabolism’, ‘plant hormone signal transduction’ and ‘RNA transport’ were the most prominent of the 20 most enriched pathways (Fig. 5, File S4).

Fig. 5.  Enrichment analysis of DEGs between ZD909G DS and ZD909Z WW. The size of a node is proportional to the number of targets in the GO category. The colour of the node represents the significance of DEG enrichment: higher significance is represented by a deeper colour.

In order to identify the DEGs involved in the metabolic pathways potentially related to pollen grains, further analysis was performed using the KEGG pathway database. In total, 44 DEGs were assigned to nine KEGG pathways (Table 2), which showed 11 genes for the most represented pathway ‘starch and sucrose metabolism’, followed by eight genes for ‘plant hormone signal transduction’ and seven genes for ‘RNA transport’. Data from Table 2 depicting nine pathways and 44 DEGs are described in further detail in the following sections.

Table 2.  KEGG pathways enriched with DEGs between ZD909G (DS) and ZD909Z (WW).

Eleven DEGs were found to be involved in starch and sucrose metabolism (Table 2). The main physiological function of starch and sucrose metabolism is to provide energy and carbon building blocks and sugar signals (Wu et al. 2013; Liao et al. 2020; Li et al. 2021). Of the 11 DEGS, seven were related to cell wall integrity, and of these seven genes, five were upregulated, including three glucan synthase genes involved in callose synthesis (Chen and Kim 2009; Seale 2020) and two encoding pectin esterase. Among the downregulated genes under drought stress, a key gene was identified encoding polygalacturonase. A gene encoding a protein kinase, likely having a regulatory role, was also downregulated, as was a gene for a pentatricopeptide repeat-containing protein implicated in mitochondrial biogenesis (Lurin et al. 2004). Four genes were classified under pentose and glucuronate interconversions, two of which, encoding pectinesterases, were upregulated (Table 2).

Eight DEGs were involved in plant hormone signal transduction (Table 2). Among them, two were upregulated; these encoded carboxylesterases with a role in pollen tube development (Iaria et al. 2016). Of the six downregulated genes, two encoded proteins associated with auxin and ethylene signal transduction. Auxin response factor 18 (ARF18) is involved in auxin signalling (Xing et al. 2011), and the EIN3-binding F-box protein is involved in ethylene signalling (Gagne et al. 2004). Two of the other downregulated genes were for transcription factors and two were potentially involved in signal transduction. Of the latter two, the gene for protein serine/threonine kinase (Zm00001d051615) was also listed under ‘starch and sucrose metabolism’. The transcription factor bHLH13 (gene Zm00001d009013) has been shown to negatively regulate jasmonic acid-mediated defence and development in Arabidopsis (Song et al. 2013).

Four DEGs were associated with ubiquinone and other terpenoid-quinones and all were downregulated; three of them are linked to redox activity. Two of the genes are known to code for 1,4-dihydroxy-2-naphthoy1-CoA synthase, which is involved in vitamin K1 synthesis in the peroxisome (Babujee et al. 2010). Another unspecified gene is involved in ubiquinone biosynthesis, which is an important component of the electron transport chain of plant mitochondria (Rose and Sheahan 2012).

There were four DEGs responsive to light stimulus by KEGG pathways enrichment analysis, of which two were downregulated and two were upregulated. Genes in this category are likely to be involved in determining flowering time and responding to drought stress (Song et al. 2017). The two upregulated genes may be responsible for oxidoreductase activity (Table 2).

Seven DEGs were categorised under RNA transport pathways, with five upregulated and two downregulated. Four of the upregulated genes are associated with translation activity, and the other codes for the transcription factor NF-YC2. This was the only upregulated transcription factor as shown in Table 2. One of the downregulated genes is for translational initiation and the other for response to viruses.

Five DEGs related to nucleotide sugar metabolism were found, and among them there were two upregulated genes: one associated with redox homeostasis (DIHYDROLIPOYL DEHYDROGENASE 1) and the other with translation (Table 2). There were three downregulated genes, one of which is related to protein kinase activity, one to cellulase inhibition, and the third to carbohydrate metabolic processes, likely to be involved in pectin biosynthesis (Borg et al. 2021).

Four DEGs were associated with protein export and ABC transporters (Table 2). The MADS-box transcription factor is in the protein export pathway and is downregulated. The MADS-box family has been implicated in several reproductive development processes (Masiero et al. 2011). One ABC transporter was upregulated and two were downregulated. ABC transporters are involved in pollen-wall biosynthesis and are essential for pollen fertility (Somaratne et al. 2017; Luo et al. 2020).

The responsive gene enrichment analysis revealed the existence of a considerable number of pathways among the DEGs, related to pollen grain development under drought stress (Fig. 6). The genes Zm00001d024481 (response to light stimulus) and Zm00001d033505 (nucleotide sugar metabolism) were connected with each other. The gene Zm00001d045078 (nucleotide sugar metabolism) was connected to genes Zm00001d033637 (ubiquinone and other terpenoid-quinone biosynthesis), Zm00001d025253, Zm00001d036151 (RNA transport) and Zm00001d039248 (Fig. 6). Furthermore, the two genes Zm00001d024481 and Zm00001d033505 and the three genes Zm00001d014866, Zm00001d014865 and Zm00001d023351 were all connected to each other. Overall, the correlation network and KEGG enrichment showed that the DEGs Zm00001d033637, Zm00001d024481, Zm00001d033505, Zm00001d025253 and Zm00001d045078 were all down-regulated in the maize pollen grains developed under drought stress. In particular, our data suggest that these DEGs and the KEGG pathways in which they were involved had a degree of correlation. Therefore, our analysis provides a foundation for further studies to identify the function of these DEGs related to drought stress in the development of maize pollen grains.

Fig. 6.  Correlation network of DEGs between ZD909G (DS) and ZD909Z (WW). The size of a node is proportional to the number of DEGs in the GO category. The colour of the node represents the significance of DEG enrichment. The thickness of the line represents the intensity of the correlation between DEGs.

The accuracy of RNA-Seq results was confirmed by analysing 22 DEGs related to drought stress by qPCR for the same total RNA sample as used in RNA-Seq (Fig. 7a). Gene specific qPCR primers and gene names are listed in Table S1. Expression patterns for qPCR were consistent with the RNA-Seq expression, except for three genes (Zm00001d031286, Zm00001d003431 and Zm00001d039212). The rate of coincidence between RNA-Seq expression and qPCR expression was 86.4%. Coincidence analysis between RNA-Seq and qPCR results showed consistency (fitted by a linear regression equation y = 0.4352x − 0.1143, R² = 0.6608), indicating the reliability of RNA-Seq expression profile in this study (Fig. 7b).

Fig. 7.  DEGs confirmed by qPCR using the same total RNA sample as used in RNA-Seq: (a) relative expression level of DEGs as found for each technique; (b) coincidence analysis of log₂(fold-change) obtained by both techniques. DEGs at false discovery rate ≤0.01 and log₂(fold-change) ≥1 from RNA-Seq.

Pollen molecular biology has received little attention in endeavours to understand the mechanisms of drought-stress tolerance (Zhang et al. 2018; Wang et al. 2019). In maize, transcriptomic studies of drought stress have been performed for a range of organs but not for pollen (Kakumanu et al. 2012; Song et al. 2017; Danilevskaya et al. 2019; Wang et al. 2019; Jia et al. 2020). Given the sensitivity of the flowering period to drought stress, it is important to understand the effect of drought stress during flowering on mature pollen (Andrade et al. 1999; Song et al. 2010). Therefore, we investigated and identified the DEGs involved in metabolic pathways in response to drought stress, using the KEGG pathway database and correlation network.

Transport of photoassimilate (mainly in the form of sucrose) from source leaves to sink organs is essential for plant development (Ruan et al. 2010; Zhu and Yang 2013). Sucrose must be degraded into hexoses or their derivatives for various metabolic and biosynthetic processes once it reaches sinks such as flowers, fruits, seeds and roots (Ruan et al. 2012; Yan et al. 2013). Anther and pollen development requires carbon nutrients (Shen et al. 2019), with sugars providing both nutrients and signals for anther and pollen development (Wu et al. 2013; Shen et al. 2019). Therefore, abnormal carbohydrate metabolism may lead to pollen abortion (Ruan et al. 2010; Yang et al. 2017).

If subjected to drought stress during pollination, maize leaf photosynthesis is inhibited (Schussler and Westgate 1995) and the genes related to sugar-processing enzymes such as sucrose synthase 7 (Song et al. 2017) are downregulated. The translocation of sucrose to ovaries and pollen grains is then reduced or even arrested (Ruan et al. 2010). However, if sucrose is injected into the stems of water-stressed maize plants, the physiological activities and metabolites are restored, which may prevent or alleviate grain abortion (Ruan et al. 2010).

In our study, the DEGs in the starch and sucrose metabolism related pathways suggest that drought affected the expression of genes related to the supply of energy and carbon skeletons in pollen grains, which may result in a shortage of energy required for pollen maturation and, consequently, abnormal development of pollen grains. The transcription of the pentatricopeptide repeat gene, important in mitochondrial biogenesis (Lurin et al. 2004), was downregulated. Together with the downregulation of a gene (Zm00001d033305) likely involved in ubiquinone biosynthesis of the electron transport chain (Rose and Sheahan 2012), the findings suggest impairment of energy production. The inhibition of expression of the exopolygalacturonase gene (Zm00001d053622) may affect the ability of the pollen to germinate and penetrate through the style (Allen and Lonsdale 1993; Holmes-Davis et al. 2005; Dai et al. 2007; Gong et al. 2015). Sugar nucleotide metabolism, important in providing precursors for cell wall synthesis (Borg et al. 2021), could also be impaired with the expression of the gene encoding UDP-glucuronate 4-epimerase (Zm00001d045078) strongly downregulated.

It is intriguing that five DEGs in starch and sucrose metabolism related pathways were also connected with other pathways. Among these five DEGs, four (i.e. Zm00001d034741, Zm00001d043725, Zm00001d003431 and Zm00001d025588) are involved in pentose and glucuronate interconversion pathways, and the other DEG (Zm00001d051615) also reportedly takes part in plant hormone signal transduction. As a protein kinase gene, Zm00001d051615 is likely to have an important signalling role in pollen tube growth. Surprisingly, five DEGs were upregulated: three were callose synthases and two were pectinesterases. Callose (Dong et al. 2005; Parre and Geitmann 2005; Chen and Kim 2009; Seale 2020) and pectin (Holmes-Davis et al. 2005; Cascallares et al. 2020) are cell wall components of the pollen grains and pollen tubes. The upregulation related to callose and pectinesterases likely affects the cell wall dynamics of the pollen grains and pollen tubes, as discussed below.

Eleven DEGs in this study were associated with pollen cell wall biosynthesis. The five downregulated genes encoded a cellulose synthase, a 1,4-alpha-glucan branching enzyme, polygalacturonase, an ABC transporter and a UDP-glucuronate 4-epimerase. Given the likely reduced carbon skeletons, this was not unexpected. However, six genes were upregulated, including three callose synthase genes, raising the question of why carbon was channelled into callose. It is known that callose is essential for properly sculptured pollen exine and pollen viability (Dong et al. 2005) and that stress stimulates callose biosynthesis (Chen and Kim 2009; Piršelová and Matušíková 2013). Perhaps the maintenance or increase in callose makes pollen more resistant to drought stress and facilitates germination and early pollen tube growth. Two pectin lyase genes were also upregulated, which may also influence pollen tube growth (Holmes-Davis et al. 2005) under stress conditions. Clearly, there are transcript modifications (up and down DEGs) supporting pollen wall changes that ultimately impaired pollen viability. The ABC transporter genes, important in pollen wall development (Luo et al. 2020), also reflect this, with upregulation of G-type and downregulation of C-type (see Table 2). The DEGs influencing pollen wall development lead to pollen abnormalities and reduced germination under drought stress, as seen in Fig. 1. Proper cell wall development in pollen grains guarantees plant sexual reproduction, and the majority of pollen grains without vitality have abnormal wall development (Zhu and Yang 2013).

Phytohormones play a key role in growth and development, as well as influencing response to abiotic stress (Weiss and Ori 2007; Chen et al. 2019). In our study, DEGs related to auxin and ethylene signalling were identified, and were downregulated. Auxin has been shown to influence pollen tube germination and growth (Wu et al. 2008; Gao et al. 2019), and impairment of auxin signalling affects pollen tube germination (see Fig. 1h). In Arabidopsis, ARF17 has been implicated in pollen wall patterning (Yang et al. 2013). Ethylene has also been implicated in pollen tube growth. In Arabidopsis, ethylene has been shown to promote pollen tube growth by affecting actin filament organisation via a cGMP-dependent pathway (Jia et al. 2018). Ethylene has also been implicated in drought-stress tolerance (Hong et al. 2017) and can reduce stress during pollen germination (Jegadeesan et al. 2018); therefore, there is a need for comprehensive understanding of ethylene signalling in the pollen grain. In this context, genes for two transcription factors (Zm00001d009013 coding bHLH13, Zm00001d039243 a receptor-like serine/threonine-protein) and a protein kinase (Zm00001d051615) were downregulated in response to drought stress. The transcription factor bHLH13 negatively regulates jasmonic acid-mediated defence and development in Arabidopsis (Song et al. 2013). This suggests that jasmonic acid, implicated in drought-stress responses in wheat (Ilyas et al. 2017), may be involved in pollen viability because the downregulation of Zm00001d009013 could enhance jasmonic acid responses. The maize transcription factor ZmbHLH13 is clearly of interest in relation to drought stress adaptation.

Five DEGs encoded transcription factors, of which four were downregulated, with only one upregulated DEG. The downregulated transcription factors ARF, the receptor-like serine/threonine-protein, bHLH13 and MADS-box 27 have been considered. The F-box protein 2 which binds the EIN3 transcription factor was also down-regulated. Transcription factor ZmNF-YC2 was strongly upregulated. NF-Y subunit B has been linked to drought stress tolerance in maize (Nelson et al. 2007); importantly, transgenic maize plants overexpressing ZmNF-YB2 showed drought tolerance under water-limiting condition. Transcription factor ZmNF-YC2 is clearly worthy of further investigation with its potential in transgenics and marker-assisted breading.

Redox homeostasis is important because reactive oxygen species (ROS) are an integral part of stress responses (Liu et al. 2019). ROS can act as a signal as well as causing oxidative damage at higher concentrations (Foyer and Noctor 2005; Mittler et al. 2011). Five DEGs related to redox activity were represented in three different pathways; two were downregulated and three were upregulated. The genes Zm00001d016211 and Zm00001d016212 encoding a vitamin K1 biosynthesis protein in the peroxisome were downregulated. Vitamin K1 is involved in plasma membrane redox activity (Bérczi and Møller 2000). Zm00001d009163, encoding the flavin enzyme dihydrolipoyl dehydrogenase, was upregulated. In plants, as in other organisms, it is likely to increase ROS production generated by dihydrolipoyl dehydrogenase-containing mitochondrial enzyme systems α-ketoglutarate (2-oxoglutarate) dehydrogenase and pyruvate dehydrogenase (Rigoulet et al. 2011; Timm et al. 2015). In terms of stress adaptation, there remains a need to link ROS and redox activity to other changes that have been discussed, notably in pollen and pollen tube wall dynamics, hormone signalling (auxin, ethylene and jasmonic acid) and key transcription factors that relate to the abnormal pollen and inhibited germination phenotypes under drought stress.

Drought stress during flowering reduces pollen grain viability and pollen tube germination, thereby lowering the fertilisation rate. In this study, biological functions and metabolism pathways of key DEGs for pollen grains were identified based on KEGG enrichment and correlation network analyses. The downregulation of polygalacturonase reduced cell wall lysis in early germination, increased callose synthase transcripts and decreased cellulase transcripts, which ultimately may inhibit pollen development and initiation of pollen tube growth. Modified auxin and ethylene metabolism and translation and reduced mitochondrial function could also impede the onset of pollen tube growth. Genes with potential adaptation roles were also altered in their transcript levels. These included genes encoding the upregulated transcription factor ZmNF-YC2 (Zm00001d028433) and the downregulated transcription factor ZmbHLH13 (Zm00001d009013). Downregulation of ZmbHLH13 could enhance jasmonic acid responses. The upregulated flavin enzyme gene DIHYDROLIPOYL DEHYDROGENASE 1 (Zm00001d009163), associated with increased ROS, was of potential interest in relating redox activity to stress adaption, given the fact that ROS can have a signalling role. In particular, the correlation network and KEGG enrichment analyses suggest that the DEGs and the KEGG pathways in which they were involved showed a degree of correlation. Therefore, findings from this study will assist in identifying potential molecular targets for improving drought tolerance of pollen grains, and provide a basis for screening traits during maize breeding for flowering drought tolerance.

Supplementary material is available online.

The data that support this study are available in the article and accompanying online supplementary material.

Youhong Song is an Associate Editor of Crop and Pasture Science but was blinded from the peer-review process for this paper.

This work has been financially supported by National Key R&D Program of China (No. 2017YFD0301307, 2017YFD0300408) for YS, and the Jiangsu Collaborative Innovation Centre for Modern Crop Production for JCL.

YZ, JCL and YS designed the experiment; YZ and SS performed the experiment; YZ and JL analysed the data; YZ, JL, SS, RR, YLR and YS wrote the manuscript. All authors read and approved the final manuscript.