An ancestral stomatal patterning module revealed in the non-vascular land plant Physcomitrella patens

The patterning of stomata plays a vital role in plant development and has emerged as a paradigm for the role of peptide signals in the spatial control of cellular differentiation. Research in Arabidopsis has identified a series of epidermal patterning factors (EPFs), which interact with an array of membrane-localised receptors and associated proteins (encoded by ERECTA and TMM genes) to control stomatal density and distribution. However, although it is well-established that stomata arose very early in the evolution of land plants, until now it has been unclear whether the established angiosperm stomatal patterning system represented by the EPF/TMM/ERECTA module reflects a conserved, universal mechanism in the plant kingdom. Here, we use molecular genetics to show that the moss Physcomitrella patens has conserved homologues of angiosperm EPF, TMM and at least one ERECTA gene that function together to permit the correct patterning of stomata and that, moreover, elements of the module retain function when transferred to Arabidopsis. Our data characterise the stomatal patterning system in an evolutionarily distinct branch of plants and support the hypothesis that the EPF/TMM/ERECTA module represents an ancient patterning system. Summary: The genetic module controlling patterning of stomata in vascular plants also functions in non-vascular plants, consistent with the idea that it represents an ancestral mechanism in plant evolution.


INTRODUCTION
Stomata are microscopic pores present in the epidermis of all angiosperms and the majority of ferns and bryophytes. Evolution of stomata has proved to be an essential step in the success and diversification of land plants over the past 400 million years (Beerling, 2007). In particular, this innovation, coupled with vascular tissues and a rooting system, enabled land plants to maintain hydration by regulating the plant-soil-atmosphere water flows under fluctuating environmental conditions (Berry et al., 2010;Raven, 2002;Vatén and Bergmann, 2012). Stomatal distribution is tightly regulated, both by endogenous developmental mechanisms that influence their number and pattern in different organs of the plant, and by modulation of these controls by a host of environmental factors Geisler et al., 1998;Hunt and Gray, 2009;MacAlister et al., 2007). This spatial control of stomatal distribution, combined with the ease of scoring phenotype on the exposed epidermis, makes them an attractive system to investigate the control of patterning in plants, a major topic highlighted in the seminal work by Steeves and Sussex (1989).
Extensive molecular genetic analyses in the model flowering plant Arabidopsis have provided significant insight into the mechanisms controlling stomatal patterning and differentiation in angiosperms Engineer et al., 2014;Pillitteri and Torii, 2012;Simmons and Bergmann, 2016). In Arabidopsis, negatively and positively acting secreted peptide signals [epidermal patterning factors (EPFs) and epidermal patterning factor-like proteins (EPFLs)] function to control where and when stomata form and ensure that stomata are separated from each other by at least one intervening epidermal cell, thus optimising leaf gas exchange (Abrash and Bergmann, 2010;Hara et al., 2007Hara et al., , 2009Hunt et al., 2010;Hunt and Gray, 2009;Sugano et al., 2010). This 'one cell spacing rule' results from the stereotypical local pattern of cell divisions by which stomata form, accompanied by cross-talk between cells. The molecular mechanism enforcing the spacing rule involves EPF(L)s interacting with transmembrane receptors, including those encoded by members of the ERECTA gene family (ERECTA, ER; ERECTA-LIKE1, ERL1; and ERECTA-LIKE2, ERL2) activity of which is modulated in stomatal precursor cells by the receptor-like protein TOO MANY MOUTHS (TMM) (Lee et al., 2015(Lee et al., , 2012Shpak et al., 2005;Torii, 2012). Binding of EPF(L)s entrains a well-characterised signal transduction pathway involving a series of mitogen-activated protein kinases, which leads to the cellular events of stomatal differentiation (Torii, 2015).
Little is known of the developmental mechanisms regulating stomatal patterning in early land plants. Fossil cuticles of 400million-year-old small branching leafless vascular land plants, such as Cooksonia, indicate stomata were generally scattered more or less evenly across stem surfaces without clustering (Edwards et al., 1998) and these authors report that in the Rhynie Chert fossil plants stomata commonly occur on 'an expanded portion of the axis just below the sporangium'. These observations suggest the existence of a stomatal patterning module early in land plant evolution but we have very limited information on the nature of the genetic module controlling this process. However, homologues of key genes regulating vascular land plant stomatal differentiation are present in the genome and are expressed during sporophyte development in the moss Physcomitrella patens O'Donoghue et al., 2013;Ortiz-Ramírez et al., 2016;Vatén and Bergmann, 2012), a basal non-vascular land plant lineage with stomata. This suggests that genetic components involved in regulating stomatal spacing have been conserved between mosses and vascular plants. This notion is further supported by complementation work performed in Arabidopsis showing that Physcomitrella patens group 1A basic helix-loop-helix transcription factors can at least partially fulfil the function of their angiosperm counterparts in the regulation of stomatal development (MacAlister and Bergmann, 2011).
Here, we use molecular genetics to compare stomatal patterning systems in a bryophyte (Physcomitrella patens) and an angiosperm (Arabidopsis thaliana). We show that P. patens has an EPF/TMM/ ERECTA module required for stomatal patterning fundamentally similar to that found in angiosperms and that elements of the module retain function when transferred to Arabidopsis. Our data characterise the stomatal patterning system in moss and are consistent with the hypothesis that the EPF/TMM/ERECTA module represents an ancient patterning system in plants.

RESULTS
To identify potential orthologues of angiosperm genes implemented in stomatal patterning in P. patens, we performed a bioinformatic analysis. As shown in Fig. 1A and Fig. S1A, a single homologue of Arabidopsis EPF1 and EPF2 exists in P. patens, PpEPF1 (see also Takata et al., 2013). Similarly, the stomatal patterning protein TMM (which is encoded by a single gene in Arabidopsis) is homologous to a single gene in P. patens, termed PpTMM (Peterson et al., 2010) ( Fig. 1C; Fig. S1B). The situation with the ERECTA genes is more complicated as six potential orthologues are found in the genome of P. patens (Villagarcia et al., 2012) (Fig. 1E; Fig. S1C).
To identify genes potentially involved in stomatal patterning, we first interrogated a microarray database (O'Donoghue et al., 2013) to ascertain which PpERECTA genes showed upregulation of expression in the developing sporophyte. All the PpERECTA genes were expressed to some level in the sporophyte but only PpERECTA1 was upregulated relative to protonemal tissue (Fig. S2A), and qRT-PCR analysis confirmed that PpERECTA1 expression was significantly upregulated in the sporophyte (Fig. S2B). This was further indicated by the analysis of two other transcriptomic data sets accessible via phytozome V11 and the eFP browser at bar.utoronto.ca (see Table S1 for accession numbers), which showed a relatively high level of PpERECTA1 expression in the sporophyte (Fig. S2C,D) (Goodstein et al., 2012;Ortiz-Ramírez et al., 2016;Winter et al., 2007). Taken together, the data suggested that PpERECTA1 expression was increased in the sporophyte and, thus, might be involved in stomatal patterning. As shown in Fig. 1B,D,F, analysis of eFP Browser data for PpEPF1, PpTMM and PpERECTA1 indicated an accumulation of the relevant transcripts in young sporophyte tissue.
Having identified genes encoding homologues for each of the components of the core EPF/TMM/ERECTA module involved in angiosperm stomatal patterning, we undertook a functional analysis in P. patens by creating a series of gene knockouts and analysing stomatal patterning in the sporophytes of the transgenic plants.
Interruption of the targeted locus in transgenic plants was confirmed by genomic PCR (Fig. S3). As shown in Fig. 2A-F, loss of PpEPF1 function led to an increase in the number of stomata per capsule. The extra stomata formed at the appropriate location at the base of the sporophyte ( Fig. 2A,B), i.e. they did not extend ectopically into the flanks of the spore capsule. As a consequence, stomata in ppepf1 knockout mutant capsules frequently occurred in clusters that were not apparent in wild-type (WT) sporophytes, where most stomata are separated from each other by at least one neighbouring epidermal cell (Fig. 2C,D). Quantification confirmed an increased number of stomata per capsule in the sporophytes of three independently generated ppepf1 knockout lines (Fig. 2E). Expression analysis confirmed the absence (lines ppepf1-2, ppepf1-3) or greatly decreased transcript level ( ppepf1-1) for PpEPF1 in these plants (Fig. 2F). Interruption of the targeted locus in transgenic plants was verified by genomic PCR (Fig. S3).
We also characterised the outcome of increased expression of PpEPF1 on stomatal formation by creating lines of transgenic P. patens in which the PpEPF1 coding sequence was constitutively overexpressed via the rice actin promoter (Fig. 2L). Sporophytes of the transgenic plants displayed a phenotype with a greatly reduced number of stomata. At the base of the sporophyte stomata were sporadic (Fig. 2G,H) and the number of stomata per capsule significantly decreased in three independent lines overexpressing PpEPF1 (Fig. 2K). Although the number of mature stomata was clearly decreased in the plants overexpressing PpEPF1, analysis of the epidermis at the base of the sporophytes of the transgenic plants indicated occasional division patterns suggestive of the formation of stomatal precursors that had failed to undergo further differentiation into the stomatal lineage (compare Fig. 2I and 2J).
To investigate the role of the PpTMM receptor, we generated independent knockout lines. Examples of the range of phenotypes observed are shown in Fig. 3B-D for comparison with the WT pattern (shown in Fig. 3A). Some capsules had exceptionally few stomata (Fig. 3C) whereas others developed numerous stomata, many of which occurred in clusters (Fig. 3D). This variation was consistently observed across all three independent pptmm knockout lines. Again, as with the ppepf1 knockout and WT lines, stomata formation remained restricted to the base of the capsule. Quantification of the transgenic sporophytes revealed that the number of stomata per capsule tended to be lower in the pptmm knockout lines than in the WT control, although this was statistically significant only in the line pptmm-3 (Fig. 3E). When the proportion of stomata forming in clusters (defined as stomata forming in pairs or higher order adjacent complexes) was measured, it was apparent that the pptmm knockout lines had a higher number of stomata in clusters than WT (Fig. 3F). Interruption of the targeted locus in transgenic plants was verified by genomic PCR (Fig. S3) and expression analysis confirmed that the three pptmm knockout lines contained no detectable PpTMM transcript (Fig. 3G).
We further investigated the role of PpTMM by analysing transgenic P. patens in which the PpTMM sequence was constitutively overexpressed (Fig. 3M). For this part of the investigation we were only able to identify a single transgenic line but analysis suggested that an increased level of PpTMM transcripts had little effect on stomatal patterning. There was a slight increase in the number of stomata per capsule (Fig. 3H,I) but quantification indicated that this was not statistically significant (Fig. 3L). The extent of stomatal clustering was similar to that observed in WT sporophytes (Fig. 3J,K).
To ascertain whether P. patens requires ERECTA gene functioning during stomatal development, we next targeted the PpERECTA1 gene. Only a single PpERECTA1 knockout line was identified and, as shown in Fig. 4A,B, stomata formed in the appropriate position at the base of the sporophyte with no obvious difference in stomatal differentiation (Fig. 4C,D) and no effect on stomatal number per capsule (Fig. 4E). Loss of PpERECTA1 gene expression in this line was confirmed by RT-PCR (Fig. 4F), as was interruption of the targeted locus by genomic PCR (Fig. S3). Because analysis of the pperecta1 knockout was unable to establish a conclusive role for this component in stomatal development, further experiments were carried out. To understand whether the PpTMM and PpEPF1 genes were acting in the same pathway as PpERECTA1 during stomatal development, a series of double knockout mutants were produced. Analysis of ppepf1-erecta1 double knockouts indicated a diminished ppepf1 phenotype. Thus, although more stomata per capsule developed compared with WT the increase was less than in the ppepf1 mutant (Fig. 4G). An even more dramatic effect was observed when pptmm-epf1 double knockouts were generated. In this situation, the phenotype of increased stomata per capsule observed in the ppepf1 knockout was found to be entirely dependent on the presence of a functional PpTMM gene (Fig. 4H). Finally, a pptmm-pperecta1 double knockout displayed a greater decrease in stomata per capsule than observed in the single pptmm and pperecta1 mutants (Fig. 4I). Analysis of epidermal regions of capsules from the different knockout combinations (Fig. 4J-O) suggested that, in addition to the differences in stomata number, loss of some EPF/TMM/ERECTA gene combinations influenced the positioning/form of stomata and the general pattern of cell division in the epidermis. For example, although loss of PpERECTA1 in a ppepf1 background led to a (A,C,E) Phylogenetic trees constructed using amino acid sequences of selected Arabidopsis EPF1 (A), TMM (C) and ERECTA (E) gene family members based on Phytozome V11 (Goodstein et al., 2012), using the neighbour-joining method (Saitou and Nei, 1987;Takata et al., 2013) on MEGA6 (Tamura et al., 2013). The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches (Felsenstein, 1985). Amino acid sequences from P. patens (Pp), Selaginella moellendorffii (Sm), Zea mays (Zm), Symphytum tuberosum (St), Medicago truncatula (Mt) and A. thaliana (At) were used to generate trees, except for ERECTA, for which S. moellendorffii and S. tuberosum gene family members were omitted, owing to the large overall number of genes in the ERECTA family. decrease in stomatal number per capsule, there was often an apparent disruption to the epidermal cell patterning in the vicinity of the stomata that were formed (Fig. 4J,M). In the case of the ppepf1pptmm double knockouts (which restored stomata number per capsule to wild-type levels), there were also alterations to the epidermal cell division planes from the patterns observed in the ppepf1 or pptmm capsules (Fig. 4J,K,N). Furthermore, guard cells that formed at the boundaries appeared stretched, taking on a shape akin to neighbouring epidermal pavement cells. This elongated guard cell conformation was also seen in the pptmm-pperecta1 line (Fig. 4L,O). Analysis of the various double knockout combinations described above confirmed the absence of the relevant transcripts ( Fig. 4P-R) and interruption of the targeted loci (Fig. S3).
In addition to suppression of stomata formation, some EPF-like peptides in angiosperms have evolved to inhibit EPF action competitively (Lee et al., 2015;Ohki et al., 2011). Most notably, AtEPFL9 (STOMAGEN) has been shown to enhance stomata formation in Arabidopsis and overexpression of STOMAGEN leads  to increased stomatal number in Arabidopsis (Hunt et al., 2010;Sugano et al., 2010). As indicated in Fig. 1 and Fig. S2A, bioinformatic analysis indicates that the P. patens genome encodes a peptide similar to EPF2 (which in Arabidopsis inhibits stomatal development), but has no apparent equivalent to STOMAGEN (which in Arabidopsis antagonises the activity of EPF2 and stimulates stomatal development). To investigate whether the stomatal patterning system in P. patens could be disrupted by overexpression of the evolutionarily distinct antagonistic peptide STOMAGEN (Fig. S4A), we overexpressed the Arabidopsis STOMAGEN gene in a WT P. patens background. Our results indicated that although the STOMAGEN transcripts accumulated to a high level (Fig. S4G), there was no apparent phenotype in terms of altered numbers of stomata per capsule (Fig. S4B-D) although abnormal epidermal cells and aberrant guard cells were occasionally observed at the base of the transgenic capsules (Fig. S4E,F).
Because our data suggested that PpEPF1, PpTMM and PpERECTA1 all play a role in stomatal patterning in P. patens, we investigated whether they might represent conserved functions by introducing the P. patens genes into the appropriate Arabidopsis genetic background, i.e. could they complement the cognate angiosperm gene function in stomatal patterning? For this experiment, we focused on the putative EPF and TMM orthologues because the respective mutants in Arabidopsis have clear phenotypes with respect to stomatal density and patterning. In leaves, loss of AtEPF1 or AtEPF2 results in increased stomatal density, with stomatal clustering being especially pronounced in atepf1 (Hara et al., 2007). In atepf2, increased density is the result of increased entry of cells to the stomatal lineage, which causes not only more stomata but also more small epidermal stomatal precursor cells (Hara et al., 2009;Hunt and Gray, 2009). In Arabidopsis tmm, stomatal phenotype varies depending on the organ. For example, in leaves stomatal density and clustering is markedly increased, in tmm inflorescence stems no stomata are found, and in the flower pedicel a gradient of stomatal density is observed. Thus, at the base of tmm pedicels there are no stomata, in the middle region a few stomata form, and at the apical region of the pedicel stomatal density exceeds that of wild type and clustering is common (Bhave et al., 2009;Geisler et al., 1998;Yang and Sack, 1995).
When PpEPF1 was constitutively overexpressed in atepf1 we found that the mutant phenotype was partially rescued, with leaves having stomatal densities that were lower than in the atepf1 background and which approached wild-type values (Fig. 5A). When PpEPF1 was overexpressed in the atepf2 mutant background only a slight recovery of stomatal density occurred (Fig. 5A) and epidermal cell density was essentially unchanged relative to that observed in atepf2 (Fig. S5). With respect to PpTMM, when this sequence was expressed in the Arabidopsis attmm mutant under control of an endogenous AtTMM promoter there was no overt restoration of stomatal density to wild-type values in leaves and stomatal clustering was still observed (Fig. 5B). However, when the pedicel was examined there was a partial rescue of the attmm phenotype. This was most obvious in the middle region where stomatal density was restored towards wild-type values whereas in the basal and apical regions stomatal densities were similar to those observed in the tmm mutant (Fig. 5C).

DISCUSSION
The control of patterning is core to development and the EPF/TMM/ ERECTA module has emerged as a paradigm for peptide signalling in plants to control the distribution of essential cellular complexes on the epidermis, the stomata. Although it is well-established that stomata arose very early in the evolution of land plants, until now it has been unclear whether the angiosperm stomatal patterning system represents an ancient, universal mechanism in the plant kingdom. Our data indicate that an essentially similar system functions in the moss P. patens, providing strong evidence that the EPF/TMM/ERECTA module represents an ancestral patterning system for stomata. Mosses and flowering plants last shared a common ancestor over 400 million years ago (Ruszala et al., 2011), suggesting that the leafless sporophytes of early vascular land plants may have deployed a patterning module comprising genes closely related to the EPF/TMM/ERECTA suite identified here.
Our data establish, firstly, that the genome of an extant bryophyte, P. patens, contains sequences homologous to all three components of the EPF/TMM/ERECTA module present in the angiosperm A. thaliana and that they are expressed at an appropriate time in development to play a role in stomatal patterning. Two of these components (PpEPF1 and PpTMM) are present as single-copy genes, consistent with them representing relatively ancient ancestral forms. The situation with the ERECTA gene family was more complicated, but our expression analysis, including the analysis of staged, dissected sporophyte tissue, allowed us to identify one member of the ERECTA family in P. patens (PpERECTA1) that was expressed at the appropriate time and place to play a role in stomata formation and which was therefore selected for further investigation (O'Donoghue et al., 2013). Functional analysis of these three genes (PpEPF1, PpTMM, PpERECTA1) indicated that they are all involved in stomatal patterning in P. patens with roles not dissimilar to those played by their putative orthologues in Arabidopsis (Geisler et al., 1998;Hara et al., 2007Hara et al., , 2009Hunt and Gray, 2009;Shpak et al., 2005;Yang and Sack, 1995). This was clearest with PpEPF1. Loss of this peptide led to an increase in stomatal clustering and in number of stomata per capsule whereas overexpression led to a decrease in stomata per capsule. This indicates a function directly comparable to that observed for AtEPF1 and, to a lesser extent, AtEPF2 in Arabidopsis where loss of function leads to an increase in leaf stomatal density and stomatal clustering (Hara et al., 2007(Hara et al., , 2009Hunt and Gray, 2009).
With respect to PpTMM, a more complicated picture emerged, consistent with the context-dependence of the tmm phenotype reported in Arabidopsis (Geisler et al., 1998;Yang and Sack, 1995). For example, in mature leaves of the Arabidopsis attmm mutant stomatal clustering is apparent and stomatal density is higher than that of WT, whereas at the base of the pedicels no stomata form, in the middle region some stomatal formation occurs (with some clustering), and at the top of the pedicel ectopic stomata form, leading to increased density and clustering relative to the WT (Bhave et al., 2009;Geisler et al., 1998). In the sporophyte of P. patens (where stomata only form at the base of the spore capsule), we observed an overall trend for a decrease in stomatal density and increase in clustering in the pptmm lines. However, these average values obscure significant spatial variation even within single spore capsules, so that on a given capsule it was not uncommon to observe both stomatal clustering and adjacent areas devoid of stomata, i.e. the phenotype encompassed elements observed on leaves and pedicels in Arabidopsis. The mechanistic basis of this variation awaits elucidation but the data indicate an overall conservation of sequence and function for TMM in P. patens and Arabidopsis and support an important ancestral role for this protein in the modulation of stomatal patterning in leafless early land plants. It is possible that the regulation of stomatal stochasticity by TMM in early land plants enabled or facilitated the later evolution of distinct stomatal patterns in different parts of the plant. One possibility is that other peptides in P. patens, encoded by genes similar to Arabidopsis EPFL6 (CHALLAH), EPFL5 (CHALLAH-LIKE 1) and EPFL4 (CHALLAH-LIKE 2), play a role in inhibiting stomatal formation in the absence of PpTMM, as is the case in Arabidopsis (Abrash and Bergmann, 2010;Abrash et al., 2011). A recent bioinformatics study has identified nine PpEPFL (CHALLAH-like) genes that are upregulated in the developing sporophyte (Ortiz-Ramírez et al., 2016;Takata et al., 2013). These genes represent a target for future work to provide a deeper understanding of peptide signalling and stomatal patterning.
In Arabidopsis, TMM modulates the activity of ERECTA proteins and the action of AtEPF2 (and possibly AtEPF1) is dependent on TMM (Hara et al., 2007;Hunt and Gray, 2009;Lee et al., 2012). To test whether PpEPF1 action requires PpTMM, we produced double mutant lines ( pptmm-epf1) and found that the ppepf1 phenotype was masked. Thus, there is an epistatic interaction between PpTMM and PpEPF1 similar to the situation reported in Arabidopsis for EPF1 or EPF2 and TMM (Hara et al., 2007(Hara et al., , 2009Hunt and Gray, 2009).
Our analysis of P. patens sporophytes lacking PpERECTA1 expression indicated no difference in stomatal number and only a slight difference in the spacing of stomata. As only one knockout line could be assessed we emphasise that this result should be interpreted with caution. However, analysis of pperecta1 in combination with either ppepf1 or pptmm indicated a more pronounced role for PpERECTA1 in stomatal patterning. For example, loss of PpERECTA1 partially rescued the phenotype shown by the ppepf1 knockout mutant. The available data indicate that there are at least five other closely related PpERECTA genes expressed in the sporophyte ( Fig. 2C; Fig. S1C) (Villagarcia et al., 2012), so the lack of phenotype in the single PpERECTA1 knockout mutant may reflect a degree of genetic redundancy in a manner similar to the redundant activity of this receptor family in Arabidopsis (Shpak et al., 2005). Further analysis of these PpERECTA genes in the context of pptmm and ppepf1 mutants may improve our insight into the role of these genes in stomatal development.
Our experiments demonstrated that for both EPF1 and TMM conservation of function extends across the evolutionary distance separating bryophytes and angiosperms, with expression of PpEPF1 and PpTMM coding sequences leading to a partial rescue of the mutant phenotype in the relevant Arabidopsis genetic backgrounds. Interestingly, overexpression of PpEPF1 in the atepf1 background was sufficient to restore stomatal number to near wild-type level whereas it was less able to rescue the related atepf2 mutant phenotype. AtEPF1 has been implicated in the spacing patterning of stomata whereas AtEPF2 is thought to be more important for the earlier asymmetric divisions required for angiosperm stomatal initiation (Hara et al., 2007(Hara et al., , 2009Hunt and Gray, 2009). Our data therefore support the idea of an ancient role for an EPF peptide ligand in stomatal patterning, with the evolution of the angiosperm EPF gene family being linked to acquisition of asymmetrically dividing cells in stomatal development (Hara et al., 2009;Hunt and Gray, 2009). In Arabidopsis, the EPF(L) gene family appears to have expanded over evolutionary time so that particular combinations of different ligands and receptors function in different organs. This divergence of EPF function linked to increased plant complexity is supported by the observed inability of the Arabidopsis STOMAGEN sequence to alter stomatal patterning in P. patens. The acquisition of such novel regulators of stomata formation may reflect an evolutionary trend to more complex developmental systems, enabling a flexible control of stomatal pattern to allow plants to adapt organs to specific environments (Abrash and Bergmann, 2010;Hronková et al., 2015;Hunt and Gray, 2009;Rychel et al., 2010;Shpak et al., 2005;Takata et al., 2013).
In conclusion, our results establish that the members of an EPF/ TMM/ERECTA ligand-receptor system are conserved between bryophytes and angiosperms, both in terms of the presence and expression of the relevant genes and in the functional conservation of their role in stomatal patterning. Our data do not provide information on the conservation (or otherwise) of the molecular interactions between the components of the EPF/TMM/ERECTA module (which have only recently become well-described in the more highly studied Arabidopsis system; Lee et al., 2015Lee et al., , 2012 and this represents an area for future investigation. Finally, the acquisition of stomata is recognised as being of fundamental importance in the evolution of land plants (Beerling, 2007;Berry et al., 2010) and our data strongly support the proposition that the genetic system regulating stomatal patterning was recruited at an extremely early stage of land plant evolution, supporting the idea that extant stomata are of monophyletic origin (Beerling, 2007;Edwards et al., 1998).

Bioinformatic analysis
Sequences were aligned using the MUSCLE (Edgar, 2004) alignment tool on MEGA6 (Tamura et al., 2013). The evolutionary history was inferred using the neighbour-joining method (Saitou and Nei, 1987) on MEGA6 (Tamura et al., 2013). The optimal trees with the sum of branch length 13.422 (EPF), 9.345 (TMM) or 20.7487 (ERECTA) are shown. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches (Felsenstein, 1985). The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Poisson correction method (Zuckerkandl and Pauling, 1965) and are in the units of the number of amino acid substitutions per site. For EPF(L) representatives, the analysis involved 79 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 33 positions in the final dataset. For TMM, the analysis involved 31 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 185 positions in the final dataset. For ERECTA, the analysis involved 82 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 24 positions in the final dataset. See Table S1 for accession numbers not indicated in the trees.

P. patens gene manipulation and expression analysis
The 5′ and 3′ flanking regions of targeted genes were amplified from P. patens genomic DNA and inserted into plasmids by conventional cloning using primers detailed in Table S2. Resulting plasmids were used as PCR templates to amplify knockout constructs. PpEPF1 was blunt-end ligated into EcoRV-digested pKS-Eco, then BsoBI digested, and a hygromycin selection cassette (obtained from pMBLH6bI) blunt-end ligated between 5′ and 3′ flanking regions to produce the Ppepf knockout construct. The pptmm and pperecta1 knockout constructs were both created by blunt-end ligating 5′ flanking sequences into Ecl136II-digested pMBL5DLdelSN (a pMBL5 derivative) containing the NPTII cassette. Resulting plasmids were digested with EcoRV and 3′ flanking sequences inserted via blunt-end ligation. To target the PpEPF1 locus in the pptmm-1 background, the Ppepf knockout construct was used. To target the PpERECTA1 locus in the pptmm-1 background, the NPTII cassette in the pperecta construct was replaced with an HPH cassette at KpnI and NsiI sites to produce a hygromycin-selective pperecta knockout construct.
Gene targeting and polyethylene glycol-mediated transformation of P. patens was performed using PCR-derived templates (Kamisugi et al., 2005;Schaefer et al., 1991). Confirmation of integration at target site was performed by genomic PCR analysis (Fig. S3). Briefly, for each independent line PCR was performed targeting a fragment spanning the 5′ genomic sequence to the transgene resistance cassette (Fig. S3B,D,F,H,J) or the 3′ genomic sequence to the transgene resistance cassette (Fig. S3C,E, G,I,K). For each gene knockout, either two or three independent transgenic lines were generated and analysed, with the exception of pperecta1 and pptmm-erecta1 for which only one line was obtained showing correct gene targeting and no expressed transcript. Expression of transgenes and absence of expression of targeted knockout genes was determined by RT-PCR using single-stranded cDNA generated from extracted RNA by M-MLV Reverse Transcriptase (Fisher Scientific). RNA was extracted using Spectrum Plant Total RNA Kit (Sigma). For expression analysis in ppepf1 and pptmm lines, 120 developing sporophytes per line were harvested and used to extract RNA. For other RT-PCR analyses, gametophyte-sporophyte mix samples were collected for each line, which contained 25 gametophores and ∼15 developing sporophytes. For quantitative RT-PCR analysis of PpERECTA1 transcript, relative expression was compared between RNA extracted from protonemal versus pooled sporophyte tissue (∼300 capsules per replicate: 100 immature, 100 mid-sized and 100 fully expanded sporophytes) in triplicated experiments. RNA integrity was verified by electrophoresis and NanoDrop ND-8000 (Fisher Scientific) and 1 µg RNA used in reactions alongside three control 'housekeeping' transcripts (Le Bail et al., 2013;Wolf et al., 2010), according to Luna et al. (2014) with slight modifications. Transcript abundance was assayed using Rotor-Gene SYBR Green PCR Kit and a Corbett Rotor Gene 6000 (Qiagen).

Plant phenotyping
Fully expanded (orange to brown coloured) spore capsules were fixed in modified Carnoy's solution (2:1 ethanol: glacial acetic acid) 6 to 7 weeks after fertilisation by flooding. Capsules of a similar size were dissected to remove associated spores, mounted between a bridge of cover slides in distilled H 2 O and stomata imaged with an Olympus BX-51 microscope fitted with an Olympus DP71 camera and Olympus U-RFL-T-200 UV lamp equipped with an LP 400 nm emission filter. Multiple fields of view were stacked and colour corrected using ImageJ. Min Intensity (bright-field) or Max Intensity (fluorescence) settings were used to compile flattened images.
Fully expanded Arabidopsis leaves were collected 7 to 8 weeks after germination, abaxial epidermal impressions produced and stomatal densities taken from two to three fields of view per leaf (Hunt and Gray, 2009). Pedicels were collected from 14-week-old plants, fixed and cleared in modified Carnoy's solution, dissected longitudinally, rinsed in 0.5% diphenylboric acid-2-aminoethyl ester (DPBA) (Sigma-Aldrich) and 0.1% Triton X-100 (v/v) for 30 s, then mounted as above. Images were collected using bright-field on an Olympus BX-51 microscope with accompanying 400 nm fluorescence ( pE-2 UV, CoolLED, Andover, UK) and 455 nm emission filter to capture fluorescence and stacked using ImageJ. Stomata were counted in areas of 180.26×262.56 µm approximately 300 µm from where pedicels were excised at the base, halfway up the stem and 150 µm from the tip. T2 and homozygous T3 plants were phenotyped. Statistical tests were performed using Graphpad Prism6 and graphs were produced using SigmaPlot version 13 (Systat Software).