Co-option and neofunctionalization of stomatal executors for defence against herbivores in Brassicales

Main

During eukaryotic evolution, the co-option of existing gene regulatory networks (GRNs) controlled by transcription factors (TFs) has been associated with the acquisition of new cell types, tissues and organs^1,2,3,4. Basic helix-loop-helix (bHLH)-type TFs play key roles in cell fate determination and cell differentiation during eukaryotic development^5,6,7.

In plants, three sister bHLH TFs (SPEECHLESS [SPCH], MUTE and FAMA) promote the differentiation of stomata in epidermal tissues^{8,9,10,11,12,13}. Stomata are valves through which CO₂ passes, placing them at the centre of the global carbon cycle; stomata are composed of a pair of specialized cells, the guard cells (GCs)^{8,9,10,11,12,13}. SPCH, MUTE and FAMA form heterodimers with other bHLH-type TFs, namely, SCREAM (SCRM, also reported as INDUCER OF CBF EXPRESSION 1 (ICE1)) and SCRM2, and control the transition from protodermal cells to meristemoids, from meristemoids to guard mother cells (GMCs), or from GMCs to GCs^14,15,16,17. Recently, it was reported that MUTE–SCRM and FAMA–SCRM heterodimers have the potential for pioneering activity (that is, as pioneer TFs), binding their target genes even when ensconced in closed or open local chromatin, and they initiate remodelling of the epigenome¹⁸. Although early stomatal development is well understood, GC differentiation remains enigmatic¹⁹. Indeed, how FAMA regulates the differentiation of GCs through its direct targets is largely unknown.

Two TF genes, ETHYLENE-RESPONSE FACTOR 51 (ERF51, also reported as DEHYDRATION-RESPONSIVE ELEMENT-BINDING PROTEIN 2F (DREB2F))^20,21 and STOMATAL CARPENTER 1 (SCAP1)^20,22,23, are potential downstream factors of FAMA because the expression of ERF51 and SCAP1 is upregulated in an oestrogen-inducible FAMA overexpression line²⁰. In addition, ERF51 is expressed in the stomatal lineages and FAMA directly binds to the WSB promoter region in FAMA overexpression lines²⁰. However, a previous study reported that an erf49 erf50 erf51 erf52 quadruple mutant did not exhibit highly penetrant stomatal phenotypes²⁰. Therefore, the physiological function of ERF51 in stomatal development remains largely unknown. Independent of the study by ref. ²⁰, SCAP1 was isolated as the causal mutated gene in a mutant impaired in CO₂-dependent changes of stomatal conductance and is specifically expressed in GCs²². Although a subset of scap1 GCs showed an abnormal morphology (the ventral cell walls of GCs appeared floppy and were often irregularly curved), the stomatal phenotypes of scap1 mutants were much weaker than those of fama mutants²². SCAP1 is probably a direct target of FAMA; however, other direct targets of FAMA may synergistically promote the differentiation of GCs together with SCAP1 (ref. ²⁴).

We and another group previously reported that FAMA regulates the differentiation of another type of specialized cells called idioblast myrosin cells (MCs), which are distributed along the vasculature in inner tissues^{25,26,27,28,29}. MCs play a critical role in a Brassicales-specific chemical defence system called the myrosinase-glucosinolate system, in which MCs are thought to protect the vasculature from herbivory attacks^27,28,29. When herbivores damage plant tissues, myrosinase (also named thioglucoside glucohydrolase) and its substrate, glucosinolate, are released from MCs and S cells (a glucosinolate-rich cell type), respectively, and they react to produce the toxic compounds isothiocyanates^27,28,29. Humans also detect isothiocyanates as pungency in some foods (such as wasabi and horseradish)³⁰. Despite the importance of the myrosinase-glucosinolate system as a plant defence mechanism, details about the downstream targets of FAMA involved in the development of MCs, and how Brassicales plants acquired this defence system (MCs and S cells) during evolution are largely unknown^27,28,29.

We previously hypothesized that FAMA-mediated GRNs were co-opted from GCs to MCs during evolution, as GCs are one of the oldest innovations of land plants, whereas MCs are a Brassicales innovation^29,31,32. However, the precise components of the GRNs regulating GC and MC development via FAMA, and how each cell identity is differentially regulated are largely unknown, as are the key genes directly regulated by FAMA.

Here we describe two direct FAMA target genes that encode TFs: ERF51, which we renamed WASABI MAKER (WSB) after its functional analysis, and SCAP1, in Arabidopsis thaliana (hereafter, Arabidopsis). Genetic and biochemical analyses revealed that FAMA deploys WSB and SCAP1 for GC differentiation, and WSB specifically for MC differentiation. We propose that the FAMA–WSB–SCAP1 GRN had an established function in the differentiation of GCs before being co-opted and reduced to the FAMA–WSB GRN for the differentiation of MCs during evolution, reflecting the neofunctionalization of FAMA and WSB. This study not only advances our understanding of how specialized cells arise during evolution, but also provides insight into the co-option and neofunctionalization of GRNs during evolution.

Results

Comprehensive atlas of genes regulated by pioneer TF FAMA

To obtain a nearly complete map of transcriptional changes driven by FAMA, we employed an oestrogen-inducible FAMA overexpression line (iFAMA for induced FAMA expression), which was previously used to show that FAMA is sufficient for the differentiation of MCs and GCs^14,20,26. Oestrogen treatment was reported to suffice for GC induction^14,20. Indeed, after 8 h of oestrogen treatment, we detected ectopic expression of the MC marker gene MYR001 (ref. ³³) (consisting of the VACUOLAR SORTING RECEPTOR HOMOLOG 1 (VSR1) promoter driving the reporter gene β-GLUCURONIDASE (GUS)) in iFAMA outside the area adjacent to the vasculature where MCs normally develop. This result suggested that MC characteristics were conferred to mesophyll cells by FAMA expression (Fig. 1a and Extended Data Fig. 1a) and indicated that 8 h of FAMA induction is sufficient to trigger FAMA-mediated downstream transcriptional changes that promote MC differentiation. A time-course expression analysis of oestrogen-treated iFAMA seedlings revealed a rapid upregulation of FAMA transcript levels relative to seedlings mock treated with diluted ethanol (−oestrogen), by >25-fold within 8 h and >100-fold after 24 h (Extended Data Fig. 1b). Transcriptome deep sequencing (RNA-seq) analysis of iFAMA mock-treated or oestrogen-treated seedlings identified 161 and 248 upregulated and 418 and 189 downregulated genes in iFAMA after 8 and 24 h of oestrogen treatment, respectively (hereafter referred to as iFAMA 8 h UP or DOWN and iFAMA 24 h UP or DOWN) compared with mock-treated iFAMA seedlings (with a false discovery rate (FDR) < 0.05) (Fig. 1b, and Supplementary Data 1 and 2). At the 24-h time point, more genes were upregulated and to a much greater extent than at 8 h (Fig. 1b). The Gene Ontology (GO) categories enriched for iFAMA 8 h UP genes were overwhelmingly related to the cell wall, suggesting that FAMA regulates cell wall modifications during the differentiation of GCs and MCs (Fig. 1c purple and Supplementary Data 3). Highly enriched categories also included the term ‘stomatal complex morphogenesis’ (Fig. 1c green and Supplementary Data 3). Strikingly, two genes associated with this category, SHAVEN3 (SHV3) and SHV3–LIKE 1 (SVL1), are redundantly required for GC morphogenesis³⁴. Thus, we successfully captured the transcriptional changes driven by FAMA in iFAMA 8 h UP, which might include key genes for the differentiation of MCs and GCs.

Co-option and neofunctionalization of stomatal executors for defence against herbivores in Brassicales — **Fig. 1: The TF gene *WSB* is a direct target of FAMA based on combinational RNA-seq analyses.**

Identification of the TF gene WSB, a primary target of FAMA

To identify key genes for the differentiation of MCs and GCs from the 161 genes regulated by FAMA, we focused on SCRMs, which interact with FAMA to promote the differentiation of MCs and GCs^17,26. The fama and scrm scrm2 mutants exhibited severe dwarfism phenotypes^14,17. We reasoned that it might be difficult to collect tissue samples with similar shapes, including similar cell types, from the wild type (WT) and these mutants due to this large size difference; an RNA-seq analysis using these mutants would probably mainly reflect the differences in cell types in the collected samples rather than a specific differential gene expression caused by loss of FAMA–SCRM complexes. Therefore, we used scrm single mutants which did not exhibit severe dwarfism and showed the low accumulation levels of the myrosinase, THIOGLUCOSIDE GLUCOHYDROLASE 1 (TGG1) in leaves and stems^17,26. In addition, we performed RNA-seq analysis on stems, as they offered samples with no morphological differences between WT and scrm although scrm leaves were slightly smaller than those of WT, thus allowing us to collect samples at similar developmental stages and similar size. The loss of SCRM function was associated with the downregulation of 154 genes (scrm DOWN), including known marker genes for both GCs and MCs (Extended Data Fig. 1c and Supplementary Data 4). To identify the downstream targets of FAMA that are expressed in both GCs and MCs, we included two additional published transcriptome datasets: (1) genes upregulated in the semi-dominant mutant scrm-D (scrm-D UP), in which mature stomata are overproduced^17,35; and (2) genes upregulated in the syntaxin of plants 22 (syp22) mutant (syp22 UP), in which mature MCs are overproduced^36,37,38,39 (Supplementary Data 5).

The overlap between these four datasets comprised only two genes with iFAMA 8 h UP: FAMA itself and At3g57600 (Fig. 1d and Supplementary Data 5). We also detected only four genes overlapping with iFAMA 24 h UP: FAMA, At3g57600, TGG1 (refs. ^36,40) and HEAVY METAL ASSOCIATED ISOPRENYLATED PLANT PROTEIN 20 (HIPP20)⁴¹ (Fig. 1e and Supplementary Data 5). We selected At3g57600, which we named WASABI MAKER (WSB) after functional analysis (see below), as it was upregulated before HIPP20 and the myrosinase genes TGG1 and TGG2 (Fig. 1f), which are expressed in both GCs and MCs (Extended Data Fig. 1d)⁴⁰. Moreover, WSB was the only gene co-expressed with FAMA, TGG1 and TGG2 in the ATTED-II database⁴² (Fig. 1g); the overlap between the four iFAMA transcriptome datasets from this study and the transcriptome dataset from ref. ²⁰ consisted of only two genes, FAMA and WSB (Supplementary Data 6). WSB encodes an APETALA2/ERF (AP2/ERF)-type TF belonging to Group IV of this family and known as ERF51 and DREB2F, with no known biological function⁴³.

WSB is a direct target of the FAMA–SCRM complexes

To examine whether FAMA is required for and directly regulates WSB expression, we determined the expression levels of WSB in iFAMA in response to oestrogen treatment using reverse-transcription quantitative PCR (RT–qPCR). Consistent with the RNA-seq data, WSB expression was rapidly induced after 8 h of oestrogen treatment and reached a >7-fold increase after 24 h (Fig. 2a). Conversely, we failed to detect WSB transcripts in fama mutants (Fig. 2b) and detected no GUS signal from the leaves of fama proWSB:GUS seedlings (Extended Data Fig. 3d). These results suggest that FAMA is necessary and sufficient for WSB expression.

**Fig. 2: *WSB* is expressed in MCs and GCs.**

FAMA was reported to directly bind to the WSB promoter region in chromatin immunoprecipitation–quantitative PCR (ChIP–qPCR) assays using FAMA overexpression lines²⁰. We thus asked whether FAMA directly binds to the WSB promoter region in GCs and MCs by performing a ChIP–qPCR analysis using a complementation line, fama proFAMA:FAMA-myc⁴⁴, in which MC-specific reporters and GC-specific reporters were expressed properly (Extended Data Fig. 2a,b), with an anti-myc antibody (Fig. 2c). Indeed, we detected a significant increase in FAMA binding to a WSB promoter region near a CTCGTG motif (-672 to -667 bp relative to the ATG) that is identical to the cognate cis element of the FAMA sister protein SPCH⁴⁵ (Fig. 2c and Extended Data Fig. 2c). The value of % input (binding strength) to the WSB promoter was higher than that to the EPIDERMAL PATTERNING FACTOR 1 (EPF1) locus to which binding of FAMA was previously shown⁴⁴. These results suggest that FAMA directly binds to the WSB promoter in GCs and MCs. WSB expression was lower in scrm mutants (Extended Data Fig. 1c) and almost undetectable in the scrm scrm2 double mutant (Extended Data Fig. 2d). In addition, SCRM-GFP (a fusion of SCRM to the green fluorescent protein)¹⁷ bound to the same region in the WSB promoter (Extended Data Fig. 2e). These results confirm that FAMA directly activates WSB expression (Fig. 2 and Extended Data Fig. 3d). Moreover, we determined that SCRM also directly activates WSB expression (Extended Data Figs. 1c and 2d,e). Collectively, these results suggest that FAMA–SCRM complexes are critical transcriptional activators for WSB.

WSB is highly expressed in both MC and GC lineages

We investigated the spatial expression pattern of WSB in planta by generating transgenic lines expressing the GUS reporter gene under the control of a 1.4-kb WSB promoter fragment (proWSB:GUS). We observed GUS staining in elongated cells with horn-like extensions, that is, MCs, along the vasculature in the inner layer of leaves (Fig. 2d left); in small round cells, corresponding to GMCs; and in GCs in the leaf epidermis (Fig. 2d right), suggesting that WSB is expressed throughout leaf development in these specific cell types (Extended Data Fig. 3a). Consistent with the RNA-seq data (Fig. 1d,e), we noticed more GUS-positive cells in the epidermis of scrm-D proWSB:GUS seedlings (Extended Data Fig. 3b) and the inner layer of leaf tissue from syp22 proWSB:GUS seedlings than in the WT (Extended Data Fig. 3c).

We also investigated the accumulation of WSB with a translational fusion construct introduced in a WSB mutant (wsb-1; see below), wsb-1 gWSB-mVenus. In at least 10 independent transgenic lines, the gWSB-mVenus reporter (encoding a fusion between WSB and yellow fluorescent protein) was highly expressed in patches in leaf inner tissues and the epidermis, based on mVenus fluorescence patterns. We detected the accumulation of WSB in the nucleus of myrosin-lineage cells in inner tissues (Fig. 2e right and Extended Data Fig. 3e), as well as GMCs and young GCs in the epidermis (Fig. 2e left). We also observed weak WSB-mVenus signals in the cytoplasm of MCs and GCs.

The myrosinase gene TGG2 is a marker for mature GCs and MCs. To clarify whether WSB is expressed in differentiating GCs and MCs, we compared the expression patterns of WSB and TGG2. Accordingly, we crossed the proTGG2:Venus-2s⁴⁶ (encoding Venus targeted to vacuoles) transgenic line to wsb-1 and transformed the resulting line with the gWSB-mTurquoise2 construct (encoding WSB fused to cyan fluorescent protein⁴⁷). We detected both fluorescent proteins in developing MCs and GCs, marking sites of Venus-2sc and WSB-mTurquoise2 accumulation (Fig. 2f). In the stomatal lineage, GMCs accumulated only WSB-mTurquoise2 (Fig. 2g). GCs just after division and young GCs presented both fluorescent markers, while mature GCs only accumulated Venus-2sc driven by the TGG2 promoter. These results suggest that WSB is expressed before TGG2 and that WSB specifically localizes to the nuclei of MCs, GMCs and young GCs. These results indicate that the expression of WSB follows the same pattern as that of its upstream regulator, FAMA, except in mature GCs^14,25,26 (Extended Data Fig. 3f).

Generation of WSB knockout mutants by CRISPR/Cas9

To clarify the biological functions of WSB in vivo, we generated loss-of-function mutants for WSB via clustered regularly interspaced short palindromic repeats (CRISPR)/CRISPR-associated protein 9 (Cas9)-mediated genome editing⁴⁸, as no T-DNA insertion mutants with insertions in the exon are available. We obtained nine mutants (wsb-1 to wsb-9) with frameshifts introducing premature stop codons (Fig. 3a,b and Extended Data Fig. 4a). All wsb mutants except wsb-7 and wsb-8 lacked a part of the sequence encoding the DNA-binding domain (the AP2 domain) (Fig. 3b and Extended Data Fig. 4a). Plant sizes of wsb mutants were comparable to those of the WT (Fig. 3c and Extended Data Fig. 4b,c).

**Fig. 3: WSB is required for the development of MCs.**

WSB is required for the differentiation of MCs

To examine whether WSB, similar to FAMA, is required for the development of MCs, we characterized MC development in wsb mutants by RT–qPCR analysis for the endogenous MC markers TGG1 and TGG2. We detected no transcripts for these genes (Fig. 3d and Extended Data Fig. 5a–c) and much lower TGG1 and TGG2 abundance in wsb mutants by immunoblot analysis (Fig. 3e, Extended Data Fig. 5d–f and Supplementary Figs. 1–4). TGG1 and TGG2 transcripts were largely undetectable in the F₁ progenies of wsb-1/wsb-2 and wsb-1/wsb-3 (Extended Data Fig. 5g,h); however, the fluorescent translational fusion construct gWSB-mVenus fully rescued TGG1 and TGG2 expression levels when introduced into the wsb-1 mutant, indicating that the WSB-mVenus fusion is fully functional (Extended Data Fig. 5i,j). In addition, WSB showed a dosage-dependent effect on the expression of myrosinase genes, based on the analysis of WSB-knockdown lines expressing an artificial microRNA designed to decrease WSB transcript levels (Extended Data Fig. 5k,l). Collectively, WSB is essential for the expression of myrosinase genes. We also examined MC development in wsb-3 mutants with the MC-specific reporter MYR001:GUS^24,31. We detected barely any activity for the reporter along the vasculature of whole wsb-3 seedlings, leaves (cotyledons, true leaves and cauline leaves), stems and flowers (sepals, petals and flower stalks), whereas MCs developed normally in the WT, suggesting that MCs fail to differentiate in wsb (Fig. 3f,g and Extended Data Fig. 6a–f). The defects in MC development seen in wsb were reminiscent of those in fama^25,26, suggesting that WSB is a critical downstream target of FAMA for the development of MCs.

WSB governs transcriptional networks for both MCs and GCs

To explore the transcriptional networks regulated by WSB, we performed a comparative RNA-seq analysis of WT and wsb-1 seedlings (Fig. 3h green and Supplementary Data 7). We identified 383 downregulated genes in wsb, including TGG1 and HIPP20 (MC and GC markers) as well as TGG2 (MC marker), suggesting that WSB is required for their expression (Fig. 3h, Extended Data Fig. 6g and Supplementary Data 7). When compared to single-cell RNA-seq (scRNA-seq) data for MCs and GCs from mature leaves⁴⁹, 54 and 42 of the 383 downregulated genes were expressed specifically in MCs and GCs, respectively, with another 114 expressed in both MCs and GCs (Fig. 3h and Supplementary Data 8). These results suggest that WSB governs transcriptional networks for the differentiation of both MCs and GCs, in agreement with the WSB expression pattern.

To verify whether WSB can activate MC-specific genes in various cell types, among the 54 genes above, we focused on the TF gene WRKY23 (hereafter referred to as MYR002), as we had previously identified it as a FAMA downstream gene²⁶ and it was reported to be an MC-specific marker gene on the basis of scRNA-seq data⁵⁰. In addition, WRKY23 is one of the genes co-expressed with WSB in the ATTED-II database⁴² (Fig. 1g). We confirmed that MYR002 expression levels are much lower in wsb-1 (Fig. 3i), suggesting that WSB is required for MYR002 expression. To test the sufficiency of WSB for MYR002 expression, we generated an oestrogen-inducible WSB-mVenus line (hereafter referred to as iWSB-mVenus). We detected WSB-mVenus in iWSB-mVenus after 6 h of oestrogen treatment, with higher levels reached after 24 h of oestrogen treatment (Fig. 3j and Supplementary Figs. 5–7). The expression levels of MYR002 rapidly rose >3-fold after 2 h of continuous oestrogen treatment, reaching >60-fold higher levels at 24 h under oestrogen treatment compared with mock-treated iWSB-mVenus seedlings (Fig. 3k). We also examined GUS staining in iWSB-mVenus MYR002:GUS seedlings generated by crossing the two independent transgenic lines. In mock-treated seedlings, consistent with a previous report⁵⁰, we specifically observed GUS staining in MCs but not in GCs (Fig. 3l). In seedlings treated with oestrogen for 48 h, we detected GUS signals throughout the cotyledons, including pavement cells and mesophyll cells, as well as a strong signal in the vasculature (Fig. 3l top middle and right). Notably, we also observed GUS signals in GCs (Fig. 3l lower right). In addition, WSB induced the expression of other MC marker genes (HIPP20, At1g23170 and TETRASPANIN 4 (TET4, At5g60220)) (Fig. 3m–o and Extended Data Fig. 6h,i). Notably, WSB did not induce the expression of MYR001 (Supplementary Fig. 8). Taken together, these results suggest that WSB can activate the expression of a subset of MC marker genes in various cell types.

SCAP1 is a stomatal lineage-specific direct target of FAMA

WSB was expressed in stomatal-lineage cells (Fig. 2), and RNA-seq analysis revealed that loss of WSB function decreased the expression of 156 GC-related genes (Fig. 3h), suggesting that WSB may be involved in stomatal development. However, in our cultivation conditions, we did not see clear changes in stomatal morphology in wsb single mutants. We hypothesized that WSB and other direct FAMA target genes might function synergistically in stomatal development. We focused on the DOF-type TF gene SCAP1 on the basis of the following reasoning, although we did not identify SCAP1 as a differentially expressed gene in the RNA-seq analysis of iFAMA (Fig.1b). First, SCAP1 was downregulated in scrm (Extended Data Fig. 1c). Second, FAMA was reported to activate SCAP1 expression in different FAMA-inducible lines²⁰, and 50% of all stomata in scap1 mutants show skewed morphologies²². Third, SCAP1 is the only other co-expressed TF-encoding gene with FAMA with strong correlation (correlation value 1–5; bold lines in Fig. 1g) in the ATTED-II database⁴² together with WSB. As the inducibility of the oestrogen-inducible system varied among different transgenic lines, we generated another FAMA overexpression line (hereafter referred to as iFAMA^strong), in which we induced FAMA more strongly than in the original iFAMA line (Extended Data Fig. 7a, compare with Extended Data Fig. 1b). In iFAMA^strong, we observed FAMA overexpression phenotypes, as reported previously²⁰ (Supplementary Fig. 9). Using iFAMA^strong and consistent with the previous report²⁰, we determined that SCAP1 expression is indeed upregulated by FAMA (Fig. 4a). Conversely, SCAP1 expression levels were lower in the fama mutant (Fig. 4b). In addition, we established that FAMA directly binds to the SCAP1 promoter regions by ChIP–qPCR (Fig. 4c) using a primer pair that amplifies a SCAP1 promoter fragment containing two CTCGTG motifs (Fig. 4c and Extended Data Fig. 7b). A proSCAP1:GUS reporter line showed specific GUS staining in GCs (Fig. 4d and Extended Data Fig. 7b)²². By contrast, the proSCAP1_delta_67_bp:GUS reporter construct lacking the two CTCGTG motifs produced no GUS signal in GCs (Fig. 4d and Extended Data Fig. 7b). Consistent with these results, we noticed two cis elements, recently reported to be critical for SCAP1 expression²³, within the 67-bp DNA fragment deleted in proSCAP1_delta_67_bp:GUS (Extended Data Fig. 7b). We barely detected any SCAP1 expression in MCs, as evidenced by the very low GUS staining in proSCAP1:GUS seedlings and Venus fluorescence in gSCAP1-mVenus seedlings (Fig. 4d,f Extended Data Fig. 7b). We rarely found a signal for SCAP1-mVenus in the nucleus of large MCs (Supplementary Fig. 10). These results confirm that SCAP1 is a stomatal lineage-specific gene that is also a direct target of FAMA.

**Fig. 4: *SCAP1* is a stomatal lineage-specific direct target gene of FAMA.**

We hypothesized that WSB might inhibit SCAP1 expression in MCs because the expression of WSB persisted longer in MC lineages than in GCs (Fig. 2e–g and Extended Data Fig. 7c). Consistent with this idea, we frequently detected the ectopic expression of SCAP1 in MCs of WSB-knockdown lines (Extended Data Fig. 7d). Importantly, we did not notice any clear difference in the expression window of SCAP1 in stomatal lineages between the WT and WSB-knockdown lines (Supplementary Fig. 11). These results suggest that the sustained expression and/or high levels of WSB repressed the expression of SCAP1 specifically in MC lineages. In addition, we detected ectopic expression of SCAP1 in a subset of MCs in syp22 leaves, in which polar auxin transport (PAT) is perturbed³⁹, as indicated by a previous transcriptome analysis (Extended Data Fig. 7e)²⁶. These results suggest that WSB and PAT repress the expression of the guard cell-specific gene SCAP1 in the MC lineage, resulting in the fate specification of MCs.

WSB and SCAP1 are co-expressed in GCs just after division

We examined the accumulation of SCAP1 protein with a translational fusion construct introduced into the scap1 mutant (scap1 gSCAP1-mVenus). In at least 10 independent transgenic lines, we detected expression of the gSCAP1-mVenus reporter as patches in the epidermis, highlighting the accumulation of SCAP1 in the nuclei of GCs, but not of GMCs (Fig. 4e). In addition, we rarely detected SCAP1-mVenus in MCs (Fig. 4f). To clarify the expression window of WSB and SCAP1 in the stomatal lineage, we compared the sizes of GMCs and GCs accumulating each protein (Extended Data Fig. 7f). Stomatal-lineage cells expressing WSB-mVenus were on average smaller than those expressing SCAP1-mVenus, while a subset of GCs accumulated both proteins, suggesting that WSB has an earlier expression window than SCAP1, although they partially overlap. To explore this hypothesis in more detail, we generated two translational reporter lines, gWSB:mVenus and gSCAP1:mTurquoise2; we transformed wsb-1 gWSB-mVenus plants using Agrobacterium harbouring gSCAP1-mTurquoise2 constructs. Using this line, we visualized both WSB and SCAP1 fused to different fluorescent proteins in the same stomatal-lineage cells from GMCs to maturing GC stages (Fig. 4g,h and Supplementary Figs. 12–15) and quantified the signal intensities for WSB and SCAP1 (Fig. 4i and Supplementary Fig. 14), and cell area (Fig. 4j and Supplementary Fig. 14). From this analysis, we conclude that WSB and SCAP1 co-occur in the same GCs just after division (stages 6 and 7 in Fig. 4h–j and Supplementary Fig. 14) when guard cell area dramatically increases, suggesting that these two TFs contribute to the differentiation of GCs at the stage just after final division.

WSB and SCAP1 act synergistically to differentiate GCs

To assess whether WSB and SCAP1 act synergistically in the differentiation of GCs, we crossed previously isolated scap1 single mutants²² with wsb-1 to obtain the wsb scap1 double mutant. First, we examined stomatal development on the abaxial epidermis of the third true leaves of seedlings at 14 DAG. Unlike the WT and the single mutants, the wsb scap1 double mutant had very few normal GCs (Fig. 5a), with fewer than 3% compared with the WT (Fig. 5b,c and Extended Data Fig. 8a,b). Instead, wsb scap1 developed irregular GCs of four types: (1) GCs without a pore (unopened, blue), (2) three GCs in one stomatal unit (three cells, green), (3) a pair of stomata with direct contact (four cells, yellow) and (4) unorganized GCs, for example, four cells in a single stoma and a discontinuous circle of stomata (unorganized, pink) (Fig. 5a–c and Extended Data Fig. 8a,b). Furthermore, a fraction of unopen GCs had strong propidium iodide (PI) staining at the centre of the division plane (Extended Data Fig. 8c–f). However, in many cases, we did not find clear pores with strong accumulation of PI, suggesting that strong PI signals indicate wall thickening rather than pore formation. We counted the number of unopen GCs with/without wall thickening (Extended Data Fig. 8c–f). Over 60% of stomata in wsb scap1 were unopened GCs at this stage, with a diameter similar to that of WT GCs just after division, suggesting that unopened GCs in wsb scap1 are generated by the proper division of GMCs, then stall at the GC stage just after division, without forming a pore. Consistent with this scenario, the expression levels of two cyclin D genes, CYCD5 and CYCD7, which promote GMC division^51,52, were not clearly different and almost comparable to those in wsb scap1, the WT and the single mutants (Extended Data Fig. 8g). To rule out the possibility that stomatal phenotypes might be due to potential off-target effects during the generation of each single mutant using genome editing and chemicals, we introduced complementation constructs into the wsb scap1 double mutant: proWSB:WSB-mTurquoise2, proWSB:WSB-mVenus, proSCAP1:SCAP1-mVenus, proSCAP1:SCAP1-mCherry and proSCAP1:SCAP1-mTurquoise2. All complementation constructs restored stomatal phenotypes to WT levels (Fig. 5a–c, Extended Data Fig. 8 and Supplementary Fig. 12), validating the hypothesis that a functional WSB or SCAP1 is required to ensure normal stomatal development.

**Fig. 5: WSB and SCAP1 are synergistically required for the development of GCs.**

Next, to examine whether unopen GCs in wsb scap1 differentiate into normal GCs or abnormal GCs at later time points, we investigated stomatal development on the abaxial epidermis of third true leaves from WT and wsb scap1 seedlings at 14 and 21 DAG (Fig. 5d,e, Extended Data Fig. 8e and Supplementary Fig. 16). The frequency of unopen GCs (blue in Fig. 5d) was lower, while that of unorganized and irregular GCs (green, yellow and pink in Fig. 5d) was higher at 21 DAG than at 14 DAG. In addition, at 21 DAG, we observed more types of abnormal stomata in wsb scap1, some of which were not seen at 14 DAG (Fig. 5e). Some of the three-cell stomata and four-cell stomata seen in 21-DAG seedlings had pores (Fig. 5e,f). These results suggest that unopen GCs at 14 DAG might have the potential to differentiate, but they failed to properly differentiate into normal GCs, instead differentiating into irregular GCs at 21 DAG.

Finally, we compared GCs just after division and mature GCs between WT and wsb scap1 at a structural level by transmission electron microscopy (TEM) (Fig. 5f). Almost all organelles in GCs just after division appeared similar in the two genotypes. However, we noticed the accumulation of starch granules (bright white colour) in the chloroplasts of GCs at the stage just after division in the WT, while starch rarely accumulated in the chloroplasts of wsb scap1 GCs at the same stage, suggesting that starch biosynthesis is inhibited or that synthesized starch is quickly degraded in the chloroplasts of GCs following the loss of WSB and SCAP1 function. In addition, consistent with our confocal microscopy observations, we detected three-cell stomata and four-cell stomata in the TEM images (Fig. 5f). In addition, in some GCs of wsb scap1, the development of vacuoles was severely inhibited (Fig. 5f three-cell stomata). Overall, these results suggest that WSB and SCAP1 act synergistically in the differentiation of GCs.

Fewer GCs result in decreased CO₂ uptake and photosynthesis. Consistent with this scenario, the wsb scap1 double mutant displayed a dwarf phenotype (Supplementary Fig. 17). The degree of dwarfism was less severe than that of fama mutants (for example, wsb scap1 produced seeds, unlike fama) (Supplementary Fig. 17). This difference might be explained by the presence of a few normal stomata in wsb scap1 (Fig. 5a–c).

FAMA represses the expression of CYCD7 to inhibit unnecessary GC divisions. Therefore, fama tumours, which are undifferentiated stomatal-lineage cells with symmetric divisions in multiple places along the longest wall, are induced by the overexpression of CYCD7 (ref. ⁵²). We thus asked whether the abnormal GCs in wsb scap1 undergo additional divisions by introducing the proFAMA:CYCD7-mVenus reporter. We observed small cell clusters with random orientation for the division plane in the epidermis instead of unopen GCs, but these cell clusters were not identical to fama tumours in terms of the orientation of the cell divisions (Supplementary Fig. 18). These results suggest that CYCD7 can induce random divisions of stomatal-lineage cells in wsb scap1 instead of the symmetric divisions or division along the longest wall previously reported in fama.

Genome-wide mapping of WSB-binding sites

To identify the global binding events of WSB, we performed ChIP-seq for WSB using iWSB-mVenus and an anti-GFP antibody. We identified 3,199 WSB-binding sites corresponding to 3,034 loci (Fig. 6a and Supplementary Data 9). Consistent with the status of WSB as a TF, 66% of all WSB-binding peaks were associated with promoters, mostly within 1,000 bp upstream of the transcriptional start site (Fig. 6a,b and Supplementary Data 9). To categorize direct targets of WSB by expression pattern, we compared these data to scRNA-seq data for MCs and GCs from mature leaves⁴⁹, revealing that 267, 403 and 489 of the 3,034 WSB target genes are expressed specifically in MCs, GCs or both, respectively (Fig. 6c). To categorize WSB targets more precisely along stomatal development, we compared the lists of WSB targets to genes expressed from the GMC stage to the mature GC stage⁵³, identifying 81 WSB target genes expressed during this developmental gradient (Fig. 6d and Supplementary Data 10). These results indicate that WSB directly binds to the genomic regions of genes that are expressed in mature MCs and GCs, suggesting that WSB participates in the maturation of MCs and GCs.

We also identified SCRM and SCRM2 among the genes bound by WSB; their protein products form heterodimers with FAMA (Fig. 6e). We focused on SCRM because we previously showed that myrosinase genes are markedly less expressed in scrm, but not in scrm2 (ref. ²⁶). SCRM expression levels were lower in wsb and upregulated by WSB upon induction of iWSB-mVenus, suggesting positive feedback regulation between WSB and SCRM (Fig. 6f,g). As SCRM was required for FAMA expression (Extended Data Fig. 1c) and the FAMA–SCRM complex directly binds to the FAMA promoter¹⁸, we monitored FAMA expression in wsb mutant and WSB overexpression lines. FAMA transcript levels were lower in wsb and upregulated by WSB (Fig. 6h,i). In wsb proFAMA:GFP, we detected a GFP signal only in GCs just after division, indicating that WSB is required for continued FAMA expression in GCs (Fig. 6j left panels). We rarely detected a GFP signal in inner tissues of wsb proFAMA:GFP, suggesting that FAMA expression in MCs almost fully depends on WSB (Fig. 6j right panels). Finally, we detected the binding of WSB to the FAMA promoter under high WSB overexpression conditions (using iWSB-mVenus treated with 0.01 or 10 µM oestrogen for ChIP-seq and ChIP–qPCR for the FAMA locus, respectively) (Fig. 6k), suggesting that WSB may weakly bind to the FAMA promoter. Collectively, these results suggest the existence of a positive feedback regulation between WSB and the FAMA–SCRM complex in both MCs and GCs. This feedback is required for the sustained, strong expression of the pioneer TF gene, FAMA, in GCs and MCs.

We then searched for factors that may promote cell and nuclear expansion, two characteristics of idioblast MCs. WSB bound to the promoter of CELL CYCLE SWITCH PROTEIN 52 A1 (CCS52A1, or FIZZY-RELATED 2 (FZR2)), which promotes endoreduplication and nuclear and cell expansion (Fig. 6e)^54,55,56. In agreement with its role, CCS52A1 was upregulated by WSB (Fig. 6l) and was also strongly expressed in MCs (Fig. 6m and Extended Data Fig. 9a). Furthermore, CCS52A1 expression was downregulated in the wsb-1 mutant (Fig. 6n), suggesting that high levels of CCS52A1 increase the ploidy, nuclear size and cell size of MCs. These three characteristics were diminished in ccs52a1 mutants (Fig. 6o–q, Extended Data Fig. 9b,c and Supplementary Fig. 19). Collectively, these results suggest that the WSB–CCS52A1 module is required for MC differentiation.

WSB was co-opted for MC differentiation during evolution

Our study revealed that WSB is required for the differentiation of both MCs and GCs, and for the expression of the pioneer TF gene FAMA. GCs are one of the oldest innovations of plants as they adapted to land, as evidenced by their presence in bryophytes, mosses and hornwort^31,32. By contrast, MCs are present only in Brassicales plants²⁹. Therefore, we hypothesized that WSB might not represent a Brassicales-specific TF but might instead have first been acquired for GC differentiation before being co-opted to form the MC differentiation pathway. To explore the evolutionary path of WSB, we generated a maximum-likelihood tree of 461 Group IV AP2/ERF TFs from 62 plant species (Fig. 7a and Supplementary Fig. 20). This high-resolution analysis revealed that, rather than a single subclass, as was previously reported⁴³, Group IV AP2/ERF TFs form four distinct subclasses, I–IV. In Arabidopsis, the single-member WSB belongs to subclass II, supporting the notion that WSB has a unique biological function. Subclass II appears to be highly conserved in angiosperms, suggesting that plants acquired WSB before the emergence of Brassicales. Like other regulators of stomatal differentiation, WSB was lost secondarily in the seagrass Zostera marina, an astomatous plant (Fig. 7a and Supplementary Fig. 20)⁵⁷. By contrast, subclass I, III and IV members are conserved, even in seagrass (Fig. 7a and Supplementary Fig. 20). These results indicate that WSB orthologues are intimately linked to the evolution of stomata. Collectively, our evolutionary analysis suggests that WSB was co-opted for MC differentiation during evolution.

**Fig. 7: Hypothetical model for the co-option of the reduced transcriptional module, FAMA–*WSB*, during the evolution of Brassicales.**

The presence of MCs is a synapomorphic characteristic of Brassicales²⁷ (a characteristic derived from a common ancestor), which raises the question: When was WSB co-opted for MCs in the common ancestor of Brassicales? Synapomorphy posits that the MC-promoting function of WSB should be conserved across Brassicales plants. To assess this prediction, we compared the sequence of the DNA-binding domain in WSB and WSB-like proteins and performed complementation tests by expressing the coding sequences of WSB orthologues under the control of the Arabidopsis WSB promoter. The DNA-binding domain is indeed highly conserved among WSB orthologues from Brassicaceae (Arabidopsis and salt cress Thellungiella parvula), Cleomaceae (giant spider-flower Tarenaya hassleriana) and basal Brassicales plants, including Caricaceae (papaya Carica papaya) (Extended Data Fig. 10a). In addition, expressing TpWSB, ThWSB or CpWSB from the Arabidopsis WSB promoter rescued the loss of MCs seen in wsb (Fig. 7b), suggesting that the MC-promoting function of WSB is deeply conserved in Brassicales. We obtained the same results when driving the coding sequence for WSB orthologues in non-Brassicales from the Arabidopsis WSB promoter in wsb (Supplementary Fig. 21). Taken together, these results suggest that the co-option of WSB might have occurred in the common ancestor of Brassicales.

Next, to clarify the evolution of the FAMA–WSB module, we examined the induction of WSB by expressing the coding sequences of FAMA orthologues under the control of the Arabidopsis FAMA (AtFAMA) promoter in an Arabidopsis fama mutant. MpSeta, a bHLH Ia gene in the astomatous plant Marchantia polymorpha and the closest homologue to AtFAMA, AtMUTE and AtSPCH⁵⁸, did not induce WSB expression (Extended Data Fig. 10b). Notably, BdFAMA from purple false brome (Brachypodium distachyon) induced WSB expression in the epidermis but not in inner tissues (Fig. 7c and Extended Data Fig. 10c). These results can be explained by a previous report stating that BdFAMA driven by the AtFAMA promoter in fama is expressed in the epidermal layer but not in inner tissues⁵⁹. Moreover, the expression of AtFAMA in inner tissues requires positive feedback regulation by AtFAMA itself^25,26. BdFAMA may not activate the AtFAMA promoter, resulting in the loss of WSB in inner tissues. These results suggest that the positive feedback-dependent expression of FAMA in inner tissues was acquired during evolution.

As FAMA directly activated WSB expression (Fig. 2), we explored the extent of conservation at WSB promoter regions in Brassicales and observed high conservation of a ~500-bp fragment in the WSB promoter regions across Brassicales plants, which includes both the FAMA-binding site and CTCGTG motifs (Extended Data Fig. 10d,e)^60,61. In addition, synteny analysis suggested that the arrangement of the WSB locus in its genomic region is conserved in Brassicales (Extended Data Fig. 10f and Supplementary Data 11). Together, these results suggest that the FAMA–WSB transcriptional module may be conserved in this plant order. Overall, we propose that Brassicales plant species co-opted the FAMA–WSB module for the evolution of the myrosinase-glucosinolate defence strategy.

Discussion

Despite the importance of GCs for gas exchange and MCs for plant defence, no component critical for MC and GC differentiation has been identified under the control of FAMA. Here we rediscovered WSB and identified SCAP1 as direct targets of FAMA (Figs. 1, 2 and 4). WSB is expressed in both MCs and GCs, while SCAP1 is specifically expressed in GCs (Figs. 2 and 4). WSB is required for the expression of myrosinase genes (Fig. 3). wsb mutants lack MCs, while the wsb scap1 double mutant has few normal GCs (Figs. 3 and 5). In addition, ChIP-seq analysis revealed that WSB governs the transcriptional networks controlling the differentiation of both MCs and CGs, with WSB forming a positive feedback circuit with the FAMA–SCRM complex, while the WSB–CCS52A1 module is a key module specific to the differentiation of MCs (Fig. 6). Finally, an evolutionary analysis revealed that WSB is conserved across stomatous angiosperms (Fig. 7a). It was also previously reported that MCs are Brassicales-specific idioblast cells²⁹. From these results, we propose that the conserved and reduced FAMA–WSB transcriptional module was co-opted for MC differentiation during evolution (Fig. 7d). Below, we discuss the conserved stomatous transcriptional module, FAMA–WSB–SCAP1, and the co-option and neofunctionalization of the stomatal executors, FAMA and WSB, for defence against herbivores in Brassicales.

To our knowledge, wsb is the first mutant to be reported in which no MCs develop without being accompanied by dwarfism (Fig. 3). By contrast, fama mutants lack MCs but also showed strong dwarfism because of their absence of stomata, thus limiting CO₂ fixation. In wsb, in addition to TGGs, 168 MC-expressed genes (for example, HIPP20 and WRKY23) whose functions in MC are unknown were downregulated (Fig. 3h). Therefore, wsb is a useful tool to explore as yet undiscovered functions of MCs besides their role in plant defence. MCs differentiate from ground meristem cells, which are the stem cells in leaf primordia. As there is no good reporter gene for these ground meristem cells, it is unclear when and how the differentiation of MCs might stop or stall after the ground meristem cell stage in wsb or fama mutants. Notably, in our study, MCs in WSB-knockdown lines, in which WSB function and MCs were not completely lost, ectopically expressed SCAP1 (Extended Data Fig. 7d). The characterization of these abnormal MCs in WSB-knockdown lines may shed light on the origin, trajectory and cell fate bifurcation of MC differentiation. Recently, it was reported that MCs exist in Arabidopsis roots during secondary growth and might differentiate from phloem parenchyma cells and/or companion cells⁶². It is an open question whether SCAP1-positive cells in WSB-knockdown lines are phloem parenchyma cells (specifically, companion cells) in leaves, even though the development of vascular tissues is different between leaves and roots. In addition, scRNA-seq or single-nucleus RNA sequencing analysis of leaf primordia, including MCs and phloem lineages, may provide useful information about the cell lineages and upstream factors of FAMA in inner tissues surrounding future vascular tissues. From the viewpoint of applied science, WSB may be a good target for gene editing to generate Brassicales vegetables/crops without pungency, as gene editing of multiple homologous genes, such as the TGGs, is time-consuming.

In this study, we showed that WSB and SCAP1 are synergistically required for the differentiation of GCs. It was previously reported that a quadruple mutant defective in ERF49–52 (ERF49, ERF50, WSB (also known as ERF51 and DREB2F) and ERF52) did not exhibit highly penetrant stomatal phenotypes²⁰. Taken together, these results suggest that WSB did not function redundantly with homologous Group IV AP2/ERF TFs but rather synergistically with SCAP1, which belongs to the DOF family of TFs. Consistent with this idea, our phylogenetic analysis revealed that only WSB belongs to subclass II of the Group IV AP2/ERF TFs (Fig. 7a). These results indicate that FAMA executes GC differentiation by deploying the TFs WSB and SCAP1 (Fig. 7d).

Consistent with the observation that SCAP1 expression is sustained in mature GCs, beyond its role in cell differentiation, SCAP1 is required for stomatal movement through the direct induction of GATED OUTWARDLY-RECTIFYING K⁺ CHANNEL (GORK) and MYB60 (ref. ²²). In agreement with the narrow expression window of WSB in GMCs and GCs just after division, light- or CO₂-induced stomatal movement in wsb mutants was comparable to that in the WT (Extended Data Fig. 10g–i), suggesting functional differences between SCAP1 and WSB. Thus, whether WSB and SCAP1 share downstream targets and/or have specific targets are open questions. The transcriptional cascade starting from FAMA that includes CYCD7, WSB and SCAP1 is a good model in which to study the timely expression of developmental regulators in plants. FAMA represses CYCD7 and activates WSB and SCAP1, sequentially and with different timing (Fig. 7d). Notably, both WSB and SCAP1 are direct targets of FAMA; however, SCAP1 was slower to respond than WSB (Figs. 2 and 4). In addition, it was recently reported that ectopic induction of FAMA in embryos induced the ectopic expression of WSB but not of SCAP1 (ref. ²¹). WSB did not suppress SCAP1 expression in the stomatal lineages (Supplementary Fig. 11). These results suggest that unknown genetic and/or epigenetic barriers may exist that modulate the induction of SCAP1 by FAMA in the stomatal lineage. We do not exclude the possibility that additional transcriptional activators are required for the induction of WSB and SCAP1 in addition to FAMA. Consistent with this prediction, the repressive histone mark trimethylation of lysine 27 on histone H3 may modulate SCAP1 expression⁶³. To induce the expression of SCAP1 by FAMA, the previous dilution of repressive histone marks may be needed. Further studies are required to clarify the role of epigenetic regulation at the SCAP1 locus in the stomatal lineage. A set of 8,484 direct target genes were recently identified for FAMA¹⁸; therefore, other direct FAMA targets may play critical roles in other aspects of stomatal development and function.

In addition to promoting GC differentiation and fine-tuning GC divisions, a third role for FAMA in GC development is fate maintenance. FAMA interacts with RETINOBLASTOMA-RELATED (RBR) and HISTONE ACETYLTRANSFERASE OF THE CBP FAMILY 1 (HAC1), and the FAMA–RBR and FAMA–HAC1 complexes both maintain the irreversible fate of GCs through histone modifications^18,44,64,65. For example, when the FAMA–RBR interaction is abolished, paired GCs differentiate inside of existing GCs and lobed GCs form. We did not observe such abnormal GCs on the epidermis of wsb scap1, suggesting that WSB and SCAP1 are not primarily involved in this third FAMA role. As with Arabidopsis FAMA, FAMA in rice (Oryza sativa) promotes this GC differentiation^66,67. BdFAMA mainly promotes GC differentiation, although BdFAMA and BdMUTE share functions⁵⁹. These results suggest that the promotion of GC differentiation is the primary role of FAMA in angiosperms. A stomatal function for WSB and SCAP1 in other species remains to be determined.

Genome-wide profiling of WSB targets revealed that WSB contributes to the full expression of SCRM. Reporter analyses revealed that WSB is required for FAMA expression (Fig. 6). FAMA was clearly downregulated in wsb, while SCRM was weakly downregulated in wsb. This difference may be explained by the broad expression pattern of SCRM, which is expressed in early stomatal lineages¹⁷ and mesophyll cells²⁶, in addition to GMCs and GCs. The binding strength of FAMA to the WSB locus is comparable to that at the EPF1 locus, to which FAMA is reported to bind (Fig. 2c). The binding strength of SCRM at the WSB locus is comparable to that at the TMM locus, to which SCRM is reported to bind (Extended Data Fig. 2e). Notably, the binding of WSB to the FAMA locus is very weak because such binding can only be detected when WSB is overexpressed, and it was not detected by ChIP-seq analysis for WSB (Fig. 6k). These results suggest that the binding of WSB to the FAMA locus is transient and/or that WSB indirectly binds to the FAMA locus, that is, it binds through other protein(s). Further experiments are required to fully elucidate the detailed binding mode of WSB to the FAMA locus.

Analyses of fluorescent reporters^14,17 and scRNA-seq⁵³ indicated that the FAMA–SCRM complex exists until GCs mature, when WSB disappears (Fig. 4). The positive feedback circuit between WSB and FAMA–SCRM underpins the sustained expression of FAMA and SCRM in the stomatal lineage (Fig. 7d). We detected weak GFP signals in GCs of wsb proFAMA:GFP just after their division (Fig. 6j), suggesting that FAMA expression is not fully dependent on WSB. Combined with the results that showed that WSB is not expressed in the epidermis of fama, we conclude that FAMA is an upstream factor of WSB. In the epidermis, it was previously reported that MUTE, which is expressed in meristemoids and early GMCs, directly induces FAMA expression in stomatal-lineage cells⁵¹. Therefore, these results suggest that MUTE induces FAMA in late meristemoids and/or early GMCs, which in turn induces WSB in GMCs and GCs just after division; finally, the FAMA–WSB feedback sustains FAMA expression in young/maturing GCs. Around the time when WSB expression declines, FAMA may start to induce the expression of SCAP1, whose encoded protein subsequently induces MYB60 and GORK expression, key factors for stomatal maturation. It remains an open question how WSB disappears during stomatal differentiation even though FAMA is still expressed. One hypothesis is that unknown stomatal lineage-specific factors might degrade WSB. Alternatively, epigenetic regulators such as HAC1 may be required for the sustained expression of FAMA¹⁸. Future studies will shed light on how WSB and epigenetic regulators work together or independently to maintain sustained FAMA expression.

How did plants adjust the expression of WSB in MCs during evolution? In cells from the stomatal lineage, FAMA expression is MUTE dependent and FAMA independent⁵¹, and importantly, MUTE is not expressed in MCs²⁵. Our comparative analysis of AtFAMA and BdFAMA indicates that a positive feedback regulation of FAMA by itself is essential for FAMA and subsequent WSB expression in MCs (Fig. 7c)^25,26. How Brassicales plants established the positive feedback regulation of FAMA is an important question for further investigation. For example, the evolution of cis elements necessary for FAMA binding in the FAMA promoter and the 3D structure of the FAMA–SCRM complex bound to these cis elements should be examined. We also do not exclude the possibility that other inner tissue-specific phytohormones and/or TFs contribute to the expression of WSB with the FAMA–SCRM complex.

One remarkable difference between GCs and MCs is their cell size^29,68. Here we showed that a positive regulator of cell and nuclear expansion, CCS52A1, was directly induced by WSB and was more strongly expressed in MCs than in GCs. Indeed, loss of CCS52A1 function resulted in smaller cells and nuclei, as well as lower ploidy, for mature MCs (Fig. 6). These results suggest that the strong and/or sustained induction of CCS52A1 is a key cellular event for the differentiation of MCs (Fig. 7d). To clarify the detailed role of CCS52A1 in the differentiation of MCs, live imaging of MC development and their ploidy change using the histone H2B reporter (for example, H2B-GFP) in WT and ccs52a1 mutants is needed in future studies. In addition, comparative transcriptome studies of MCs from WT and ccs52a1 are needed to elucidate the comprehensive role of CCS52A1 in MCs using single-cell analysis approaches. Further analysis of other WSB targets will provide a detailed molecular mechanism of MC differentiation.

In addition to co-option, we also revealed the neofunctionalization of FAMA and WSB during evolution (Fig. 7d). First, the neofunctionalization of FAMA led to its essential role in its own activation in the MC lineage. We do not exclude the possibility that the positive feedback regulation of FAMA also occurs in the stomatal lineage, but it is not essential because some cells in fama-1 tumours express the FAMApro:GFP reporter¹⁴. Second, the neofunctionalization of WSB is reflected by its prolonged expression window in MC lineages. Prolonged expression of WSB strongly induces the expression of CCS52A1 for MC differentiation (increased ploidy and cellular and nuclear expansion) (Fig. 6) and repressed SCAP1 expression to prevent MCs from acquiring stomatal characters (Extended Data Fig. 7d). The stomatal-lineage expression window of SCAP1 in WSB-knockdown lines (amiRNA-WSB) is indistinguishable from that in the WT, suggesting that the inhibitory effects of WSB toward SCAP1 are MC lineage specific (Supplementary Fig. 11). For the prolonged expression of WSB, promoter activity and/or WSB stability may be critical. Indeed, WSB homologues outside Brassicales driven by the AtWSB promoter rescued the absence of MCs in the wsb mutant (Supplementary Fig. 21), suggesting that the AtWSB promoter is important for the development of MCs and neofunctionalization of WSB. Analyses of the evolution of cis regulatory elements in the AtWSB promoter are needed for future studies. We cannot address a possible neofunctionalization of WSB protein stability in Brassicales, as we used codon-optimized sequences for each WSB homologue tested, which may improve translational efficiency. These two neofunctionalization events of FAMA and WSB are connected by the transcriptional layers between FAMA and WSB. Taken together, these results suggest that the positive feedback-mediated sustained or high expression of WSB (WSB dynamics) in the MC lineage may have driven the evolution of MCs. Thus, we provide evidence here for the co-option and the neofunctionalization of stomatal executors, FAMA and WSB, for defence against herbivores in Brassicales.

Methods

Plant materials and growth conditions

The Arabidopsis (Arabidopsis thaliana) accession Columbia-0 (Col-0) was used for all lines except for MYR002:GUS (in C24). The T-DNA insertion mutants SALK_100073 (fama-1), SALK_003155 (scrm) and SAIL_808_B10 (scrm2-1) were obtained from the Arabidopsis Biological Resource Center (ABRC) at Ohio State University. The iFAMA MYR001:GUS²⁶, MYR001:GUS³³, proTGG2:VENUS-2s⁴⁶ and syp22 (ref. ³⁶) germplasms were reported previously. proFAMA:FAMA-myc⁴⁴ was provided by Dr Dominique Bergmann (Stanford University). scrm scrm2 (ref. ¹⁷), scrm-D¹⁷ and proSCRM:SCRM-GFP¹⁷ were provided by Dr Keiko Torii (University of Texas at Austin). scap1 (ref. ²²) and proSCAP1:GUS²² were provided by Dr Juntaro Negi (Kyushu University). MYR002:GUS (proWRKY23:GUS⁶⁹) was provided by Dr Tom Beeckman (Ghent University). proFAMA:GFP⁶⁴ was provided by Dr EunKyoung Lee and Dr Abel Rosado (University of British Columbia). ccs52a1-1 (refs. ^54,56) and proCCS52A1:GUS^54,56 were provided by Dr Masaaki Umeda (Nara Institute of Science and Technology). ccs52a1-5 (ref. ⁵⁵) was provided by Dr Eva Kondorosi (Institute of Plant Biology) and Dr Peter Mergaert (Institute for Integrative Biology of the Cell). fama proAtFAMA:MpSETA⁵⁸ was provided by Dr Tomoo Shimada (Kyoto University). wsb mutants were generated by CRISPR/Cas9-mediated gene editing (see below). The wsb scap1 double mutant, wsb-1/wsb-2 F₁ plants and wsb-1/wsb-3 F₁ plants were generated by crossing. Mutants and transgenic plants used in this study are listed in Supplementary Data 12. Seeds were surface sterilized with 70% (v/v) ethanol and then sown onto Murashige and Skoog (MS) medium (Wako) solidified with 0.5% (w/v) gellan gum (Wako) and containing 1% (w/v) sucrose. The seeds were stratified at 4 °C for 2–5 days to break dormancy and grown at 22 °C for 20 days under continuous light (100 µE s⁻¹ m⁻²) or long-day conditions (16 h light/8 h dark) in a plant growth chamber (BiOTRON, LPH-411SP, NIPPON MEDICAL and CHEMICAL INSTRUMENTS) with fluorescent light tubes. Humidity in the plant growth chamber was kept at 40–50%. After 14–21 days, plants were transferred to pots containing vermiculite and Metro-Mix for subsequent growth. The humidity in the plant culture room was kept at 40–50%.

Plasmid construction and transgenic Arabidopsis lines

The Gateway cloning system (Life Technologies) was used for plasmid construction. For transcriptional GUS constructs, a 1,444-bp WSB promoter fragment, a 1,703-bp SCAP1 promoter fragment (proSCAP1 67bp_del) and a 1,704-bp HIPP20 promoter fragment were individually cloned into pENTR D-TOPO (Invitrogen by Thermo, K240020). The resulting subcloned promoters were individually recombined into the binary vector pBGWFS7 or pHGWFS7 using LR reactions (LR Clonase II Enzyme mix, Invitrogen by Thermo, 11791100). For translational C-terminal fusion constructs of WSB, the coding sequence of mVenus or mTurquoise2 was inserted in-frame and downstream of a 2.2-kb WSB genomic fragment (including 1.4 kb of the promoter sequence) without the stop codon and with a GGSG linker sequence. For translational C-terminal fusion constructs of SCAP1, the coding sequence of mVenus or mCherry or mTurquoise2 was inserted in-frame and downstream of a 2.7-kb SCAP1 genomic fragment (including 1.8 kb of promoter sequence) without the stop codon and with a GGSG linker sequence. After cloning into pENTR D-TOPO, the constructs were recombined into the binary vector pGWB601 (ref. ⁷⁰) using an LR reaction. For the inducible overexpression construct of WSB-mVenus, the coding sequence of mVenus was inserted in-frame and downstream of the 0.9-kb WSB coding sequence without a stop codon and with a GGSG linker sequence. After cloning into pENTR D-TOPO, the constructs were recombined into the binary vector pMDC7 (ref. ⁷¹) using an LR reaction. For the proWSB:TpWSB, proWSB:CpWSB, proWSB:ChWSB, proWSB:AmWSB, proWSB:BdWSB and proWSB:OrWSB constructs, the coding sequences of TpWSB, CpWSB, ChWSB, AmWSB, BdWSB and OrWSB codon-optimized for Arabidopsis were synthesized (Eurofins Japan) and cloned into pENTR D-TOPO. The resulting constructs were then recombined into the binary vector R4pGWB601 with pENTR 5′-TOPO proWSB using an LR reaction. For the proFAMA:H2B-tdTomato, proFAMA:mScarlet-N7, proFAMA:mCherry-RCI2A, proFAMA:CYCD7-mVenus and proFAMA:BdFAMA-mVenus constructs, a 3,102-bp FAMA promoter fragment was cloned into pENTR 5′-TOPO, and the coding sequences of mScarlet-N7, mCherry-RCI2A, CYCD7-mVenus and BdFAMA-mVenus (BdFAMA was synthesized by Eurofins Japan) were cloned into pENTR D-TOPO, respectively. pENTR D-TOPO H2B-tdTomato was provided by Dr T. Goh (Nara Institute of Science and Technology). The resulting constructs were recombined into the binary vector R4pGWB601 (ref. ⁷⁰) using an LR reaction. For proUBQ11:mCitrine-RCI2A, the UBQ11 promoter was cloned into pENTR 5′-TOPO, and the coding sequence of mCitrine-RCI2A was cloned into pENTR D-TOPO. The resulting constructs were recombined into the binary vector R4pGWB601 (ref. ⁷⁰) using an LR reaction. Two artificial microRNAs (amiRNAs) against WSB were designed using the WMD3-Web MicroRNA Designer (http://wmd3.weigelworld.org/cgi-bin/webapp.cgi)⁷² and amplified using the following primers from the pRS300 vector⁷²: amiRNA-F, amiRNA-R, WSB_amiRNA3_1, WSB_amiRNA3_2, WSB_amiRNA3_3, WSB_amiRNA3_4, WSB_amiRNA6_1, WSB_amiRNA6_2, WSB_amiRNA6_3 and WSB_amiRNA6_4. The amplified amiRNA-WSB DNA fragment was cloned into the pENTR D-TOPO plasmid before being recombined into binary vector pFAST-R02 (ref. ⁷³) using an LR reaction to generate amiRNA-WSB. For transcriptional GUS constructs, a MYB60 promoter fragment (1,000 bp), a KAT1 promoter fragment (2,000 bp), a TET4 promoter fragment (1,000 bp) and an At1g23170 promoter fragment (1,398 bp) were cloned into pENTR D-TOPO (Invitrogen by Thermo, K240020). The resulting subcloned promoters were individually recombined into the binary vector pBGWFS7 or pHGWFS7 using LR reactions. All binary constructs were transformed individually into Agrobacterium tumefaciens strain GV3101. Plants were transformed with Agrobacterium cultures harbouring each binary vector using the floral-dip method⁷⁴. T₁ seedlings were selected on medium containing 10 mg l⁻¹ BASTA or 25–50 mg l⁻¹ hygromycin B, or the FAST-R reporter⁷³. Mutants and transgenic plants used in this study are listed in Supplementary Data 12 and the primer sets are listed in Supplementary Data 13.

Oestrogen treatment

Estradiol (Sigma-Aldrich, E8875-1G) was dissolved in 70% (v/v) ethanol to 10 mM concentration and used at 10 µM (1,000× dilution) or 0.01 µM (1,000,000× dilution, only for ChIP-seq experiments using iWSB-mVenus lines) as the final working concentration in water. For mock treatment, a 0.07% (v/v) ethanol solution was used. Seeds were sown onto MS medium (Wako) solidified with 0.5% (w/v) Gellan gum (Wako) and containing 1% (w/v) sucrose; seedlings at 4 DAG were transferred to water containing 10 µM oestrogen or water with 0.07% (v/v) ethanol, vacuumed at least three times and incubated for the indicated times.

Generation of wsb mutants by CRISPR/Cas9

To generate wsb mutants by CRISPR/Cas9, two single guide RNAs (sgRNAs) were manually designed to target the sequence encoding the DNA-binding domain of WSB. The sgRNAs fused to the Arabidopsis tRNA-Gly gene were amplified using pGTR (Addgene, plasmid 63143)⁷⁵ as PCR template and utilizing primers listed in Supplementary Data 13. The resulting PCR products harbouring the two sgRNA sequences were cloned into the pKI1.1R vector (Addgene, plasmid 85808)⁴⁸ digested with AarI using NEBuilder HiFi enzyme (NEB). pKI1.1R harbouring the two sgRNAs was transformed into Arabidopsis Col-0 by the floral-dip method using Agrobacterium strain GV3101 (ref. ⁷⁴). T₁ transformants and T₂ transgenic plants without the CRISPR/Cas9 construct were selected using the FAST-R reporter⁷³.

GUS staining and microscopy

Samples were first placed into ice-cold 90% (v/v) acetone for 15–30 min and then into β-glucuronidase (GUS) staining solution containing 0.5 mg ml⁻¹ X-Gluc (Gold Biotechnology, G1281C5), 0.1 M sodium phosphate buffer (pH 7.0), 10 mM EDTA, 0.5–5 mM potassium ferricyanide, 0.5–5 mM potassium ferrocyanide and 0.1% (v/v) Triton X-100. Samples were vacuum infiltrated with GUS staining solution before being incubated at room temperature for 2–24 h²⁶. Tissue sectioning was performed as described previously³⁹. The sections were stained with 0.01% (w/v) toluidine blue (Wako Chemicals). Representative images were captured under an AX-70 light microscope (Olympus) and an AXIO Zoom V16 (Zeiss) microscope.

Confocal laser-scanning microscopy

Fluorescence micrographs were obtained with a confocal laser-scanning microscope (Leica SP8 FALCON or STELLARIS 5) using dry objectives, a pulsed laser (440 nm), a diode laser (448 nm), a white-light laser (WLL), a HyD/HyD S detector and a gating system. The WLL allows for time-gated acquisition (gating system) to remove the autofluorescence from chlorophylls. The laser wavelengths used were mainly 440 nm (mTurquoise2), 488 nm (GFP), 514 nm (mVenus) and 552 nm (mCherry). For multicolour imaging, the sequential mode was used. The images were analysed and processed using LAS X (https://www.leica-microsystems.com/products/microscope-software/p/leica-las-x-ls/) (Leica), Fiji (https://fiji.sc/) and Photoshop 2021 (Adobe) software. To observe inner tissues, ClearSee solution (Fujifilm)⁷⁶ was used. The samples were fixed overnight in 4% (w/v) paraformaldehyde in phosphate-buffered saline. Samples were cleared in ClearSee for at least 1 week. Before observations, SR2200 was added at a 5,000× dilution (SCRI Renaissance Stain 2200, Tokyo Future Style) for the staining of cell walls, or 80 μg ml⁻¹ Hoechst 33342 solution (DOJINDO) for DNA staining. For quantification of stomatal features, cell peripheries were visualized following staining with 1 mg ml⁻¹ propidium iodide (Molecular Probes, P1304MP). Methods for the staining of cell outlines are provided in Supplementary Data 14. The same confocal laser-scanning microscope settings were used when comparing the WT and mutants.

RNA-seq analysis

Libraries for RNA‐seq were prepared using the Breath Adapter Directional sequencing method⁷⁷. The libraries were sequenced on a Next‐Seq 500 instrument (Illumina). Mapping to the Arabidopsis reference genome (TAIR10) was conducted using Bowtie with the following options ‘‐‐all ‐‐best ‐‐strata ‐‐trim5 8’. The number of reads mapped to each gene was counted⁷⁷. After normalization, the FDR and fold-change were calculated using the edgeR package in R. Differentially expressed genes were identified (FDR < 0.05). The data have been deposited at the DNA Data Bank of Japan (DRA013687).

RT–qPCR

Leaves or stems were frozen in liquid nitrogen immediately after sample collection. An RNeasy Plant Mini kit (Qiagen, 74104) was used to extract total RNA. An RNase-Free DNase set (Qiagen, 79254) was used to eliminate genomic DNA contamination in RNA samples. Reverse transcription was performed using PrimeScript RT Master Mix (Takara, RR036A). Quantitative PCR was performed as described previously using FastStart Essential DNA Green Master Mix (Roche, 6924204001)⁷⁸. Arabidopsis ACT2 was used as the internal reference transcript. Each experiment was performed at least three times. The relative expression level of each gene was calculated using the 2^−ΔΔCT method⁷⁹. Primers are listed in Supplementary Data 13.

ChIP–qPCR experiment

For ChIP–qPCR, ChIP was carried out as described previously⁸⁰. Samples were collected from seedlings at 4 DAG. For each sample, 300–600 mg of fresh seedling tissue was crosslinked in 1% (w/v) formaldehyde for 15 min. After quenching the crosslinking with 0.93% (w/v) glycine for 5 min, tissues patted dry with Kimwipes were frozen in liquid nitrogen and kept at −80 °C until use. Tissues were ground to a fine powder with an ice-cold mortar and pestle. Using nuclear extraction buffer, chromatin was isolated from the nuclear extracts. Fragmentation was conducted using a UD-201 Ultrasonic Disruptor sonicator (TOMY) or Bioruptor II (BMBio). After preclearing, the indicated antibodies were added and the mixtures were rotated overnight at 4 °C. The following antibodies were used: anti-Myc (sc-40 X, Santa Cruz Biotechnology, 3 µl per sample) and anti-GFP (SAB4301138, SIGMA, 5 µl per sample). For immunoprecipitation to capture DNA–protein complexes, Dynabeads with Protein A or G (Thermo Fisher) were used. Beads were washed in low-salt buffer and 250 mM LiCl buffer, and DNA was eluted from the beads overnight at 65 °C. The resulting DNA was purified using a QIAquick PCR Purification kit (Qiagen). DNA was quantified on a LightCycler 480 (Roche) instrument using FastStart Essential DNA Green Master Mix (Roche). The ratio of ChIP to input DNA (% input) was compared on the basis of the reaction threshold cycle for each ChIP sample compared with a dilution series of the corresponding input sample. Relative values were normalized to the negative control locus of the TA3 retrotransposon (At1g37110)⁸⁰. At least three independent experiments were performed. Primers are listed in Supplementary Data 13.

ChIP-seq experiments

ChIP-seq was performed as previously described with slight modifications⁸¹. Two grams of iWSB-mVenus seedlings at 3 DAG were bathed in 0.01 µM estradiol in water for 24 h. Samples were then rapidly frozen in liquid nitrogen and stored at −80 °C until use. The frozen tissues were thoroughly ground using a mortar and pestle to obtain a fine powder. Chromatin was then fixed and isolated using a nuclei isolation buffer (10 mM HEPES, 1 M sucrose, 5 mM KCl, 5 mM MgCl₂ and 5 mM EDTA) containing 1% (w/v) formaldehyde (Thermo Scientific). After chromatin fragmentation using an ultrasonicator (Covaris M220), immunoprecipitation of the resulting chromatin fragments was carried out using anti-GFP antibodies (SAB4301138, Sigma, diluted 1:1,000) and Dynabeads protein A (Thermo Scientific) at 4 °C. The validity of the antibodies was confirmed by the suppliers and by immunoblot analysis in the laboratory before use. Following immunoprecipitation, DNA was purified using a DNA Cleanup kit (New England Biolabs). The resulting DNA was used as a template to generate a sequencing library using a ThruPLEX DNA-seq kit (Rubicon Genomics), following manufacturer instructions. The immunoprecipitated fraction was analysed using a NovaSeq 6000 instrument (Illumina). The resulting FASTQ file underwent quality assessment using FastQC (v.0.11.7) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/); the raw reads were subjected to trimming using Trimmomatic (v.0.38)⁸². The reads were mapped to the Arabidopsis TAIR10 genome using Bowtie2 (v.2.3.4.2)⁸³. Peak calling was performed with MACS2 (v.2.2.6)⁸⁴ using the callpeak command with the parameter ‘–nomodel -q 0.01’. The read counts were calculated using featureCounts (v.1.6.3). Motif analysis was conducted using HOMER (v.4.1) (http://homer.ucsd.edu/homer/). Heat maps were generated using deeptools (v.3.2.1)⁸⁵. Binding peaks were visualized in the Integrative Genomics Viewer (v.2.8.13)⁸⁶. The data have been deposited at the DNA Data Bank of Japan (DRA016932).

SDS–PAGE and immunoblot analysis

SDS–PAGE and immunoblot analysis were performed as described previously³⁶. The antibodies used in this analysis were anti-GFP (1,000× dilution; SAB4301138, Sigma), anti-TGG1 (5,000× dilution)³⁶, anti-TGG2 (5,000× dilution)³⁶, anti-ACTIN (2,000× dilution; A0480, Sigma-Aldrich) and anti-histone H3 (1,000× dilution; AB1791, Abcam). Membranes were exposed to Immobilon Western HRP Substrate (Millipore, WBKL0500), and luminescence signals were detected with a CCD imager (ImageQuant LAS 4000, GE Life Sciences). Stripping of membranes was performed with Restore PLUS Western Blot Stripping Buffer (Thermo Fisher, 46430). SimplyBlue SafeStain (Invitrogen, 465034) was used for Coomassie brilliant blue staining of membranes to ensure equal loading. As molecular weight standards, SeeBlue Plus2 Pre-stained standard (Thermo Fisher) and MagicMark XP Western Protein (Thermo Fisher) were used for SDS–PAGE and immunoblot analysis, respectively. Full images of immunoblots and Coomassie brilliant blue-stained membranes are shown in Supplementary Figs. 1–7.

Quantification of stomatal features

For confocal microscopy, cell peripheries of the third true leaves at 14 and 21 DAG were visualized following staining with 1 mg ml⁻¹ propidium iodide (Molecular Probes, P1304MP). Images were captured using a confocal laser-scanning microscope (Leica SP8 FALCON) with a WLL, a HyD detector and a gating system. Approximately 40 sequential confocal slices were used for the construction of Z-stack images covering a 0.308-mm² field. Brightness and contrast were uniformly adjusted using Photoshop 2021 (Adobe). The number of each category of stomata was scored in the 0.308-mm² field (WT n = 9, wsb n = 8, scap1 n = 8, wsb scap1 n = 15, gWSB-mTurquoise2 n = 24, gWSB-mVenus n = 9, gSCAP1-mVenus n = 21, gSCAP1-mCherry n = 13).

Transmission electron microscopy

The third true leaves of WT and wsb scap1 seedlings were cut and immediately fixed with 2% (w/v) paraformaldehyde and 2% (w/v) glutaraldehyde in 0.05 M cacodylate buffer (pH 7.4) at 4 °C overnight. After fixation, the samples were washed 3 times with 0.05 M cacodylate buffer for 30 min each and postfixed with 2% (w/v) osmium tetroxide in 0.05 M cacodylate buffer (pH 7.4) at 4°C for 3 h. The samples were infiltrated with propylene oxide (PO) twice for 30 min each, placed into a 1:1 (v/v) mixture of PO and resin (Quetol-651, Nisshin EM) for 3 h, and transferred to 100% resin for 3 h. The samples were allowed to polymerize at 60 °C for 48 h. The polymerized resin blocks were ultrathin sectioned to 80 nm with a diamond knife using an ultramicrotome (Ultracut UCT, Leica) and the sections were mounted onto copper grids. The sections were stained with 2% (w/v) uranyl acetate at room temperature for 15 min and washed with distilled water, followed by secondary staining with 1× lead stain solution (Sigma-Aldrich) at room temperature for 3 min. The grids were observed with a transmission electron microscope (JEM-1400Plus, JEOL) at an acceleration voltage of 100 kV. Digital images (3,296 × 2,472 pixels) were taken with a CCD camera (EM-14830RUBY2, JEOL).

Phylogenetic analysis

Protein sequences related to Group IV members of the Arabidopsis AP2/ERF TF family⁴³ were retrieved from the proteome databases in Phytozome 13 (https://phytozome-next.jgi.doe.gov/) and MarpolBase (https://marchantia.info/)⁸⁷ using BLASTP searches. The amino acid sequences of WSB (At3g57600), At1g75490, At2g38340, At2g40220, At2g40340, At2g40350, At3g11020, At5g05410 and At5g18450 were used as queries for BLASTP searches with a threshold e-value of 12. A total of 461 protein sequences were collected from the following 62 plant species (gene name prefixes are given in parentheses): Amaranthus hypochondriacus (AH), Amborella trichopoda (evm_27.TU.AmTr_v1.0_scaffold), Anacardium occidentale (Anaoc.), Ananas comosus (Aco), Anthoceros agrestis (AagrOXF_evm.model.utg), Anthoceros angustus (AANG), Anthoceros punctatus (Apun_evm.model.utg), Aquilegia coerulea (Aqcoe), Arabidopsis halleri (Araha.), Arabidopsis lyrata (AL), Arabidopsis thaliana (AT), Asparagus officinalis (evm.TU.AsparagusV1_), Beta vulgaris (EL10Ac), Brachypodium distachyon (Bradi), Brassica oleracea var. capitata (Bol), Brassica rapa (Brara.), Capsella rubella (Carub.), Carica papaya (evm.TU.supercontig), Ceratodon purpureus (CepurGG1), Ceratopteris richardii (Ceric), Cinnamomum kanehirae (CKAN), Citrus sinensis (orange1.), Coffea arabica (evm.TU.Scaffold), Corymbia citriodora (Cocit.), Dioscorea alata (Dioal.), Eucalyptus grandis (Eucgr.), Glycine max (GlymaFiskIII.), Gossypium hirsutum (Gohir.), Helianthus annuus (HanXRQChr), Hordeum vulgare (HORVU.MOREX.r3.), Kalanchoe laxiflora (Kalaxd.), Lactuca sativa (Lsat_1_v5_gn), Linum usitatissimum (Lus), Lotus japonicus (Lj), Manihot esculenta (Manes.), Marchantia polymorpha (Mp), Medicago truncatula (Medtr), Nymphaea colorata (Nycol.), Oryza sativa (Os), Phaseolus vulgaris (Pv5-593.), Physcomitrium patens (Phpat.), Picea abies (MA), Poncirus trifoliata (Ptrif.), Populus trichocarpa (Potri.), Ricinus communis (gene name starts with number), Salix purpurea (Sapur), Salvinia cucullata (Sacu_v1.1_s), Setaria italica (Seita.), Solanum lycopersicum (Solyc), Solanum tuberosum (Soltu.), Sorghum bicolor (SbiSC187.), Sphagnum fallax (Sphfalx), Sphagnum magellanicum (Sphmag), Spinacia oleracea (Spov3_), Spirodela polyrhiza (Spipo), Theobroma cacao (Thecc.), Thuja plicata (Thupl.), Triticum aestivum (Traes), Vitis vinifera (VIT_), Vigna unguiculata (Vigun), Zea mays (Zm) and Zostera marina (Zosma). A multiple sequence alignment was created using MAFFT (v.7.505) with the ‘–auto’ option (https://mafft.cbrc.jp/alignment/software/)⁸⁸, trimmed with trimAL (v.1.4.rev15) with the ‘-gt 1 -cons 10’ options (http://trimal.cgenomics.org/trimal)⁸⁹ and used for phylogenetic analysis in RAxML (v.8.2.12) with the ‘-f a -m PROTGAMMAAUTO’ options with 1,000 bootstrap replicates (https://cme.h-its.org/exelixis/web/software/raxml/)⁹⁰. The tree was visualized in R (v.4.2.1) with the ggtree package (v.3.4.2)⁹¹.

Measurement of stomatal responses to light and CO₂

Stomatal conductance in intact leaves was measured using a gas-exchange system (LI-6400, Li-Cor) as described previously⁹². Leaves from dark-adapted plants were illuminated with red light (300 µmol m⁻² s⁻¹) for 1 h, and then blue light (10 µmol m⁻² s⁻¹) was superimposed onto red light for 20 min. For the measurement of CO₂-induced stomatal responses, the ambient CO₂ concentration was lowered from 350 ppm to 100 ppm and then increased to 800 ppm in the dark. For stomatal aperture measurements, epidermal strips prepared from dark-adapted leaves were incubated in 5 mM MES-bistrispropane (pH 6.5), 50 mM KCl and 0.1 mM CaCl₂ in the dark or under red light (50 µmol m⁻² s⁻¹) and blue light (10 µmol m⁻² s⁻¹) for 2 h at 24 °C. Stomatal apertures on the abaxial epidermis were measured using an inverted microscope (Eclipse TS100, Nikon).

Statistics and reproducibility

Statistical analyses were performed using Microsoft Excel, GraphPad Prism 9 and GraphPad Prism 10. Student’s t-tests were used to compare two samples. For comparison of three or more samples, one-way analysis of variance (ANOVA) was used to test for significant interactions. When significant, Tukey–Kramer tests were used for pairwise comparisons. For experiments conducted without statistical analysis (for example, CLSM observation and GUS staining), a minimum of three independent biological replicates were performed to confirm reproducibility. The source data behind the graphs in the paper and all statistical analyses including exact P values are provided in Supplementary Data 15.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.