Fay-Wei Li

Paul Brouwer

Lorenzo Carretero-Paulet

Shifeng Cheng

Jan De Vries

Pierre-Marc Delaux

Ariana Eily

Nils Koppers

Li-Yaung Kuo

Zheng Li

Mathew Simenc

Ian Small

Eric Wafula

Stephany Angarita

Michael S Barker

Andrea Bräutigam

Claude Depamphilis

Sven Gould

Prashant S Hosmani

Yao-Moan Huang

Bruno Huettel

Yoichiro Kato

Xin Liu

Steven Maere

Rose Mcdowell

Lukas A Mueller

Klaas G J Nierop

Stefan A Rensing

Tanner Robison

Carl J Rothfels

Erin M Sigel

Yue Song

Prakash R Timilsena

Yves Van De Peer

Hongli Wang

Per K. I Wilhelmsson

Paul G Wolf

Xun Xu

Joshua P Der

Henriette Schluepmann

-S Wong

Kathleen M Pryer

Fern genomes elucidate land plant evolution and cyanobacterial symbioses

Ferns are the closest sister group to all seed plants, yet little is known about their genomes other than that they are generally colossal. Here, we report on the genomes of Azolla filiculoides and Salvinia cucullata (Salviniales) and present evidence for episodic whole-genome duplication in ferns-one at the base of 'core leptosporangiates' and one specific to Azolla. One fernspecific gene that we identified, recently shown to confer high insect resistance, seems to have been derived from bacteria through horizontal gene transfer. Azolla coexists in a unique symbiosis with N 2 -fixing cyanobacteria, and we demonstrate a clear pattern of cospeciation between the two partners. Furthermore, the Azolla genome lacks genes that are common to arbuscular mycorrhizal and root nodule symbioses, and we identify several putative transporter genes specific to Azolla-cyanobacterial symbiosis. These genomic resources will help in exploring the biotechnological potential of Azolla and address fundamental questions in the evolution of plant life.

T he advent of land plants ~474-515 Myr ago 1 led to complex vegetational innovations that shaped emerging terrestrial and freshwater ecosystems. Bryophytes, lycophytes, ferns and gymnosperms dominated the global landscape before the ecological radiation of flowering plants 90 Myr ago [START_REF] Schneider | Ferns diversified in the shadow of angiosperms[END_REF] . The first complete plant genome sequence (Arabidopsis thaliana) was published in 2000 [START_REF]The Arabidopsis Genome Initiative Analysis of the genome sequence of the flowering plant Arabidopsis thaliana[END_REF] , followed by reference genomes for all other major lineages of green plants, except ferns. A dearth of genomic information from this entire lineage has limited not only our knowledge of fern biology but also the processes that govern the evolution of land plants.

The relatively small genome (0.75 Gb) [START_REF] Obermayer | Nuclear DNA C-values in 30 species double the familial representation in pteridophytes[END_REF] of Azolla is exceptional among ferns, a group that is notorious for genomes as large as

Articles

NaTure PlaNTS Asia [START_REF] Lumpkin | Azolla: botany, physiology, and use as a green manure[END_REF] . The Azolla symbiosis is unique among plant-bacterial endosymbiotic associations because the cyanobiont is associated with the fern throughout its life cycle, being vertically transmitted during sexual reproduction to subsequent generations [START_REF] Zheng | Cellular responses in the cyanobacterial symbiont during its vertical transfer between plant generations in the Azolla microphylla symbiosis[END_REF] . In all other land plant symbiotic associations, the relationship must be renewed each generation. The Nostoc symbiont is not capable of autonomous growth and its genome shows clear signs of reduction, with several housekeeping genes lost or pseudogenized [START_REF] Ran | Genome erosion in a nitrogen-fixing vertically transmitted endosymbiotic multicellular cyanobacterium[END_REF] . With a fossil record that extends back to the mid-Cretaceous period, Azolla probably shares a ~100-Myr-old co-evolutionary relationship with Nostoc [START_REF] Hall | Studies on fossil Azolla: Azolla montana, a Cretaceous megaspore with many small floats[END_REF] .

To better understand genome size evolution in Azolla and its closely related lineages, we obtained genome size estimates for all five genera of Salviniales (Supplementary Table 1). We found them to be at least an order of magnitude smaller than any other fern species (Fig. 1a), and, most notably, the genome of Salvinia cucullata, which belongs to the sister genus to Azolla, is only 0.26 Gb, the smallest genome size ever reported in ferns. This unanticipated discovery afforded us the opportunity to include a second fern genome for comparison.

Results

Genome assembly and annotation. To gain insight into fern genome evolution, as well as plant-cyanobacterial symbioses, we sequenced the genomes of A. filiculoides (Fig. 1b) and S. cucullata (Fig. 1c) using Illumina and PacBio technologies. The assembled Azolla and Salvinia genomes have N50 contig size of 964.7 Kb and 719.8 Kb, respectively. The BUSCO (Benchmarking Universal Single-Copy Orthologs) assessment and Illumina read-mapping results indicate high completeness for both assemblies (Supplementary Fig. 1 and Supplementary Table 2). We identified 20,201 and 19,914 high-confidence gene models in Azolla and Salvinia, respectively, that are supported by transcript evidence or had significant similarity to other known plant proteins (Supplementary Figs. 123, Supplementary Table 3 and Supplementary Discussion). Salvinia genes are much more compact, with a mean intron length half of that in Azolla (Supplementary Fig. 1). In addition to introns, differences in the repetitive content explain some of the nearly threefold difference in genome size. Azolla has more of every major category of repeat, but 191 Mb of the 233-Mb difference in the total

Articles

NaTure PlaNTS repetitive content are made of retroelements, especially Gypsy and Copia long terminal repeat retrotransposons (LTR-RTs; Fig. 1d and Supplementary Fig. 4). DNA transposon profiles are similar for the two ferns except that Azolla has substantially more SOLA elements than does Salvinia (Fig. 1e).

Insights into gene family evolution in land plants. The genomes of Azolla and Salvinia offer a new opportunity to examine the evolution of plant genes and gene families across all Viridiplantae (land plants plus green algae). We classified genes into orthogroups from 23 genomes (12 angiosperms, 2 gymnosperms, 2 ferns, 1 lycophyte, 2 mosses, 2 liverworts, 1 charophyte and 1 chlorophyte; Supplementary Table 5) and reconstructed the gene family evolution-gain, loss, expansion and contraction-across the green tree of life (Supplementary Fig. 5 and Supplementary Table 5). To investigate the origin of genes linked to seed development, we examined orthogroups containing 48 transcription factors that express exclusively in Arabidopsis seeds [START_REF] Le | Global analysis of gene activity during Arabidopsis seed development and identification of seed-specific transcription factors[END_REF] . Homologues of 39 of them were detected in ferns or other seed-free plants, indicating that many seed transcription factors were present before the origin of seeds (Supplementary Table 6). Similarly, only a handful of transcription factor families arose along the branch that led to seed plants (Supplementary Table 7); rather than relying on entirely novel transcription factors, it seems instead that an expansion of pre-existing transcription factor families had a greater role in seed plant evolution [START_REF] Wilhelmsson | Comprehensive genome-wide classification reveals that many plant-specific transcription factors evolved in streptophyte algae[END_REF] . Indeed, ancestral gene number reconstructions of MADSintervening keratin-like and C-terminal (MIKC)-type MADS box genes (orthogroup 23) show that these important developmental regulators more than doubled in number from 15 in the ancestral vascular plant to 31 in the ancestral euphyllophyte (here, Salviniales plus seed plants; Supplementary Table 5).

In a recent study on the evolution of plant transcription-associated proteins, which include transcription factors and transcriptional regulators [START_REF] Wilhelmsson | Comprehensive genome-wide classification reveals that many plant-specific transcription factors evolved in streptophyte algae[END_REF] , ferns were exclusively represented by the Pteridium aquilinum transcriptome. The finding that the transcriptional regulator Polycomb group EZ (PcG_EZ) was lost in ferns is corroborated here by our whole-genome data (Supplementary Table 8). Conversely, the transcription factor ULTRAPETALA, which originated at the base of euphyllophytes and is present in P. aquilinum, was apparently secondarily lost in Salviniales (Supplementary Table 8). YABBY, an important transcription factor that patterns leaf polarity in flowering plants, is absent in our fern genomes and in the genome of the lycophyte Selaginella moellendorffii [START_REF] Banks | The Selaginella genome identifies genetic changes associated with the evolution of vascular plants[END_REF] (Supplementary Table 8). Interestingly, a YABBY homologue was recently identified in a separate lycophyte species-Huperzia selago [START_REF] Evkaikina | The Huperzia selago shoot tip transcriptome sheds new light on the evolution of leaves[END_REF] -suggesting that YABBY has been lost at least twice in land plant evolution (in Selaginella and in ferns). How the differential retention of YABBY shaped the evolution of the vascular plant body plan requires further studies.

Among the orthogroups specific to seed plants, 1-aminocyclopropane-1-carboxylic acid (ACC) oxidase is of special interest because it converts ACC to ethylene-the last step in the ethylene biosynthetic pathway (Fig. 2). Ethylene is a critical plant hormone that controls various important physiological responses (for example, fruit ripening, flowering time, seed germination and internode elongation). Because ethylene responses are known in bryophytes, lycophytes and ferns [START_REF] Van De Poel | Ethylene in Plants[END_REF] , it is intriguing to find that ACC oxidase only evolved within seed plants, a result confirming that seed-free plants must possess an alternative ethylene-forming mechanism [START_REF] Osborne | Evidence for a non-ACC ethylene biosynthesis pathway in lower plants[END_REF] . Two other mechanisms, found in bacteria and fungi, result in ethylene formation: one via the 2-oxoglutarate-dependent ethylene-forming enzyme and the other through the non-enzymatic conversion of 2-keto-4-methylthiobutyric acid (KMBA) into ethylene [START_REF] Van De Poel | Ethylene in Plants[END_REF] . We did not identify ethylene-forming enzyme in seed-free plant genomes, suggesting the absence of the ethylene-forming enzyme-based biosynthetic pathway. Seed-free plants may possibly synthesize ethylene non-enzymatically via KMBA; however, further biochemical studies are needed to test this hypothesis. Interestingly, ACC synthase (upstream of ACC oxidase) is present in seed-free plants, albeit in a lower copy number (< 3) compared to seed plants, which average 9.3 copies (Fig. 2 and Supplementary Fig. 6). Paralogues of ACC synthase in seed plants are differentially regulated in response to varying developmental or environmental stimuli [START_REF] Tsuchisaka | A combinatorial interplay among the 1-aminocyclopropane-1-carboxylate isoforms regulates ethylene biosynthesis in Arabidopsis thaliana[END_REF] . Thus, it is plausible that the expansion of the ACC synthase family was coupled with the origin of ACC oxidase in seed plants to create a regulated ethylene biosynthetic pathway.

The history of whole-genome duplication in ferns. Our Multi-tAxon Paleopolyploidy Search (MAPS) [START_REF] Li | Early genome duplications in conifers and other seed plants[END_REF] phylogenomic analyses of the Azolla and Salvinia genomes (Fig. 3a), together with all available transcriptome data from other ferns, support two whole-genome duplication (WGD) events: a recent WGD event occurring in Azolla following its divergence from Salvinia and an earlier WGD predating the origin of 'core leptosporangiates' (sensu Pryer at al. [START_REF] Pryer | Phylogeny and evolution of ferns (monilophytes) with a focus on the early leptosporangiate divergences[END_REF] ), a large clade comprising the heterosporous, tree and polypod ferns. The observed peaks of duplication associated with the inferred WGDs exceeded the 95% confidence intervals of our birth and death simulations for gene family evolution in the absence of WGDs. This high number of shared gene duplications is readily explained by a significant episodic birth event, such as a WGD. The discovery that Azolla experienced a genome duplication independent of other heterosporous ferns is not entirely surprising because Azolla has nearly twice the number of chromosomes of other heterosporous ferns, including Salvinia and Pilularia [START_REF] Rice | The Chromosome Counts Database (CCDB)-a community resource of plant chromosome numbers[END_REF][START_REF] Wood | The frequency of polyploid speciation in vascular plants[END_REF] (Fig. 1a).

To further substantiate the two WGD events identified by MAPS, we examined the distribution of synonymous distances (K s ) between syntenic paralogues within each of the genomes, as well as syntenic orthologues between Azolla and Salvinia. In the Azolla genome, we detected 242 syntenic blocks comprising 988 syntelog pairs. By contrast, only 83 syntenic blocks with 254 syntelog pairs could be found in Salvinia. Between Azolla and Salvinia, 3,587 pairs of syntenic orthologues were detected, clustering into 356 syntenic genomic blocks. We fit Gaussian mixture models to identify peaks in the K s distributions (Fig. 3b and Supplementary Fig. 7). The main peak for Azolla-Salvinia orthologue pairs centres at ~1.0, which marks the species divergence between the two genera. To the left of this peak is the major Azolla intragenomic peak (~0.8), whose position confirms the Azolla-specific WGD event (Fig. 3b). To the right of the Azolla-Salvinia divergence peak is the Salvinia intragenomic K s peak (~1.2-1.3), which matches a minor Azolla intragenomic peak, consistent with the proposed pre-core leptosporangiates WGD (Fig. 3b). Moreover, despite the antiquity of the WGDs and species divergence (Fig. 1a), we were still able to detect Azolla-Salvinia syntenic regions in a 2:1 or 2:2 syntenic relationship (Fig. 3c), respectively, corroborating the Azolla-specific and the older WGD events. The confirmation of these two WGDs in ferns further allows us to characterize patterns of gene retention following WGD. We found that Azolla syntenic paralogues are enriched for transcription-related genes (Supplementary Table 9), similar to what was observed in Arabidopsis and other angiosperms [START_REF] Maere | Modeling gene and genome duplications in eukaryotes[END_REF] . Likewise, protein kinases, another functional category commonly retained after WGD in seed plants, are significantly enriched in Salvinia syntenic paralogues (Supplementary Table 9). Additional genomic data are needed to better characterize the distribution of WGD events across the fern tree of life and to compare patterns of post-WGD gene fractionation with those documented in seed plants.

The pentatricopeptide repeat family and RNA editing. The pentatricopeptide repeat (PPR) family is the largest gene family found in the Azolla and Salvinia genomes, with the Azolla genome encoding over 2,000 PPR proteins and the Salvinia genome over 1,700 PPR proteins. PPRs are implicated in organellar RNA processing [START_REF] Barkan | Pentatricopeptide repeat proteins in plants[END_REF] , and the large repertoire of PPRs correlates well with the extensive RNA editing we observed in the organellar genomes of Salviniales: 1,710 sites in Azolla organelles and 1,221 sites in Salvinia (Supplementary Table 10). These editing events include both C-to-U conversions (~70%) and U-to-C conversions (~30%). The number of PPR genes and the degree of RNA editing greatly exceed that found in seed plants and most bryophytes [START_REF] Takenaka | RNA editing in plants and its evolution[END_REF] . Of the sequenced plant genomes, only that of S. moellendorffii [START_REF] Banks | The Selaginella genome identifies genetic changes associated with the evolution of vascular plants[END_REF] has more PPR genes [START_REF] Cheng | Redefining the structural motifs that determine RNA binding and RNA editing by pentatricopeptide repeat proteins in land plants[END_REF] , correlating with the hyperediting seen in lycophytes [START_REF] Oldenkott | Chloroplast RNA editing going extreme: more than 3400 events of C-to-U editing in the chloroplast transcriptome of the lycophyte Selaginella uncinata[END_REF] . However, there are no U-to-C editing events in Selaginella, making the Azolla and Salvinia genome sequences a novel and valuable resource for identifying the unknown factors catalysing these events.

More than half of the plastid transcripts and two-thirds of the mitochondrial transcripts in Azolla and Salvinia require start codon creation by C-to-U editing or stop codon removal by U-to-C editing before translation is possible. Most stop codon edits (76%) and start codon edits (62%) are shared between Azolla and Salvinia plastomes (as opposed to only 19% in internal ACG codons; Supplementary Fig. 10). This persistence of start and stop codon edits suggests that their loss is selected against, that is, creating the translatable sequence by RNA editing has an advantage over having it encoded by the genome. This argues that these particular RNAediting events are not selectively neutral [START_REF] Gray | Evolutionary origin of RNA editing[END_REF] and supports editing as a control mechanism for gene expression in fern organelles.

Only ~55-60% of PPR proteins (1,220 in Azolla and 930 in Salvinia) contain domains associated with RNA editing in other plants. Although sufficient to account for the number of editing events observed (assuming each protein can specify one or a few sites as in other plants), this leaves a very large number of PPR proteins (~700 in Azolla and ~600 in Salvinia) with unknown functions. By comparison, flowering plants contain only 200-250 PPR proteins that lack editing domains.

Origin and evolution of a fern insecticidal protein. Ferns are remarkable for their high levels of insect resistance compared to flowering plants [START_REF] Hendrix | An evolutionary and ecological perspective of the insect fauna of ferns[END_REF] . Recently, Shukla et al. [START_REF] Shukla | Expression of an insecticidal fern protein in cotton protects against whitefly[END_REF] isolated a novel insecticidal protein, Tma12, from the fern Tectaria macrodonta. Transgenic cottons carrying Tma12 exhibit outstanding resistance to whitefly, yet show no decrease in yields, demonstrating tremendous agricultural potential. Tma12 has a high similarity to chitinbinding proteins (Pfam PF03067), but its evolutionary origin is unknown. Here, we found a Tma12 homologue to be present in the Salvinia genome (henceforth ScTma12), as well as in a few 1,000 Plants (1KP) [START_REF] Matasci | Data access for the 1,000 Plants (1KP) project[END_REF] fern transcriptomes, but not in Azolla or any other publicly available plant genomes. Phylogenetic analyses position the fern Tma12 sequences together with bacterial sequences, and are most closely related to the chitin-binding proteins from Chloroflexi (Fig. 4). We investigated whether this insecticidal protein was more likely a result of horizontal gene transfer (HGT) from bacteria to ferns or produced by fern-associated microorganisms. ScTma12 is in a 646,687-bp scaffold (Sacu_v1.1_s0099) and has an 247-bp intron. The genes upstream and downstream of ScTma12 are all clearly plant genes, and we found no abnormality in read-mapping quality, nor an abrupt change in read coverage (Supplementary Fig. 9), which together speak against the sequence being a contamination from a bacterial source. It has been argued that differential loss of genes in eukaryotes is the rule and gene acquisition by HGT rather rare [START_REF] Ku | Endosymbiotic origin and differential loss of eukaryotic genes[END_REF] . The concerted loss of Tma12 in each of the other Viridiplantae lineages is unlikely but cannot entirely be ruled out. However, functional HGT into eukaryotes does occur [START_REF] Li | Horizontal transfer of an adaptive chimeric photoreceptor from bryophytes to ferns[END_REF][START_REF] Husnik | Functional horizontal gene transfer from bacteria to eukaryotes[END_REF] and ScTma12 might represent such a case that contributed to the welldocumented resistance of ferns against phytophagous insects.

Azolla-cyanobacterial symbiosis.

To explore the co-evolutionary history of the Azolla-Nostoc symbiosis, we resequenced five other Azolla species and assembled each of their cyanobiont genomes. We then compared the cyanobiont phylogeny to the host species phylogeny and found a clear cospeciation pattern, with just one exception (the placement of Azolla caroliniana; Fig. 5a). Although such a pattern has been hinted at before [START_REF] Van Coppenolle | Genetic diversity and phylogeny analysis of Anabaena azollae based on RFLPs detected in Azolla-Anabaena azollae DNA complexes using nif gene probes[END_REF][START_REF] Zheng | Genetic diversity and classification of cyanobacteria in different Azolla species by the use of PCR fingerprinting[END_REF] , we provide unequivocal The ethylene-forming pathway involves the Yang cycle, where ACC is synthesized from S-adenosyl-methionine (SAM; also known as AdoMet) by ACC synthase. ACC oxidase then catalyses the conversion of ACC to ethylene. We found that ACC oxidase is unique to seed plants (green) and its origin probably drove the expansion of the ACC synthase gene family (orange; Supplementary Fig. 6) to create a regulated ethylene biosynthetic mechanism.

evidence from whole-genome data. The genetic basis for this persistent symbiosis is undetermined. In plants, two other mutualistic associations-the arbuscular mycorrhizal (AM) and the nitrogenfixing root nodule (RN) symbioses-have been well characterized.

Whereas the AM symbiosis is formed between almost all land plants and a single fungal clade (Glomeromycota) [START_REF] Parniske | Arbuscular mycorrhiza: the mother of plant root endosymbioses[END_REF] , the RN symbiosis is restricted to a few angiosperm lineages (mostly legumes) that associate with various nitrogen-fixing bacterial symbionts (for example, Rhizobium and Frankia). Despite these distinct differences, both symbioses require that a common symbiosis pathway (CSP) be established [START_REF] Parniske | Arbuscular mycorrhiza: the mother of plant root endosymbioses[END_REF] . This pathway is highly conserved in all land plants [START_REF] Delaux | Algal ancestor of land plants was preadapted for symbiosis[END_REF] , except for those that have lost the AM symbiosis [START_REF] Delaux | Comparative phylogenomics uncovers the impact of symbiotic associations on host genome evolution[END_REF][START_REF] Bravo | Genes conserved for arbuscular mycorrhizal symbiosis identified through phylogenomics[END_REF] , such as A. thaliana and three aquatic angiosperms [START_REF] Delaux | Comparative phylogenomics uncovers the impact of symbiotic associations on host genome evolution[END_REF][START_REF] Bravo | Genes conserved for arbuscular mycorrhizal symbiosis identified through phylogenomics[END_REF] . We investigated whether the CSP might have been co-opted during the evolution of the Azolla-Nostoc symbiosis by searching for six essential CSP genes in the Azolla and Salvinia genomes, as well as in transcriptomic data from other ferns in the 1KP data set [START_REF] Matasci | Data access for the 1,000 Plants (1KP) project[END_REF] (Supplementary Table 11). Although DMI2 (also known as SYMRK), DMI3 (also known as CCaMK), IPD3 (also known as CYCLOPS) and VAPYRIN were found in other ferns, the Azolla and Salvinia genomes completely lacked orthologues (Fig. 5b). IPD3 and VAPYRIN do not belong to multigene families [START_REF] Delaux | Algal ancestor of land plants was preadapted for symbiosis[END_REF] and homologues were not detected. Although homologues of DMI2 and DMI3 were identified, phylogenetic analyses confirmed that they are not orthologous to the symbiotic genes (Supplementary Data). In addition, for DMI3, we searched the Azolla and Salvinia homologues for two motifs (threonine 271 and the calmodulin-binding domain) that are critical for symbiosis. Both motifs are missing from these sequences, confirming the absence of DMI3. CASTOR and POLLUX are paralogues resulting from a gene duplication event in the ancestor of seed plants, and although pre-duplicated homologues are present in Salvinia and other seed-free plants, they are absent in Azolla (Fig. 5b). The co-elimination of the CSP genes suggests the lack of AM symbiosis in Azolla and Salvinia and that the nitrogenfixing Azolla-Nostoc symbiosis does not rely on this pathway.

To identify genes important for the Azolla-Nostoc symbiosis, we treated A. filiculoides with erythromycin to remove the cyanobiont (AzCy-) and compared its gene expression patterns with the wild type (AzCy+ ). Experiments were carried out in conditions where the nitrogen nutrient (ammonium nitrate) was either supplied (N+ ) or withheld (N-) from the growth media. Results from nifH real-time PCR confirmed the complete absence of cyanobacteria in AzCy-and showed that the addition of the nitrogen nutrient suppresses symbiotic N 2 fixation in AzCy+ (Supplementary Fig. 10), consistent with an earlier study [START_REF] Brouwer | Metabolic adaptation, a specialized leaf organ structure and vascular responses to diurnal N 2 fixation by Nostoc azollae sustain the astonishing productivity of Azolla ferns without nitrogen fertilizer[END_REF] . A large portion of the transcriptome is affected by the presence or absence of cyanobionts, with 6,210 and 2,125 genes being differentially transcribed under N-and N+ conditions, respectively (Fig. 5c and Supplementary Discussion). Of these, over 33% have at least a twofold expression difference. In response to nitrogen starvation, the Azolla transcriptomes remained moderately stable when the cyanobiont was present, but shifted drastically once it was absent (Fig. 5d). This finding suggests that the presence of the cyanobiont buffers the transcriptomic profile of Azolla from fluctuations in environmental nitrogen availability. We focused primarily on those genes that are differentially expressed between the nitrogen treatments when the cyanobiont is present, and to a lesser extent on when the cyanobiont is absent (Fig. 5e and Supplementary Discussion). Because the cyanobacterial N 2 -fixation rate is strongly induced in the N-condition, we expect these genes to be candidates involved in nutrient exchange or in communication with the cyanobiont to promote N 2 fixation. A total of 88 upregulated and 72 downregulated genes were identified (Fig. 5e). Among the upregulated genes is a paralogue of the ammonium transporter 2 subfamily (AfAMT2-4; Azfi_s0034. g025227; Fig. 5e and Supplementary Fig. 11) that is probably dedicated to ammonium uptake from the Azolla leaf cavity where the cyanobiont resides; homologous ammonium transporters have been implicated to participate in the AM and RN symbioses [START_REF] Breuillin-Sessoms | Suppression of arbuscule degeneration in Medicago truncatula phosphate transporter4 mutants is dependent on the ammonium transporter 2 family protein AMT2;3[END_REF][START_REF] D' Apuzzo | Characterization of three functional high-affinity ammonium transporters in Lotus japonicus with differential transcriptional regulation and spatial expression[END_REF] . There is also a paralogue of the molybdate transporter gene family (AfMOT1; Azfi_s0167.g054529) that is most likely specialized for supplying molybdenum, a required co-factor for nitrogenase, to the cyanobiont. One of the legume MOT1 genes was recently found to facilitate nitrogenase activity in RN symbiosis [START_REF] Tejada-Jiménez | Medicago truncatula molybdate transporter type 1 (MtMOT1.3) is a plasma membrane molybdenum transporter required for nitrogenase activity in root nodules under molybdenum deficiency[END_REF] . In addition to these two transporters, we identified a chalcone synthase paralogue in this candidate gene set. Chalcone synthase catalyses the production of naringenin chalcone and is the first committed step in flavonoid biosynthesis. Interestingly, naringenin and naringin both have significant effects on promoting cyanobacterial growth [START_REF] Żyszka | Modulation of the growth and metabolic response of cyanobacteria by the multifaceted activity of naringenin[END_REF] and differentiation [START_REF] Cohen | Flavonoid-induced expression of a symbiosisrelated gene in the cyanobacterium Nostoc punctiforme[END_REF] . Naringin is also a hormogonium-repressing factor [START_REF] Cohen | Flavonoid-induced expression of a symbiosisrelated gene in the cyanobacterium Nostoc punctiforme[END_REF] . Because hormogonia lack heterocysts and cannot fix nitrogen, naringin (or related flavonoids) could act as a plant signal to boost N 2 fixation in the cyanobiont (Supplementary Discussion).

Although the ancient and intimate nature of the Azolla-Nostoc relationship suggests that gene transfer from Nostoc to the Azolla nuclear genome may have occurred over time, a thorough homology search found no evidence of Nostoc-to-Azolla HGT. However, we did discover a cyanobacteria-derived gene in the Azolla genome, but one that is shared with other ferns. This gene encodes a squalene-hopene cyclase (SHC), which mediates the cyclization of squalene into hopene, and is thought to be the evolutionary progenitor of many classes of eukaryotic and prokaryotic sterol cyclases. We found SHC homologues in both the Azolla and the Salvinia genomes, as well as in 40 fern 1KP transcriptomes. Our reconstructed gene phylogeny clearly shows that the fern SHCs are nested among cyanobacteria sequences (Supplementary Fig. 12). Although no homologue can be found in seed plants or in green algae, the SHC is also present in bryophyte (moss and liverwort) genomes and transcriptomes. Interestingly, these bryophyte SHCs are not related to those of ferns but are embedded in other bacterial SHC lineages (the monophyly of land plant SHCs is rejected by the Swofford-Olsen-Waddell-Hillis test [START_REF] Church | Automation and evaluation of the SOWH test with SOWHAT[END_REF] , P < 0.005). This finding implies a complex evolutionary history for SHCs in land plants, possibly featuring independent transfers of SHC from different prokaryotic lineages to mosses, liverworts and ferns. We are confident that these SHC genes are not from contaminants because the gene phylogeny largely mirrors the species phylogeny; furthermore, the SHC genes were not assembled into stray scaffolds in the genomes of Azolla, Salvinia, Physcomitrella [START_REF] Rensing | The Physcomitrella genome reveals evolutionary insights into the conquest of land by plants[END_REF] or Marchantia [START_REF] Bowman | Insights into land plant evolution garnered from the Marchantia polymorpha genome[END_REF] . In addition, we detected the triterpene products of SHC, hop-22(29)-ene, diplopterol and tetrahymanol, in S. cucullata biomass, providing direct evidence for SHC activity in this fern (Supplementary Fig. 13). Similar observations of SHC-synthesized triterpenes have been made in polypod ferns [START_REF] Ageta | A fern constituent, fernene, a triterpenoid hydrocarbon of a new type[END_REF][START_REF] Shinozaki | Squalene cyclase and oxidosqualene cyclase from a fern[END_REF] and mosses [START_REF] Marsili | Triterpenes from mosses[END_REF] . Because hopenes have an important role in plasma membrane stability in prokaryotes (similar to steroids in eukaryotes) and have been shown to confer low-temperature adaptation and stress tolerance [START_REF] Sáenz | Hopanoids as functional analogues of cholesterol in bacterial membranes[END_REF] , it is plausible that the convergent evolution of hopene biosynthesis in seed-free plants, through independent HGTs from bacteria, might have contributed to the early adaptations of land plants to diverse and adverse environments. Functional studies are needed to confirm this hypothesis.

We anticipate that the availability of the first genomic data from ferns will continue to lead to vital insights into the processes that govern the evolution of plant genes and gene families. The implementation of fern data into the existing comparative genomic framework will enhance our understanding the plant tree of life.

Methods

Flow cytometry and genome size estimation. To estimate the genome sizes of S. cucullata, P. americana, Regnellidium diphyllum and Marsilea minuta (Supplementary Table 1), we used the Beckman chopping buffer to extract nuclei from fresh leaves, following the protocol of Kuo and Huang [START_REF] Kuo | Determining genome size from spores of seedless vascular plants[END_REF] . The nuclei extractions were mixed with those from standards, stained with 1/50 volume of propidium iodide solution (2.04 mg ml -1 ) and incubated at 4 °C in darkness for 1 h. For each species, three replicates were performed on the BD FACSCan system. For S. cucullata, we used A. thaliana (0.165 pg per C) [START_REF] Praça-Fontes | C-value reassessment of plant standards: an image cytometry approach[END_REF] as the standard, and for all other samples, we used Zea mays 'CE-777' (2.785 pg per C) [START_REF] Praça-Fontes | C-value reassessment of plant standards: an image cytometry approach[END_REF] . For each peak (in both standard and sample), over 1,000 nuclei were collected with cross-validation values lower than 5%, except for those of A. thaliana 2n nuclei peaks, which ranged from 5.5% to 5.9%. To calculate the 2C-value of S. cucullata, we used a formula of: (0.66 pg × (F -S2n) + 0.33 pg × (S4n -F))/(S4n -S2n). For all other samples, we used: 5.57 pg × F/S2n, where 0.66, 0.33 and 5.57 pg are the 4C-value of A. thaliana, the 2C-value of A. thaliana and the 2C-value of Z. mays 'CE-777' , respectively. S2n, S4n and F are the relative fluorescence amount (that is, the peak mean value) of the standard 2n nuclei, standard 4n nuclei and the sample 2n nuclei, respectively.

Genome and transcriptome sequencing. A. filiculoides was collected from the Galgenwaard ditch in Utrecht, the Netherlands, and propagated directly or sterilized as described in Dijkhuizen et al. [START_REF] Dijkhuizen | Is there foul play in the leaf pocket? The metagenome of floating fern Azolla reveals endophytes that do not fix N 2 but may denitrify[END_REF] . A. filiculoides (sterilized without cyanobiont) DNA was extracted, then sequenced on PacBio RSII at 51× coverage [START_REF] Dijkhuizen | Is there foul play in the leaf pocket? The metagenome of floating fern Azolla reveals endophytes that do not fix N 2 but may denitrify[END_REF] and Illumina HiSeq2000 (100 bp paired-end; ~86× coverage; Supplementary 1) and DNA was extracted by a modified cetyltrimethylammonium bromide (CTAB) procedure [START_REF] Beck | Does hybridization drive the transition to asexuality in diploid Boechera?[END_REF] . Illumina libraries with a 500-bp insert size were prepared and sequenced on Illumina HiSeq2000 (100 bp paired-end; ~50× coverage; Supplementary Table 12). S. cucullata was originally collected from Bangladesh and subsequently cultured at Taiwan Forestry Research Institute, Dr. Cecilia Koo Botanic Conservation Center and Duke University (Supplementary Table 1). Genomic DNA was purified using a modified CTAB procedure [START_REF] Beck | Does hybridization drive the transition to asexuality in diploid Boechera?[END_REF] and sequenced on both PacBio RSII (10 SMRT cells; 46× coverage) and Illumina HiSeq2000 platforms (1 lane of 125 bp pairedend; 215× coverage; Supplementary Table 12). S. cucullata RNA from the floating and submerged leaves was separately extracted using the Sigma Spectrum Plant Total RNA kit, each with three biological replicates. To examine patterns of RNA editing, one library per leaf type was prepared by the Illumina Ribozero Plant kit (that is, not poly-A enriched), whereas the other two were done by the Kapa Stranded mRNA-seq kit. These six RNA libraries were pooled and sequenced in one lane of Illumina HiSeq2000 (125 bp paired-end).

Genome assembly. We assembled the PacBio reads from A. filiculoides and S. cucullata genomes using PBcR [START_REF] Berlin | Assembling large genomes with single-molecule sequencing and locality-sensitive hashing[END_REF] , and the resulting drafts were then polished by Quiver [START_REF] Chin | finished microbial genome assemblies from long-read SMRT sequencing data[END_REF] (A. filiculoides) or Pilon [START_REF] Walker | Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement[END_REF] (S. cucullata). Plastid genomes were separately assembled using Mitobim [START_REF] Hahn | Reconstructing mitochondrial genomes directly from genomic next-generation sequencing reads-a baiting and iterative mapping approach[END_REF] and annotated in Geneious [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF] with manual adjustments. The PBcR contigs were filtered to remove plastome fragments. Although the A. filiculoides strain we sequenced was surface sterilized and treated with antibiotics to remove its associated microbiome, other endophytes could still persist, as shown by Dijkhuizen et al. [START_REF] Dijkhuizen | Is there foul play in the leaf pocket? The metagenome of floating fern Azolla reveals endophytes that do not fix N 2 but may denitrify[END_REF] . Thus, we thoroughly assessed the A. filiculoides and S. cucullata assemblies to filter out all possible non-plant scaffolds. We used BlobTools 66 in combination with SILVA [START_REF] Quast | The SILVA ribosomal RNA gene database project: improved data processing and web-based tools[END_REF] and UniProt 68 databases to infer the taxonomy for each scaffold. We removed all scaffolds that were classified as bacteria or fungi and also those that had a skewed GC content and read coverage. The completeness of each final assembly was assessed by BUSCO [START_REF] Simão | BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs[END_REF] with the Plants set, and by using BWA [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] and HISAT2 [START_REF] Kim | HISAT: a fast spliced aligner with low memory requirements[END_REF] to map Illumina reads to the assemblies (Supplementary Table 2).

Repeat annotation. RepeatModeler 72 was used to generate species-specific repeat libraries for masking and annotation. Consensus repeat sequences with homology to known plant genes were removed from the repeat libraries. Homology was defined as having a significant (E-value < 1 × 10 -5 ) blastx [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] hit to a subset of the PlantTribes 74 v1.1 database that does not contain transposable element-related terms. Filtered RepeatModeler libraries were annotated with the name of the highest-scoring significant Repbase 75 v22.04 full database sequence (tblastx [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] , E-value < 1 × 10 -5 ) and the highest-scoring significant Dfam 76 v2.0 profile hidden Markov model (HMM) (hmmsearch [START_REF] Eddy | Accelerated profile HMM searches[END_REF] , E-value < 1 × 10 -5 ).

LTR-RTs were discovered using structural criteria by the GenomeTools 78 program LTRHarvest [START_REF] Ellinghaus | an efficient and flexible software for de novo detection of LTR retrotransposons[END_REF] with the following modifications to the default settings: a LTR similarity threshold of 0.01, an allowed LTR length range of 100-6,000 bp, an allowed distance between LTRs of a single element range of 1,000-25,000 bp and the number of bases outside LTR boundaries to search for target-site duplications set to 10. The GenomeTools program LTRDigest 80 was used with a set of 138 transposable element-related Pfam profile HMMs to annotate protein-coding domains in the internal regions of LTR-RTs.

We used 38 previously published A. filiculoides RNA-seq libraries and 6 S. cucullata libraries (see above) to assemble transcriptomes for facilitating gene model predictions. Reads from A. filiculoides and S. cucullata libraries were processed using a combination of Scythe [START_REF] Buffalo | Scythe-a Bayesian adapter trimmer v.0.994 beta[END_REF] and Sickle [START_REF] Joshi | Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files[END_REF] or SOAPnuke [START_REF] Chen | SOAPnuke: a MapReduce acceleration supported software for integrated quality control and preprocessing of high-throughput sequencing data[END_REF] , with adapter and contaminating sequences discovered using FastQC 84 (v0.11.5). Approximately 627 million (A. filiculoides) and 259 million (S. cucullata) cleaned paired reads went into the assemblies. Libraries from experimental replicates were combined and assembled de novo by Trinity 85 (v2.0.6) and in a reference-guided manner using HISAT2 71 (v2.0.4) and StringTie 86 (v1.2.2), except for nine libraries published in de Vries et al. [START_REF] De Vries | Cytokinin-induced promotion of root meristem size in the fern Azolla supports a shoot-like origin of euphyllophyte roots[END_REF] for which only a reference-guided approach was used. All programs used default parameters, and Trinity was run with the additional --trimmomatic option. StringTie results were merged using StringTie --merge, combined with the Trinity output, and were purged of redundant sequences using the GenomeTools sequniq utility [START_REF] Gremme | GenomeTools: a comprehensive software library for efficient processing of structured genome annotations[END_REF] .

Putative centromere sequences were first identified by searching the genome assemblies with Tandem Repeat Finder 87 to identify very high-copy (> 100 repeats) tandem repeats with a motif length in the range of 185-195 bp. Motif sequences were extracted from the Tandem Repeat Finder output and clustered using USEARCH [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF] . A single major cluster was identified for each species and the sequences were separately aligned using MAFFT [START_REF] Katoh | MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform[END_REF] . Multiple sequence alignments for each species were used to generate a profile HMM representing the putative centromere sequences. We next used hmmsearch [START_REF] Eddy | Accelerated profile HMM searches[END_REF] to search the genome assemblies again to identify all regions with similarity to the centromere profile HMMs. Genomic regions with significant HMM matches were identified and these regions were annotated in a GFF3 format.

Gene prediction. Protein-coding genes were predicted using MAKER-P 90 (v2.31.8), and three MAKER-P iterations were performed: (1) repeat masking and creation of initial gene models from transcript and homologous protein evidence;

(2) refinement of initial models with SNAP [START_REF] Korf | Gene finding in novel genomes[END_REF] ab initio gene predictor trained on initial models; and (3) final models generated using SNAP and the ab initio gene predictor AUGUSTUS [START_REF] Stanke | AUGUSTUS: a web server for gene prediction in eukaryotes that allows user-defined constraints[END_REF] trained on gene models from the second iteration.

Masking was performed by RepeatMasker 93 (v4.0.5) using the previously described species-specific repeat libraries and the full Repbase v22.04 database. After masking, gene models were inferred from transcripts and homologous protein sequences by first aligning to the genomes using BLAST+ [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] (v2.2.31) blastn or blastp, and then refined using the functions est2genome and protein2genome from the splice-site aware alignment program Exonerate 94 (v2.2.0). We included the previously described A. filiculoides or S. cucullata transcriptomes and the set of protein sequences consisting of the full Swiss-prot database (downloaded 18 June 2016), Amborella trichopoda v1.0 proteins, A. thaliana TAIR10 proteins, Chlamydomonas reinhardtii v5.5 proteins, Oryza sativa v7.0 proteins and Physcomitrella patens v3.3 proteins (from Phytozome [START_REF] Goodstein | Phytozome: a comparative platform for green plant genomics[END_REF] ). Gene models with an annotation edit distance (AED) score of < 0.2 were used to train SNAP, which was used during the second iteration of MAKER-P. SNAP was retrained for the final iteration using gene models from the second iteration with an AED score of < 0.2 and a translated protein length of > 200 amino acids. Prior to training AUGUSTUS [START_REF] Stanke | AUGUSTUS: a web server for gene prediction in eukaryotes that allows user-defined constraints[END_REF] , redundant sequences, defined as those sharing ≥ 70% sequence similarity in significant (E-value < 1 × 10 -5 ) HSPs from an all-by-all blastn alignment, were removed from the training set. Final non-redundant sets of 5,013 (A. filiculoides) or 6,475 (S. cucullata) gene models were used to train AUGUSTUS [START_REF] Stanke | AUGUSTUS: a web server for gene prediction in eukaryotes that allows user-defined constraints[END_REF] .

Phylogenomic inference and placement of WGDs from nuclear gene trees.

To infer ancient WGDs, we used a gene-tree sorting and counting algorithm, implemented in the MAPS tool [START_REF] Li | Early genome duplications in conifers and other seed plants[END_REF] . We selected four species of heterosporous ferns (two Azolla, one Salvinia and one Pilularia) and representatives of three other leptosporangiate lineages (Blechnum, Lygodium and Dipteris). The MAPS algorithm uses a given species tree to filter collections of nuclear gene trees for subtrees consistent with relationships at each node in the species tree. Using this filtered set of subtrees, MAPS identifies and records nodes with a gene duplication shared by descendant taxa. To infer and locate a potential WGD, we plotted the percentage of gene duplications shared by descendant taxa by node: a WGD will produce a large burst of shared duplications, appearing as an increase in the percentage of shared gene duplications [START_REF] Li | Early genome duplications in conifers and other seed plants[END_REF] .

We circumscribed and constructed nuclear gene family phylogenies from multiple species for each MAPS analysis. We translated each transcriptome into amino acid sequences using the TransPipe pipeline [START_REF] Barker | net: bioinformatic tools for ecological and evolutionary genomics[END_REF] . Using these translations, we performed reciprocal protein BLAST (blastp) searches among data sets for each MAPS analysis using an E-value cut-off of 10 -5 . We clustered gene families from these BLAST results using OrthoFinder with the default parameters [START_REF] Emms | OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy[END_REF] and only kept gene families that contained at least one gene copy from each taxon in a given MAPS analysis. We discarded the remaining OrthoFinder clusters. We used PASTA [START_REF] Mirarab | PASTA: ultra-large multiple sequence alignment for nucleotide and amino-acid sequences[END_REF] for automatic alignment and phylogenetic reconstruction of gene families, employing MAFFT [START_REF] Katoh | MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform[END_REF] for constructing alignments, MUSCLE [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF] for merging alignments and RAxML 100 for tree estimation. The parameters for each software package were the default options for PASTA. For each gene family phylogeny, we ran PASTA until we reached three iterations without an improvement in the likelihood score using a centroid breaking strategy. We used the best-scoring PASTA tree for each multi-species nuclear gene family to infer and locate WGDs using MAPS.

For the null simulations, we first estimated the mean background gene duplication rate (λ) and the gene loss rate (μ) with WGDgc 101 . Gene count data were obtained from OrthoFinder clusters associated with each species tree. λ = 0.0031 and μ = 0.0039 were estimated using only gene clusters that spanned the root of their respective species trees, which has been shown to reduce biases in the maximum likelihood estimates of λ and μ 101 . We chose a maximum gene family PPR annotation. We conducted a targeted annotation for PPR genes because they are generally only weakly expressed and thus often lack transcriptome support. Open reading frames from the nuclear genome assemblies were translated into amino acid sequences using the "getorf " tool from the EMBOSS (v.6.5.7) package 113 with a minimum size restriction of 300 nucleotides. These open reading frames were searched for PPR motifs using the hmmsearch tool from the HMMER3 package [START_REF] Eddy | Accelerated profile HMM searches[END_REF] . The PPR motif models and parameters used follow those of Cheng et al. [START_REF] Cheng | Redefining the structural motifs that determine RNA binding and RNA editing by pentatricopeptide repeat proteins in land plants[END_REF] . Motifs were assembled into full PPR tracts and the best model for each PPR was determined [START_REF] Cheng | Redefining the structural motifs that determine RNA binding and RNA editing by pentatricopeptide repeat proteins in land plants[END_REF] .

To study the prevalence and location of RNA editing, non-poly(A)-enriched RNA-seq data were filtered to remove adapters, low-quality reads and reads with ≥ 5% Ns. Clean reads were aligned against the assembled plastid and mitochondrial genome assemblies using TopHat 2.0 (ref. 114 ). One of the inverted repeat regions in the plastid genomes was removed before mapping. Only uniquely mapped reads were retained as input for SAMtools 115 to call mismatches between RNA and the corresponding DNA. Differences between corresponding RNA and DNA sequences were identified as the putative RNA-editing sites. The RNA-editing level was defined as the number of altered reads divided by the total mapped reads for each site.

Phylogeny of the insecticidal protein Tma12. We used BLASTp [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] to search for Tma12 (Genbank accession: JQ438776) homologues in Phytozome [START_REF] Goodstein | Phytozome: a comparative platform for green plant genomics[END_REF] , 1KP transcriptomes [START_REF] Matasci | Data access for the 1,000 Plants (1KP) project[END_REF] and the NCBI Genbank non-redundant protein database. Although Tma12 homologues are present in fern transcriptomes and in the S. cucullata genome, no significant hit was found in any other plant genomes or transcriptomes. In addition, the majority of the Tma12 protein is composed of a chitin-binding domain that belongs to the PF03067 Pfam family. This family does not contain any plant genes but is predominantly represented in the genomes of Actinobacteria, insects and fungi. To trace the origin of fern Tma12 genes, we downloaded representative sequences containing PF03067 and PF08329 (as the outgroup) from UniProt and Genbank and reconstructed the phylogeny using IQ-TREE 116 . We then used this preliminary phylogeny (Supplementary Data) to construct a more focused data set to narrow down the phylogenetic affinity of Tma12. PartitionFinder 117 was used to infer the optimal codon partition scheme and substitution models, and RAxML 100 was used for maximum likelihood phylogeny inference and to calculate bootstrap branch support. Azolla phylogeny. From the resequencing data (Supplementary Table 12), we compiled both plastome and nuclear phylogenomic data sets to infer the Azolla species phylogeny. S. cucullata was used as the outgroup. For the plastome phylogeny, we concatenated nucleotide alignments from 83 protein-coding genes and used PartitionFinder 117 to identify the optimal data partition scheme and the associated nucleotide substitution models. RAxML 100 was used for maximum likelihood phylogeny inference and to calculate bootstrap branch support. For the nuclear data set, we focused on genes that, based on the gene family classification results, are single copy in both A. filiculoides and S. cucullata genomes. We used HybPiper 118 to extract the exon sequences from each of the resequenced species. The 'bwa' option was used in HybPiper instead of the 'blastx' default. We then filtered out genes with more than two species missing or having an average sequence length shorter than 75% of the one in A. filiculoides. This resulted in a final data set of 2,108 genes. Sequence alignments and gene tree inferences were done in PASTA [START_REF] Mirarab | PASTA: ultra-large multiple sequence alignment for nucleotide and amino-acid sequences[END_REF] , with the default setting, except that RAxML 100 was used to estimate the best tree on the final alignment. To infer the species tree from these gene trees, we used the multi-species coalescent method implemented in ASTRAL-III (v5.6.1) 119 . The tree topology from the plastome and nuclear data sets were identical, and all nodes received bootstrap support of 100 and a local posterior probability of 1.0. Cyanobiont phylogeny. To compare the host and symbiont phylogenies, we assembled the cyanobiont genomes from five additional Azolla species (Supplementary Table 12) using the resequencing data generated from total DNAs, including sequences derived from both the host and the cyanobiont. To extract the cyanobiont genomes from each of the Azolla species, we first filtered out chloroplast sequences by using BWA [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] (default parameters) to map the total clean DNA reads against each chloroplast genome reference. In this step, ~3-4% of the reads were filtered out, which is necessary to remove plastid ribosomal RNAs that are highly similar to ones in the cyanobionts. For each of the five Azolla species, we then mapped the filtered reads to the published cyanobiont reference (N. azollae 0708 isolated from A. filiculoides [START_REF] Ran | Genome erosion in a nitrogen-fixing vertically transmitted endosymbiotic multicellular cyanobacterium[END_REF] ; Genbank accession: NC_014248) using BLAST [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] (alignment criteria: E-value ≤ 1 × 10 -5 , sequence identity of ≥ 90% and an aligned coverage of ≥ 90%). Only the aligned reads were assembled by Mitobim 64 (iterations = 5) using N. azollae 0708 (ref. [START_REF] Ran | Genome erosion in a nitrogen-fixing vertically transmitted endosymbiotic multicellular cyanobacterium[END_REF] ) as a reference. Gene prediction for each assembled cyanobiont was performed by the Prodigal program 120 . Transfer RNAs were predicted by tRNAscan-SE 121

Statistical parameters

When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main text, or Methods section).

n/a Confirmed

The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.

A description of all covariates tested

A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated

Clearly defined error bars

State explicitly what error bars represent (e.g. SD, SE, CI)

Our web collection on statistics for biologists may be useful.

Software and code

Policy information about availability of computer code

Data collection

Provide a description of all commercial, open source and custom code used to collect the data in this study, specifying the version used OR state that no software was used.

Data analysis

We employed a number of software for data analyses in this study, which was described in detail in the materials and methods, including the parameters used, versions (if applicable), and citations.

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information. The axis labels state the marker and fluorochrome used (e.g. CD4-FITC).

nature research | reporting summary

The axis scales are clearly visible. Include numbers along axes only for bottom left plot of group (a 'group' is an analysis of identical markers).

All plots are contour plots with outliers or pseudocolor plots.

A numerical value for number of cells or percentage (with statistics) is provided.

Methodology Sample preparation

We used flow cytometry to estimate the genome sizes of Pilularia americana, Regnellidium diphyllum, Marsilea minuta, Salvinia cucullata.

Instrument

BD FACSCan system (BD Biosciences, Franklin Lake, NJ, USA) Software BD FACSCan system (BD Biosciences, Franklin Lake, NJ, USA)

Cell population abundance Pilularia americana:

Replicate 1: sample peak particle number = 1514, standard1 peak particle number = 1154. Replicate 2: sample peak particle number = 1834, standard1 peak particle number = 1371. Replicate 3: sample peak particle number = 1450, standard1 peak particle number = 1036.

Regnellidium diphyllum: Replicate 1: sample peak particle number = 1222, standard1 peak particle number = 1737. Replicate 2: sample peak particle number = 1180, standard1 peak particle number = 1613. Replicate 3: sample peak particle number = 1137, standard1 peak particle number = 1759.

Marsilea minuta: Replicate 1: sample peak particle number = 1892, standard1 peak particle number = 1118. Replicate 2: sample peak particle number = 1850, standard1 peak particle number = 1209. Replicate 3: sample peak particle number = 1892, standard1 peak particle number = 1227.

Salvinia cucullata: Replicate 1: sample peak particle number = 1084, standard1 peak particle number = 1484, standard2 peak particle number = 1170. Replicate 2: sample peak particle number = 1129, standard1 peak particle number = 1552, standard2 peak particle number = 1253. Replicate 3: sample peak particle number = 1229, standard1 peak particle number = 1584, standard2 peak particle number = 1500.

Gating strategy

For particle acquisition, we set a threshold of FL2-H = 52 for the samples of Pilularia americana, Regnellidium diphyllum, and Marsilea minuta. For Salvinia cucullata, a threshold of FL2-H = 100 is applied.

Tick this box to confirm that a figure exemplifying the gating strategy is provided in the Supplementary Information.

Fig. 1 |

 1 Fig. 1 | Genome size evolution in Salviniales. a, Members of Salviniales have smaller genome sizes than other ferns (averaging 1C = 12 Gb) 6 . Two wholegenome duplication (WGD) events identified in this study were mapped onto the phylogeny, with divergence time estimates obtained from Testo and Sundue 129 . b,c, Whole genomes were assembled from A. filiculoides (b) and S. cucullata (c). d,e, The genome of S. cucullata has substantially reduced levels of RNA (d) and DNA (e) transposons compared to A. filiculoides. Image in panel c courtesy of P.-F. Lu.

Fig. 2 |

 2 Fig. 2 | Evolution of ethylene biosynthesis.The ethylene-forming pathway involves the Yang cycle, where ACC is synthesized from S-adenosyl-methionine (SAM; also known as AdoMet) by ACC synthase. ACC oxidase then catalyses the conversion of ACC to ethylene. We found that ACC oxidase is unique to seed plants (green) and its origin probably drove the expansion of the ACC synthase gene family (orange; Supplementary Fig.6) to create a regulated ethylene biosynthetic mechanism.

Fig. 3 |

 3 Fig.3| The history of WGD in Azolla and Salvinia. a, MAPS analysis identified two WGD events: one specific to Azolla (orange circle) and one predating the core leptosporangiates (green circle). The blue line illustrates the percentage of subtrees indicative of a gene duplication shared by the descendants at each node. The grey lines display the gene birth-death simulation results without WGD. The species divergence dates are from Testo and Sundue 129 . b, Density plots from fitting Gaussian mixture models to K s distributions estimated from pairs of syntenic paralogues within the Azolla and Salvinia genomes, as well as of syntenic orthologues between Azolla and Salvinia. c, Examples of synteny between Azolla and Salvinia genomic regions. The left and right panels display a 2:1 and 2:2 syntenic relationship between Azolla and Salvinia regions, respectively. Each subpanel represents a genomic region in Azolla or Salvinia, with gene models on both strands shown above and below the dashed line. High-scoring sequence pairs (HSPs) in protein-coding sequences are marked by short vertical bars above the gene models. Selected HSP links between genomic regions are depicted as coloured lines crossing the subpanels, whereas others (for example, the HSP links between the two Azolla genomic regions in the left panel) are left out for clarity. Collinear series of HSPs across genomic regions indicates a syntenic relationship between the regions concerned. Genomic regions conserved in duplicate after the WGD that occurred prior to the divergence between Azolla and Salvinia should show a 2:2 syntenic relationship, whereas regions conserved in duplicate after the Azolla-specific WGD should show a 2:1 syntenic relationship with Salvinia regions. The left and right panels can be regenerated at https://genomevolution. org/r/ujll and https://genomevolution.org/r/ukys, respectively.

Fig. 4 |

 4 Fig. 4 | Origin of a fern insecticidal protein. Phylogenetic analysis of the chitin-binding domain Pfam PF03067 shows that the fern Tma12 insecticidal protein was probably derived from bacteria through an ancient HGT event. The numbers above the branches are bootstrap (BS) support values (BS = 100 is omitted), and the thickened branches indicate BS > 70. The tree is rooted based on the result from a broader phylogenetic analysis of PF03067 and PF08329 (Supplementary Data). The pink star denotes the sequence from the S. cucullata genome.

Fig. 5 |

 5 Fig.5| Azolla-cyanobacterial symbiosis. a, The cyanobiont phylogeny largely mirrors the host species phylogeny, indicating a convincing cospeciation pattern between the two partners. All nodes received a maximum likelihood bootstrap support of 100%, and for the host phylogeny, all nodes also received a local posterior probability of 1.0 from the ASTRAL 119 analysis. Both the nuclear and the plastome data sets gave the same topology for the host, and the branch lengths shown here were from the plastome tree. Scale bars represent 0.01 substitutions per site. b, The CSP genes were lost in the Azolla and Salvinia genomes (empty boxes), whereas orthologues can be found in other fern transcriptomes (red boxes). *Arabidopsis lacks the CSP genes and does not have AM symbiosis. c, Cyanobionts have a large effect on the Azolla transcriptome. d, The Azolla transcriptome responds to nitrogen starvation more significantly when cyanobionts are absent than when they are present. PC, principal component. e, Candidate genes involved in nutrient transport and communication with cyanobionts.

  using a bacterial tRNA gene structure model. The presence of rRNA sequences (gene number and structure) for each cyanobiont 1 nature research | reporting summary April 2018 Corresponding author(s): Fay-Wei Li Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see Authors & Referees and the Editorial Policy Checklist.

Table 12 )

 12 with library insert sizes of 175 bp and 340 bp. RNA sequencing (RNA-seq) data from A. filiculoides of the Galgenwaard ditch used for annotation included the following RNA profiles: (1) at four time points during the diel cycle of fern sporophytes grown with or without 2 mM ammonium nitrate for 1 week 42 ; (2) of different reproductive stages comparing fern sporophytes, microsporocarps and megasporocarps collected at noon 58 ; (3) of roots treated with cytokinin, indole-3acetic acid (IAA) or none[START_REF] De Vries | Cytokinin-induced promotion of root meristem size in the fern Azolla supports a shoot-like origin of euphyllophyte roots[END_REF] ; and (4) of sporophytes with or without cyanobacterial symbionts grown with or without ammonium nitrate for 2 weeks then collected at noon. Plant materials of A. caroliniana, Azolla mexicana, Azolla microphylla, Azolla nilotica and Azolla rubra were obtained from the International Rice Research Institute (Supplementary Table

  1. Prepare buffer for use. a. Allocate appropriate amount of Backmen stock buffer to a 50-ml tube based on an estimation of 1-1.5 ml per sample. b. Add 0.04 g PVP-40, 5 μl 2-mercaptoethanol, 1 μl RNase per ml of buffer. 2. Extract sample and standard nuclei by chopping leaf tissue a. Add 500 μl of buffer to a glass Petri dish. b. Add a (~400 mm2) piece of young leaf to the Petri dish, and chop it with a razor on ice until most tissue slices are less than 1 mm in size. c. Filter the chopped sample and standard into a 2.0-ml tube through a 30-μm nylon mesh. d. Add additional buffer to the sample, and ensure that the filtered leaf nuclei solution is greater than 500 μl in volume or more depending on need. 3. Staining nuclei solutions a. Mix sample nuclei and standard leaf nuclei solutions into a 500-μl volume in 2.0-ml tubes. b. Add 10 μl PI solution (2.04 mg/ml ) into each of mixed nuclei solutions. c. Incubate in the dark at 4 °C for 1 h for staining.

	Recipes
	Backmen stock buffer
	1.0% Triton X-100
	50 mM Na2SO3
	50 mM Tris-HCl (pH 7.5)
	ddH2O (the solvent)
	Note: Store at 4 °C up to 1 year.

Gb[START_REF] Hidalgo | Is there an upper limit to genome size?[END_REF] and averaging 12 Gb[START_REF] Sessa | Evolutionary genomics of ferns and lycophytes[END_REF] . Azolla is one of the fastest-growing plants on the planet, with demonstrated potential to be a significant carbon sink. Data from the Arctic Ocean show that, ~50 Myr ago, in middle-Eocene sediments, an abundance of fossilized Azolla characterizes an ~800,000-year interval known as the 'Azolla event'[START_REF] Brinkhuis | Episodic fresh surface waters in the Eocene Arctic Ocean[END_REF] . This period coincides with the shift from the early Eocene greenhouse world towards our present icehouse climate, suggesting that Azolla had a role in abrupt global cooling by sequestering atmospheric carbon dioxide[START_REF] Speelman | The Eocene Arctic Azolla bloom: environmental conditions, productivity and carbon drawdown[END_REF] . Azolla is also remarkable in harbouring an obligate, N 2 -fixing cyanobacterium, Nostoc azollae, within specialized leaf cavities. Because of this capability, Azolla has been used as 'green manure' for over 1,000 years to bolster rice productivity in Southeast

NATuRE PLANTS | VOL 4 | JULY 2018 | 460-472 | www.nature.com/natureplants

Acknowledgements

We are grateful to 123 backers from Experiment.com who crowdfunded the initial work. We thank Z. Fei for providing comments and suggestions, M. Harrison for discussion on Azolla transporter genes, J. Shaw and D. Weston for providing access to the Sphagnum genome and T. Nishiyama for the Jungermannia genome, and P.-F. Lu for providing the image for Fig. 1c. This project was partly supported by the Shenzhen Municipal Government of China (no. JCYJ20150529150409546), the National Science Foundation Doctoral Dissertation Improvement Grant DEB-1407158 (to K.M.P. and F.-W.L.) and the German Research Foundation Research Fellowship VR132/1-1 (to J.d.V.). Computational support was provided by the Duke Compute Cluster and the Centre for Information and Media Technology at University of Düsseldorf.

Articles

NaTure PlaNTS size of 100 for parameter estimation, which was necessary to provide an upper bound for numerical integration of node states 101 . We provided a prior probability distribution of 1.5 on the number of genes at the root of each species tree, such that ancestral gene family sizes followed a shifted geometric distribution with a mean equal to the average number of genes per gene family across species.

Gene trees were then simulated within each MAPS species tree using the GuestTreeGen program from GenPhyloData 102 . We developed ultrametric species trees from the topological relationships inferred by the 1KP Consortium analyses and median branch lengths from TimeTree 103 . For each species tree, we simulated 4,000 gene trees with at least one tip per species: 1,000 gene trees at the λ and μ maximum likelihood estimates, 1,000 gene trees at half the estimated λ and μ, 1,000 trees at three times λ and μ, and 1,000 trees at five times λ and μ.

Classification of syntenic duplicates and microsynteny analysis. To distinguish gene duplicates as syntenic or tandem, we used the SynMap 104 tool from the CoGe 105 platform, with default parameters and the Quota Align algorithm to merge syntenic blocks. Sets of syntenic paralogues or orthologues (defined by a collinear series of putative homologous genes) were extracted using the DAGChainer algorithm, whereas duplicates within ten genes apart in the same genomic region were identified as tandem duplicates (Supplementary Discussion). Results for within Azolla and Salvinia genome comparisons, as well as between Azolla and Salvinia, can be regenerated using the links https://genomevolution. org/r/tozk, https://genomevolution.org/r/toz7 and https://genomevolution.org/r/ toyy, respectively. Microsynteny analyses were performed using the GEvo tool from CoGe 105 . We used the default setting to define the minimum number of collinear genes for two regions to be called syntenic. Non-coding regions were masked in the two genomes to include only the protein-coding sequences. The two example microsyntenies shown in Fig. 5c can be regenerated at https://genomevolution. org/r/ujll and https://genomevolution.org/r/ukys.

Gaussian mixture model analysis of K s distributions. Estimates of K s were obtained for all pairs of syntenic paralogous and orthologous genes using the CODEML program 106 in the PAML package (v4.8) 107 on the basis of codon sequence alignments. We used the GY model with stationary codon frequencies empirically estimated by the F3 × 4 model. Codon sequences were aligned with PRANK (version 100701) using the empirical codon model 108 (setting -codon) to align coding DNA, always skipping insertions (-F). Only gene pairs with K s values in the range of 0.05-5 were considered for further analyses. Gaussian mixture models were fitted to the resulting frequency distributions of K s values by means of the densityMclust function in the R mclust version 5.3 package 109 . The Bayesian information criterion was used to determine the best-fitting model for the data, including the optimal number of Gaussian components to a maximum of nine. For each component, several parameters were computed including the mean and the variance, as well as the density mixing probabilities and the total number of gene pairs. Gene family classification and ancestral reconstruction. The OrthoFinder [START_REF] Emms | OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy[END_REF] clustering method was used to classify complete proteomes of 23 sequenced green plant genomes, including A. filiculoides and S. cucullata (Supplementary Table 5), into orthologous gene lineages (that is, orthogroups). We selected taxa that represented all of the major land plant and green algal lineages, including six core eudicots (A. thaliana, Lotus japonicus, Populus trichocarpa, Solanum lycopersicum, Erythranthe guttata and Vitis vinifera), five monocots (O. sativa, Sorghum bicolor, Musa acuminata, Zostera marina and Spirodella polyrhiza), one basal angiosperm (A. trichopoda), two gymnosperms (Pinus taeda and Picea abies), two ferns (A. filiculoides and S. cucullata), one lycophyte (S. moellendorffii), four bryophytes (Sphagnum fallax, P. patens, Marchantia polymorpha and Jungermannia infusca) and two green algae (Klebsormidium flaccidum and C. reinhardtii). In total, 16,817 orthogroups containing at least two genes were circumscribed, 8,680 of which contain at least one gene from either A. filiculoides or S. cucullata. Of the 20,203 annotated A. filiculoides genes and the 19,780 annotated S. cucullata genes, 17,941 (89%) and 16,807 (84%) were classified into orthogroups, respectively. The details for each orthogroup, including gene counts, secondary clustering of orthogroups (that is, super-orthogroups) 110 and functional annotations, are reported in Supplementary Table 5.

We used Wagner parsimony implemented in the program Count 111 with a weighted gene gain penalty of 1.2 to reconstruct the ancestral gene content at key nodes in the phylogeny of the 23 land plants and green algae species (Supplementary Table 5). The ancestral gene content dynamics-gains, losses, expansions and contractions-are depicted in Supplementary Fig. 5. Complete details of orthogroup dynamics for the key ancestral nodes that include seed plants, such as Salviniaceae, euphyllophytes and vascular plants, are reported in Supplementary Table 5.

Transcription-associated protein characterization.

Transcription-associated proteins comprise transcription factors that bind in a sequence-specific manner to cis-regulatory DNA elements and transcriptional regulators that act via protein-protein interaction or chromatin modification. We conducted genomewide, domain-based annotation of transcription-associated proteins according to previous studies [START_REF] Wilhelmsson | Comprehensive genome-wide classification reveals that many plant-specific transcription factors evolved in streptophyte algae[END_REF]112 . A total of 1,206 (6%, Azolla) and 983 (7%, Salvinia) proteins were sorted into families; this amount is similar to Selaginella but less than in gymnosperms or angiosperms (Supplementary Table 8).

Articles

NaTure PlaNTS

was confirmed by mapping the rRNAs of N. azollae 0708 against each assembled genome cyanobiont sequence using BLAST. We used mugsy 122 to generate the whole-genome alignment, which resulted in a nucleotide matrix of 5,354,840 characters. IQ-TREE 116 was used for model testing and maximum likelihood tree inference. Because the N. azollae genome is reduced in size and is significantly diverged from other cyanobacteria, we could not find an appropriate outgroup to root the cyanobiont tree. To overcome this, we used STRIDE 123 implemented in OrthoFinder [START_REF] Emms | OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy[END_REF] to locate the root by reconciling gene trees. STRIDE was run with the default setting, except that MAFFT [START_REF] Katoh | MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform[END_REF] was used for alignment and RAxML 100 for tree inference. The root was found to be the node placing the A. nilotica cyanobiont as sister to a clade comprising all other cyanobionts. The reconciled species tree is identical to the tree reconstructed from the whole-genome alignment.

Identification of the CSP genes. The Medicago truncatula DMI2, DMI3, IPD3, CASTER/POLLUX and VAPYRIN sequences were used as queries, as in a previous study [START_REF] Delaux | Algal ancestor of land plants was preadapted for symbiosis[END_REF] , to search against the genomes and transcriptomes from species listed in Supplementary Table 11 using tBLASTn [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] . For liverworts and ferns from the 1KP data set [START_REF] Matasci | Data access for the 1,000 Plants (1KP) project[END_REF] , non-annotated transcriptomes were used as targets, with the longest open reading frame of each contig extracted and translated. For A. filiculoides and S. cucullata, both the annotated gene models and the unannotated scaffolds were used. All hits that matched already annotated gene models were discarded prior to subsequent analyses. No homologues were identified in the two fern genomes for IPD3 and VAPYRIN. Protein sequences for DMI2/SYMRK, DMI3/CCaMK and CASTOR/POLLUX were aligned using MAFFT [START_REF] Katoh | MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform[END_REF] . The best substitution model for each alignment (JTT for all alignments) was determined using MEGA6 (ref. 124 ). Phylogenetic trees were generated using RAxML 100 on the CIPRES platform 125 , and node support was assessed with 100 rapid bootstrap pseudoreplicates.

Quantitative real-time PCR of nifH. Quantitative real-time PCR for the N. azollae nifH gene was conducted using total RNA extracted from A. filiculoides. Primers were derived from Brouwer et al. [START_REF] Brouwer | Azolla domestication towards a biobased economy?[END_REF] . ThermoFisher Superscript IV was used to generate complementary DNA from the RNA. The cDNA was then used for quantitative PCR with the Roche SYBR Green Master Mix on a Chromo4 realtime PCR machine with the Opticon platform. The relative gene expression was calculated using the 2 ΔC(t) method, with the cyanobacteria present/nitrogen absent condition as the reference.

Azolla symbiosis transcriptome analysis. We used RNA-seq to compare gene expression patterns of AzCy+ and AzCy-individuals grown with or without ammonium nitrate. Each condition and treatment combination has three biological replicates. RNA-seq reads were mapped to the A. filiculoides genome by HISAT2 [START_REF] Kim | HISAT: a fast spliced aligner with low memory requirements[END_REF] , and read counts for each gene were calculated using the HTSeq software package 126 . We used the rlog function in the DESeq2 package 127 for data normalization and carried out differential expression analysis in DESeq2 to identify upregulated and downregulated genes with an adjusted P value of 0.005. Distance clustering and principal component analysis were used to examine the relatedness of samples and conditions as a quality-control measure.

Azolla-cyanobacteria HGT. To identify cyanobiont-derived genes in the A. filiculoides genome, we first investigated a potential orthologous relationship between any Azolla genes and cyanobacteria. For this, we used the Azolla genome assembly as a query for a DIAMOND BLASTx 128 against a protein data set of 11 cyanobacterial genomes. This resulted in 30,312 Azolla genome contigs hitting 8,779 different cyanobacterial proteins that were used as a query in a tBLASTn [START_REF] Camacho | BLAST+ : architecture and applications[END_REF] against the Azolla genome; 340 Azolla contigs had reciprocal hits. To investigate whether these represent possible Nostoc-to-Azolla transfers or just examples of plastid-to-nucleus transfers, we used the 340 Azolla contigs for another BLASTx against the cyanobacteria and extracted all 51,743 BLASTx-aligned Azolla sequences. These highly redundant protein sequences were used for a DIAMOND BLASTp against the non-redundant database of NCBI. Almost all of the sequences had streptophyte proteins as the top hit, and when not, phylogenetic analysis clearly placed them within streptophytes.

Phylogeny of SHC.

Homologues of SHC and oxidosqualene cyclase were obtained by searching against Phytozome [START_REF] Goodstein | Phytozome: a comparative platform for green plant genomics[END_REF] , 1KP transcriptomes [START_REF] Matasci | Data access for the 1,000 Plants (1KP) project[END_REF] and the NCBI Genbank non-redundant protein database. Protein alignment was done in MUSCLE [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF] . We used IQ-TREE 116 to find the best-fitting amino acid substitution model and infer the phylogeny using maximum likelihood. Bootstrap support was assessed with 1,000 pseudoreplicates. To test whether the monophyly of fern, lycophyte, moss and liverwort SHC could be rejected, we conducted a Swofford-Olsen-Waddell-Hillis test using SOWHAT [START_REF] Church | Automation and evaluation of the SOWH test with SOWHAT[END_REF] . We compared the best maximum likelihood topology against the topology with all land plant SHC constrained to be monophyletic. SOWHAT was run with 1,000 replicates.

Detection of SHC-synthesized triterpenes.

Freeze-dried S. cucullata biomass was Soxhlet extracted in a 9:1 DCM:MeOH mixture for 24 h. The total lipid extracts obtained were dried over Na 2 SO 4 followed by evaporation of the solvent by a gentle stream of N 2 . Aliquots of the total lipid extracts were methylated with diazomethane to convert the acid groups into corresponding methyl esters, purified over a SiO 2 column and silylated using bis(trimethylsilyl)trifluoracetamide (BSTFA) in pyridine at 60 °C for 20 min to convert the hydroxy groups into the corresponding trimethylsilyl ethers. The total lipid extracts were on-column injected on a Thermo Trace GC Ultra Trace DSQ gas chromatography mass spectrometry (GC-MS) onto a CP-sil 5CB-fused silica column (30 m × 0.32 mm internal diameter, film thickness: 0.10 μ m). The GC-MS was operated at a constant flow of 1.0 ml min -1 . The GC oven was programmed starting at 70 °C to rise to 130 °C at a rate of 20 °C per min and then to 320 °C at a rate of 4 °C per min, followed by an isothermal hold for 20 min.

Reporting Summary. Further information on experimental design is available in the Nature Research Reporting Summary linked to this article. Data availability. The genome assemblies and annotations can be found at www.fernbase.org. The raw genomic and transcriptomic reads generated in this study were deposited in the NCBI SRA under the BioProject PRJNA430527 and PRJNA430459. The sequence alignments and tree files can be found in the Supplementary Data.

Author contributions

Competing interests

The authors declare no competing interests.

Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/ s41477-018-0188-8.

Reprints and permissions information is available at www.nature.com/reprints.

Correspondence and requests for materials should be addressed to F.-W.L.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Articles

NaTure PlaNTS

Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Open Access

Life sciences study design

All studies must disclose on these points even when the disclosure is negative.

Sample size

For phylogenetic analyses, bootstrapping datasets were sampled between 100 to 1000 times, which is the field standard.

Data exclusions Gene models without transcript or homology supports were excluded (see the supplementary discussion).

Replication

The RNA-seq experiments were done with three biological replicates per treatment.

Randomization The plant cultures for RNA-seq were placed on the same growth chamber shelf, but the positions were randomized in terms of nutrient treatments and symbiont types.

Blinding

Blinding is not applicable in this study.

Reporting for specific materials, systems and methods Obtaining unique materials The plant materials used in this study are available upon request (to F.-W. Li or H. Schluepmann)