Loss or duplication of key regulatory genes coincides with environmental adaptation of the stomatal complex in Nymphaea colorata and Kalanchoe laxiflora

The stomatal complex is critical for gas and water exchange between plants and the atmosphere. Originating over 400 million years ago, the structure of the stomata has evolved to facilitate the adaptation of plants to various environments. Although the molecular mechanism of stomatal development in Arabidopsis has been widely studied, the evolution of stomatal structure and its molecular regulators in different species remains to be answered. In this study, we examined stomatal development and the orthologues of Arabidopsis stomatal genes in a basal angiosperm plant, Nymphaea colorata, and a member of the eudicot CAM family, Kalanchoe laxiflora, which represent the adaptation to aquatic and drought environments, respectively. Our results showed that despite the conservation of core stomatal regulators, a number of critical genes were lost in the N. colorata genome, including EPF2, MPK6, and AP2C3 and the polarity regulators BASL and POLAR. Interestingly, this is coincident with the loss of asymmetric divisions during the stomatal development of N. colorata. In addition, we found that the guard cell in K. laxiflora is surrounded by three or four small subsidiary cells in adaxial leaf surfaces. This type of stomatal complex is formed via repeated asymmetric cell divisions and cell state transitions. This may result from the doubled or quadrupled key genes controlling stomatal development in K. laxiflora. Our results show that loss or duplication of key regulatory genes is associated with environmental adaptation of the stomatal complex.


Introduction
Stomata are a pore-like structure in multiple organs, including leaves and stems, which facilitates gas and water exchange. When environmental conditions are unfavourable, plants can regulate water evapotranspiration and reduce CO 2 uptake by opening and closing the stomata. For instance, Crassulacean acid metabolism (CAM) plants are adapted to arid conditions 1 . The stomata in CAM plants remain closed during the day to reduce evapotranspiration while staying open at night to absorb CO 2 . These physiological traits make CAM plants resistant to diverse stresses, including strong irradiance and drought 2 .
Stomatal structure is highly conserved across land plants. The basic core structure with two guard cells surrounding the stomatal pore has remained unchanged during evolution 3 . However, the patterning of the mature stomatal structure differs among plant groups and can be generally summarized by three classes: anomocytic, stephanocytic, and paracytic 4 . The widely used model plant Arabidopsis thaliana exhibits anomocytic stomata. However, there are a few species (for example, CAM families) among the eudicots with paracytic stomata 5 .
Most grass species have paracytic mature stomata 6 . Amborella trichopoda in ANITA possesses stephanocytic stomata 7 . The diverse architecture of mature stomatal structures may suggest the evolution of their different developmental regulations and their adaption to different environments.
In A. thaliana, meristem mother cells (MMCs) undergo up to three asymmetrical divisions to form guard mother cells (GMCs). In grasses, meristemoids divide asymmetrically to form GMCs, and the lateral neighbouring axial cell lineage surrounding the GMC undergoes asymmetric division to give rise to lateral subsidiary cells (LSCs) 8 . In A. trichopoda, however, protodermal cells can directly become GMCs or divide asymmetrically to produce a GMC 9 . Hence, the regulation of stomatal development is highly diverse in different groups of land plants.
In the past, A. thaliana and Oryza sativa were often used as model systems to study stomatal patterning and development. Based on those studies, we now have a good understanding of the basic molecular network behind stomatal development. In A. thaliana, a complex signalling cascade of several genes has been identified to promote stomatal development. The secreted peptides of the EPIDERMAL PATTERNING FACTOR (EPF)/EPF-LIKE (EPFL) family act with a mitogen-activated protein kinase (MAPK) cascade to regulate the activity of basic-helixloop-helix (bHLH) transcription factors 10 . EPF1 and EPF2 specifically bind to leucine-rich repeat receptor (LRR) kinase complexes that include members of TOO MANY MOUTHS receptor-like protein (TMM) and the ERECTA family (ER). EPF1 is expressed in late-stage meristemoids, GMCs and young guard cells, whereas EPF2 is expressed in early-stage protodermal cells 11,12 . In the downstream pathway, a number of mitogen-activated protein (MAP) kinases, including MAPKKK YODA, MPKK4/5, MPKK7/9, and MAPK MPK3/6, were found to transduce the signalling for stomatal development 13 . Five bHLH transcription factors positively regulate the stomatal-lineage transition and differentiation. For example, SPEECHLESS (SPCH), MUTE, and FAMA act sequentially to promote the cellular transition in a stagespecific manner. SPCH regulates asymmetric divisions in MMC and MUTE involved in GMC differentiation 14,15 . FAMA promotes the last step to form GCs 16 . Two additional bHLH proteins, SCREAM/ICE1 and SCREAM2, act redundantly to heterodimerize SPCH, MUTE, and FAMA to coordinate the regulation 17 .
Polarity information is critical in stomatal development and directs asymmetric cell division and possibly cell fate determination. In A. thaliana, two unique polarity proteins, POLAR LOCALIZATION DURING ASYM-METRIC DIVISION AND REDISTRIBUTION (POLAR) and BREAKING OF ASYMMETRY IN THE STOMA-TAL LINEAGE (BASL), show mostly overlapping localization during asymmetric stomatal divisions 18,19 . In the grass, the asymmetric division taking place in the lateral neighbouring cell to produce the subsidiary cell relies on two LRR receptor-like kinases, PANGLOSS1 (PAN1) and PAN2 20,21 . PAN proteins are located at the poles in SMCs at the site of contact with GMCs, which precedes the polar accumulation of small GTPases (ROPs) and F-actin 22 . Interestingly, recent observations in Brachypodium distachyon found that BdMUTE regulates subsidiary cells through cell-to-cell movement 23 . In contrast, the MUTE homologue in A. thaliana is immobile 23,24 .
Although stomata morphologies across land plants have been widely examined, questions on the early evolution of angiosperms and the adaptation of stomata to diverse environments remain to be answered. It is not clear how molecular regulation of stomatal development evolved and how that relates to the diverse stomata morphologies among the land plants. Immediately above the root node of angiosperm evolution is the ANITA grade (basal angiosperms), which includes Amborella, Nymphaeales, Illiciaceae, Trimeniaceae and Austrobaileyaceae 7 . In this study, we took advantage of the newly sequenced genome of Nymphaea colorata (not released yet), a typical base angiosperm, to examine stomata regulation in early angiosperm evolution.
The structure and function of stomata are important for environmental adaptation. In some species, stomata underwent radical modifications to facilitate habituation to a particular environment. A recent study indicated that Z. marina lost all the genes involved in stomatal differentiation, which is coincident with its marine habituation. Nymphaea colorata is also an aquatic plant, so it is interesting to know if its stomata-related genes also changed during evolution. By contrast, Kalanchoe laxiflora, a CAM species, has adapted to drought conditions and has evolved specialized stomata functions. To understand how the evolution of the molecular regulation of stomatal development is associated with environmental adaptation, we analysed stomatal morphologies and related regulatory cascades in both Nymphaea colorata and K. laxiflora. Our analysis showed that although generally conserved, loss or duplication of key genes could be associated with structural and physiological renovations required for individual adaptation of plants to local environments.

Plant materials and growth condition
A. thaliana Columbia seeds were germinated and grown on 1/2 MS medium with 1% agar, 1% sucrose and 0.05% (wt/vol) morpholinoethansulfonic acid monohydrate (pH 5.7) under a 16/8-h light/dark cycle at 23°C. Plants were imaged 3-4 days after planting. O. sativa and K. laxiflora were grown at 28°C with a 16/8-h light/dark photoperiod. N. colorata were cultivated in water at 23°C in the greenhouse. Leaves of Spirodela polyrhiza were collected in winter 2017 at the Fujian Agriculture and Forestry University.

Microscopy and image processing
For Differential Interference Contrast (DIC) imaging, the protocol was modified slightly according to Raissig et al. 23,25 . Samples from the mid-regions of leaves were cut into small squares and cleared using a solution (ethanol: acetic acid glacial, in proportions 4:1 by volume) to remove chlorophyll; then, samples were subjected to a basic solution (a mixture of 7% NaOH in 60% ethanol). Finally, samples were washed briefly with 40% ethanol and mounted in water for visualization and microscopy analysis. Samples were examined using a Nikon ECLIPSE Ni-U microscope fitted with a Nikon DS-Ri 2 digital camera. Images were processed using ImageJ.

Phylogenetic analysis
We surveyed a number of genomes, such as A. thaliana, K. laxiflora, Sorghum bicolor, O. sativa, Zea mays, Ananas comosus, S. polyrhiza, and A. trichopoda, from Phytozome v12. Nelumbo nucifera and Phalaenopsis equestris were retrieved from ftp://ftp.ncbi.nih.gov/genomes/. Ginkgo biloba was found from GigaDB (http://gigadb.org/). N. colorata was recently sequenced by Liangsheng Zhang's Lab in Fujian Agriculture and Forestry University, and sequences were available in the water lily genome database (eplant.org). To obtain probable orthologous genes, we performed BLASTp (protein query-proteins database) and tBLASTn (protein query-nucleic acid database) c-e Micrograph of stomata at different developmental stages in adaxial leaf surfaces. c Squared patterning, a protodermal cell. d Large round cells are putative GMCs (orange arrow). e Stage with maturing stomata (red arrow). Schematic diagram of stomatal development. A protodermal cell (pale blue) that differentiated directly into a guard mother cell (orange); then, the GMC divided into GCs (red) searches to selectively look for similar protein sequences from these genomes 26 . A MAFFT (Multiple Sequence Alignment program) was chosen to produce an alignment of all amino-acid sequences with a BLAST score of at least 60 against A. thaliana 27 . The phylogenetic tree was reconstructed using the maximum likelihood (ML) method in FastTree2 28 .
Protein domains were identified using the National Center for Biotechnology Information conserved domain search tool. PEST domains were identified using emboss. bioinformatics.nl/cgi-bin/emboss/epestfind.

Loss of stomatal development genes in N. colorata
It was reported that different stomatal development patterns occur in plants of the ANITA grade. A. trichopoda possesses mostly perigenous and mesoperigenous stomata 9 . In this species, protodermal cells can directly become GMCs or divide asymmetrically to produce GMCs and stomatal lineage ground cells 9 . However, in Nymphaea, protodermal cells seemed to skip asymmetric divisions and directly gave rise to GMCs 7,9 . It is still to be determined whether asymmetric division is an ancestral stomata-forming step during evolution.
To gain a deeper understanding of the ancestral development of stomatal structure, we performed anatomic observation of the stomatal structure in N. colorata. We found that N. colorata stomata are only present on the adaxial surface of the floating leaf, with each stoma surrounded by 4-8 neighbouring cells (Fig. 1a). On the abaxial surface of N. colorata, we only found hydropote complexes with lens-shaped cells and bowl-shaped cells, which appeared to be surrounded by specialized rosettes of epidermal cells (Fig. 1b). It was hypothesized that the hydropote in Nymphaea colorata is homologous to stomatal complexes, and its functions and morphologies are highly associated with aquatic habitats 29 . Similarly, another floating plant, S. polyrhiza, has lost stomata on  NcSPCH shares the bHLH domain (orange) and C-terminal SMF domain (light blue) with AtSPCH but has no protein degradation-associated PEST domain (grey) and has a shorter MAPK target domain (yellow). Both NcMUTE and AtMUTE genes have a unique conserved region (MUTE unique, dark blue) and lack some residues preceding the bHLH domain that are present in all the other bHLH Ia members with various lengths. Both NcFAMA and AtFAMA genes have high AA sequence similarity and harbour three unique domains (FAMA unique 1, red; FAMA unique 2, blue; Ia extension, brown). Both NcICE-like and AtICE1/2 have highly conserved bHLH domains, potential PEST domains and ACT domains (green)  Figure S1). These results reveal that floating plants tend to lose stomata or create special stomata-like structures to adapt to the aquatic environment. It can also be exemplified by seagrass, Zostera marina, in which no stomata are present on leaves, and coincidently, entire stomatal genes are lost to adapt to the marine lifestyle 30 . Although anatomical descriptions of stomatal development have been reported for many taxa, little is known about the evolution of the molecular machine of stomatal formation across land plants.
One way to understand the evolution of these essential regulators of stomatal development is to analyse their phylogenies. This is currently feasible based on the genome sequences for many species, including the eudicots A. thaliana and K. laxiflora; the monocot plants O. sativa and Z. mays. To facilitate our understanding of the early evolution of these regulators, we included basal angiosperms A. trichopoda, and we recently sequenced the genome of an early-divergent angiosperm N. colorata (see Materials and methods for information on genome data) (Fig. 2a). To understand some special features of stomata formation in N. colorata, we analysed the potential orthologues of A. thaliana genes involved in stomatal formation using the unique unpublished genome data of water lily. In line with A. thaliana, we found high conservation of the core genes required for stomatal formation in N. colorata, including an orthologue of an SPCH-like gene, NcSPCH (Fig. 2b); orthologue of a MUTE-like gene, NcMUTE (Fig. 2c); orthologue of a FAMA-like gene, NcFAMA (Fig. 2d), and two orthologues of an ICE/SCRM-like gene, NcICE1 and NcSCRM2 (Fig. 2e). We further analysed the conservation of the homologous domain of these proteins and found a high degree of domain conservation (Fig. 3). However, we also found a number of genes missing from the N. colorata genome, including the peptide ligands EPF2, MPK6, and AP2C3 and the polarity controllers BASL and POLAR (Fig. 4). Interestingly, the function of lost genes seems to be highly specific to the asymmetric stomatal development stages.

Stomatal development gene duplications in K. laxiflora
Whole-genome duplications (WGDs) are a common phenomenon during evolution, and the resulting gene duplications (GDs) provide redundant functions or specified novel functions [31][32][33][34] . WGDs are the source of functional diversity or novelty in the genome for adaption to environmental changes 35 . It has been suggested that two distinct WGDs occur in the K. laxiflora lineage and generate four gene copies across the genome 36 .

Novel formation of subsidiary cells in K. laxiflora
CAM increases water-use efficiency and drought resistance in plants, which is characterized by nocturnal opening and diurnal closing of the stomata 36 . Therefore, stomatal control in the leaves is particularly important for this type of plant to reduce evapotranspiration in the daytime and increase carbon dioxide (CO 2 ) collection at night 2,36 . The physiological traits probably improve the resistance of CAM plants to diverse environmental stresses, including drought 1,2 .
To gain a better understanding of the stomatal complex in CAM plants, we performed anatomical observation of K. laxiflora, a member of the eudicot CAM family. In K. laxiflora, stomata are surrounded by three to four small subsidiary cells in adaxial leaf surfaces (Fig. 6a). Similarly, we found that the stomata of Phalaenopsis equestris, another CAM monocot species, is also surrounded by approximately four subsidiary cells ( Figure S3). This innovation of stomatal architecture could derive from differential regulation of stomatal formation. We found that in K. laxiflora, stomata formed via a series of asymmetric cell divisions and cell state transitions: protodermal cells entered the stomatal lineage and took on a MMC identity; the MMC underwent three or four asymmetrical divisions to form GMC and Stomatal lineage ground cell (SLGC) (Fig. 6d-g). The GMC underwent a symmetric division to form a pair of guard cells, and SLGCs eventually became subsidiary cells surrounding the guard cell (Fig. 6b, c). NF not found It is widely accepted that different stomatal patternings reflect the asymmetric division of precursor cells and lateral divisions of neighbouring cells 37 . For example, in anomocytic stomata occurring in the eudicot A. thaliana (Fig. 7a, b), the MMC underwent three asymmetric divisions to give rise to a GMC and SLGCs, which was followed by a transition from SLGCs to pavement cells (Fig. 7c). Although both A. thaliana and K. laxiflora are eudicots, K. laxiflora possesses stephanocytic stomata (Fig. 7d, e). Developmentally, there is a similarity between these two types of stomata: meristemoids undergo a series of asymmetric divisions to produce SLGCs surrounding guard cells (Fig. 7f), and different cell fate choices of SLGCs finally give rise to different stomatal complexes ( Figure S4). In monocot species such as O. sativa, the type of mature stomata is named the paracytic type, in which the guard cell is surrounded by two subsidiary cells (Fig. 7g, h). In this type, the stomatal meristemoid divides asymmetrically to form a larger SLGC and a smaller meristemoid that directly forms the GMC. Before the GMC divides, it induces neighbouring cell files to adopt an SMC identity, which subsequently forms SCs via asymmetric divisions. The GMC then undergoes symmetric mitosis to eventually form guard cells (Fig. 7i). Therefore, subsidiary cells can develop through different ways: one is through asymmetric division in O. sativa, and the other is through SLGC differentiation in K. laxiflora. In K. laxiflora, subsidiary cells are noticeably visible, but little is known about the factors defining subsidiary cell identity. In Brachypodium distachyon, subsidiary cells are formed through asymmetric divisions. BdMUTE is an orthologue of A. thaliana MUTE that has been identified as sufficient for SC formation based on its acquisition of cell-to-cell mobility 23 . In A. thaliana, AtMUTE, which is associated with GMC identity, is nonmobile. The question is whether the KalaxMUTE could also specify SC identity by being mobile. To address this, we compared MUTE orthologues of the representative species with B. distachyon, A. thaliana and K. laxiflora to test potential mobility motifs in K. laxiflora (Fig. 8). Our results show high conservation in the bHLH functional domain. The differences in potential mobility residues of KalaxMUTE

Discussion
Stomatal patterning is diverse among different land plants. In Physcomitrella patens, stomata exhibit partial or complete division to form a single GC or paired GCs, respectively 38 . Moss does not have genes encoding MUTE or SPCH and uses genes encoding two bHLH proteins, PpSMF1 and PpSCRM1, to promote stomatal formation 39 . In A. thaliana, the stomata are surrounded by two kidney-shaped guard cells, and polar localization of BASL is required for a series of asymmetric divisions to form the stomatal structure 40 . In O. sativa, polar localization of PAN protein is responsible for subsidiary cell asymmetry in the stomatal complex 10 . In B. distachyon, BdMUTE is necessary and sufficient for SC formation. However, AtMUTE in A. thaliana defines GC precursor fate 23 .
Overall, it appears that the function of most genes is conserved during stomatal formation across plant evolution, but there are novel genes recruited to regulate unique aspects of stomatal patterning in some species.
The regulatory machine of stomata development appeared to be flexible and adaptable during evolution. The adaptation pressure could quickly change the division and differentiation pattern during stomata formation. For example, all the genes involved in stomatal differentiation are lost in seagrass Zostera to enhance its adaptation to marine lifestyle 30 . Plants of the ANITA grade form specialized structures in the epidermal cells to adapt to its habitat 29 . Similarly, N. colorata has lost genes, which could be associated with its unique stomatal development. However, further molecular and genetic manipulations are needed for functional verification.
Compared with our understanding of stomatal development in model systems, little is known about the molecular evolution of stomatal morphology, particularly i Diagrams illustrating stomatal development for the stomatal complex. Cell protoderm files (pale blue) asymmetrically divide to create a meristemoid (green), and the meristemoid differentiates into a GMC (orange). Then, neighbouring cell files (SMC, pale purple) divide asymmetrically to form SCs (blue). Finally, the GMC divides once symmetrically to form GCs (red), and the GCs and SCs terminally differentiate and form mature dumbbell-shaped stomata. Key: protodermal cell that will give rise to the stomatal lineage, pale blue; MMC (meristemoid mother cell), pale green; meristemoid, green; SLGCs (stomatal-lineage ground cell), white; GMC (guard mother cell), orange; GCs (guard cells), red; SMC (subsidiary mother cell), pale purple; SCs (subsidiary cells), blue in basal angiosperms. Alongside the completion of the genome, we are beginning to find the comparative molecular basis of the evolution of stomatal development and identify orthologues of stomatal regulator genes in a selected range of phylogenetic taxa. However, it is still technically difficult to analyse the function of orthologues. In the N. colorata genome, we found that a number of the genes that are highly specific to the stomatal asymmetric division were missing. Taken together, these results suggest that most core regulators of stomata formation remain conserved during evolution, whereas some gene loss events can occur to modify stomata formation processes, such as asymmetric division. These changes at the genetic and morphological levels of individual species may result from adaptation to inhabitant environments rather than evolutionary changes.
Recent studies have indicated that WGD events are ubiquitous in the evolution of angiosperms, and WGDs tend to retain multiple family duplications to increase the frequency of multiplication and the function of genes 41 . Thus, WGDs are widely thought to provide genomic novelties and complexities to promote plant adaptation to environments 42 . Large-scale GDs involved in stomata development through WGDs in K. laxiflora have been identified 36 .
Analysis of the genes involved in stomata formation showed that the protein sequences of the core genes Fig. 8 Alignment of grass and eudicot MUTE orthologues to identify potential mobility residues. MUTE orthologues of the representative grass species Brachypodium (BdMUTE-Bradi1g18400) and rice (OsMUTE-LOC_Os05g51820) were aligned with the MUTE orthologues of the representative eudicot species Arabidopsis (AtMUTE-AT3G06120) and Kalax.0004s0103/Kalax.0418s0025/Kalax.0268s0032/Kalax.0539s0032 using ClustalW (http://www.genome.jp/tools-bin/clustalw). The bHLH domain spans the first 50 amino acids and is indicated. Green shaded amino acids represent high similarity, whereas yellow shaded amino acids represent intermediate similarity. Candidate amino acids that are either consistently different between grasses and eudicots or are conserved among grasses but not in eudicots, or vice versa, are marked with a red asterisk and represent potential mobility motifs required to instigate and pattern stomata are conserved in K. laxiflora (Table 1). It is unclear whether the expression or protein modification of these regulators is different in K. laxiflora compared with that in A. thaliana. Indeed, the duplication of stomata regulator genes appears to be a common theme in K. laxiflora, but the extent to which this represents a divergence in gene function requires further studies.
It seemed that genes encoding critical developmental regulators were more likely to be retained during evolution 43,44 . For stomatal development, subsidiary cells can occur from an adjacent cell file or the same cell as the guard cells. Based on sequence conservation, the mobility of KalaxMUTE could be similar to its homologue in Arabidopsis. Thus, it is less likely that the modification of KalaxMUTE leads to featured stomatal subsidiary cells in K. laxiflora. Further work is needed to investigate whether the gene gains in K. laxiflora are associated with subsidiary cell establishment.