Azolla ferns testify: seed plants and ferns share a common ancestor for leucoanthocyanidin reductase enzymes.

Questions about in vivo substrates for proanthocyanidin (PA) biosynthesis and condensation have not been resolved and wide gaps in the understanding of transport and biogenesis in 'tannosomes' persist. Here we examined the evolution of PA biosynthesis in ferns not previously reported, asking what PAs are synthesised and how. Chemical and gene-expression analyses were combined to characterise PA biosynthesis, leveraging genome annotation from the floating fern Azolla filiculoides. In vitro assay and phylogenomics of PIP-dehydrogenases served to infer the evolution of leucoanthocyanidin reductase (LAR). Sporophyte-synthesised (epi)catechin polymers, averaging only seven subunits, accumulated to 5.3% in A. filiculoides, and 8% in A. pinnata biomass dry weight. Consistently, a LAR active in vitro was highly expressed in A. filiculoides. LAR, and paralogous fern WLAR-enzymes with differing substrate binding sites, represent an evolutionary innovation of the common ancestor of fern and seed plants. The specific ecological niche of Azolla ferns, a floating plant-microbe mat massively fixing CO2 and N2 , shaped their metabolism in which PA biosynthesis predominates and employs novel fern LAR enzymes. Characterisation of in vivo substrates of these LAR, will help to shed light on the recently assigned and surprising dual catalysis of LAR from seed plants.

Introduction

Azolla is a genus of floating ferns known for the Eocene Azolla event when floating fern mats covered parts of the North Pole Sea and adjacent fresh water bodies: the mats captured substantial amounts of CO2 as they sunk to the seafloor, presumably contributing to climate cooling (Brinkhuis et al., 2006; Neville et al., 2019). The mats are deemed undesirable by ecologists as they cause anoxia in the water column and outcompete plants that cannot rely on atmospheric dinitrogen (N2) in freshwater ecosystems (Wagner, 1997). There is renewed interest in Azolla fern mats for high‐yield protein‐rich biomass production in the complete absence of nitrogen fertiliser and, possibly, to build soils in subsiding delta regions (Brouwer et al., 2018; Brouwer et al., 2019). Azolla mats, moreover, represent a unique model in which ferns construct a floating ecological niche that includes several symbiotic bacteria (Dijkhuizen et al., 2018). The most prominent symbiont is the obligate filamentous cyanobacterium Nostoc azollae, with an eroded genome that fixes sufficient N2 for the symbiosis, even at the highest growth rates (Ran et al., 2010; Brouwer et al., 2017).

Niche construction was shown to influence the evolution of proanthocyanidin (PA; Fig. 1) production by plant–soil feedback (Arnoldi et al., 2020). The PA content of Azolla filiculoides and Azolla pinnata was determined by colorimetric assay at 4.4% and a staggering 19.9% of the DW, respectively (Brouwer et al., 2019). PA, therefore, could similarly define the niche of fern mats that contain little to no lignin (Nierop et al., 2011). Fern PA biosynthesis pathways, however, have not been documented. This is unfortunate as PAs constitute a key step in evolution and have not been reported in the other, more ancient, plant lineages such as the bryophytes and the lycophytes (Markham, 1988; Alonso‐Amelot et al., 2004).

Fig. 1

General structure of flavonoids. Flavonoids are made of the A, B, and C rings and key features of their structures are hydroxyl groups at the numbered carbon atoms within their structure (a). Flavan‐3‐ols can polymerise into proanthocyanidins, with the C4–C8 (green) being the most widespread carbon–carbon bond (b).

Symbiotic cyanobacteria in Azolla could, in principle, contribute to PA biosynthesis because cyanobacteria contain the flavonoid biosynthesis pathway (FBP) that, in angiosperms, provides the substrates for PA biosynthesis: cyanobacteria have recently been found to synthesise flavonols such as dihydroquercetin and quercetin (Żyszka‐Haberecht et al., 2018). By contrast, dihydroflavonol 4‐reductase (DFR) may have evolved in the plant lineage (Campanella et al., 2014). Enzymes shared between PA and anthocyanidin biosynthesis, including DFR and anthocyanidin synthase (ANS), have already evolved in mosses and spike mosses, and have been identified already in ferns but not tested (Campanella et al., 2014).

Evidence for a role for PA and biosynthesis intermediates in the Azolla symbiosis interaction is sketchy. The flavanone naringenin, the precursor of most flavonoid classes including PA, was isolated from the symbiosis and shown to stimulate growth of Nostoc punctiforme, a free‐living facultative plant symbiont related to N. azollae; its rhamnoglycoside, naringin, is known to induce inhibition of hormogonia (Cohen et al., 2002; Żyszka et al., 2017). In addition, simple trichomes of the leaf cavities from Azolla ferns, known to mediate the interactions between symbiotic bacteria and fern host, contained PA (Pereira & Carrapiço, 2007). To begin understanding the role of PA in Azolla symbioses, chemical characterisation and knowledge of the biosynthesis pathway are required.

PA can range from soluble oligomers to large, insoluble polymers with a degree of polymerisation (DP) of over 50 flavan‐3‐ol subunits; PAs of flavan‐3‐ol subunits contain both or either catechin and epicatechin (Khanbabaee & van Ree, 2001; Dixon et al., 2004). The flavanol subunits with their characteristic aromatic A‐ and B‐rings (Fig. 1), and defining hydroxyl groups, are produced by enzymes of the FBP (Yonekura‐Sakakibara et al., 2019), but enzymes synthesising the recently discovered cysteinyl‐linked substrate have not been identified (Liu et al., 2016). Moreover, enzymes catalysing subunit polymerisation into PA have not been identified; polymerisation involving flavonoid carbocations without enzymes occurs at neutral pH unlike the acidic pH typical of the vacuole (Liu et al., 2016; Wang et al., 2020). PAs, however, are generally trapped in vacuoles, thus preventing any interaction with cytosolic proteins (Zhao et al., 2010; Brillouet et al., 2013). Uncertainty about the in vivo substrates in PA biosynthesis persists.

Recent advances have revealed a dual role of leucoanthocyanidin reductase (LAR) from an array of seed plants beyond the synthesis of (+)‐catechin from leucocyanidin: the enzymes also accept 4β‐(S‐cysteinyl)‐epicatechin (Cys‐EP) which they cleave to Cys and (−)‐epicatechin, the polymer starter subunit (Yu et al., 2019). Cys‐EC has been detected in several plant extracts and both genetic and biochemical evidence supports a role for LAR in determining DP (Liu et al., 2016, 2018). Characterising FBP and PA polymerisation in the more ancient lineage of ferns, specifically those accumulating large amounts of PA, may help to close gaps in our current understanding of PA biosynthesis. The improved transcriptome assemblies and orthogroup classifications available from the 1000 plants transcriptomes initiative (1KP), published genome assemblies annotations and transcript profiles of ferns including A. filiculoides render the task possible (Brouwer et al., 2017; Li et al., 2018; Wong et al., 2019).

This work, therefore, integrates PA chemical analysis with genomics to reconstruct the biosynthetic pathway in a fern. Using A. filiculoides and A. pinnata, we first determined flavonoid intermediates in biomass, and subunits in purified PA, using liquid chromatography coupled to mass spectrometry (LC‐MS), and thermally assisted hydrolysis and methylation combined with gas chromatography and mass spectrometry (THM‐GC‐MS). Secondly, we reconstructed a pathway for PA biosynthesis in A. filiculoides with the knowledge from the expression of isoenzymes in diel transcript expression profiles, leveraging information from the first Azolla genome annotation (Li et al., 2018). We finally focused on the gene of secondary metabolism most prominently expressed in A. filiculoides investigating its phylogenetic relationship to known LAR enzymes and the activity of its recombinant protein.

Materials and Methods

Chemicals

Phenolic standards (chlorogenic acid, catechin, naringin, tannic acid and quercetin) and substrates for the enzyme assay (NADPH and dihydroquercetin) were purchased from Sigma‐Aldrich. Trace element analysis grade ethyl acetate, dichloromethane and methanol were purchased from VWR International (Amsterdam, the Netherlands), anhydrous acetone from Merck and acetonitrile from Biosolve.

Plant materials and growth conditions

Azolla filiculoides originated from a ditch in Utrecht (the Netherlands; 52.075°N 5.149°E). Azolla pinnata was provided by the International Rice Research Institute and originated from Sigiriya in Sri Lanka (IRRI code PI 535 in Watanabe et al., 1992). The plants were grown in 16 h light with an intensity of 100 µmol m−2 s−1, at 21°C and 70% humidity. The growth medium contained 0.65 mM NaH2PO4·H2O, 1.02 mM K2SO4, 1 mM CaCl2·2H2O, 1.65 mM MgSO4·7H2O, 17.9 µM Fe‐EDTA, 9.1 µM MnCl2·4 H2O, 1.6 µM Na2MoO4·2H2O, 18.4 µM H3BO3, 0.8 µM ZnSO4·7H2O, 0.3 µM CuSO4·5H2O, 0.2 µM CoCl2·6H2O. The plants were kept in the linear growth phase by continuous harvesting (Supporting Information Fig. S1) after which the biomass was freeze dried immediately. Biomass of both species was collected over a period of 5 wk and pooled. Under these growth conditions, the sporophytes were entirely green.

Purification of proanthocyanidins

Freeze‐dried biomass was first extracted in a Soxhlet extractor with dichloromethane : methanol, 9 : 1 for lipid removal. The residual pellet was air dried and extracted three times with aqueous acetone (70% v/v) shaking at 200 rpm for 24 h and pooled, then the mixture was centrifuged at 5000  and the supernatant was collected from three replicate extractions. The acetone was evaporated in a rotary evaporator. The remaining aqueous extract was extracted with ethyl acetate yielding the ethyl acetate fraction (EAF). The remaining aqueous extract was then mixed with an equal amount of methanol and loaded onto a column with Sephadex LH‐20 resin and the flow‐through and subsequent aqueous methanol washes were collected as the flow‐through (FT) fraction. The column was washed with aqueous methanol (50% v/v) until the FT was colourless. Finally, elution was performed with aqueous acetone (70% v/v) yielding the SF fraction.

Determination of dry weight, C and N‐content

Freeze‐dried biomass was weighed as such, while the pellets and the phenolic fractions were first dried for 24 h at 105°C and cooled for 1 h in a desiccator for dry weight determination. Here, 1–2.5 mg of dried material was used for total carbon and nitrogen content determination using an elemental analyser (Fisons NA 1500 CNS).

Colorimetric assay to estimate total phenol content

Freeze‐dried biomass, the dried pellets and phenolic fractions were extracted with aqueous acetone (70% v/v) for 24 h shaking at 200 rpm. Folin–Ciocalteu reagent was added to the extract to determine the phenol content, as described by Waterman & Mole (1994). Tannic acid was used as a standard.

Colorimetric assay to estimate proanthocyanidin content

Proanthocyanidin content was estimated using the acid butanol assay: 1.2 ml of butanol‐HCl solution (5% v/v) and 20 μl of 2% NH4Fe(SO4)2 in 2 N HCl was added to the freeze‐dried biomass, the dried pellets and phenolic fractions followed by vigorous mixing and incubation for 50 min at 95°C. The absorbance was measured at 550 nm. Proanthocyanidins from the black spruce Picea mariana were used as a standard (Nierop et al., 2005) because this PA is also very rich in (epi)catechin units and, therefore, resembles the reactivity of Azolla PA towards the assay.

Estimate of average degree of polymerisation

A vanillin assay as described by Butler et al. (1982) was performed on purified proanthocyanidins: the dried sample was dissolved in either methanol or glacial acetic acid, both containing 4% concentrated HCl and 0.5% vanillin. Samples were incubated for 20 min at 30°C, after which the absorbance was measured at 500 nm using a spectrophotometer. The average DP was estimated using the ratio of absorbance in methanol over glacial acetic acid.

Extraction of phenolic compounds in methanol and their analysis using liquid chromatography/mass spectrometry (LC‐MS)

Methanolic extracts for LC‐MS analysis were obtained by extracting freeze‐dried biomass with aqueous methanol (75% v/v) acidified with formic acid (0.1% v/v) as described by De Vos et al., 2007. Methanolic extracts and fractions obtained from proanthocyanidin purification were analysed by LC‐MS as described by van der Hooft et al. (2012). Briefly, separation of extracted compounds using a Waters Acquity (U)HPLC was over a Luna C18 column (150 mm × 2.0 mm; 3 μm particles; Phenomenex) with a water–acetonitrile gradient in 0.1% v/v formic acid. Detection was by both a Waters Photodiode Array (200–600nm) and a LTQ Orbitrap FT‐MS hybrid mass spectrometer (Thermo Fisher Scientific, Nieuwegein, the Netherlands) with negative electrospray ionisation, a mass resolution of 70 000 FWHM and a mass range of m/z 100–1300.

Thermally assisted hydrolysis and methylation combined with gas chromatography/mass spectrometry (THM‐GC‐MS) analysis of phenolic compounds

The method was as described by Nierop et al. (2005). Briefly, freeze‐dried biomass was pressed onto Curie‐Point wires together with a drop of aqueous tetramethylammonium hydroxide (25% v/v) and heated for 5 s at 600°C in a Curie‐Point Pyrolyzer (Horizon Instruments, Heathfield, UK). Separation of the products was performed on‐line using gas chromatography (Carlo Erba GC8060) on a fused silica column coated with CP‐Sil 5 (Varian; 25 m × 0.32 mm; 0.40 μm). The carrier gas was helium. Oven temperature was programmed as 1 min at 40°C, heating to 180°C at a rate of 7°C min−1 and a final heating to 320°C at a rate of 20°C min−1 which then was held for 15 min. Detection was using the Fisons MD800 mass spectrometer (70 eV electron ionisation, mass range m/z 45‐650, cycle time 0.7 s). The reference standards are listed in Fig. S2.

Identification of candidate genes of phenylpropanoid pathway (PBP) and FBP

The work was initiated with transcript assembly accessions from Brouwer et al., 2017 and then revised using the first assembly annotation of the A. filiculoides Galgenwaard genome (Li et al., 2018). Proteins predicted were accessed at ftp://ftp.fernbase.org/Azolla_filiculoides/Azolla_asm_v1.1/. Fasta files of both the high‐ and low‐confidence gene models were used for local Blast (https://www.ncbi.nlm.nih.gov/books/NBK279690/) with sequences of enzymes in the literature experimentally validated. In addition, predictions were used from the Kegg database and Mercator tool (https://plabipd.de/portal/mercator‐sequence‐annotation; Lohse et al., 2014; Kanehisa et al., 2016). For each enzyme of the FBP, proteins identified were then aligned using Mafft software (https://mafft.cbrc.jp/alignment/software/, Katoh & Standley, 2013) to experimentally validated enzymes, GBlocks was used to remove poorly aligned regions (Talavera & Castresana, 2007) and phylogenetic trees were constructed using PhyML (Guindon et al., 2010; Katoh & Standley, 2013). Identified candidate genes were named according to the most closely related known genes of the PBP and FBP.

Phylogenetic analyses of the LAR‐like genes

Using the Vitis vinifera LAR accession number (GSVIVT01011958001), a precomputed ‘LAR‐orthogroup’ was extracted from the 1KP database v.2 consisting of 6500 protein sequences ( http://jlmwiki.plantbio.uga.edu/onekp/v2/, Wong et al., 2019), then re‐sampled to include all sequences from key model plants representing land‐plant lineages, Azolla and Salvinia as well as a few other ferns and known PIP enzymes, then Azolla LAR and WLAR sequences were added. The final, resulting 808 protein sequences were aligned using Mafft‐linsi (Katoh & Standley, 2013). The alignment was filtered limiting column and sequence gaps using trimAll (Capella‐Gutiérrez et al., 2009), leaving 584 unique sequences with 305 parsimony informative sites. From the trimmed alignment, a maximum likelihood phylogenetic tree was computed using Iq‐Tree (Nguyen et al., 2015; Trifinopoulos et al., 2016) using extended model fitting (best‐fit model LG + R7; Kalyaanamoorthy et al., 2017) and both ultra‐fast bootstraps and SH‐aLRT support (Hoang et al., 2018). Trees were visualised in iTol (Letunic & Bork, 2019). All details concerning the phylogenetic analyses are available on the GitHub page: github.com/lauralwd/LAR_phylogeny_gungor‐et‐al‐2020.

Analysis of the RNA‐sequencing data to reveal transcript accumulation in ferns and Arabidopsis

Azolla RNA‐sequencing raw data were for A. filiculoides Galgenwaard grown without or with nitrogen (2 mM NH4NO3) in the medium and sampled in three biological replicates at 2, 8 and 14 h into the 16 h light period, then at 4 h in the 8 h darkness period of the diel cycle (Brouwer et al., 2017). Arabidopsis RNA‐sequencing raw data were from 3‐wk‐old, unstressed Arabidopsis thaliana plants collected at four time points spanning the diel cycle (Coolen et al., 2016). Quality, trimmed, paired‐end reads were aligned to the gene models of each plant and only paired matches aligning to only one gene model were counted (CLC Genomics Workbench), as described by Brouwer et al. (2017). We realise that for isoenzymes with high similarities some reads were lost as multimappers, leading to a possible underestimation of the absolute read counts for some of these isoenzymes, but this handicap does not change interpretation of differential expression when comparing differing time points. Read counts were normalised to total mapped reads per sample as reads per million (RPM).

Cloning, expression and purification of recombinant AzfiLAR is detailed in Methods S1.

Assay of the leucoanthocyanidin reductase

The enzyme reaction contained 5 mM DTT, 0.5 mM dihydroquercetin, 0.5 mM NADPH, 5 μM VvDFR and 5 μM AzfiLAR in a total volume of 100 μl with the buffer at pH 6.5, as described by Maugé et al., 2010. The reaction was performed for 2 h at 30°C and stopped by extraction with ethyl acetate and drying overnight. The pellet was stored at −20°C and dissolved in aqueous methanol (50% v/v) before LC‐MS analysis, essentially as described by Maugé et al. (2010), except that separation was on a C18 column (Varian; 150 × 3.0 mm; 3 μm particles) with an ultra‐high performance (UHP) LC instrument coupled to a 6560 Ion Mobility Q‐TOF mass spectrometer (Agilent Technologies, Amstelveen, the Netherlands).

Accession numbers

Novel DNA sequences were deposited at EBI's ENA Browser under accession number PRJEB39515.

Results

Ferns from either sections of the Azolla genus share a similar phenolic profile

To draw a more general conclusion about PA production in the Azolla genus, we first investigated which phenolic compounds could be detected in two Azolla species: A. filiculoides (section Azolla) thriving in temperate regions and A. pinnata (section Rhizosperma) thriving in (sub)tropical regions (Heckman, 1997; Metzgar et al., 2007).

LC‐MS profiles of aqueous methanol extracts from both species were dominated by phenylpropanoid‐type caffeoylquinic acids (e.g. peaks 2, 12, 14 and 16; Fig. 2; Table 1). The profiles further contained small amounts of the flavonol glycosides quercetin‐3‐O‐rutinoside, quercetin‐3‐O‐glucoside and kaempferol‐3‐O‐rutinoside (peaks 8–10; Fig. 2; Table 1); compared with A. filiculoides these were more abundant in A. pinnata, much like in A. imbricata that is also a species from the Rhizosperma section (Qian et al., 2020). Azolla pinnata furthermore contained more (epi)catechin trimer (peak 7; Fig. 2; Table 1). By contrast, THM‐GC‐MS of the whole biomass yielded a large peak of 1,3,5‐trimethoxybenzene characteristic of flavonoids and PA in both A. filiculoides and A. pinnata (peak 8; Fig. S3; Table S1), and confirming the prevalence of PA reported earlier in Azolla ferns (Nierop et al., 2011; Brouwer et al., 2019). In addition, THM‐GC‐MS profiles of whole biomass contained the THM products typical of caffeoylquinic acids in both species (peaks 5 and 12; Fig. S3; Table S1).

Fig. 2

Representative total ion chromatograms obtained during liquid chromatography‐mass spectrometry (LC‐MS) analyses of methanolic extracts from Azolla filiculoides and Azolla pinnata biomass. Base peak intensity (ion counts per second; y‐axis) of compounds eluting from the column at indicated time after injection (min; x‐axis). Peak numbers correspond to the compounds listed in Table 1. Upper, extract from A. filiculoides; lower, extract from A. pinnata.

Table 1

Overview of the (poly)phenolic compounds, their retention time (RT) and their m/z as ([M‐H]‐), that could be assigned to main peaks detected in the LC‐MS chromatograms of the Azolla filiculoides and Azolla pinnata extracts analysed in this study (Figs 2, 4).

Peak number	RT (min)	[M‐H]−	Compound
1	8.6	353.0874	1‐Caffeoylquinic acid
2	12.1	353.0875	5‐Caffeoylquinic acid
3	12.9	353.0875	3‐Caffeoylquinic acid
4	13.6	577.1346	(Epi)catechin dimer
5	15.1	353.0875	4‐Caffeoylquinic acid
6	15.3	289.0715	Epicatechin monomer
7	16.6	865.1876	(Epi)catechin trimer
8	22.1	609.1455	Quercetin‐3‐O‐rutinoside
9	23.2	463.0877	Quercetin‐3‐O‐glucoside
10	25.2	593.1504	Kaempferol‐3‐O‐rutinoside
11	25.7	515.1192	Dicaffeoylquinic acid I
12	26.3	515.1192	Dicaffeoylquinic acid II
13	27.1	515.1192	Dicaffeoylquinic acid III
14	28.5	515.1192	Dicaffeoylquinic acid IV
15	32.6	515.1192	Dicaffeoylquinic acid V
16	37.4	677.1503	Tricaffeoylquinic acid I
17	38.8	677.1503	Tricaffeoylquinic acid II
18	39.5	677.1503	Tricaffeoylquinic acid III

Fig. 4

Liquid chromatography‐mass spectrometry (LC‐MS) analyses of the three fractions obtained from polyphenol extractions for Azolla filiculoides and Azolla pinnata. Representative chromatograms of the ethyl acetate (EAF), the Sephadex LH‐20 flow‐through (FT) and the Sephadex LH‐20 eluted (SF) fractions for A. filiculoides (upper) and A. pinnata (lower). Peak numbers correspond to the compounds listed in Table 1. The large series of peaks in the SF of both species could be assigned to (epi)catechin polymers, thus proanthocyanidin (PA).

Phenolic profiles are similar amongst species of the Azolla genus, we conclude. In A. pinnata, a member of the Rhizosperma section, flavonol glycosides, tricaffeyolquinic acids and (epi)catechin make up a higher proportion of the LC‐MS profile than in A. filiculoides. PAs were not in the range m/z 100–1300 detected by the LC‐MS profiles of whole‐fern extracts. The strong flavonoid signals observed by THM‐GC‐MS in both species, likely therefore originated from PA with a molecular weight > 1300. To research this hypothesis, we set out to purify PAs from the polyphenol extracts and study these in more detail.

Azolla PA consists mainly of (epi)catechin subunits with an average degree of polymerisation of about seven

To purify PA, A. filiculoides and A. pinnata whole biomass was extracted according to the scheme in Fig. 3(a); this resulted in three fractions: the EAF, the flow‐through of Sephadex LH‐20 including washes (FT), and the eluate from Sephadex LH‐20 with 70% aqueous acetone (SF). Firstly, the dry weight and the nitrogen and carbon content was determined for each fraction. The FT contained most of the initial biomass, while the SF had the highest carbon content (Fig. 3b), close to that of pure PA from the black spruce Picea mariana (59–66%) (Nierop et al., 2005).

Fig. 3

Polyphenol extraction from whole Azolla biomass and key features of the fractions obtained from it. (a) Overview of the extraction process, fractions obtained included the lipid extracted residue X, the (poly)phenol extracted residue Y, the ethyl acetate fraction (EAF), the flow‐through (FT) and the fraction bound to Sephadex LH‐20 then eluted with 70% acetone (SF). (b) Chemical analysis of whole Azolla filiculoides and Azolla pinnata biomass and the fractions obtained from (a). Total phenolic and proanthocyanidin (PA) content are respectively in percentage tannic acid and pure black spruce PA equivalents. Portion stands for the percentage of the initial biomass in each fraction. Values represent means of three biological replicates with SD. na, not available.

Secondly, total phenols and PA content of each fraction were determined using colorimetric assays and expressed in equivalents of tannic acid (polygalloyl esters, hydrolyzable tannins) and black spruce tannin (poly‐flavan‐3‐ols, proanthocyanidins), respectively. In line with the carbon content, the SF were highly enriched in phenols and especially in PA (Fig. 3b). The EAF were also highly enriched in phenols, with much lower levels of PA except for the EAF from A. pinnata that contained a significant amount of PA (Fig. 3b). The FT contained relatively low levels of phenols and PA.

Subsequently, the content of each fraction was profiled by both LC‐MS and THM‐GC‐MS. LC‐MS showed that the EAF and the FT contained mostly caffeoylquinic acids (Fig. 4; Table 1), also present in the original methanol extract (Fig. 2). THM‐GC‐MS of both fractions revealed THM products characteristic of caffeoylquinic acids (Fig. S4; Table S1). The caffeoylquinic acids (peaks 16, 17 and 18; Fig. 4; Table 1) were enriched in EAF and FT fractions compared with the original methanol extract (Fig. 3; Table 1). The (epi)catechin trimer (peak 7; Fig. 4; Table 1) in the EAF of A. pinnata corresponded to the PA content of this fraction and, therefore, suggested that a proportion of A. pinnata PA was extracted with ethyl acetate (Fig. 3b).

The SF of both Azolla species were massively enriched in PA, yielding a large series of PA peaks that were eventually restricted by the LC‐MS detection limit of m/z 1300 (Fig. 4; Table 1). The THM products derived from flavonoid A‐ and B‐rings, especially the highly abundant 1,3,5‐trimethoxybenzene and 1‐methyl‐2,4,6‐trimethoxybenzene (peaks 8 and 10, respectively; Fig. S4; Table S1), were most likely to have been derived from the flavan‐3‐ol subunits of the PA (Fig. S4; Table S1). Considering in addition the carbon‐to‐nitrogen ratios and the PA contents determined by colorimetric assay in these fractions, we concluded that the SF of both species contained highly purified PA. The SF furthermore contained large proportions, that is 81% and 86%, of PA found in the initial biomass of A. filiculoides and A. pinnata respectively (Fig. 3b). Finally, the mass balance and total phenolic content of the fractions, indicated that PAs were the most abundant, soluble, phenolic compounds in these Azolla species.

Most flavonoid B‐rings detected in the SF of both species with THM‐GC‐MS had two hydroxyl groups (Fig. S4; Table S1). We therefore concluded that both A. filiculoides and A. pinnata PAs mainly consisted of catechin and epicatechin subunits. An additional vanillin assay on the SF of both species indicated an average DP of about 7 (data not shown). Although A. filiculoides and A. pinnata differed slightly in their phenol content, the characteristics of their PAs were the same. We therefore assumed that PA biosynthesis pathway was most likely to be shared amongst species of the Azolla genus and focused on the genome of A. filiculoides for subsequent identification of enzymes from the FBP.

Azolla filiculoides invests more in transcripts of the FBP than Arabidopsis thaliana

The first genome annotation of A. filiculoides (Li et al., 2018) served to identify genes involved in PA production. Automatic functional annotation of each candidate gene was systematically verified or corrected with previously identified genes of the phenylpropanoid and the FBP from A. thaliana and Vitis vinifera (FBP including LAR). Transcript profiles at four time points during the diel cycle of A. filiculoides sporophytes (Brouwer et al., 2017) were then compared with those from 3‐wk‐old and nonstressed control A. thaliana plants (Coolen et al., 2016).

Firstly, seven PAL‐like, two C4H‐like and eight 4CL‐like genes were identified (Table S2). All these genes had a higher read count, normalised per million reads and averaged over four time points in the diel cycle, in the ferns compared with Arabidopsis plants and, therefore, a higher transcript accumulation in A. filiculoides compared with A. thaliana (Fig. 5). Secondly, four CHS‐like, one CHI‐like, two F3H‐like, eight F3’H‐like, two DFR‐like, one FNSII‐like, one LAR‐like and one FLS/ANS‐like genes were identified (Table S2). The latter having a shared phylogeny with both FLS and ANS. Most of these genes had again a higher transcript accumulation in ferns compared with A. thaliana, exceptions were FNSII‐like and FLS/ANS‐like (Fig. 5).

Fig. 5

Transcript accumulation of enzymes involved in the phenylpropanoid and the flavonoid biosynthesis pathways in Azolla filiculoides sporophytes and Arabidopsis thaliana 3‐wk‐old plants. The average transcript accumulation, normalised as reads per million (RPM), over four time points during the day was summed up for each isoenzyme of a catalytic step to reflect the transcriptional investment of each plant in each conversion. 4CL, 4‐coumarate‐CoA ligase; ANR, anthocyanidin reductase; ANS, anthocyanidin synthase; C4H, cinnamic acid 4‐hydroxylase; CHI, chalcone isomerase; CHS, chalcone synthase; DFR, dihydrokaempferol 4‐reductase; F3H, flavanone 3‐hydroxylase; F3′H, flavonoid 3′‐hydroxylase; F3′5′H, flavonoid 3′5′‐hydroxylase; FLS, flavonol synthase; FNS, flavone synthase; IFS, isoflavone synthase; LAR, leucoanthocyanidin reductase; PAL, phenylalanine ammonia‐lyase.

A. filiculoides, we concluded, generally accumulates more transcripts of the phenylpropanoid pathway and the FBP compared with A. thaliana. Markedly, A. filiculoides LAR‐like was the 44th most accumulating transcript of all genes annotated and one of the highest accumulating transcripts of genes encoding enzymes involved in secondary metabolism. Given the role of LAR enzymes in PA biosynthesis in seed plants, the observation was consistent with the large accumulation of PA in A. filiculoides measured in Fig. 3. We next wondered whether the expression of LAR‐like correlated with that of other genes encoding enzymes for PA biosynthesis.

Candidate genes involved in PA production in Azolla filiculoides are co‐regulated

As PAs are so abundant in Azolla, enzymes for their biosynthesis may dominate the FBP. We therefore analysed the transcript accumulation of PBP and FBP candidate genes at four different time points during the diel cycle in the RNA‐sequencing data following Brouwer et al., 2017 and using the genome annotation from Li et al., 2018. Given that many isoenzymes were found for each of the conversions of both the PBP and FBP, and given that we did not know a priori to which final product they would contribute, we added up the read counts of the isoenzymes to estimate transcriptional investment over the diel cycle in each enzymatic step. The CHS‐like, CHI‐like, DFR‐like and LAR‐like genes shared similar transcript accumulation patterns in ferns under different conditions: transcripts peaked at 08:00, 2 h into the 16 h light period, but did not peak when nitrogen was supplied (Fig. 6a ). Generally, transcripts of only one or two isoenzyme genes dominated the read counts for every enzymatic step (Table S2). The shared expression pattern observed when summing up read counts indicated co‐regulation of enzymes for PA biosynthesis that dominated the transcriptional investment in secondary metabolism pathways in the ferns. R2R3MYB transcription factors from class VIII‐E that generally control FBP and PBP in seed plants were present in A. filiculoides (Jiang & Rao, 2020), their read counts were low compared with those from class VIII‐D (Fig. 6b,c). Class VIII‐D R2R3MYB from the clade containing AtNOECK(MYB106) were the most highly expressed (Fig. 6d, Azfi_s0001.g000839 and Azfi_s0016.g014344). Correlated changes in read counts and expressed potential R2R3MYB regulators strengthened the hypothesis that A. filiculoides LAR‐like (AzfiLAR‐like from this point forward) may be involved in PA production. We therefore investigated this gene in more detail.

Fig. 6

Changes in transcripts of proanthocyanidin (PA) biosynthesis enzymes and potential R2R3 MYB‐regulators during the diel cycle. Raw read data from (Brouwer et al., 2017) were mapped to gene models predicted in the Azolla filiculoides genome assembly v.1 (Azfi_1; Li et al., 2018). Data were available from sporophytes grown without (−N) and with nitrogen (2 mM NH4NO3, +N), collected in triplicates at 2 h, 8 h, 14 h and 20 h with the light period of 16 h commencing at 6 h. Transcript accumulation is given as reads per million reads mapping (RPM; Supporting Information Table S2). Changes in transcript accumulation of PA biosynthesis enzymes are relative to 8 h time point in sporophytes without or with nitrogen (a). Transcript accumulation of R2R3‐MYB assigned to the VIII classes D (b) and E (c) by Jiang & Rao (2020) are averages from three biological replicates, with SD.

AzfiLAR‐like phylogeny places it at the base of other characterised LARs and its recombinant proteins synthesises catechin in vitro

LAR belongs to the PIP protein family, including the enzymes pinoresinol‐lariciresinol reductase (PLR), isoflavone reductase (IFR) and phenylcoumaran benzylic ether reductase (PCBER) (Maugé et al., 2010). A phylogenetic tree with AzfiLAR‐like and other well characterised members of the PIP protein family placed AzfiLAR‐like at the base of the branch splitting into gymnosperm and angiosperm LAR (Fig. S5; Table S3). AzfiLAR‐like protein had only 42.6% amino acid identity with V. vinifera LAR (VvLAR). 3D‐modelling of the proteins revealed that their structures largely overlapped, with minor changes in the active site (data not shown). We therefore proceeded with testing whether recombinant AzfiLAR‐like protein was able to convert leucoanthocyanidins into trans‐flavan‐3‐ols in vitro.

AzfiLAR‐like was cloned, and the protein produced via recombinant protein expression in E. coli and purified (Fig. S6). Enzymatic assay was performed in a double reaction with VvDFR because leucoanthocyanidins are highly unstable and not commercially available: VvDFR can produce leucoanthocyanidins from dihydroflavonols such as dihydroquercetin (DHQ) (Petit et al., 2007). Products of the enzymatic assay were analysed using LC‐MS. Incubation of DHQ with only AzfiLAR‐like did not yield any new peak (Fig. 7a), while catechin was produced when both VvDFR and AzfiLAR‐like were present (Fig. 7b). The control incubation of DHQ, with only VvDFR, did not result in the production of catechin (data not shown). Catechin production could therefore only be ascribed to LAR activity by AzfiLAR‐like (AzfiLAR from this point forward).

Fig. 7

Leucoanthocyanidin reductase activity of recombinant AzfiLAR. (a) Dihydroquercetin (DHQ) was incubated with only recombinant AzfiLAR in the presence of NADPH. (b) DHQ was incubated with both recombinant AzfiLAR and Vitis vinifera DFR (VvDFR) in the presence of NADPH. Representative chromatograms of the catechin and DHQ, respectively. Main ion m/z 291.08 and m/z 305.06 are shown.

AzfiLAR orthologues are found in many ferns and LAR and PCBER evolved from a common ancestor in the last common ancestor of seed plants and ferns

We next wondered whether LAR enzymes were found in further ferns and other lineages of seed‐free land plants. The 1KP orthogroup database v.2 (Wong et al., 2019) contained an orthogroup of the Vitis vinifera LAR (NP_001267887) which we used for phylogenetic analyses. The orthogroup sequences were subsampled by including all transcripts for each species chosen, and concentrating on species with known genomes and PIP proteins. We monitored changes in trees computed after including firstly more algae species and secondly more species from lycophytes and bryophytes. Typically, green algae IFR/PLR‐like enzymes clustered and we, therefore, rooted the phylogenetic tree as shown in Fig. 8 on this algal clade. The orthogroup contained sequences of all known PIP‐family dehydrogenases. We observed that the PLR clade, well supported by bootstrapping, evolved early during land‐plant evolution and that the IGS&EGS enzymes evolved at least twice (Fig. 8).

Fig. 8

Phylogeny of LAR and other PIP‐family enzymes across land‐plant lineages. Protein sequences were retrieved as a single orthogroup containing the Vitis vinifera LAR from the 1KP orthogroup database_2 (Wong et al., 2019), subsampled and supplemented with guide and Azolla filiculoides LAR‐like sequences. The 808 sequences were aligned with Mafft‐linsi (Katoh et al., 2019), then trimmed using trimAL (Capella‐Gutiérrez et al., 2009) and, using iQ‐Tree (Nguyen et al., 2015), the phylogenetic tree computed with a resulting 305 parsimony informative sites, a best‐fit substitution model LG + R7 and 1000 bootstrap determined by SH‐aLRT (Kalyaanamoorthy et al., 2017). The tree was annotated in iTol (Letunic & Bork, 2019; details of all sequences are given in Supporting Information Fig. S8). Nodes with bootstrap support equal or > 80% SH‐aLRT are indicated by circles. EGS and IGS clustered in two groups: IGS/EGS‐I and IGS/EGS‐II. Abbreviations: EGS, eugenol synthase; IFR, isoflavone reductase; IGS, isoeugenol synthase; LAR, leucoanthocyanidin reductase; PCBER, phenylcoumaran benzylic ether reductase; PLR, pinoresinol‐lariciresinol reductase; WLAR, fern‐specific LAR‐like.

The clade containing the PCBER/IFR/group_2_IGS&EGS and LAR enzymes was separated from the PLR clade with over 93.8% aLRT bootstrap support and did not contain sequences from lycophytes, bryophytes and algae (Fig. 8). The well supported PCBER/IFR/group_2_IGS&EGS clade may have originated from LAR‐like enzymes in the common ancestor of seed plants and ferns: a cluster of branches with many fern sequences, including AzfiLAR and its orthologue in A. caroliniana, was found at the base of the well supported seed‐plant LAR clade. The AzfiLAR clade was well supported with 99.9% aLRT bootstrap and contained sequences from Adiantum tenereum, Ceraptorpteris thalictoides as well as many other ferns, indicating that LAR activity is more widespread in ferns. Additional fern LAR‐like genes were reproducibly found in two separate clades in all trees we computed; we named these clades WannabeLAR_1 (WLAR1) and WannabeLAR_2 (WLAR2) (Fig. 8). Inspection of the alignment of the proteins encoded in A. caroliniana transcripts 2016335 and 2016336, AzfiLAR and VvLAR, showed that the proteins differed in key amino acids of the active site compared with the known LAR enzymes but not in the NADPH binding site. WLAR1 and WLAR2 enzymes were found in A. filiculoides and other ferns (Fig. S7) and are, therefore, not an artefact of 1KP transcript assembly. Dipteris conjugata, for example, had transcripts from LAR, WLAR1 and WLAR2 clades (Fig. S7).

Taken together, phylogenetic analysis of the 1KP orthogroup revealed that LAR enzymes are more generally found in ferns and that they evolved in the common ancestor of ferns and seed plants. In addition ferns contain paralogues of the LAR‐genes of unknown function.

Discussion

Soluble procyanidins dominated the secondary metabolite content of the biomass from the sporophyte symbiosis. Candidate genes encoding enzymes for PA biosynthesis found in A. filiculoides were upregulated in sporophytes and included a LAR, but FLS/ANS‐like was not. Recombinant AzfiLAR was active in vitro and AzfiLAR orthologues were commonly found in fern species but were not detected in other seed‐free plant lineages. Additional WLAR were found with differing substrate binding sites from LAR.

The procyanidin type PA accumulating in Azolla ferns have an optimum degree of polymerisation for protein aggregation and possibly, therefore, a role as antifeedant

PA could be purified with high yield and characterised from both A. filiculoides and A. pinnata, and their structures were similar: procyanidins with an average DP estimated at seven (Figs 3b, 4; Table 1). Assay of a large diversity of PA for protein aggregation revealed that the DP is a key determinant for the capacity of procyanidin type PA to aggregate proteins in solution, with stiff globular proteins, such as bovine serum albumin, requiring more PA than the more flexible and less compacted gelatin (Ropiak et al., 2017). Moreover, the smallest PA with largest effect on protein aggregation had a DP of seven which happens to also be the average DP of Azolla PA. We concluded, therefore, that Azolla PA may have evolved to optimally aggregate proteins and thus may have an ecological role as antifeedant. Consistently, PA formed precipitating PA–protein aggregates during protein extraction from Azolla biomass which reduced the protein yields from extractions (Leizerovich et al., 1988; Zeller et al., 2015; Brouwer et al., 2019). This also explained why Azolla biomass with about 20% DW crude protein may only be included in animal diets in limited proportions without affecting growth (Leterme et al., 2009; Acharya et al., 2015; Brouwer et al., 2019).

Staining with vanillin–hydrochloric acid confirmed localisation of PA in vacuoles of specialised trichomes located inside the leaf cavity where N. azollae fix nitrogen (Pereira & Carrapiço, 2007). The PAs therefore may have additional roles in the interaction with symbionts and as redox‐active compounds important in the ‘after‐life’ processes in the Azolla mats. More research is needed, manipulating levels of individual compounds in the pathways, to evaluate physiological and ecological roles of PA and PA biosynthesis intermediates in Azolla. Such manipulations may best be done by specifically altering key biosynthetic enzymes in vivo.

The common ancestor of ferns and angiosperms already contained all enzymes of the angiosperm PA biosynthesis pathway

Genes of all enzymes of PA biosynthesis characterised from angiosperms were found in the genome of A. filiculoides except for F3′5′H, in line with the lack of flavan‐3‐ols with three hydroxyl groups attached to the flavonoid B‐ring, and ANR (Fig. 6). Presence of either F3′5′H, or ANR in Azolla cannot be completely excluded, however, due to the incomplete nature of our first fern assembly and annotation reported in Li et al., 2018. Results are consistent with our chemical analyses, and earlier examples that detected flavone and flavonol derivatives in various Azolla species (Teixeira et al., 2001; Qian et al., 2020). Azolla was shown, in addition, to synthesise red‐coloured 3‐deoxyanthocyanins correlating with snail and tadpole deterrence (Cohen et al., 2002a, 2002); the fern turned completely red upon 3‐deoxyanthocyanin accumulation, for example when exposed to cadmium that led to CHS and DFR‐like transcript accumulation in Azolla imbricata (Dai et al., 2012). Steps leading to final 3‐deoxyanthocyanin formation remain uncharacterised, final steps might be followed by DFR‐like reduction or ANS‐like dehydrogenation/dehydroxylation reactions, and subsequent glycosylation (Liu et al., 2010).

The phenylpropanoid, flavonol and lignin biosynthesis pathways are widely distributed in seed‐free plant lineages (Tohge et al., 2013; De Vries et al., 2017). Specifically, the products of PLR, the lignans, have been described in bryophytes and lycophytes, as well as other vascularised plants (Koeduka, 2018; Markulin et al., 2019) consistent with results in Fig. 8. By contrast, PA appear to be absent in bryophytes and lycophytes (Davies et al., 2020). LAR, therefore, could have evolved from a PLR.

Recent analyses aimed at evaluating the production of antioxidants in cyanobacteria, revealed the accumulation of catechin and epicatechin in several Nostoc and Anabaena species, with catechin present in all filamentous cyanobacteria tested (Blagojević et al., 2018). Moreover, a DFR was annotated from a N. punctiforme strain associated with lichen (Tabuchi Yagui., NCBI‐accession RCJ41608). In our hands, searches for DFR identified the N. punctiforme gene reproducibly, but returned none from the N. azollae genome, searches for LAR were fruitless. We concluded, therefore, that enzymes until DFR were likely to have already evolved in cyanobacteria, but that N. azollae had lost the capability.

Interestingly, synthesis of some of the flavonols was repressed when cyanobacteria were grown under conditions in which they needed to fix nitrogen (Blagojević et al., 2018), suggesting that they may be important for the N2‐fixation process, possibly as antioxidants. Synthesis of flavonols for N. azollae that are locked into nitrogen fixation in the leaf cavities, may thus have been taken over by the fern host in the Azolla symbiosis. Some fungi synthesise leucoanthocyanidins, suggesting that these metabolites are not restricted to cyanobacteria and plant lineages (Bu et al., 2020). A thorough analysis of cyanobacterial FBP enzymes will reveal whether enzymes such as DFR appeared in the plant, fungal and cyanobacterial lineages by convergent evolution or horizontal gene transfer.

AzfiLAR evolved in the common ancestor of ferns and seed plants, but do fern LAR also control the DP of PA?

Although ferns were known to accumulate PA, neither LAR nor ANR had been identified in this lineage before (Alonso‐Amelot et al., 2004; Nierop et al., 2011). The AzfiLAR recombinant protein exhibited LAR activity in vitro and phylogenetic analyses placed it at the base of the LAR leaf containing characterised LAR from seed plants (Figs 8, S5). Apparent absence of PA in bryophytes and lycophytes, and our phylogenetic analyses suggested that LAR could be a key evolutionary innovation in the last common ancestor (LCA) of ferns and seed plants. This was suggested in spite of the fact that enzyme activity of PIP proteins from a single clade may not be predicted with certainty (Koeduka, 2018; Markulin et al., 2019).

Additionally, two paralogous clades of the AzfiLAR containing WLAR1 and WLAR2 were identified in ferns. Given their differing substrate binding site from genuine LAR and PCBER (Fig. S7), more work is needed to understand the function of these novel fern‐specific reductases. These enzymes are particularly interesting, however, because seed‐plant LAR enzymes catalyse two different reactions, possibly catalysed by distinct enzymes in lineages that evolved earlier in land‐plant evolution. VvLAR and MtLAR convert Cys‐EP into epicatechin as well as leucoanthocyanidins into trans‐flavan‐3‐ols (Liu et al., 2016; Yu et al., 2019).

In vitro, epicatechin produced by LAR from Cys‐EP is required as starter unit to which epicatechin carbocations condense sequentially, bringing about polymerisation. The LAR regulate the DP of PA, therefore, by influencing the amount of starter. AzfiLAR may similarly encode an enzyme cleaving Cys‐EP. Consistent with the hypothesis is the low DP of PA and the absence of a highly expressed candidate gene for ANR in Azolla. We concluded that the high AzfiLAR expression in A. filiculoides may maintain the DP optimum for antifeedant function and will need to test this hypothesis in the near future. Sadly, we have yet to establish mutagenesis or genome editing protocols for Azolla to verify the roles of specific enzymes and regulators of the PA pathway in vivo.

Coordinate expression of the FBP enzymes is also found in seed plants and suggests a conserved mode of regulation

The coordinate expression observed during the diel cycle of isoenzymes of the FBP (Fig. 6) would be consistent with a regulating transcription factor complex controlling anthocyanin accumulation, as it exists in dicotyledonous plants (Albert et al., 2014; Albert, 2015). Light regulation of anthocyanin biosynthesis in A. thaliana is well known and mediated by regulating the stability of the activating R2R3‐MYB transcription factors PAP1 and PAP2 that are part of the MYB‐bHLH‐WD40 (MBW) complexes more generally regulating differentiation of trichomes, as well as expression of genes in secondary metabolism (Bulgakov et al., 2017). Class VIII‐E R2R3MYB were, however, much less expressed than class VIII‐D, the MIXTA1‐type known to specify glandular trichome development (Shi et al., 2018; Jiang & Rao, 2020). Glandular trichomes are abundant and critical to attract N. azollae in the specialised structures of the Azolla fern host (Peters & Perkins, 1993).

That transcripts of the FBP enzymes begin accumulating in the dark when nitrogen is supplied, or soon after the start of the light period when N. azollae uses light energy to fixate N2 is consistent with observations in cyanobacteria in which flavonol biosynthesis is upregulated when nitrogen is available (Brouwer et al., 2017; Blagojević et al., 2018). Identification of the MYB and bHLH transcription factors that control the FBP pathway will be of particular use to engineer ferns with altered levels of PA in order to study the specific relation of PA characteristics and quantities with symbiont development, ecological fitness and biomass digestibility. Such studies have begun successfully in symbiotic seed‐plant species, such as species of clover (Liu et al., 2018).

Author contributions

PB grew plants, purified tannins and analysed the fractions with the help of KN, RdV and IvdM. PB, EG and HS analysed genes of the condensed tannin biosynthesis pathways, DCS subcloned AzfiLAR, EG expressed and tested recombinant AzfiLAR, helped by JST for MS analyses of the products. LWD, EG and HS computed PIP protein phylogenies. EG, PB and HS wrote the manuscript with input from all other authors; all authors agreed to the list of authors and the identified contributions of those authors.

Fig. S1 Photographs of typical Azolla sporophytes and their continuous culturing setup.

Fig. S2 Partial THM‐GC‐MS chromatograms of the main (poly)phenolic standards.

Fig. S3 Partial THM‐GC‐MS chromatograms of whole A. filiculoides and A. pinnata biomass.

Fig. S4 THM‐GC‐MS of fractions from polyphenol extractions for A. filiculoides and A. pinnata.

Fig. S5 Phylogenetic relationship of A. filiculoides LAR and characterised PIP enzymes.

Fig. S6 SDS‐PAGE of the proteins used in the enzymatic assay of Fig. 7.

Fig. S7 Alignments of LAR, PCBER and related sequences (WLAR) from ferns and V. vinifera.

Fig. S8 Phylogeny of LAR and other PIP‐family enzymes across land‐plant lineages.

Methods S1 Cloning, expression and purification of recombinant AzfiLAR.

Table S1 Compounds detected with THM‐GC‐MS in Azolla biomass during this study.

Table S2 Transcript accumulation of phenylpropanoid/flavonoid isoenzymes and relevant R2R3‐MYB transcription factors in the A. filiculoides genome.

Table S3 PIP enzymes used to generate the phylogenetic tree in Fig. S4.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Click here for additional data file.