Transcriptional profiling reveals conserved and species-specific plant defense responses during the interaction of the early divergent plant Physcomitrium patens with Botrytis cinerea

Bryophytes were among the first plants that colonized earth and they evolved key defense mechanisms to counteract microbial pathogens present in the new environment. Although great advances have been made on pathogen perception and subsequent defense activation in angiosperms, limited information is available in early divergent land plants. In this study, a transcriptomic approach uncovered the molecular mechanisms underlying the defense response of the bryophyte Physcomitrium patens against the important plant pathogen Botrytis cinerea. A total of 3.072 differentially expressed genes were significantly affected during B. cinerea infection, including genes encoding proteins with known function in angiosperm immunity and involved in pathogen perception, signaling, transcription, hormonal signaling, metabolic pathways such as shikimate and phenylpropanoid, and proteins with diverse role in defense against biotic stress. Similarly as in other plants, B. cinerea infection leads to downregulation of genes involved in photosynthesis and cell cycle progression. These results highlight the existence of evolutionary conserved defense responses to pathogens throughout the green plant lineage, suggesting that they were probably present in the common ancestors of land plants. Moreover, several genes acquired by horizontal transfer from prokaryotes and fungi, and a high number of P. patens-specific orphan genes were differentially expressed during B. cinerea infection, indicating that they are part of the moss immune response and probably played an ancestral role related to effective adaptation mechanisms to cope with pathogen invasion during the conquest of land. Key Message Evolutionary conserved defense mechanisms present in extant bryophytes and angiosperms, as well as moss-specific defenses are part of the immune response of the early divergent land plant Physcomitrium patens.


Introduction
Bryophytes, including mosses, liverworts, and hornworts, are small non-vascular plants that were among the first plants that colonized land approximately 450 million years ago (Mya).
Due to the transition from an aquatic to a terrestrial environment, bryophytes acquired adaptation mechanisms to cope with different kinds of abiotic stresses, including desiccation stress, variations in temperature, and UV-B radiation, as well as defense mechanisms against microorganisms present in the air and soil (Rensing et al. 2008;Ponce de León and Montesano 2017). Bryophytes diverged before vascular plants appeared and represent therefore excellent models to reveal ancient defense mechanisms against pathogens.
However, in contrast to the vast information available in angiosperms-pathogen interactions, only few evidences show the recognition of pathogens and subsequent activation of defense responses in early land plants such as mosses and liverworts.
Plants are in permanent contact with microbial pathogens, including fungi, bacteria, viruses and oomycetes. In order to detect the presence of the pathogen, angiosperms have evolved complex signaling and perception pathways. Microorganisms have carbohydrate-and protein-based signals that are essential for microbial survival, for example flagellin or chitin, classified as pathogen-associated molecular patterns (PAMPs) (Boller and Felix 2009).
Based on their conservation, and the fact that they are not produced in plant cells, plants synthesize different plasma membrane localized pattern recognition receptors (PRRs), including receptor-like kinases (RLKs), that recognize PAMPs to control plant immunity. In response to PAMPs, plants trigger a defense response called PAMP-triggered immunity (PTI) or basal resistance, which is the first level of defense that restricts pathogen infection in most plant species (Jones and Dangl 2006). Pathogens adapted to their host plants have evolved strategies to interfere and inhibit plant defense by the action of pathogen-secreted virulence factors known as effectors that target key PTI components (Boller and Felix 2009).
Plants have a second layer of immune system consisting of receptors encoded by resistance (R) genes to detect directly or indirectly the effector proteins leading to effector-triggered immunity (ETI) (Jones and Dangl 2006). Both PTI and ETI result in a rapid burst of extracellular reactive oxygen species (ROS), activation of mitogen-activated protein kinases (MPKs), increase in hormone synthesis, and induction of genes with different roles in plant defense (Jones and Dangl 2006). ETI is often accompanied by a hypersensitive response (HR), a type of programmed cell death that restricts the pathogen to the site of infection.
In bryophytes, most of the knowledge related to defense activation has been generated during the interaction of microbial pathogens with the moss Physcomitrella patens (Ponce de León and Montesano 2017), which has been recently renamed as Physcomitrium patens (P. patens) (Rensing et al. 2020). Pathogenic bacteria, fungi and oomycetes infect P. patens tissues leading to maceration, necrosis and cell death (Ponce de León 2011). P. patens senses the presence of fungal pathogens by perceiving chitin through the receptor CERK1 and consequently activates a MPKs signaling cascade, which is necessary for plant immunity (Bressendorff et al. 2015). Hormonal levels increases after pathogen colonization in moss tissues, including salicylic acid (SA), abscisic acid (ABA), auxin and the precursor of jasmonic acid (JA), cis-oxophytodienoic acid (OPDA), since JA is not synthesized in bryophytes (Ponce de León et al. 2012;Mittag et al. 2015). Expression levels of several genes 5 P. patens Gransden wild type colonies were cultivated on solid BCDAT medium (Ashton and Cove 1977) under standard long-day conditions (22°C, 16-h light/8-h dark regime under 60-80 μmol m 2 s -1 white light) for 3 weeks before spray inoculation with a B. cinerea 2x10 5 spores/mL suspension or water (mock). Three time points corresponding to spore germination (4 hours post inoculation; hpi), germ tubes elongation (8 hpi), and moss cells colonization (24 hpi) were analyzed. Three independent biological replicates of tissue were harvested at each time point and treatment for RNA extraction, immediately frozen in liquid nitrogen, and stored at -80˚C.

RNA extraction, cDNA library preparation and sequencing
Frozen samples were pulverized with a mortar and pestle and total RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Germany), including a RNase-Free DNase treatment in column (Qiagen, Germany), according to manufacturer's protocol. RNA quality control, library preparation, and sequencing were performed at Macrogen Inc. (Seoul, Korea). RNA integrity was checked before library preparation using an Agilent Technologies 2100  (Li et al., 2009), for further analysis.

Differential expression analysis
The reads were counted by the FeatureCounts software ver. 1.6.4 (Liao et al., 2013).
Additionally, for default options, the following parameters were set: Stranded (Reverse), Count fragments instead of reads -p, Allow read to contribute to multiple features True, Count multi-mapping reads/fragments -M and Reference sequence file P. patens v3.3.
Differential expression analyses were performed using EdgeR software ver. 3.24.1 (Robinson et al., 2010; with p-value adjusted threshold 0.05, p-value adjusted method Benjamini and Hochberg (1995)  The expression level of twelve selected genes was analyzed to validate RNA-Seq results via quantitative reverse transcription PCR (RT-qPCR). cDNA was generated from 1 μg of RNA using RevertAid Reverse transcriptase (Thermo Scientific) and oligo (dT) according to the manufacturer's protocol. RT-qPCR was performed in an Applied Biosystems QuantStudio 3 thermocycler using the QuantiNova Probe SYBR Green PCR Kit (Qiagen, Germany); mix proportions and cycling parameters were used as described in manufacturer's instructions.
Relative expression of each gene was normalized to the quantity of constitutively expressed Ubi2 gene (Le Bail et al., 2013), using the 2 -ΔΔCt method (Livak and Schmittgen 2001). Gene expression of P. patens inoculated tissues was expressed relative to the corresponding watertreated samples at the indicated time points, with its expression level set to one. Each data point is the mean value of three biological replicates. Student's t-test was performed to determine the significance for quantitative gene expression analysis using GraphPad Prism software ver. 8.0.2. P-values <0.01 were considered statistically significant. Primer pairs used for qPCR analyses are provided in SupplementaryTable S1, in all cases amplification efficiencies was greater than 95%.
Additionally, to default options, for every step, the following parameters were set for Blast:  (Wilson el al., 2005).

Global gene expression analysis of P. patens plants during B. cinerea infection
In previous studies, we have shown that B. cinerea infects P. patens tissues and that fungal biomass starts to increase after 8 hours post-inoculation (hpi), reaching high levels at 24 hpi (Ponce de León et al. 2012). To have more insights into P. patens defense response to B.
cinerea infection, RNAseq profiling was performed at three different stages of the colonization event; spore germination (4 hpi), germ tubes elongation (8 hpi), and hyphae proliferation and tissue colonization (24 hpi) ( Fig. 1). At 8 hpi, some plant cells start to respond to B. cinerea, evidenced by cell wall modifications detected with solophenyl flavine, and initial penetration attempts of some cells by infectious hyphae were visible ( Fig. 1d-e).
At 24 hpi, hyphae colonized P. patens cells and fungal structures could be clearly distinguished within P. patens cells (Fig. 1f). In our working conditions, disease symptoms were visible at 24 hpi evidenced by darkening of the tissues and heavy maceration at 48 hpi ( Fig. 1g-j), which are typical symptoms caused by this pathogen in angiosperms.
The transcriptomes of water-treated and B. cinerea-inoculated P. patens tissues were analyzed in three biological replicates. A total of 86.8-96.5% of the reads in the libraries mapped successfully to the genomes of P. patens (nuclear, chloroplast and mitochondria).
Reads mapped uniquely to the P. patens nuclear genome (approximately 208 million reads) were used for further analyses (Supplementary and LR, and included metabolic processes such as secondary metabolic process, phenylpropanoid, erythrose 4-phosphate/phosphoenolpyruvate family amino acid, Lphenylalanine and aromatic amino acid family metabolic process, and response to fungus, defense response and response to chemical (Fig. 3). Other upregulated defense-related GO terms present at MR and LR included response to abiotic stimulus, external stimulus, wounding, UV-B, and salt, as well as regulation of hormone levels and jasmonic acid biosynthetic process. Moreover, at LR, cell wall metabolic processes were enriched, including cell wall organization or biogenesis and regulation of cell wall macromolecule.
Collectively, biotic, abiotic related processes and phenylpropanoids metabolism, accounted for more than 70% and 50% of the GO biological process (BP) terms identified in the upregulated genes at MR and LR, respectively.
Significant enriched GO terms for downregulated genes were only identified at LR and were related to photosynthesis, including photosynthesis light reaction, generation of precursor metabolites and energy, protein-chromophore linkage, response to light stimulus, carbon fixation, photorespiration and chlorophyll metabolic process. Overrepresented downregulated DEGs included genes encoding chlorophyll binding proteins, proteins related to photosystem I and II, as well as many other proteins related to photosynthesis. Other significantly enriched GO terms were cell cycle process, cell cycle, nuclear division, cell division, chromosome segregation, and anaphase-promoting complex-dependent catabolic process. Within these overrepresented downregulated DEGs, genes encoded cyclins, cyclindependent kinases, centromeric proteins, cell division control proteins, kinesins, as well as other proteins involved in cell cycle progression and cell division (Supplementary Table   S5, S6). Taken together, these results indicate that during MR, P. patens activates a general respond that encompasses several defense pathways that are maintained during LR, while at LR important processes such as photosynthesis, cell division and cell cycle progression are downregulated.   ROS accumulation can lead to an HR response, which is a type of programmed cell death (PCD) (Lamb and Dixon 1997 cinerea colonization (Fig. 1) . 6; Supplementary Table S12). Most of these B. cinerearesponsive species-specific genes (92%) are described as with unknown functions in the Phytozome database. Of the 599 orphan genes, 376 showed no blast hits, meaning that they do not have any homologs and were de novo created genes. In addition, 137 orphan genes belong to multimember gene families, and 12 were not expressed in the P. patens Atlas project (Perroud et al. 2018) (Fig. 6). Two B. cinerea inducible orphan genes were secreted from moss tissues after chitosan treatment (Lehtonen et al. 2014), including Pp3c19_19780, Pp3c18_6260, suggesting that they play a role in moss defense. Taken together, the acquisition of horizontal transferred gene as well as orphan genes could have provided adaptive advantages to encounter microbial colonization in basal land plants. cinerea, including AP2/ERF, WRKY, MYB, bHLH, PLATZ and NAC (Mengiste et al. 2003;Berrocal-Lobo et al. 2002;Zhao et al. 2012;AbuQamar et al. 2017), were also differentially expressed in fungal-infected moss tissues. Compared with algae, the numbers of members of TFs families including MYB, WRKY, and AP2/ERF have expanded significantly in land plants and might be associated to terrestrial adaptation (Rensing et al. 2008;Rinerson et al. 2015;Pu et al. 2020). The high number of AP2/ERF TFs, as well as WRKYs and associated regulatory proteins with VQ motif, upregulated during B. cinerea infection, suggest their involvement in defenses to counteract pathogens among these adaptation mechanisms.

Discussion
Previous work also showed induced expression of AP2/ERF family members by other biotic stresses such as elicitors of P. carotovorum, chitosan and Phytophthora capsici (Alvarez et al. 2016, Bressendorff et al. 2016Overdijk et al. 2016). This result is consistent with a central regulatory role played by the AP2/ERF family during different stress conditions in P. patens (Hiss et al. 2014). Interestingly, three B. cinerea-inducible NAC TFs (Pp3c13_20650, Pp3c3_12890 and Pp3c4_20430), that are involved in P. patens water-conducting cell development, showed cell wall thickening and PCD when overexpressed in P. patens (Xu et al. 2014). Collectively, these findings indicate that defense signaling pathways were already  (Peng et al. 2017). The SA pathway was probably present in the common ancestor of land plants since all genes of this pathway are present in bryophytes but not in algae (Wang et al. 2015). However, the recent finding that P. patens NPR1-like receptor did not complement Arabidopsis npr3-4 double mutant suggests that fine-tuning of SA-induced defense responses probably evolved later in angiosperms (Peng et al. 2017 are not present in P. patens and moss cells do not respond to BR (Prigge et al. 2010;Wang et al. 2015). In addition, P. patens lacks a GA biosynthetic pathway downstream of entkaurenoic acid (KA), indicating that KA metabolites instead of GA may play physiological roles (Miyazaki et al. 2018). The shikimate and phenylpropanoid pathways played a crucial role in adaptation mechanisms that have evolved in early land plants to cope with abiotic stresses such as UV radiation and desiccation, as well as microbial attack (Emiliani et al. 2009). Our data highlight the importance of these pathways in moss defense against pathogens since a very high number of biosynthetic and regulatory genes were upregulated during B. cinerea infection. Similar results were obtained during angiosperms-B. cinerea interactions (Windram et al. 2012;Haile et al. 2019Haile et al. , 2020. The contribution of phenylpropanoids to bryophyte defense against pathogens was evidenced in several other P. patens-pathogen interactions (Oliver et al. 2009;Reboledo et al. 2015: Alvarez et al. 2016, and during infection of the liverwort M. polymorpha with the oomycete P. palmivora (Carella et al. 2019). P. patens genome duplicated 45 Mya and several gene families expanded, including PALs and CHSs (Koduri et al. 2010;Wolf et al. 2010). The presence of different CHSs activities in P. patens , Cinnamic acid levels increases after P. carotovorum elicitor treatment (Alvarez et al. 2016), and caffeic acid, coumaric acid, 4-hydroxybenzoic acid and caffeoyl quinic acid have been detected in P. patens tissues, some of which have antimicrobial properties (Erxleben et al. 2012;Richter et al. 2012). Several B. cinerea-responsive DEGs were annotated as CHI, F3H and FLS. However, unlike angiosperms, P. patens does not seem to contain enzymes required for the synthesis of flavonones and flavones. Both CHIs are enhancers of flavonoid production genes encoding type IV CHI proteins with no CHI activity (Ngaki et al. 2012;Waki et al. 2020), and there are no evident orthologs for F3H and FLS (Koduri et al. 2010;Kawai et al. 2014). Thus, P. patens probably lacks the later steps in flavonoid biosynthesis, in contrast to M. polymorpha, which accumulates different types of flavonoids (Albert et al. 2018). P. patens has a pre-lignin pathway instead of true lignin (Weng and Chapple, 2010), and we showed previously that polyphenolic compounds reinforce the cell wall upon B.  Renault et al. (2017) showed that the pre-lignin pathway of P. patens involves the formation of caffeic acid coupled to the hexose-derived ascorbate pathway, rather than the shikimate pathway, as is reported for angiosperms. This leads to the formation of soluble caffeoyl-threonic acids, which are essential for the synthesis of the moss cuticle.
In addition, transcript levels of several B. cinerea-responsive MYB, belonging to R2R3 MYB (Pu et al., 2020), and bHLB TFs are involved in phenylpropanoid biosynthesis in angiosperms (Xu et al. 2013 (Remy et al., 1994), and gene acquisition following fungal interaction might have represented adaptive benefits for resistance to pathogen colonization.
In addition, P. patens has 13 % orphan genes, which are species-specific genes with no orthologs in other plants (Zimmer et al. 2013). Our findings show that a large proportion of B. cinerea-responsive P. patens DEGs (599 genes: 19%) correspond to orphan genes, from which 376 were de novo created genes, suggesting that they could represent innovative adaptive strategies to biotic stress. The involvement of orphan genes in cold acclimation (Beike et al. 2015), indicate that they also participate in abiotic stress tolerance. Moreover, twelve B. cinerea-responsive P. patens orphan genes were not expressed in any condition present in the large-scale RNA-seq P. patens Altlas project (Perroud et al. 2018), suggesting that they could represent pathogen-specific orphans. Deciphering the roles played by orphan genes, most of which have no annotations, during biotic and abiotic stress in P. patens, will help to understand how they contributed to stress adaptation of early-diverging land plants.
In conclusion, our results reveal that conserved defense mechanisms between extant bryophytes and angiosperms, as well as moss-specific defenses are part of the P. patens immune system, which could have been pivotal for land colonization by plants. Further studies on P. patens-pathogen interactions will contribute to uncover the molecular 26 mechanisms underlying moss-specific defenses and their involvement during coevolution of ancient land plants and pathogens.

Availability of supporting data
The sequencing raw data from the RNA-Seq libraries were deposited on the Sequence Read