The maize single-nucleus transcriptome comprehensively describes signaling networks governing movement and development of grass stomata

Abstract

The unique morphology of grass stomata enables rapid responses to environmental changes. Deciphering the basis for these responses is critical for improving food security. We have developed a planta platform of single-nucleus RNA-sequencing by combined fluorescence-activated nuclei flow sorting, and used it to identify cell types in mature and developing stomata from 33,098 nuclei of the maize epidermis-enriched tissues. Guard cells (GCs) and subsidiary cells (SCs) displayed differential expression of genes, besides those encoding transporters, involved in the abscisic acid, CO₂, Ca²⁺, starch metabolism, and blue light signaling pathways, implicating coordinated signal integration in speedy stomatal responses, and of genes affecting cell wall plasticity, implying a more sophisticated relationship between GCs and SCs in stomatal development and dumbbell-shaped guard cell formation. The trajectory of stomatal development identified in young tissues, and by comparison to the bulk RNA-seq data of the MUTE defective mutant in stomatal development, confirmed known features, and shed light on key participants in stomatal development. Our study provides a valuable, comprehensive, and fundamental foundation for further insights into grass stomatal function.

IN A NUTSHELL

Background: As the outmost layer of cells, the leaf epidermis controls gas and water exchange between the plant and its environment, and provides protection from various biotic and abiotic stresses. Within the epidermis, stomata determine photosynthetic and water-use efficiencies, responding rapidly to changing environmental conditions. Unlike dicot stomata, grass stomata contain two subsidiary cells (SCs) flanking two dumbbell-shaped guard cells (GCs). Coordinated regulation of SCs and GCs allows higher speeds of stomatal movement in response to environmental changes, and results in higher water-use efficiencies.

Question: Can we obtain a comprehensive understanding concerning the molecular basis underlying the fast stomatal response and stomatal development in grass? What are the key players in stomata movement and development, such as membrane channels, components of key signaling pathways, and GC or SC cell wall-specific proteins?

Findings: We firstly isolated propidium iodide-stained nuclei from maize epidermal peels by fluorescence-activated sorting, and developed a high-throughput platform of single nucleus RNA sequencing. Assisted by reporter lines involving transgenic maize and rice, we acquired the transcriptional profiles of the different cell types in maize peels, including GCs, SCs, epidermal cells, the attached mesophyll, and vascular-related cells. We found GC- and SC-specific genes encoding previously reported and new transporters, components of abscisic acid, CO₂, Ca²⁺, starch metabolism, and of blue light signaling pathways identified in other species. Using the same technical approach, we analyzed the maize leaf base sample, obtaining new candidates involved in grass stomatal development. Among these candidates, we identified cell wall-related genes that may play important roles in dumbbell-shaped morphogenesis of developing and mature stomata.

Next steps: Functional experiments are needed to test the hypotheses emerging from the single nucleus transcriptome data. Transcriptomic comparison of maize stomata with other monocot crops would provide deeper insights into the evolutionary innovation of the dumbbell-shaped stomata.

Introduction

Stomata, generally found in the outmost cellular layer of the aerial parts of plants, control gas and water exchange with the environment, thereby determining photosynthetic and water-use efficiencies (Hetherington and Woodward, 2003; Lawson and Blatt, 2014). In grasses, they incorporate two lateral subsidiary cells (SCs) flanking a pair of elongated, dumbbell-shaped guard cells (GCs), to form the unique “graminoid” morphology of the stomatal complex (Hepworth et al., 2018; Nunes et al., 2020). The coordinated regulation of SCs and GCs is perceived as an essential evolutionary innovation, allowing higher speeds of stomatal movement in response to rapid changes in the environment, leading to higher water-use efficiencies and considerably greater success in niche adaptation (Raissig et al., 2017; Nunes et al., 2020). In maize (Zea mays) and other monocots, genes encoding the transporters of cations, anions, abscisic acid (ABA), and sugar, as well as signaling pathway components, have been identified that sense and transmit a variety of environmental signals, including light, CO₂ concentration, and temperature, in order to regulate osmotic potential and turgor during stomatal movement (Kollist et al., 2014; Gray et al., 2020), Many of these genes are expressed in a GC- or SC-specific manner (Chen et al., 2017; Nunes et al., 2020). A complete understanding of the molecular basis underlying stomatal movement, particularly regarding how SCs cooperate with GCs to achieve rapid responses to environmental changes, has nevertheless remained elusive.

Stomata also provide an excellent model system to study development, and key monocot players are beginning to emerge. These include the closely-related basic helix–loop–helix (bHLH) domain transcription factors SPEECHLESS (SPCH) (Raissig et al., 2016), MUTE (Raissig et al., 2017; Wang et al., 2019a), and FAMA (Liu et al., 2009; Wu et al., 2019), other important genes regulating Stomatal Density and Distribution1 (SDD1), Epidermal Patterning Factor 2 (EPF2), EPF9, EPF1, Pangloss1 (PAN1), and Actin-Related Proteins 2/3 (ARP2/3), which function at different stages of stomatal development (Cartwright et al., 2009; Facette et al., 2015; Apostolakos et al., 2018; Hepworth et al., 2018; McKown and Bergmann, 2020; Yu et al., 2020). Although these prior studies provide an excellent basis for our understanding of important events in stomatal development, the way in which these cell states are progressively regulated along the developmental trajectory is largely unknown.

A systematic and thorough identification of the compendium of different genes that are expressed during graminoid stomatal movement, and of the inventory of gene products accumulated at the various developmental stages, would greatly advance our understanding of this important plant feature. A comprehensive approach to analyze transcriptional regulation requires quantitative examination of the contributions of the different cell types present in the tissue samples. For plants, this can involve the prior preparation of protoplasts, but this approach involves unavoidable perturbation of the biological system associated with enzymatic removal of the cell wall under hypertonic conditions. One protoplast production, thereby avoiding perturbation of the studied biological system whilst providing a comprehensive sampling of its cells (Liang et al., 2019), is the use of isolated nuclei as the source of transcriptional information as a surrogate for the cell (Macas et al., 1998; Zhang et al., 2008; Grindberg et al., 2013). Single-nucleus RNA sequencing (snRNA-seq) has subsequently been used in the analysis of fresh and preserved samples of tissues and organs, and has proved to be a powerful and popular platform for discovering the priority of critical regulators in animals (Liang et al., 2019; Kalish et al., 2020). So far, snRNA-seq has been successfully applied in plants to dissect the cell atlas in Arabidopsis thaliana roots, seedlings, buds, and other tissues (Farmer et al., 2021; Sunaga-Franze et al., 2021).

In this report, we described and validated snRNA-seq as an independent platform to characterize a cell atlas for maize leaf peels, following isolation of propidium iodide (PI)-stained nuclei by fluorescence-activated sorting (FAS). Based on molecular, genetic, and functional evidence, we further identified signature regulators of different cell types in peels, revealing the molecular basis of stomatal movement by analysis of GC- and SC-specific marker genes. These data allowed us to map onto GCs and SCs the specific signaling components/pathways related to photosynthesis, turgor regulation, ion transport, and ABA signaling. Based on these observations, we propose a dynamic model whereby stomata respond quickly and efficiently to intrinsic ABA and environmental signals via coordination between GCs and SCs.

Results

Data generation, clustering, and assessment of cell types in maize mature peels

Generally, peels consist of five cell types in characteristic proportions: a dominant majority of pavement cells (PCs), equal representation of GCs and SCs but at lower numbers, and trace representation of mesophyll cells (MCs) and of bundle sheath-related cells (BSCs) (Supplemental Figure S1A). To capture these cell types from the epidermal peels that represent an experimentally challenging tissue for protoplast production, we focused on the use of nuclei as the source of information concerning cellular transcription. To adapt this procedure for the maize epidermis, we optimized the production of nuclear suspensions with regards to epidermal excision, the homogenization procedure, the required amounts of tissue, and the choice of buffer. We further optimized the method of purification of the nuclei via flow cytometry and sorting, in terms of the sorting rate, the sample concentration, and the sample and sheath liquid composition, to enable fast sample processing, whilst maintaining the highest degree of nuclear integrity as well as compatibility with the 10x Genomics Chromium System (Supplemental Figure S1, B–D; Galbraith et al., 2021; Galbraith and Sun, 2021). Three abaxial mature peel samples were processed (samples Pv2, Pv3a, and Pv3b), using versions 2 and 3 of the 3' single-cell RNA-seq kit from 10x Genomics. This generated a total of 24,131 nuclei: 4,999 nuclei in sample Pv2, 10,229 in Pv3a, and 8,903 in Pv3b (Supplemental Table S1). The corresponding median numbers of genes detected per nucleus were 552 (sample Pv2), 1,332 (Pv3a), and 1,363 (Pv3b), consistent with the manufacturer’s claim that version 3 is optimized for trace mRNA discovery (Supplemental Table S1). The total numbers of genes detected in samples Pv2, Pv3a, and Pv3b were 21,406, 25,034, and 25,788, respectively accounting for 54.42%, 63.66%, and 65.58% of the total predicted genes in the maize genome (Supplemental Table S1).

To identify the contributions specific to different cell types, we first integrated and clustered the results from three peel samples using Seurat (Butler et al., 2018; Stuart et al., 2019) after assessing the organelle-encoded genes and the removal of batch effects (Supplemental Figure S2, A and B). Dimensionality reduction visualization was performed using Uniform Manifold Approximation and Projection (UMAP; Becht et al., 2019). The clusters were further re-grouped by analyzing the population relationships under different resolution values, and by examining the presence of known marker genes from different cell types, such as Mildew resistance Locus O for PCs (Acevedo-Garcia et al., 2014), C4 characterized genes for MCs and BSCs (Matsuoka et al., 2001), the SC-specific sugar transporter Closed Stomata1 (CST1; Wang et al., 2019b) and Z. mays K⁺ channel 1 (ZMK1; Büchsenschütz et al., 2005), and the GC-specific KAT1-like Shaker potassium channel 2 (KZM2; Gao et al., 2017; Figure 1, A and B). Surprisingly, when inspecting the structural and regulatory genes involved in anthocyanin biosynthesis, we identified one sub-population of PCs expressing a complete set of these genes (Supplemental Data Set 1; Supplemental Figure S3).

Figure 1

Cell type identification in snRNA-seq data of the maize peels. A, Visualization of nuclear classification after integrating three maize peel replicates using UMAP. Dots and colors represent individual nuclei and different cell types. PCs-N, PCs without active anthocyanin biosynthesis; PCs-A, PCs with active anthocyanin biosynthesis. B, Dot plot of known marker genes used for cell-type identification. The proportion and average expression levels are, respectively, denoted by the sizes and colors of the circles. The full names and gene IDs are provided in Supplemental Data Set 2. C and D, Expression of ZmGRP6 (Zm00001d045877) and ZmSCS1 (Zm00001d014325) in the epidermis. The expression patterns were analyzed by UMAP plots and promoter-driven ZmHTA6-YFP reporter lines. The YFP fluorescence signals of ZmGRP6 and ZmSCS1 are respectively specific to SCs in transgenic maize (C), and to GCs in transgenic rice (D). The cell outlines are revealed by staining with PI. The cluster assignment in the UMAP plots is shown in Figure 1A. The insets with boxed regions and dashed arrows show the GC cluster at higher resolution. Scale bar = 20 μm. The non-specific expression of ZmGRP6 in GC cluster and ZmSCS1 in SC cluster might be partially due to ambient RNA during nuclei purification.

Overall, we obtained six stable clusters from five cell types: PCs-A (PCs with active anthocyanin biosynthesis), PCs-N (PCs lacking active anthocyanin biosynthesis), MCs, BSCs, SCs, and GCs (Figure 1A). A strong correlation was found in the differentially expressed genes (DEGs) from the corresponding cell types by independent analyses (Supplemental Figure S4), indicating high reproducibility of the method. Similar proportions of cell types and UMI/gene detection numbers were also observed when comparing the three samples (Supplemental Figure S2, C–E), except for an absence of GCs and PCs-A in Pv2; this may be due to the relatively low gene detection number produced by the version 2 kit, and different lighting conditions in the growth environment, respectively. A total of 1,513 nuclei were obtained from SCs, and 44 from GCs. The low discovery rate of GCs might be due to their smaller cell size and more irregular dumbbell or elliptical nuclear shape (Brown and Johnson, 1962; Galatis, 1980). Distinct numbers of DEGs were predicted from GCs and SCs (289 and 1,186 genes, respectively), which might be due to the imbalance in numbers of nuclei from GCs and SCs, and result in our missing more of the GC-specific or -enriched genes. The DEG analysis predicted 251 marker genes for the GCs and 170 for the SCs, with only 20 genes being shared by the two cell types (Supplemental Data Sets 2 and 3).

For further confirmation of the specific clustering, the promoters of a SC-specific cell wall structure gene Glycine‐rich RNA‐binding Protein 6 (ZmGRP6) and a GC-specific SnRK2-interacting calcium sensor gene (ZmSCS1) were used to construct HTA6-eYFP-based reporters for transformation into maize and rice (Oryza sativa), respectively, with HISTONE H2A 6 (ZmHTA6) providing topogenic localization to the nucleus. Rice was included for transformation experiments based on its more efficient and rapid transformation system, its relatively close evolutionary relationship to maize, and the likelihood of conserved regulatory mechanisms of stomatal movement and development across monocots. Although the degree to which expression patterns in transgenic rice recapitulate those in transgenic maize has not been formally described in the scientific literature, the cell-type specificity of the YFP fluorescence within the transgenic plant epidermis did support the clustering results and the cell type assignments (Figure 1, C and D; Supplemental Figure S5, A and B).

Molecular mapping of GC osmoregulation related to photosynthesis, starch metabolism, and transporter activity in stomatal responses

GCs possess complex signal transduction networks and modified metabolic pathways for turgor modulation (Schroeder et al., 2001). Many documents emphasize the importance of GC starch degradation and of chloroplast activity and transporters in modulating stomatal closure and opening (Lawson and Blatt, 2014; Daloso et al., 2017). The way that GCs has comprehensively integrated the components/pathways remains an unresolved question. The cell atlas now acquired from maize peels provides an opportunity to develop a comprehensive model of the key interactions in the regulation of stomatal movement.

First, chloroplasts are a characteristic feature of GCs that are involved in osmotic regulation (Willmer and Fricker, 1996; Daloso et al., 2017). Three distinctive photosynthetic cell types are typically found in leaves: MCs, BSCs, and GCs (Lawson and Blatt, 2014). The functions of the GC chloroplasts are still debated, especially regarding what kind of roles are played by photosynthesis and starch metabolism in stomatal movement (Lawson and Blatt, 2014; Daloso et al., 2017). Our snRNA-seq study provides a detailed analysis of chloroplast functions in the three cell types (Figure 2A;  Supplemental Figure S6; Supplemental Data Set 4). In MCs, we identified an increased level of transcripts for marker genes functioning in the light reaction, including those encoding components of Light-Harvesting Complexes I and II (LHCI and LHCII), reaction centers in photosynthesis I and II (PSI and PSII), and linear and cyclic electron transport components. In contrast, only PSI-LHCI and cyclic electron transport transcripts were overrepresented in BSCs. Transcripts of genes involved in C4-associated cycles were observed in MCs and BSCs, respectively. β-Carbonic Anhydrase 1 (β-CA1), Phosphoenolpyruvatec Carboxylase (PEPC), and Pyruvate, pi Dikinase (PPDK) were highly expressed in MCs, whereas transcripts of genes encoding NADP-Malic Enzyme (NADP-ME) and Phosphoenolpyruvate Carboxykinase (PCK), the activities of which, respectively, mediate the release of CO₂ from malic acid and oxaloacetate, were mainly detected in BSCs. Transcripts of genes encoding Calvin cycle enzymes were represented only in BSCs.

Figure 2

Transcriptome dynamics and signaling pathway in GCs and SCs. A, Dot plot of photosynthesis-related genes within the different cell types, predominantly in MCs, BSCs, and GCs. B, Dot plots of known genes or orthologs involved in stomatal movement. The proportions and average expression values are, respectively, denoted by the sizes and colors of circles. The gene names are inferred from Arabidopsis and rice orthologs if the cell-type specific or enriched expression not reported in maize. C, Expression of GC and SC marker genes in the epidermis. The expression patterns were analyzed and displayed using UMAP plots, and by the images of the cellular locations of fluorescence from promoter-driven ZmHTA6-YFP reporters in transgenic rice (OST1, Zm00001d033339; ZIFL1, Zm00001d005002; CHX20, Zm00001d035631; CNGC15, Zm00001d017281; and PGX3, Zm00001d040725) and in transgenic maize and rice (ZmLecRLK1, Zm00001d046802). The cell outlines are revealed by staining with PI. The cluster assignment in the UMAP plots is shown in Figure 1A. The insets with boxed regions and dashed arrows show the GC cluster at higher resolution. Scale bar = 20 μm.

Compared with the situation in MCs and BSCs, GCs exhibit a relatively low level of expression of transcripts encoding components of the PSII-LHCII and PSI-LHCI complexes, and the initial carboxylation phase of the C4 pathway. Only a few of these were predicted as DEGs (Figure 2A;  Supplemental Data Set 4). Transcripts of ferredoxin–NADP⁺ oxidoreductase (FNR), involved in linear electron transport, showed low expression, whereas one component of NAD(P)H dehydrogenase (NDH)-dependent cyclic electron transport, photosynthetic NDH subunit of subcomplex B 4 (PnsB4), was highly expressed (Figure 2A), suggesting that, within GCs, the light reactions mainly produce ATP. These conditions might result in changes in the redox state of GC chloroplasts when the capacity for CO₂ assimilation in GCs is very low (Hampp et al., 1982; Reckmann et al., 1990; Zhang et al., 2001a, 2001b; Lawson and Blatt, 2014). In line with this, transcripts of one orthologous gene encoding plasma membrane-localized NADPH oxidase (NOX) of A. thaliana, which produces apoplastic reactive oxygen species and likely initiates extracellular signaling events (Xie et al., 2014), were found to be enriched in GCs (Supplemental Data Sets 3 and 6).

We then focused our attention on the analysis of transcript profiles associated with starch metabolism in GCs, based on the starch–sugar hypothesis predicted to play a role in the regulation of the osmotic potential of the cell during stomatal closure (Azoulay-Shemer et al., 2016; Lawson and Matthews, 2020). We found that transcripts of genes encoding proteins involved in starch metabolism were highly expressed (Figure 2A), including those encoding ADP-glucose pyrophosphorylase (AGPase), Starch Synthases 3 and 4 (SS3 and SS4), starch-branching enzyme, starch degradation glucan-water dikinase (GWD), and β-Amylase 2 (BAM2). Since expression of different suites of starch metabolism genes was observed in GCs and BSCs (transcripts of BAM3 were accumulated in BSCs, and GBE in GCs), this implies the physicochemical properties of the starch granules may differ between BSCs and GCs.

The mechanism through which osmoregulation signals control stomatal movement mainly involves regulation of transporters (Lawson and Matthews, 2020; Nunes et al., 2020). As expected, marker genes related to transporters/channels for turgor regulation linking GCs and SCs were captured in our snRNA-seq data. These included transcripts of genes or orthologous genes encoding K⁺ transporters KZM2(Gao et al., 2017) and High‐Affinity K⁺ Transporter 1 (HAK1; Qin et al., 2019), Ca²⁺ transporter Cyclic Nucleotide-Gated Channel 15 (CNGC15; Wang et al., 2021), malate efflux transporter Aluminum-activated Malate Transporter 1 (ALMT12; Meyer et al., 2010; Luu et al., 2019), high-affinity chloride transporter ZmNPF6.4 (Wen et al., 2017), and an iron transporter (Vacuolar Iron Transporter 1; VIT1; Zhang et al., 2012a), which were all specifically expressed in GCs (Figure 2B;  Supplemental Data Set 5).

In comparison, K⁺ transporters, including ZMK1 (Büchsenschütz et al., 2005), High‐Affinity K⁺ Transporter 2 (ZmHKT2; Cao et al., 2019b), and orthologs of rice High-Affinity K⁺ 2/10 (OsHAK2/10; Zhang et al., 2012b), were identified in the list of SC marker genes (Figure 2B). Previous research has shown that Z. mays Outward-Rectifying K⁺ channel (ZORK, Zm00001d037289) is expressed in both GCs and SCs (Büchsenschütz et al., 2005), whereas ZORK and its closest paralog (Zm00001d044717), both of which clustered as the ortholog of OsK5.2 (Nguyen et al., 2017), are present in the set of marker genes for GCs and SCs, respectively. Sulfate Transporter 1.2 (SULTR1.2), responsible for sulfate transport, Iron and manganese transporters Natural Resistance-Associated Macrophage Protein 6 (OsNRAMP6), and Metal Tolerance Protein 11 (OsMTP11; Yoshimoto et al., 2002; Peris-Peris et al., 2017; Ma et al., 2018) each have one ortholog that is specifically expressed in maize SCs. In contrast, transcripts of one ortholog of the vacuolar Na⁺ or K⁺/H⁺ antiporter 1/2 (NHX1/2), which in Arabidopsis regulates cell turgor and functions in the stomatal movement (Barragán et al., 2012), are highly enriched in both GCs and SCs.

Transcriptome dynamics and signaling pathways linking GCs and SCs

GCs work as multisensory valves in response to endogenous and environmental factors such as phytohormones, light, and intracellular CO₂ (Schroeder et al., 2001; Lawson and Matthews, 2020). How these signaling pathways are integrated into the well-defined stomatal responses involving the GCs and SCs of the grass stomatal complex remains unknown. For example, ABA is found to promote stomatal closure mediated through kinases/phosphatases, secondary messengers, and ion channel regulation in GCs (Schroeder et al., 2001; Murata et al., 2015). Some of the components in ABA signaling have orthologs identified in maize GCs (Figure 2B;  Supplemental Data Set 5). Transcripts of orthologous genes involved in ABA responses, including Open Stomata 1 (OST1) from SNF1-Related Protein Kinase 2 (SnRK2) group 3 (Mustilli et al., 2002), MYB60/96 (Cominelli et al., 2005), maize Slow Anion Channel-associated 1 (ZmSLAC1; Qi et al., 2018), and SLAC1 Homolog 3 (SLAH3; Zhang et al., 2016), were specific to GCs. In contrast, transcripts of orthologs of AtSnRK2 group 1, independent of ABA regulation (Kulik et al., 2011), and ABA-responsive Kinase Substrate 1/3 (AKS1/3), repressed by ABA through phosphorylation (Takahashi et al., 2013), were present as SC marker genes. Of the two OST1 orthologs, which originated from a maize-specific gene duplication event and were paralogs of previously reported ZmOST1 (Vilela et al., 2013), one was confirmed by rice transformation to be uniquely expressed in GCs (Figure 2C;  Supplemental Figure S5C). Transcripts of the orthologs of Arabidopsis NITRATE TRANSPORTER 1.2 (ZmNRT1.2; Kanno et al., 2012), which showed high levels of ABA and nitrate transport activities when expressed in Xenopus laevis oocytes (Figure 3, A and B), were highly enriched in SC cluster, but only slightly enriched in GC cluster (Figure 3C).

Figure 3

Integrated signaling pathways of ABA or CO₂ between GCs and SCs that regulate stomatal responses. A, ABA contents in X. laevis oocytes expressing injected transgenic transcripts and the corresponding water-injected control, preloaded with 46 nL of 1 mM (12.1 ng) ABA. The data are presented as means ± se, n = 4, 3 oocytes per replicate. P < 0.01. B, Current–Voltage relationships from TEVC recordings of whole X. laevis oocytes expressing ZmNRT1.2- and the water-injected control. The data are presented as means ± se, n ≥ 12, P < 0.01 by t test. C, Violin plots of ZmCYP707A4 (Zm00001d020717), ZmNRT1.2 (Zm00001d046383), Zmβ-CA1 (Zm00001d044099), and ZmMPK12 (Zm00001d024568). Each point represents the expression value (UMI) in a single cell. D, The volume distributions of GCs and SCs. A total 25 stomatal complexes were measured, comprising 50 GCs and 50 SCs. The data are presented as means ± se, n = 50. The border of the boxes, horizontal lines, and whiskers indicate the first and third quartiles, the median values, the minimum and maximum values, respectively. E, The ABA concentrations in GCs and SCs. ABA concentrations were obtained by ABA content/cell divided by the mean volumes of GCs and SCs. ABA contents were measured for 571–848 SCs and 3,160–4,370 GCs. The border of the boxes, horizontal lines, and whiskers indicate the first and third quartiles, the median values, the minimum and maximum values, respectively. F, Current–voltage relationship from TEVC recordings of whole X. laevis oocytes expressing cRNA combined from OST1+AtSLAC1, AtRHC1 + AtHT1 + OST1 + AtSLAC1, and Zmβ-CA1 + AtRHC1 + AtHT1 + OST1 +  AtSLAC1. OST1 activates SLAC1 and induces anion currents, HCO3⁻ produced by β-CA1 enhances the interaction of RHC1 with HT1 and relieves inhibition of OST1 by HT1, and thus increases anion currents by SLAC1. G, Water-loss measurement of the detached leaves from WT and zmmpk12 mutant maize. The data are presented as means ± se, n ≥ 3, P < 0.01. H, Stomatal conductance of the zmmpk12 mutant under different concentration of CO₂. Control, 400 p.p.m; High, 1,200 p.p.m; Low, 50 p.p.m. The data are presented as means ± se, n ≥ 3, P < 0.01 by t test.

Since information is limited about stomatal movements related to other plant hormones, such as auxin, ethylene, and salicylic acid (Murata et al., 2015), we assessed the relative expression of the genes related to signaling of these hormones in GCs and SCs. An additional novel finding is that orthologs of PIN-formed protein 1 (PIN1), involved in auxin efflux, and of GC-specific Zinc-Induced Facilitator-Like 1 (ZIFL1), which indirectly modulates auxin efflux, likely by regulating PIN2 abundance to regulate stomatal closure (Remy et al., 2013), are both present as GC marker genes (Figure 2B;  Supplemental Data Set 5). ZmZIFL1 displayed highly enriched expression in GCs in the rice transformation experiments (Figure 2C;  Supplemental Figure 5D), whereas the transcripts encoding Auxin-Responsive Protein 32 (SAUR32) were enriched in SCs. We observed that SCs exhibited higher levels of expression of transcripts for the gene encoding orthologs to 1-aminocyclopropane-1-carboxylate oxidase 2 (AtACO2), which catalyzes the final step of ethylene formation (Sun et al., 2017). Indeed, 1-aminocyclopropane-1-carboxylate (ACC) is detected in both GCs and SCs (Supplemental Table S2).

To identify genes that are differentially expressed in GCs and SCs in response to blue light signaling, we compared single-nucleus transcription profiles for cell type-specific enrichment of expression of a blue light signal sensor and its targets. Transcripts encoding maize Phototropin 1 (ZmPHOT1) (Suzuki et al., 2014), the ortholog of AtPHOT1 in blue light-mediated stomatal movement (Kinoshita et al., 2001), were highly accumulated in SCs. In contrast, transcripts encoding other components of blue light signaling were elevated in GCs, including orthologs of Convergence of Blue light and CO₂ 1/2 (CBC1/2), which interact with PHOT1 (Hayashi et al., 2020); BLUS1, which is directly phosphorylated by PHOT1 (Takemiya et al., 2013); and long Hypocotyle5 (HY5), which links blue light signaling to the circadian clock (Hajdu et al., 2018), to plasma membrane H⁺-ATPase 5 (AHA5) (Ueno et al., 2005), and to Cation/H⁺ exchanger 20 (CHX20) in light signaling pathways (Padmanaban et al., 2007). Of these, the transcript encoding ZmCHX20 displayed GC-specific expression in our rice transformation experiments (Figure 2C;  Supplemental Figure S5E).

Calcium ions are important second messengers in stomatal signaling (Kudla et al., 2010; Murata et al., 2015). Intriguingly, transcripts of homologs of the Arabidopsis GC-specific Ca²⁺ channel OSCA1 (Yuan et al., 2014; 67.3% identity) were elevated in both GCs and SCs, whereas, in particular, transcripts of orthologs of genes encoding the cGMP-activated nonselective Ca²⁺-permeable channel CNGC15 (Wang et al., 2021) and SCS (Tarnowski et al., 2020) were specifically expressed in GCs. The specificity of expression of ZmSCS1 and ZmCNGC15 was supported by the results of rice transformation (Figures 1, D and 2, C; Supplemental Figure S5, A and F). Similarly, one ortholog of the Arabidopsis Calcium-dependent Protein Kinase 13, a calcium-insensitive CPK implicated in stomatal movement (Ronzier et al., 2014), showed elevated expression in SCs.

Integrated ABA and CO₂ signaling pathways between GCs and SCs that regulate stomatal responses

To understand how ABA and changes in levels of CO₂ are integrated into the complex four-celled stomatal structure of grasses, we reconstructed these paths focusing on GC- and SC-specific expressed components, using a number of different approaches. For example, expression of transcripts of ABA 8'-hydroxylase (ZmCYP707A4), one member of the CYP707 family, of which AtCYP707A1 is involved in the ABA degradation by GCs (Okamoto et al., 2009), was specific to SCs (Figure 3C). This implies that the fast dynamic response for ABA metabolism for coordinating stomatal closure and opening takes place between GCs and SCs.

To provide more evidence on this, we measured the concentrations of ABA and its catabolites in GCs and SCs from light-adapted leaves. We separately collected GCs and SCs from the freeze-dried peels, using laser capture microdissection (LCM), and then measured the volumes of the GCs and SCs by quantitative confocal microscopy to allow calculation of intracellular ABA concentrations (Figure 3D;  Supplemental Figure S7; Jonkman et al., 2020). Two ABA catabolites, phaseic acid (PA) and dihydrophaseic acid (DPA), were undetectable, whereas ABA was detected in both GCs and SCs. The concentration of ABA in SCs was found to be 75% of that in GCs, given that SC volumes were estimated to be about 7.86 times of those of GCs for our LCM samples (Figure 3, D and E; Supplemental Tables S2 and S3). Consistent with this, transcripts for an ortholog of the ATP-Binding Cassette transporter G 22 (ABCG22), likely responsible for ABA transport (Kuromori et al., 2011), also displayed elevated expression in GCs and SCs (Figure 2B). The different ABA homoeostasis features recorded for GCs and SCs might be associated with coordination of stomatal closure, for which the GCs need to be actively supplied with ABA from SCs and vice versa.

Questions as to whether the CO₂ signal is sensed directly by GCs or by other cell types (MCs or SCs), and as to how these signals converge on stomatal closure signaling pathways are also matters of debate (Frommer, 2010; Hu et al., 2010). In our data, high leaf temperature 1 (HT1), which controls stomatal aperture in response to CO₂ in Arabidopsis (Hashimoto et al., 2006), and two HT1 orthologs, derived from another maize-specific gene duplication, were found to be specifically expressed in GCs. The maize ortholog of CO₂ sensor β-CA1 (Kolbe et al., 2018) showed elevated expression in SCs as well as in MCs (Figure 3C). To assess whether SC-enriched transcripts of Zmβ-CA1 can reconstitute the OST1 activation of the SLAC1 pathway, we co-injected Arabidopsis cRNAs of SLAC1, OST1, Arabidopsis RESISTANT TO HIGH CO₂ (RHC1), and HT1 with or without Zmβ-CA1 into oocytes. We found that OST1-activated SLAC1 was significantly inhibited by HT1 (Figure 3F). When Zmβ-CA1, RHC1, HT1, OST1, and SLAC1 were co-expressed, HT1-inhibited SLAC1 activity was partially recovered (Figure 3F), indicating that Zmβ-CA1 could convert CO₂ into HCO₃⁻, and that SCs have the potential for involvement in a CO₂-HCO₃⁻ sink and in CO₂ sensing. To further confirm this synergistic effect, ZmMPK12, a marker gene co-expressed both in GCs and SCs (Figure 3C) and one ortholog of AtMAPK12 in the CO₂ signaling pathway (Jakobson et al., 2016), was knocked out using a reverse genetics approach. One zmmpk12 EMS mutant carrying a premature stop codon displayed a lower fresh weight than the wild-type (WT) under well-watered conditions (Supplemental Figure S8, A–D). A more prominent wilting phenotype in the mutant was observed following drought treatment, consistent with the results provided by measurement of water-loss rate and transpiration rate (Figure 3G;  Supplemental Figure S8E). Under conditions of high CO₂ levels, the zmmpk12 EMS mutant displayed impaired stomatal closure (Figure 3H). This strongly suggests that ZmMPK12 is involved in the regulation of maize stomatal movement. Further experiments are underway to refine this distinctive role for the ZmMPK12 gene in GCs and SCs, and should help elucidate its importance in stomatal responses. It seems likely that the evolutionary segregation of CO₂ sensing and anion channel regulation into distinct cellular locations would represent an efficient way for stomatal turgor control.

Developmental tracking of stomata in the maize leaf base

To obtain insights into the molecular basis of GC and SC development, we captured and analyzed 8,967 nuclei from a leaf base sample (Supplemental Table S1; Supplemental Figure S9, A and B), obtaining a total of 14 clusters and their DEGs (Clusters 0–13, Figure 4A, Supplemental Data Set 6). To annotate the clusters, we first integrated our data with those reported for two scRNA-seq datasets derived from maize shoot apical meristems (SAMs) including the six most recently initiated leaf primordia (SAM + P6) based on 10× Genomics platform (Satterlee et al., 2020). A high level of concordance could be observed between the snRNA-seq and scRNA-seq datasets (Supplemental Figure S9, C and D). The cell types/states were primarily assigned according to the cluster annotation from scRNA-seq by Satterlee et al. (2020), except for Clusters 0, 2, and 5 (Supplemental Figure S9E): Cluster 1 was overrepresented by cells in leaf primordia; Cluster 7 was enriched in cells of the G2/M and S phases; and four clusters (3, 10, 11, and 13) and five clusters (4, 6, 8, 9, and 12) showed high association with epidermis and vasculature cells, respectively. Next, we employed previously reported marker genes and those discovered in our mature peel samples to further associate the epidermal cells in the leaf base with specific cell types (Figure 4B). Two marker genes found in PCs of the peels, encoding GDSL esterase/lipase 3 (GELP3, Zm00001d011661) and receptor-like protein 7, were highly enriched in Cluster 3, suggesting that this cluster represents developing PCs. The enrichment of transcripts for Liguleless1 (Moreno et al., 1997), a good marker gene for Cluster 10, implies this cluster represents ligules and/or auricles. Two previously reported marker genes in maize stomata showed specific expression in Cluster 11, GC-specific FAMA and SC-specific CST1 (Wang et al., 2019b). In addition, 66 marker genes from the GCs and SCs of mature peels were also expressed as marker genes in the leaf base sample, 37.8% of which (25 genes) displayed unique expression in Cluster 11 (Supplemental Figure S10). Taken together, this information suggests that Cluster 11 represents developing stomata. In comparison, the two SAM + P6 datasets had few cells clustered with the specific expression of GC and SC marker genes or of genes known to be involved in stomatal development.

Figure 4

Developmental trajectories of maize stomata. A, UMAP visualization of the nuclear classification of the leaf base. Colors represent their respective Louvain components, and curved lines represent cell populations. Cluster 4 lies adjacent to the other four vasculature clusters in the 3D-UMAP plot. B, Dot plot of known marker genes used for the identification of cell types in the developing epidermis. The proportions and average expression levels are, respectively, denoted by the sizes and colors of the circles. C, Developmental stages of stomata defined in this study and visualization of genes along the pseudo-time trajectory of stomatal development produced by Monocle. Colors indicate different cell types/states, and the arrows show the developmental decisions. D–H, Microscope observation of YFP fluorescence in maize and rice transformation lines. D, CST1; (E), ZmGELP1 (Zm00001d040239); (F), ZmGELP2 (Zm00001d018474); G, ZmMPK12 (Zm00001d024568); H, ZmMUM2 (Zm00001d042656). The cell outlines are stained using PI. Scale bar = 20 μm. I, Model of stomatal development according to previous studies and our data. Six main developmental stages are shown: (1) Specification of the stomatal file; (2) GMC formation; (3) Formation and polarization of SMC; (4) Asymmetric division of SMCs to generate SCs; (5) Symmetric division of GMCs producing paired GCs; (6) Differentiation and morphogenesis of the four-celled stomatal complex.

To obtain detailed insight into the networks regulating stomatal development, a pseudo-time trajectory was constructed using all of the 116 nuclei of Cluster 11 (Figure 4C). Three genes encoding SPCH are found in the maize genome and, in this cluster, only one was identified by expression as being detected in a single nucleus; expression of MUTE was only detected in three nuclei (Supplemental Figure S11). From these, we conclude that Cluster 11 lacks protodermal cells (PDCs), most likely due to the blending step during the sample preparation, which, therefore, includes stages only from a very early developmental stage of GMCs to the stage representing young GCs (YGCs) and young SCs (YSCs). Consequently, only guard mother cells (GMCs), subsidiary mother cells (SMCs), YGCs, and YSCs were defined along the developmental trajectory through examination of the expression of previously reported players and orthologous genes (Figure 4C;  Supplemental Figure S11). SMCs were identified by expression of PAN1 and ARP2/3 (Gray et al., 2020; Nunes et al., 2020), and YSCs by CST1 expression, whose specific expression was also observed in YSCs as well as in the mature SCs (Wang et al., 2019b). GMCs–YGCs were identified by expression of FAMA, with GMCs being demarcated by expression of SCREAM (SCRM/ICE1) and EPF2 (Lu et al., 2019). The remaining nuclei were defined as being derived from YGCs. Transitions of cell states were largely verified by analyzing the transcriptional dynamics with a steady-state deterministic model using scVelo (Bergen et al., 2020) based on the ratio of unspliced and spliced RNA transcripts (Supplemental Figure S12).

As expected, for transgenic rice, when the ZmMPK12 promoter was employed to regulate the expression of a nuclear YFP marker, enrichment of fluorescence was observed in young and mature stomata (Figure 4G). This is consistent with ZmMPK12 expression being a stomatal marker for the developing YGCs–YSCs, as well as the mature GCs–SCs (Figures 3, C and 4, C). Similarly, transgenic maize plants expressing the Pro_CST1:CST1-eYFP:Ter_CST1 transgene also displayed strong fluorescence signals in YSCs (Figure 4D). We selected two of the marker genes encoding GELPs (ZmGELP1 and ZmGELP2) for maize transformation to further verify the identification and positioning of YGCs and YSCs along the stomatal developmental trajectory. As seen in Figure 4, E and F, ZmGELP1 expression was specific to YGCs-GCs, and ZmGELP2 specific to young and mature stomata. Importantly, one β-galactosidase 16 gene (ZmMUM2, named in this study), orthologous to mucilage-modified2 (MUM2) in Arabidopsis (Dean et al., 2007), was specifically expressed in GMCs, and the observations of fluorescence within transgenic maize corroborated the predicted transcription patterns (Figure 4H). These transgenic experiments therefore provided more evidence for the assignments of cell states along the developmental trajectory, as well as the new marker gene predictions.

Further inspection of the expression of orthologous genes known to control stomatal development revealed a nearly complete recapitulation of the model previously proposed for stomatal development (Figure 4I;  Supplemental Figure S13). The expression of the gene encoding too many mouths (TMMs), a membrane receptor-like protein playing a multifaceted role during stomatal formation (de Marcos et al., 2016), is specifically detected in GMCs but is absent from SMCs, YSCs, and YGCs (Supplemental Figure S11). The signals perceived by TMM come from the EPF family for stomata production (Hughes et al., 2017; Dunn et al., 2019). EPF1 transcripts are detected only at the middle time point of GMC development. In contrast, EPF2 mRNA transcript levels are highly enriched across the entire GMC stage of development (Figure 4C). Transcription of genes encoding SDD1, a subtilisin-like serine protease that may activate EPF family peptides as a negative regulator in stomatal development (Von Groll et al., 2002; Fanourakis et al., 2020), was specific to GMCs (Supplemental Figure S11). The signal from TMM-ERf is further transduced to repress the expression of SPCH, MUTE, and FAMA by the MAPK cascade, thereby ensuring normal stomatal development (Chen et al., 2020). ZmMAPKKK38 (Liu et al., 2015) and ZmMPK12 (Liu et al., 2013) transcripts were found across all stages of the trajectory (Figure 4C;  Supplemental Figure S11). FAMA was highly expressed in GMCs and YGCs (Figure 4C), as well as within mature GCs. The model recapitulated from our data, therefore, provides a detailed framework outlining the regulatory network in stomatal development from GMCs to GCs. Experimental evidence is currently lacking for the SMC developmental stage, and further work is needed to confirm this aspect.

The potential roles of cell wall-related genes in dumbbell-shaped morphogenesis of developing and mature stomata

The cell wall undergoes multiple dynamic changes during stomatal development, including the formation of a new cell wall during the symmetrical division of GMCs, the degradation of the cell wall between two YGCs to form pores, the thickening of the GC wall at stomatal maturation, and the dynamic deformation accompanying stomatal movement (Galatis and Apostolakos, 2004; Rui et al., 2018; Spiegelhalder and Raissig, 2021). A most interesting question is how wall polysaccharides, lipids, and proteins are coordinately synthesized, deposited, reorganized, modified, and degraded in GCs over the functional lifetime of stomata. Our data indicated that multiple genes having orthologs reported to be involved in the modulation of cell wall plasticity are expressed in mature SCs and GCs (Figure 2A;  Supplemental Table S4). In the mature SCs, this includes genes encoding ATP-Binding Cassette G type 29 (ABCG29), involved in cell wall lignin biosynthesis (Alejandro et al., 2012), Endo-1,4-β-glucanase 11 (glycosyl hydrolase 9B8, GH9B8), involved in cell wall loosening (Xu et al., 2014), a non-specific Lipid Transfer Protein (LTP2), involved in the interface integrity of cell wall with cuticle (Jacq et al., 2017), glycosyltransferase (GT43), associated with secondary cell wall synthesis (Lee et al., 2014), and xyloglucan endotransglucosylase/hydrolase (OsXTH25), involved in cell wall modification and loosening (Yokoyama et al., 2004). In addition, a gene encoding an ortholog of AtMPK18, involved in microtubule organization (Walia et al., 2009), is present as an SC marker gene. One gene encoding Wall-Associated receptor Kinase 1 (ZmWAK1, named in this study), homologous to WAK1 which, in Arabidopsis, encodes a cell wall protein that binds to pectin and interacts with glycine-rich proteins (Park et al., 2001), was specifically expressed in SCs. The SC-specific expression of cell wall glycine-rich protein ZmGRP6 was confirmed by our transformation experiments (Figure 1C;  Supplemental Figure S5B), implying conservation of the GRP/WAK signaling pathway in cell wall expansion (Mangeon et al., 2017). Homologous transcripts encoding L-type lectin RLK, potentially functioning in sensing cell wall-plasma membrane connections (Gouget et al., 2006), and one G-type lectin RLK (ZmlecRLK1, named in this study), were overrepresented in SCs. Of these, ZmlecRLK1 was specifically expressed in the SCs of our transgenic maize and rice plants (Figure 2C;  Supplemental Figure S5, G and H). For GCs, reversible elongation and contraction may be mediated by the activities of orthologous genes encoding galacturonosyltransferase 7 in pectin biosynthesis (Atmodjo et al., 2011), polygalacturonase involved in expansion3 (PGX3) in pectin degradation (Rui et al., 2017), alpha expansin 7 (OsEXPA7) in cell wall extension (Jadamba et al., 2020), XTH28 in controlling mechanical properties of the cell wall by restructuring the xyloglucan cross-links (Kurasawa et al., 2009), and GH9B1 in cell wall loosening (Xu et al., 2014). Of these, rice transformation with the Pro_PGX3:HTA6-YFP construct displayed highly GC-enriched expression (Figure 2C;  Supplemental Figure S5H). Consistent with this, pectin homogalacturonan nanofilament expansion has been shown to drive shaping the cells in plant epidermal cells (Haas et al., 2020). Together these provide candidates for further investigation of the mechanical modulation of the wall properties of SCs and GCs associated with stomatal movement.

To identify the key genes in stomatal cell wall development, we examined the lists of cell-type-specific marker genes, and found a large number having orthologs involved in cell wall metabolism, particularly of pectin (Figures 4, I and 5, A; Supplemental Figure S14; Supplemental Table S4). Orthologs of two genes involved in pectin metabolism, Pectin Methylesterase 53 and Pectin Acetylesterase 12, are mainly expressed at the GMC and SMC stages (Francis et al., 2006; Gou et al., 2012). Orthologs of genes involved in the modification and degradation of cellulose and lignin were expressed across all the four cell types: ZmMUM2 (orthologous to AtMUM2; Dean et al., 2007) and endoglucanase 24 (orthologous to Cellulase 6; He et al., 2018) dominantly expressed in the GMCs, of which the specific expression of ZmMUM2 was verified by our maize transformation experiments, β-glucosidase 47 (BGLU47, orthologous to AtBGLU45/46; Chapelle et al., 2012) was mainly expressed in YGCs, and Peroxidase 20 (PRX20, orthologous to AtPRX12 and rice PRX20; Irshad et al., 2008) in YSCs. Orthologs of Arabidopsis GPI-anchored non-specific LTPs 1 and 2 (ZmLTPG12 and ZmLTPG22; Lee et al., 2009; Kim et al., 2012), involved in the lipid transfer associated with cuticular metabolism, were preferentially expressed in the SMCs and GMCs. LTPG10, expressed at the YSC and YGC stages, may regulate lipid transfer during stomata maturation. Members of GELPs (ZmGELP1 and ZmGELP2) were mainly expressed following the symmetric division of GMCs, and may play roles in wax or cutin metabolism of the stomatal cell wall, similar to those observed for their Arabidopsis and rice homologs Occlusion of Stomatal Pore 1 (Tang et al., 2020), Cutin Synthase2 (CUS2; Hong et al., 2017), and Wilted Dwarf and Lethal 1 (Park et al., 2010). Further exploration appears warranted into the way in which the fine temporal and spatial changes of the activities of genes, through regulating the formation and degradation of cell walls, are involved in the development and three-dimensional shaping of stomata in maize.

Figure 5

The heatmaps of genes involved in maize stomatal development of different cell types/states and of candidate genes regulated by ZmMUTE. A, The heatmap of marker genes specifically expressed in different cell types/states of maize stomata. B, Venn diagram of the marker genes from bulk RNA-seq data (left) and Cluster 11 of the leaf base snRNA-seq sample (right). C, Expression heatmap of the shared marker genes in the bulk RNA-seq data from ZmMUTE mutants. bzu2, the ZmMUTE mutant. Scale bar = 20 μm.

ZmMUTE/BZU2 controls not only asymmetrical division in the early GMC stage, but also cell wall formation in stomata

Maize Zmmute/bzu2 mutants lack normal stomatal complexes, and we have previously reported that maize ZmMUTE/BZU2 functions as a mobile signal in the establishment and differentiation of SMCs (Wang et al., 2019a). To further examine the role(s) played by MUTE in maize stomatal development, we generated bulk RNA-seq data from samples of the WT and ZmMUTE mutant plants at three stages of early stomatal development, respectively, enriched with PDCs, GMCs, and YGCs–YSCs. A total of 999 DEGs were identified when comparing the WT and bzu2-1 mutants at each of the three developmental stages (adjusted P < 0.05; fold change > 2) (Supplemental Figure S15A; Supplemental Data Set 7). The intersection of the up-regulated DEGs from the 14 clusters of snRNA-seq data and the DEGs of the three developmental stage comparisons from the bulk RNA-seq data comprised 137 genes (Supplemental Figure S15B), of which 28 genes were previously identified as marker genes for Cluster 11 (Figure 5B;  Supplemental Data Set 8). All 28 DEGs are from the GMCs and YGCs–YSCs developmental stages in the bulk RNA-seq, which is consistent with the expression pattern of ZmMUTE promoter–reporter lines (Wang et al., 2019a). Hierarchical clustering based on the expression profiles of the 28 genes in the bulk RNA-seq data identified 18 that showed increased transcript levels in GMCs and YGCs–YSCs of the WT, but that were greatly downregulated in the corresponding stages of the mutant (Figure 5C). The minimum E-box consensus binding sequence characteristic of bHLH family proteins including MUTE (Wang et al., 2019a) was detected in 24 of these 28 genes (Supplemental Data Set 8), which further narrowed down the candidate targets of ZmMUTE. The 28 DEGs included homologs to (1) genes known as key stomatal developmental regulators, including EPF2, controlling stomatal density and the differentiation of GMCs (Lu et al., 2019), and TMM (de Marcos et al., 2016), (2) ZmCST1 (Wang et al., 2019b) and a gene encoding a monosaccharide transporter, involved in sugar transport, (3) genes involved in cuticle and wax biosynthesis, including 3-ketoacyl-CoA synthase, GELPs, and ECERIFERUM 3; genes involved in cellulose hydrolysis, β-glucosidase; genes involved in cell expansion, a rapid alkalinization factor family gene, and (4) genes involved in pectin dynamic metabolism that show increased expression in bzu2-1 mutants and enrichment in GMCs: pectin methylesterase38 and a polygalacturonase ortholog of AtPGLR (Hocq et al., 2020). Other important regulatory components directly or indirectly impacted by ZmMUTE were also identified, including genes encoding ZmSCS1, ICP4 domain-containing protein, two NAC domain-containing proteins, a protein phosphatase 2C family protein, and a protein kinase (Figure 5C;  Supplemental Data Set 8). Together, these observations provide a detailed overview of the regulatory role of ZmMUTE in stomatal development, and a valuable guide for the design of further experiments.

Discussion

We optimized and established the snRNA-seq procedure for analysis of the cell types in maize mature peels and leaf base samples based on the 10× Chromium platform. The marker genes demonstrating specific expression of GCs and SCs were verified by transgenic experiments, which indicated the established procedures were reliable and accurate. Notably, the corroboration by maize and rice transformation of the bioinformatics results and predictions also implies conservation of gene functions in monocots, as seen for the transgenic plants regulated by the ZmlecRLK1 promoter, which displayed SC-specific expression in both maize and rice. Using our methods, the gene detection number is comparable to the latest snRNA-seq studies using Arabidopsis young roots (Farmer et al., 2021) and maize young above-ground parts (Marand et al., 2021). The lower gene detection number, relative to that seen for scRNA-seq technology applied to protoplasts (Denyer et al., 2019; Jean-Baptiste et al., 2019; Ryu et al., 2019; Seyfferth et al., 2021), had been considered to limit the use of snRNA-seq as a stand-alone technique for a complete discovery of cell-type heterogeneity in complex tissues and organs (Shaw et al., 2021). However, our results conclusively demonstrate that snRNA-seq can produce sufficient signature genes for detailed comparative analyses.

The stomatal cell states and the stomatal cell fate-determining factors have been characterized based on scRNA-seq technology in A. thaliana (Lopez-Anido et al., 2021; Zhang et al., 2021). Whereas the trajectories of maize stomatal development have been examined by integrating data from a single-nucleus Assay for Transposase-Accessible Chromatin using sequencing and snRNA-seq methods, the biological features of stomata were not investigated (Marand et al., 2021). In our study, we have generated snRNA-seq data from maize leaf bases with developing stomata. Known orthologous key players and novel genes were identified from the developmental trajectory of the stomatal cells, highlighting the active metabolism of the cell wall and the different roles of ZmMUTE. Our study, therefore, provides an invaluable resource for in-depth functional research on stomatal movement and development in monocots, and also guarantees that comparative studies with dicots are worthwhile. A model based on cell type specifically expressed genes has emerged that provides deeper understanding of the molecular mechanisms underlying stomatal movement (Figure 6).

Figure 6

Model of cell type-specific genes that are potentially involved in maize stomatal movement. The known genes or orthologs involved in stomatal movement are displayed. Transporters of non-signaling cations/anions and signaling molecules are indicated in blue and orange, respectively. Genes in cell wall plasticity are shown in gray. The signaling modules of ABA, CO₂, Ca²⁺, and blue light in GCs and SCs are indicated via different background colors.

How do SCs and GCs co-operate to regulate stomata movement in responding to the changes in the external environment? Prior to the availability of the single-cell transcriptomic approach, it was not possible to get a comprehensive view of the molecular basis as to how SCs might serve as an ion reservoir and/or provide mechanical assistance to GC movements (Hetherington and Woodward, 2003), facilitating rapid opening and closing (Figure 6). In terms of ABA signaling, a model for ABA action has been proposed and validated, in which the soluble PYR/PYL/RCAR receptors function at the apex of a negative regulatory pathway to directly regulate PP2C phosphatases, which in turn regulate SnRK2 kinases (Park et al., 2009; Santiago et al., 2009). Our results indicate that transcripts of the maize orthologs of the core ABA signaling component OST1 and of key components in the downstream signaling transduction, MYB60/96, MPK18, and SLAC1, are all found in GCs. However, the transcripts of genes encoding ABA transporters and hydroxylases, such as ZmNRT1.2 and ZmCYP707A4, were identified in SCs (Figure 6). Therefore, it becomes reasonable to propose that the closure of stomata might be affected by the ABA signaling pathway in GCs combined with ABA degradation in SCs. This suggests the important role of ABA signaling between GCs and SCs for coordinating speedy stomatal movement. These results also imply an essential function of SCs in providing a dynamic model that allows stomata to respond quickly and efficiently to intrinsic ABA and environmental signals via SC-mediated ABA degradation.

Another aspect implicating GC and SC coordination in stomatal modulation involves CO₂ and blue light signal sensing and action. For example, the maize Zmβ-CA1 ortholog showed elevated expression in SCs as well as in MCs and, therefore, could enhance the activity of AtSLAC1 in the reconstituted CO₂-sensing pathway, whereas two HT1 orthologs, derived from another maize-specific gene duplication, were detected as specifically expressed in GCs. This suggests that, in maize, CO₂/${\text{HCO}}_3^{-}$⁠, acting as an intracellular messenger in SCs, is transported to the GCs which then initiates Ca²⁺ signaling and modulates SLAC1 activity via OST1 or HT1 (Figure 6). This suggests these dynamic spatial and temporal changes for modulation of GC signaling are important facets of stomatal behavior, since the timely, rapid, and efficient regulation of stomatal closure and opening benefits plants either encountering variations in environmental conditions over the day, or sudden exposure to different stresses. The presence of GCs and SCs in grasses would, therefore, provide an ideal and convenient location for sensory or regulatory mechanisms for superior dynamic performance in drier and more rapidly fluctuating environments (Pillitteri and Torii, 2012; Chen et al., 2017).

There are at least four functional mechanisms whereby GC chloroplasts might contribute to stomatal opening, including those Calvin cycle activities that contribute to stomatal regulation, to the supply of ATP, to BLUS, and to starch metabolism (Lawson, 2009; Lawson and Blatt, 2014). Our transcriptome data reconstituted the full suite of key components associated with these chloroplast contributions in GCs, BSCs, and MCs. For instance, sucrose transport from the apoplast has been long regarded as the main carbohydrate source for GCs (Lawson and Matthews, 2020). CST1 (Sweet1b) has been reported to transport glucose specifically in SCs and may regulate stomatal movement (Wang et al., 2019b). An examination of all maize orthologs of Arabidopsis and rice plasma membrane sugar transporters identified only sugar transport protein1 (STP1) and STP5 orthologs in the set of GC marker genes (Figure 6). AtSTP1 has been reported as being responsible for glucose transport from MCs to GCs (Flutsch et al., 2020). The limited capacity of GCs for carbon fixation, coupled to the specific expression of STP1/5 in GCs, implies that glucose is the main carbohydrate imported into GCs, where it may originate from neighboring SCs, PCs, or apoplastic connections to the MCs. It could also be that glucose, acting as a diffusible mesophyll signal, is involved in the active control mechanism governing stomatal opening.

In our data, we detected proteins previously identified as being involved in maize stomatal development, including SDD1, FAMA, and EPFs, as well as proteins more recently reported as being involved in this process, such as ZmPAN1, ARP2/3, and a leucine-rich receptor-like protein (McKown and Bergmann, 2020; Nunes et al., 2020). Interestingly, transcripts for genes related to the sophisticated morphology, mechanics, and cell wall structure of grass stomata were enriched in the developing stomata. Our previous work indicated that BZU2/MUTE is required for early events in SMC initiation and differentiation, as well as for the last symmetrical division of GMCs that produces the two YGCs (Wang et al., 2019a). For the MUTE-targeted genes, 12 of these 28 DEGs were involved in cell wall plasticity modification, biosynthesis, and sugar transport steps in the unique dumbbell-shaped GCs of grass stomata (Supplemental Figure S13). Importantly, besides receiving mobile MUTE signal from GCs for determination of the antecedent SMC fate (Raissig et al., 2017; Wang et al., 2019a) and coordinating GC turgor regulation, we note that SCs could contribute to dumbbell-shaped GC formation (Nunes et al., 2020). This implies a more sophisticated relationship between GCs and SCs associated with stomatal movement and formation. Our data, therefore, help to refocus the transcriptomic framework for interpreting the role of MUTE in the determination of GMC fate and the initiation of intercellular signaling in the establishment of the four-celled stomatal structure of grasses.

Materials and methods

Plant materials and growth conditions

Zea mays ssp. mays var. B73 was used in snRNA-seq experiments. The seeds of the Pv2 sample were grown at Henan University in a greenhouse, with a 20–28°C 16-h day/8-h night (300 µmol m⁻²s⁻¹, sunlight complemented with high-pressure mercury lamps) and ambient humidity. For all the other samples, seeds were sown in a mix of nutrient soil and vermiculite (1.5:1), and grown under controlled environmental conditions in a step-in chamber (Percival PR1010L, Henan University, China) with a 25°C/22°C 16-h day/8-h night light cycle (300 µmol m⁻²s⁻¹) at 65% relative humidity.

Sample preparation, nuclei sorting, and sequencing of snRNA-seq libraries using the 10× Genomics Chromium platform

For abaxial epidermis preparation, peels of the first and second leaves from around 100 seedlings (10 DAG) were stripped mechanically using micro-forceps. For the leaf base sample, basal 1-cm sections of the second leaf of 10-day-old seedlings were collected. The accumulated peels and the leaf base sample were kept under ice water prior to individual processing for nuclear sorting as previously described (Galbraith et al., 2021; Galbraith and Sun, 2021). Briefly, blot-dried materials (0.3 g fresh weight) were chopped 200–300 times using a new double-edged shaving razor blade in 2 mL Galbraith buffer (Galbraith et al., 1983) supplemented with 2% bovine serum albumin (BSA) and the addition of 20 µL dithiothreitol (DTT) (1 M DTT in DEPC-treated diH₂O), 20 µL Ambiont^TM RNase inhibitor (40 U/μL, Invitrogen, Waltham, MA, USA), and 10 µL RNAsin (40 U/μL, Promega, Madison, WI, USA) in a 5-cm diameter disposable sterile plastic petri dish. The homogenates were separately filtered through 30-µm CellTrics filters (Sysmex, USA), and PI was added to a final concentration of 50 µg/mL (Shanghai Yuanye Biotech Ltd., China). Nuclei were flow sorted from the filtered homogenates using a BD FACSAria II or a Beckman Coulter MoFlo XDP, operated in the single-droplet sorting mode, using 70 or 75-µm diameter flow tips, with 0.01-M phosphate buffered saline dissolved in DEPC-treated water as the sheath buffer, and sorting into a 2-mL protein LoBind tube (Eppendorf, Hamburg, Germany) prerinsed with 2% BSA. A total of 40,000 nuclei were typically recovered within 5–10 min. The integrity of the sorted nuclei was estimated in terms of the nuclear morphology under phase contrast microscopy (Zeiss Axio Imager M2 or Nikon TS100-F). The concentrations of the sorted nuclei were confirmed using an automatic cell counter (Luna classic or BioRad TC10) or a hemocytometer.

Individual samples of the sorted nuclei (20,000 nuclei within a total volume not exceeding 33.8 µL for the version 2 kit and 46.6 µL for the version 3 kit) were employed for 10× Chromium processing. Second-strand cDNA production, and sample amplification steps were done as described by the manufacturer. The peels from greenhouse-grown samples were processed using the 10× Genomics single-cell kit v2, to produce sample Pv2, and the other three samples (Pv3a, Pv3b, and leaf base sample) were produced with the updated single-cell kit version 3.0. Peel preparation was slightly different for the three samples: peels of Pv2 were stripped with water, Pv3a and Pv3b were stripped with D-sorbitol medium (310 mM) with a difference that the leaves of sample Pv3b were totally immersed in sorbitol medium (310 mM) during stripping such that the stripped peels had no opportunity to contact the air. Rapid processing and storing samples on ice was emphasized to obtain high-quality, intact nuclei, and to limit RNA degradation. Centrifugation was not employed, in order to avoid non-specific nuclear adhesion. Sequencing libraries were constructed according to the Illumina Nextera protocol and sequenced using a NovaSeq 6000.

Nucleus clustering and DEG prediction of snRNA-Seq

The barcode-gene-UMI matrices were obtained by mapping the Illumina sequencing data to the maize B73 reference genome (B73 RefGen_V4) using Cell Ranger Single-Cell Software (verson 1.3.1) with default parameter settings, followed by removal of organelle-encoded genes. Considering that some mRNAs may not be processed completely to remove introns, all the coding regions were changed to exons in the genome General Transfer Format (GTF) annotation file. Before clustering, barcodes of low quality were filtered from the matrices as follows: for Pv3a and Pv3b, barcodes with abnormal expression levels were identified and were removed by application of the Isolation Forest algorithm (sklearn.ensemble.Isolation Forest Python module) using default settings. For the Pv2 and leaf base samples, the barcodes with detected gene numbers less than 300 and 800 were, respectively, discarded. Next, the remaining barcodes in each dataset were independently clustered based on 4,500 highly variable genes and the top 70 Principal Components (PCs) (resolution = 0.8) using the Seurat R package version 3.0. For the three peel samples, integration was performed using the following functions in Seurat (Butler et al., 2018): NormalizeData, to normalize the snRNA-seq datasets, FindIntegrationAnchors, to discover integration anchors across three samples using the first 55 dimensions, IntegrateData, to generate the integrated assay, ScaleData, to scale the integrated assay, RunPCA, to run PCA analysis (55 PCs), FindNeighbors, to construct the shared nearest neighbor graph, and FindClusters, to cluster the nuclei based on Louvain. RunUMAP was used to visualize the data, with clustering analysis being performed on the integrated assay at resolutions of 0.02–0.8. The final clusters corresponding to different cell types were obtained by manual assignment based on the known marker genes of all the cell types in the anatomic observation of the peels.

DEGs within the final clusters were identified using Findmarkers in Seurat (Butler et al., 2018) with the default non-parametric Wilcoxon rank sum test. Marker genes were selected from the upregulated DEGs with logFC > 0.7 and adjusted P < 0.05. Gene ontology (GO) enrichment was performed by singular enrichment analysis implemented in AgriGO (Tian et al., 2017). The featured pathways were checked manually for gene members in each family. Pseudotime developmental trajectories were inferred by Seurat (Butler et al., 2018) and monocle 3 (Cao et al., 2019a) using all 116 stomatal nuclei of Cluster 11. The nuclei were partitioned and organized into a trajectory using the learnGraph function in monocle 3 (Cao et al., 2019a) with default parameters.

Sample collection, sequencing, and analysis of bulk RNA-seq from the ZmMUTE mutants

According to our understanding of the spatiotemporal process of stomatal development in maize plants and our morphological observations under the microscope of the bzu2-1 (ZmMUTE) mutants (Wang et al., 2019a), 0–0.5, 0.5–1, and 1–1.5 cm segments of the basal leaves of eight DAG seedlings, enriched with PDCs, GMCs, and YGCs–YSCs, respectively, were collected from B73 WT and bzu2-1 mutants, and were processed for bulk RNA-Seq sequencing. Each sample was made from three seedlings and three biological replicates were generated, which generated 16 RNA-seq libraries and datasets. RNA samples were prepared using the TRIzol reagent (Life Technologies, Carlsbad, CA, USA), the quality of each sample being assessed with a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). cDNA libraries were constructed according to the instructions provided by the TruSeq RNA Library Preparation kit, and were sequenced using an Illumina HiSeq2500. Quantification of gene expression levels and prediction of DEGs were done using B73 RefGen_V4 as previously described (Qin et al., 2020). Genes with a fold change >2.0 and an adjusted P < 0.05 for each developmental stage were considered as representing significant DEGs.

Identification of maize orthologs of Arabidopsis and rice genes, and prediction of subcellular localization

Orthologous genes were generally obtained by searching the ortholog and paralog predictions with high confidence in the Ensemble plants website (Bolser et al., 2016). A phylogenetic pipeline was also constructed to determine the orthologs across other representative species, including Klebsormidium nitens, Marchantia polymorpha, Amborella trichopoda, A. thaliana, Vitis vinifera, Populus trichocarpa, Ananas comosus, O. sativa, and Z. mays. In brief, the gene families of the nine species were first clustered by OrthoFinder (Emms and Kelly, 2019), then multiple sequence alignments were performed on the ortho-groups using mafft (Katoh and Standley, 2013), and trimmed by trimal (Capella-Gutierrez et al., 2009). Maximum-likelihood phylogenetic trees were produced using IQ-TREE (Nguyen et al., 2015) with 1,000 ultrafast bootstraps, and were visualized with the R package ggtree (Yu et al., 2017). Orthologous genes were checked manually within the resulting trees.

The subcellular location information was mainly based on the literature, with other information coming from cropPAL2020 (Hooper et al., 2020); https://crop-pal.org) and SUBA4 (Hooper et al., 2017); https://suba.live) datasets, by searching the maize genes or corresponding Arabidopsis orthologs.

Transformation and observation of the predicted cell-type specific genes in maize and rice

Reporter constructs of Pro_gene:HTA6-YFP were transformed into Z. mays inbred line B104 or Oryza sativa L. japonica cv Zhonghua 11, using Agrobacterium tumefaciens strain EHA105 as previously described (Wang et al., 2019a). The maize transformation vector pWMB110, containing the bar gene as a selection marker, was cleaved by restriction endonuclease BamHI and HindIII. The rice transformation vector pBWA(V)H-ccdb-yfp (from Wuhan BioRun Biosciences Co., Ltd), containing hyg as a selection marker, was cleaved by Eco311. ZmHTA6 (Zm00001d047787), an ortholog of A. thaliana HTA6 (AT5G02560.2) was cloned from B73 leaf cDNA. Six alanines (GCT GCC GCA GCG GCT GCC) were added as a linker between the ZmHTA6 and YFP coding sequences to ensure appropriate folding of the chimeric ZmHTA6-YFP protein. The ZmHTA6-YFP was fusion into pWMB110 and pBWA(V)H-ccdb-yfp, respectively, by homologous recombinase (Vazyme, ClonExpress II One Step Cloning Kit).

The promoter regions (1.5-kb upstream of the start codon) of the target genes were amplified from B73 genomic DNA using the Phanta Max Super-Fidelity DNA Polymerase (Vazyme, China) and ligated upstream of a chimeric translational fusion ZmHTA6-YFP in the vector by homologous recombinase (Vazyme, ClonExpress II One Step Cloning Kit). The vectors, sequenced to confirm the construction from positive clones, were transformed into A. tumefaciens strain EHA105. The sequences of all primers are provided in Supplemental Table S5. The leaves of two or three independent lines in the T0 or T1 generation were observed for the gene expression pattern localizations based on YFP fluorescence emission. Cell outlines were revealed by staining with PI. Images were acquired using a Nikon A1 confocal microscope with 488 and 561-nm argon laser excitation.

Genotyping, phenotyping under drought-stress, and CO₂ experiments of the ZmMPK12 mutant

The Ethyl Methanesulfonate (EMS) mutant ZmMPK12 (Zm00001d024568, stock ID EMS4-0e151f) was ordered from the MEMD project (http://www.elabcaas.cn/memd/). The mutant Single-Nucleotide Polymorphism (SNP) allele was identified by polymerase chain reaction (PCR) sequencing with the primer pair (F: TCAATGCATCATGTTCATGCAAG/R: CTTCAACGACGTCTACATCGTCTAC). One homozygous mpk12 allele (−/−) plant, and one heterozygous MAPK12 allele (−/+) plant were acquired from the progeny of EMS4-0e151f. To avoid influences of other mutant sites, we planted the progeny of the heterozygous allele (−/+) and harvested the seeds. Individuals representing the F₂ segregation population were grown in plastic pots (6 cm in diameter). For drought-stress experiments, each pot contained the same quantity of sieved soil, and was planted with four seeds, with the genotype of each plant being verified by PCR. We stopped watering after the seedlings of the ZmMPK12 mutants and the WT grew out of the soil surface, and observed the phenotypes after 10 days. The fresh weight of the whole plant was recorded for both the WT and the mutant under well-watered and drought conditions. The detached third leaves of 16-day-old seedlings were employed for measurement of water loss. Water loss was recorded every 20 min for 240 min at 22°C at 65% humidity using three biological replicates, each comprising two individuals. The third leaves of 16-day-old seedlings were taken for measurement of transpiration rate and stomatal conductance under different CO₂ concentration using a portable photosynthesis system LI-6800 (LI-COR Inc., Lincoln, NE, USA). For CO₂ experiments, maize seedling leaves were stabilized at 400 p.p.m CO₂, then stomatal conductance was recorded for an additional 40 min at 400 p.p.m, then at 1,500 p.p.m for 40 min, then recovered to 400 p.p.m for 40 min, and finally 50 p.p.m for 40 min.

TEVC recording and ABA content measurement in X. laevis oocytes

The method and solution for NO₃⁻ current recording followed (Huang et al., 1999). The method and solution for Zmβ-CA1 activated anion current followed (Tian et al., 2015). For ABA content measurements, each oocyte was injected with 12.1-ng ABA. After incubation for 5 h, the oocytes were collected, and the ABA content was measured via ultra-high-performance liquid chromatography–tandem mass spectrometry (UPLC–MS/MS).

Collection of GCs and SCs by LCM and measurement of ABA and its catabolites, and of ACC

Peels were prepared as for those used in snRNA-seq and were spread out flat on a polyethylene tetraphthalate membrane-covered slide (MembraneSlide 1.0, Zeiss), dried by blotting paper, frozen in liquid nitrogen, and immediately desiccated in a ScanVac CoolSafe freeze drier. The specimen slide carrying the freeze-dried peels was used for LCM with a Zeiss PALM CombiSystem. SCs and the pairs of GCs were collected separately. The cell numbers were calculated from fluorescence microscope images based on the autofluorescence of the cell wall after suspension in 75% methanol.

Around 3,000–4,000 GCs and 600–800 SCs were obtained for each sample, and were suspended in 2-mL tubes with 500 mL of extraction solution (75% methanol in water). After ultrasonic homogenization for 2 min and vortexing for 10 min, the samples were kept in dark for 16 h, dried under vacuum, resuspended in 50 μL 40% methanol in water, and centrifuged (13,523 g, 4°C, for 50 min). ABA and its catabolites, PA and DPA, and also ACC were quantified by UPLC–MS/MS using an Xevo TQ-XS system (Waters, Milford, MA, USA) equipped with an ESI ion source. Chromatographic separation employed an ACQUITY UPLC HSS T3 column (2.1 × 100 mm, 1.8 µm) maintained at 40°C. The flow rate of the mobile phase, composed of solvent A (0.1% formic acid, water) and solvent B (methanol), was set to 0.3 mL min⁻¹. The linear-gradient system was set as following: 0–2 min, 2% B; 2–10 min, to 80% B; 10–12 min, 80% B; 12–13 min, to 2% B; 13–15 min, 2% B. The autosampler temperature was set to 4°C, and the sample injection volume was 10 μL. The data were collected under the negative ion mode, with multiple reaction monitoring. Precursor and fragment ions are: ABA (m/z 263.16–153.01), d₆-ABA (m/z 269.20–159.10), PA (m/z 279.10–139.01), d₃-PA (m/z 282.10–142.01), DPA (m/z 281.00–171.00), d₃-DPA (m/z 284.00–174.00), ACC (m/z 102.10–56.10), d₄-ACC (m/z 106.07–60.10). Data were analyzed using the spectrometer software (Masslynx version 4.2).

Measurement of cell volume and calculation of ABA concentrations in SCs and GCs

A total of 25 stomata were imaged in peels using a multiphoton Confocal Microscope A1 MP (Nikon). Cell wall autofluorescence was excited (2P) at 750 nm, with emission being monitored at 405 nm. The cell volume was inferred by 3D reconstruction from z-stacks with a manual definition of the cell border, using Imaris (Oxford Instruments; Jonkman et al., 2020). ABA concentrations in SCs and GCs were calculated based on the total cell number and the volume ratio of the median SC volume to the median GC volume (7.86).

Data availability

To advance the data mining in functional studies, we have constructed a web server to allow users to monitor the cell type-specific expression patterns of target genes in the peel and leaf base samples (http://songlab.henu.edu.cn:3838/2021maizestomata/). The feature, violin, and dot plots can be freely browsed and downloaded in a user-friendly way. The scRNA-seq data and bulk RNA-seq data have been deposited in the GenBank SRA database as BioProject PRJNA740934. All flow cytometric data (FCS files) and maximum-likelihood phylogenetic trees with alignment files are available on Github (https://github.com/MingzXia/mays).

Accession numbers

All the gene information mentioned in the article is available in the Supplemental Data Sets 4, 5, and 8.

Supplemental data

The following materials are available in the online version of this article.

Supplemental Figure S1. Protocols of sample collection and processing for snRNA-seq based on the 10× Genomics system in maize (Supports Figure 1).

Supplemental Figure S2. Evaluation and clustering of the three peel samples (Supports Figure 1).

Supplemental Figure S3. The biosynthetic pathway of anthocyanins in peels (Supports Figure 1).

Supplemental Figure S4. Independent clustering and association analysis of the three peel samples (Supports Figure 1).

Supplemental Figure S5. Fluorescence observation of the marker genes in the mature stomata of transgenic rice and maize (Supports Figure 2).

Supplemental Figure S6. C4 pathway and starch metabolism in MCs, BSCs, and SCs from peels (Supports Figure 2).

Supplemental Figure S7. The 3D image of a stomatal complex and capture of GCs and SCs via laser microdissection (Supports Figure 3).

Supplemental Figure S8. Genotyping and phenotyping of zmmpk12 mutants under drought stress (Supports Figure 3).

Supplemental Figure S9. Leaf base sample preparation for snRNA-seq and integration analysis with the scRNA-seq datasets from maize SAMs plus the six most recently initiated leaf primordia (SAM + P6; Supports Figure 4).

Supplemental Figure S10. Feature plot of marker genes shared by Cluster 11 of the leaf base sample and GCs/SCs of the peel samples (Supports Figure 4).

Supplemental Figure S11. Feature plot of marker genes in maize developing stomata from the leaf base snRNA-seq sample (Supports Figure 4).

Supplemental Figure S12. The models of transition in cell states or types inferred by scVelo (Supports Figure 4).

Supplemental Figure S13. Regulatory network during the maize GMC-YGC development (Supports Figure 5).

Supplemental Figure S14. The gene model of cell wall-related marker genes during stomatal development (Supports Figure 5).

Supplemental Figure S15. Marker genes shared by bulk RNA-seq data and the leaf base snRNA-seq sample (Supports Figure 5).

Supplemental Table S1. Quality control measures and statistics of three maize peel samples from snRNA-seq.

Supplemental Table S2. Measurement of ABA and ACC levels in GCs and SCs collected by LCM.

Supplemental Table S3. Volume calculations for SCs and GCs from 25 stomatal complexes using Imaris.

Supplemental Table S4. Marker genes likely involved in cell wall metabolism of stomatal development from the leaf base sample.

Supplemental Table S5. Information of marker genes used in transformation, the EMS mutant, and electrophysiological experiments.

Supplemental Data Set 1. DEGs involved in the biosynthesis and regulation of anthocyanin in the anthocyanin-accumulating cluster (PCs-A).

Supplemental Data Set 2. DEGs present in all six clusters of the three peel samples by snRNA-seq.

Supplemental Data Set 3. Marker genes of the SCs and GCs from the peel samples and developmental stomata from the leaf base sample.

Supplemental Data Set 4. Marker genes involved in photosynthesis, starch, and disaccharide metabolism in MCs, BSCs, and GCs of the peels.

Supplemental Data Set 5. Marker genes are ikely involved in stomatal movement in the GCs and SCs from the peel samples.

Supplemental Data Set 6. DEGs present in the 14 clusters of the leaf base sample by snRNA-seq.

Supplemental Data Set 7. DEGs from pairwise comparisons of the WTs and ZmMUTE mutants by bulk RNA-seq analysis.

Supplemental Data Set 8. DEGs shared by the cluster representing developing stomata in snRNA-seq and the three pairwise comparisons of WT to the bzu2-1 mutant in bulk RNA-seq.

These authors contributed equally (G.S., M.X., J.L., and W.M.).

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (https://academic.oup.com/plcell) is: Chun-Peng Song ([email protected]).

Acknowledgments

We thank Weidong Chen for assistance in nuclear sorting, Lianlian Wang for technical assistance in generating the confocal images, Yu Zhang and Pengyu Wang for suggestion as to how to capture SCs and GCs by LCM, Changqing Zhang for valuable suggestions on molecular experiments, Yun Feng and Chen Qi from the Center for Biological Imaging (CBI), Institute of Biophysics, CAS, for their help in analyzing the confocal images using Imaris, and Rainer Hedrich at the University of Wuerzburg for kindly providing the AtOST1 and AtSLAC1 vectors, the Public Technology Service Center of Xishuangbanna Tropical Botanical Garden (XTBG), CAS for providing the HPC platform.

Funding

This work is supported in part by the National Science Foundation of China (NSFC) grants to GS (31771414) and to CPS (U21A20206), the Program for Innovative Research Team (in Science and Technology) in University of Henan Province (21IRTSTHN019), the Department of Science and Technology of Henan Province to GS (18HASTIT041), and the 111 Project of China (D16014). DG also acknowledges on-going support from the USDA through the University of Arizona College of Agriculture and Life Sciences.

Conflict of interest statement. The authors do not declare any competing interests.

Author notes

†

Senior authors.