The moss P. patens (Hedw.) Bruch & Schimp accession Gransden was used in this study. Plant material was axenically cultivated in liquid minimal medium (250 mg/L KH₂PO₄, 250 mg/L KCl, 250 mg/L MgSO₄ 7H₂O, 1,000 mg/L Ca(NO₃)₂ 4H₂O, and 12.5 mg/L FeSO₄ 7H₂O) under standard growth conditions as described previously (Frank et al., 2005a). After 4 days protonemal tissue was disrupted for 10 s with an Ultra-Turrax device (IKA, Staufen, Germany) and transferred into fresh liquid medium.

(+)-cis, trans-abscisic acid (ABA; Duchefa) was added to freshly homogenized protonemal liquid cultures that were adjusted to 100 mg/L dry weight. ABA was dissolved in 100 μM KOH and mock treatments were performed with 100 μM KOH only.

Light and epifluorescence microscopy was performed either with an Olympus BX41 microscope or with a Zeiss Axioplan and images were taken with a Canon EOS D30 camera. Confocal microscopy was performed with a Zeiss LSM 510 UV microscope, excitation 543 nm; emission spectra for chlorophyll and propidium iodide were separated by linear unmixing. For cell wall staining, propidium iodide (Fluka) was added to a final concentration of 4% and slides were washed with medium after 5 min incubation.

Cryo-SEM analysis was performed using a Philips XL30 ESEM with Cryo Preparation Unit Gatan Alto 2500. Plant material was applied to a specimen holder with freeze hardening glue and biological samples were preserved by fast-freezing in liquid nitrogen. Afterward the specimen holder was inserted into the sputter chamber and coated with gold.

For TEM, after removal of supernatant, 2.5% glutaraldehyde in 50 mM cacodylate buffer (pH 5.7) was added to protonemal culture and incubated for 5 h. Samples were washed five times for 10 min with cacodylate buffer, embedded in low melting grade agarose (Sigma) and fixed with 1% OsO₄ for 4 h at 4°C. After washing with cacodylate and 2× water, samples were dehydrated using an ethanol series and finally washed two times each with ethanol and xylol and imbibed in xylene-epoxy resin 1:1 overnight. Samples were embedded in Epon 812 resin (Sigma) by incubation for 2 days at 60°C. Ultrathin slices (70–100 μm) were cut with a Leica microtome and observed in an electron microscope CM 10 (Philips).

RNA isolation, processing of RNA, microarray hybridization, washing and scanning as well as feature extraction and normalization were carried out as described previously (Wolf et al., 2010). Briefly, differentially expressed genes (DEGs) were called using Cyber-T (Baldi and Long, 2001) with false discovery rate (FDR) correction according to Benjamini and Hochberg (1995), using a q-value cutoff of 0.05. For comparison of different lists, Venn diagrams were created with the online tool Venny (Oliveros, 2007) while UpSet plots were generated as described previously (Khan and Mathelier, 2017).

For cDNA synthesis, 4 μg total RNA was treated with DNase I (Fermentas, Sankt Leon-Rot, Germany) for 1 h at 37°C in order to eliminate remaining genomic DNA. The enzyme was inactivated by heating at 65°C for 10 min. Reverse transcription was performed according to the TaqMan Applied Biosystems (Roche) manufacturer’s recommendations. qRT-PCR reactions with 90 ng of cDNA were performed with the SensiMix^TM SYBR No-Rox Kit (Bioline) according to the manufacturer’s recommendations. All qRT-PCRs were performed in three biological and three technical replicates. The qRT-PCR program was adjusted to initial denaturation and hot start at 95°C for 10 min followed by 40 cycles of amplification with 96°C for 15 s, 60°C for 30 s, and 72°C for 20 s. The SYBR Green signals were measured at each cycle and melting curves were calculated to prove primer specificities. Gene expression values were normalized to the PpEF1α control gene and the relative quantifications were calculated based on Advance Relative Quantification provided by the LightCycler^® 480 software release 1.5.0. Quantifications were based on the delta-delta Ct method.

The gene ontology (GO) enrichment analyses for P. patens and A. thaliana genes were performed with PANTHER (Mi et al., 2010) and the PANTHER database. PANTHER’s tool accesses a comprehensive list of GO annotations from the GO Consortium and is updated monthly to ensure the most current annotation data. To provide an overview on over-/under-represented GO terms Fisher’s exact test with subsequent Bonferroni correction was performed and q-values < 0.05 were considered as significant. For the GO clouds in Supplementary Figure S1, GO bias analyses were performed as in Widiez et al. (2014) based on the v1.6 GO association and visualized using Wordle¹. For Supplementary Figure S3 we used the data from this study as well as from a previous study employing the same microarray (Hiss et al., 2014). Genes with highest expression in either protonema treated 180 min with ABA (250) or brown sporophytes (55) were identified, so that maximum expression of all replicates over the complete array dataset had to have a lower expression than the minimum expression of the ABA treated/brown sporophyte replicates.

Global transcriptome analysis was calculated as previously described (Busch et al., 2013). The rationale behind the analysis is the idea that cells change their transcriptomes after stimulation in concerted fashion, i.e., there are hundreds to thousands of genes that are up- or down-regulated together, either immediately or in delayed fashion. In order to detect those genes that considerably contribute to the concerted change of the transcriptome, we first ranked all genes according to their variance in expression after stimulation starting/ending with those genes that have the strongest/weakest response to the stimulation. We then grouped the ranked list of genes into subsets of 50–500 genes and compared the change in mutual information (MI) and Pearson correlation (PC) after stimulation of these subsets with the same change of the whole transcriptome. Gene subsets that considerably contribute to the overall change in gene expression should behave similar to the whole transcriptome, as quantified by a minimal difference of the MI and PC between the respective subsets and the transcriptome. Thus, we consider all genes that show the maximal temporal response as well as those that are most correlated with the overall transcriptome change as important for the ABA response, which is depicted in Figure 6. In brief, the method calculates the Euclidean distance of the PC and the MI between the whole transcriptome and subsets of genes that have been ordered from large to small temporal response. The idea is to determine those genes that are (i) strongly regulated over time or (ii) whose temporal dynamics correlate well with the whole transcriptome (computed between the global and ordered gene subset response for subsets of different sizes). In an alternative approach we calculated the significance of differential regulation over time for each gene by fitting a reduced model to all time points using a 3rd order polynomial or a full model by fitting the ABA response and control separately. The p-value was then calculated from an analysis of variance between the models, with a small p-value indicating a differential dynamic response under treatment. P-values were FDR corrected using the R/Bioconductor package qvalue².

Studies on Funaria hygrometrica (belonging to the Funariaceae as P. patens) showed that exogenous application of ABA to protonemal filaments induces the formation of vegetative diaspores (brachycytes or brood cells) and preformed filament breaking points (tmema cells) (Schnepf and Reinhard, 1997). Brachycytes are highly drought tolerant and are still able to grow after 3 years of dormancy (Bopp and Werner, 1993). The number of formed brachycytes depends on the concentration and duration of the ABA application. Funaria brachycytes are characterized by thickened cell walls, small vacuoles, plasmatic lipid droplets and reduced plastidal starch (Schnepf and Reinhard, 1997). Tmema cells undergo programmed cell death (PCD) and support fractionation of the filaments so that the brood cells can be propagated (Bopp et al., 1991).

Although the action of ABA on P. patens has been considered to be similar to that in Funaria (Decker et al., 2006) and also induces brachycytes (Busch et al., 2013), concomitant ultrastructural changes have not been studied in detail. We thus applied ABA to protonemata of P. patens in order to analyze whether they respond in a similar fashion as in Funaria and to characterize cellular changes. Although the intracellular ABA concentration in P. patens is ∼1 μM, increasing to ∼10 μM under abiotic stress (Minami et al., 2003; Beike et al., 2015), ABA is usually applied in the 25–100 μM range in order to achieve developmental changes (Richardt et al., 2010; Busch et al., 2013). Here, we find that exogenous application of 100 μM ABA to protonema leads to the formation of rounded brachycytes and empty tmema cells within 2 weeks (Figure 1A–C). Intracellular changes induced by ABA include the disintegration of the central vacuole, an increase in plastidal starch granules and cytoplasmic oleosomes, and structural changes to the cell wall with an increase in thickness (Figure 1E–H). After 1 week of ABA treatment the cell wall was on average 43% thicker (53 ± 11 μm) than in the control (37 ± 8 μm), which is a significant increase (t-test, p < 0.01, Figure 1D,I,J). These analyses indicate that brachycytes in P. patens form similarly to those in Funaria, albeit with the notable difference that more starch granules are formed inside the plastids (Figure 1I,J). P. patens brachycytes develop from chloroplast-rich chloronema filaments are characterized by round shape, thick cell wall, loss of the large central vacuole, and increase of starch as well as lipid storage (Figure 1E,I,J). They are thus well suited as vegetative diaspores and represent the end point of the developmental progression triggered by the ABA signaling pathway.

We performed differential gene expression analysis upon ABA application using an established microarray system representing 27,828 out of 35,307 v1.2 genes (Wolf et al., 2010). Protonema grown in liquid culture was exposed to 10 μM ABA, since it was previously shown that 10 μM ABA is sufficient to induce a molecular response but is too low to cause phenotypical changes (Richardt et al., 2010). Transcriptional changes were analyzed at different time points after ABA treatment. We detected a rapid induction of 202 differentially expressed genes (DEGs) after 30 min, 663 genes after 60 min, and 936 genes after 180 min, while no down-regulated genes were detected after 30 and 60 min, and only eight genes were down-regulated after 180 min (Figure 2A and Supplementary Table S1). Five DEGs detected after 30 min were also detected at 60 min but not at 180 min, whereas 197 genes were induced at 30 min and maintained as up-regulated until 180 min of ABA treatment. 81 genes were transiently up-regulated after 60 min of treatment while 397 genes were identified as late ABA-responsive genes, detected only after 180 min of the treatment (Figure 2B). To evaluate the reliability of the microarray results, qRT-PCR was performed on five randomly selected genes that are members of transcription associated genes and represent different expression groups (early and late up-regulation, transient up-regulation, and down-regulation) and which confirmed the microarray data (Figure 2C).

Since we observed cell wall thickening in P. patens protonemal cells in response to ABA we searched for genes that might have roles in ABA-mediated cell wall modification. For this purpose, we prepared a list of genes that are either involved in cell wall biosynthesis or related to cell wall regulation from previously published P. patens studies. This list contains 71 pectin-related genes (McCarthy et al., 2014), 12 callose synthase genes (Schuette et al., 2009), 11 cellulose synthase genes (Goss et al., 2012), 9 arabinogalactan genes (Fu et al., 2007), 10 cutin synthase genes (Yeats et al., 2014), 30 expansins (Carey and Cosgrove, 2007), and 32 xyloglucan endotransglucosylase/hydrolase genes (Yokoyama et al., 2010; Supplementary Table S2). The list was compared to the identified DEGs and we found that six genes, Phypa_120256 (galactan galactosyltransferase), Phypa_207105 (Pectin lyase-like superfamily protein), Phypa_92683 (glycosyl transferase), Phypa_114674 (xyloglucan endotransglycosylase), Phypa_206446 (alpha-expansin), and Phypa_51702 (arabinogalactan protein) are differentially regulated upon ABA treatment (Table 1). This suggests that the proteins encoded by these genes might be involved in the observed cell wall thickening.

We analyzed whether there is any effect of ABA application on the expression of genes encoding enzymes involved in ABA biosynthesis or proteins acting in ABA signaling. P. patens ABA biosynthetic, “core” regulatory components and ABA signaling genes have been identified previously (Stevenson et al., 2016). We calculated fold changes of all genes at 30, 60, and 180 min of ABA application. From this analysis it is obvious that genes encoding enzymes of the ABA biosynthesis pathways are not affected, except two ABA-induced NCED genes that are known to be regulated by ABA in other plant species and encode the enzyme that catalyzes the rate-limiting step of ABA biosynthesis (Sussmilch et al., 2017). Furthermore, almost all genes known to be involved in ABA-dependent signaling and regulation except those encoding the putative ABA receptors and ABI4 were found to be upregulated (Table 2).

To examine in which biological process, cellular component and molecular function the ABA responsive genes are involved, GO analysis was performed. Functional categories with FDR corrected p-values (q-values) < 0.05 are shown in Supplementary Table S3. Enriched terms in “biological process” include response to water deprivation, negative regulation of proteolysis, endopeptidase, peptidase, hydrolase and catalytic activity, monocarboxylic acid metabolic process, lipid metabolic process, and oxidation-reduction process. In the GO “molecular function” category, the most enriched terms are endopeptidase inhibitor activity, enzyme inhibitor activity, oxidoreductase activity and catalytic activity. These results show that plant responses that are required for stress adaptation are enhanced after ABA treatment. In the GO “cellular component” category, the most abundant terms are extracellular space, endoplasmic reticulum, vacuole, organelle membrane, integral and intrinsic component of membrane and cytoplasm. The GO category “protein class” includes protease inhibitor, membrane traffic protein, oxidoreductase, enzyme modulator and transferase. These terms highlight the importance of ABA regulated genes in diverse processes.

To provide evidence that our analysis of ABA-induced differential gene expression is robust and novel for the identification of early induced genes we compared our data with two previously published data sets for P. patens (Komatsu et al., 2013; Stevenson et al., 2016). Stevenson et al. (2016) performed RNA sequencing analysis using the same ABA concentration as in our study (10 μM) and analyzed transcriptional changes after 60 min of ABA treatment whereas Komatsu et al. (2013) applied a lower ABA concentration of 1 μM and analyzed microarray gene expression after 180 min of ABA treatment. In our study we performed time series analysis of differentially expressed genes in response to 10 μM ABA for 30, 60, and 180 min. The comparison with both data sets indicates that there are 413 commonly up-regulated and 36 commonly down-regulated DEGs (Figure 3A). When we compare all our identified DEGs with both published data sets, we find 571 and 411 commonly regulated DEGs with Komatsu et al. and Stevenson et al., respectively (Figure 3B). Each study detects a set of unique DEGs, with our study detecting the highest number, i.e., 417 (Figure 3B). In addition, we have uniquely identified 16 early ABA inducible genes (Figure 3C and Table 3). Interestingly, seven of those genes encode proteins of unknown function, while five of the annotated genes might be involved in abiotic stress responses, i.e., Phypa_188559 is an F-box protein family whose homolog is heat-inducible in A. thaliana (Lim et al., 2006), Phypa_51481 and Phypa_38595 are encoding histidine kinase-related proteins, Phypa_167873 encodes a chloroplast targeted DnaJ chaperone and Phypa_90223 a CBS domain protein. These 16 ‘early’ genes were still up-regulated at 60 min of ABA treatment but at 180 min the expression of two genes (Phypa_117954 and Phypa_188559) was decreased, while the remaining 14 genes remained up-regulated (Table 3).

In order to assess the conservation of ABA-induced changes in gene expression we searched for orthologs of all ABA-regulated P. patens genes in the A. thaliana genome using the “G:Profiler orthology search” tool (Reimand et al., 2016) with maximum stringency only allowing a single best hit prediction per query gene. Out of a total of 1,030 differentially expressed P. patens genes, we found putative orthologs for 621 genes in the A. thaliana genome. Some P. patens DEGs were found to be paralogs since they were orthologous to the same A. thaliana gene leading to a final set of 549 unique A. thaliana orthologs of ABA-regulated P. patens genes (Supplementary Table S4).

To compare the ABA-dependent regulation of the identified A. thaliana orthologs with the ABA-regulated P. patens DEGs we made use of four large-scale studies that analyzed ABA-responsive genes in A. thaliana (Matsui et al., 2008; Wang et al., 2011; Liu et al., 2013; Zhan et al., 2015). These studies were selected because they have performed expression analysis with different ABA concentrations, durations of treatment and ages of the plants. In the first selected study 4-weeks-old A. thaliana seedlings were treated with 10 μM ABA for 6 h (Liu et al., 2013) and 406 up-regulated and 381 down-regulated genes were detected in response to ABA. In the second study, transcriptional changes were observed from 2-weeks-old A. thaliana plants that were subjected to 100 μM ABA treatment for 2 and 10 h (Matsui et al., 2008). In total, this study identified 3,137 up-regulated and 2,619 down-regulated genes. In the third study, 5-weeks-old plants were subjected to 50 μM ABA for 3 h (Wang et al., 2011) and 596 genes were identified as up-regulated and 441 genes were down-regulated. In the fourth study, 2 weeks old plants were subjected to 100 μM ABA for 6 h (Zhan et al., 2015) leading to the identification of 2,380 up-regulated and 2,359 down-regulated genes (Figure 4 and Table 4).

These studies revealed a total number of 8,743 non-redundant ABA-regulated A. thaliana genes and their transcriptional overlap inferred from the four studies is depicted in Table 4, Supplementary Table S5, and Figure 4. In total, 103 genes were found to be commonly regulated in all four studies. Liu et al. (2013) had identified the lowest number of genes regulated by ABA. Among the remaining three studies (Matsui et al., 2008; Wang et al., 2011; Zhan et al., 2015) 262 genes are commonly identified to be ABA regulated. In total, 2,255 genes are commonly regulated in at least two studies (Table 4).

Out of 549 unique A. thaliana orthologs of ABA-regulated P. patens genes, 238 genes were found to be ABA-regulated in at least one of the above mentioned A. thaliana studies (Figure 4, Table 4, and Supplementary Table S5). There are six genes (rubber elongation factor protein, protein phosphatase 2C family protein, NAD (P)-binding Rossmann-fold superfamily protein, NON-YELLOWING1, heme oxygenase-like and ABI1) that were found to be differentially expressed in response to ABA in all studies. Among our study and the four A. thaliana studies by Matsui et al. (2008), Wang et al. (2011), Liu et al. (2013), and Zhan et al. (2015) are 22, 180, 42, and 123 commonly regulated genes, respectively (Figure 4, Table 4, and Supplementary Table S5). This points to a considerable conservation of ABA-mediated regulation in P. patens and A. thaliana.

Gene ontology enrichment analysis of the 238 orthologs in which ABA-regulation was conserved in P. patens and A. thaliana indicates that in the “molecular function” category the majority of genes (110 genes) are involved in catalytic activity, 13 genes have the GO term coenzyme binding and 8 genes are involved in nucleic acid binding. In the “biological process” category the top enrichment GO terms are response to stimulus, response to stress, response to chemical and response to abiotic stimulus. In “cellular component” the highest fold enrichment was observed for GO terms chloroplast membrane, plastid membrane, chloroplast envelope, plastid envelope, vacuolar membrane, whole membrane and vacuole. These terms suggest evolutionary conservation of primary responses in different plastid processes as well as transcriptional regulation by ABA. The complete list of GO-enriched terms is presented in Supplementary Table S3.

In total 621 genes regulated by ABA in P. patens are non-conserved between A. thaliana and P. patens (Supplementary Table S8). Many of those might be lineage specific, since only 117 of these have GO terms associated with them in the most recent gene annotation version (Lang et al., 2018). There are 311 P. patens DEGs that have putative orthologs in A. thaliana but are not ABA regulated in A. thaliana, while there are 310 genes that are both conserved and regulated by ABA in both organisms. The ABA regulated P. patens genes that are conserved with A. thaliana do not show a dominant over-representation of GO terms (Supplementary Figure S1), regardless of whether the genes are also regulated by ABA in A. thaliana or not. However, many metabolic processes enabled by conserved genes are halted upon ABA action, as can be seen from many depleted terms that are associated with genes controlling macromolecule metabolism. The 409 ABA regulated P. patens genes not conserved with A. thaliana do not show many dominant depletion terms, but many enriched terms that reflect responses to abiotic stimuli/stresses, among them UV, cold and ROS (Supplementary Figure S1).

In order to detect DEGs under ABA stimulation over time, we measured the transcriptome response at 30, 60, and 180 min after ABA stimulation. To learn which genes control the developmental decision toward brachycytes, we performed a global transcriptome time course analysis as in Busch et al. (2013). A principal component analysis (PCA) on the transcriptomes shows a clear separation along the first principal component (PC 1) of the ABA-stimulated gene response vs. the control time points (Figure 5). Furthermore, there is a transient gene response at the 30 and 60 min time points of ABA treatment, which are separated along the second principal component (PC 2). Next, we analyzed which genes respond differently over time in the ABA-stimulated vs. the control cells by fitting a 3rd order polynomial to response of each gene under treatment and control separately and together and comparing the goodness of fit via analysis of variance (ANOVA) (Supplementary Table S6 and Figure 6). Using this analysis, we found 1,792 DEGs regulated over time at a FDR adjusted p-value threshold < 0.001 (as compared to 1,030 in the pairwise analyses). Among these there are 63 out of 1,380 annotated genes that encode transcription-associated proteins (TAPs). The analysis confirms the significant differential regulation of the five TAPs tested (Figure 2C), Phypa_141045, Phypa_172697, Phypa_163821, Phypa_72483, and Phypa_28324 (Supplementary Figure S2).

In total, 40 differentially expressed TAPs were found in all pairwise comparisons and 35 of these were also detected as differentially expressed according to the time series analysis (Supplementary Table S7). Only one TAP encoding a member of the AP2/EREBP TF family was down-regulated (Phypa_163821) whereas the other 39 TAPs were up-regulated in response to ABA. Most of the differentially expressed TAPs belong to the group of transcription factors (TF, 35) whereas only a few are transcriptional regulators (TR, 5). The latter are members of the Argonaute family mainly involved in small RNA-directed regulation of nuclear-encoded transcripts (Phypa_172642 and Phypa_141045), or they belong to the Sigma 70-like protein family that act as cofactors of the plastid-encoded RNA polymerase to control plastidic gene expression (Phypa_86427, Phypa_160930, and Phypa_230140). The five differentially expressed TAPs that were up-regulated already after 30 min stayed activated after 60 and 180 min. Half of the differentially expressed TAPs were up-regulated after 60 min and remain up-regulated even after 180 min (Figure 7A).

Out of 1,792 genes that were differentially regulated over time according to our time series analysis, 328 genes were classified as ABA-responsive according to the criteria defined by Timmerhaus et al. (2011). Out of the 40 differentially expressed TAPs from the pairwise comparisons an ABA responsive element was predicted for ten of the differentially expressed TAP genes (Timmerhaus et al., 2011) and an activation by ABA or salt was shown for 25 using a TAP-specific microarray (Richardt et al., 2010). The list of ABA-regulated TAPs found in our study thus fully incorporates both sets indicating the high sensitivity and reliability of this study (Figure 7B).

Since ABA plays an important role in the stress response of plants, we compared the set of differentially expressed TAPs identified in our analysis to other studies in P. patens that analyzed transcriptional changes in response to different abiotic stresses. Strikingly, we found that all TAPs that were previously identified to be differentially expressed during cold stress, dehydration and UV-B radiation are also differentially regulated in response to exogenously applied ABA (Figure 7C), supporting a common role of ABA in diverse molecular stress responses. We found 13 TAPs activated by ABA and cold stress (Beike et al., 2015), five activated by ABA and UV-B (Wolf et al., 2010) and nine activated by ABA and dehydration (Hiss et al., 2014). A single TAP gene (Phypa_126548) that is activated by ABA and all other stresses encodes an AP2/EREBP family protein previously shown to be prominently involved in the P. patens stress response (Richardt et al., 2010; Hiss et al., 2014). The cold-responsive TAPs that were also identified in our data set show an activation by cold primarily after 8–24 h. In our ABA time course most of them show an activation at the two later time points (60 and 180 min), suggesting that their activation upon cold requires ABA biosynthesis and/or ABA-release from glucosyl esters (Figure 7D).