Should we hail the Red King? Evolutionary consequences of a mutualistic lifestyle in genomes of lichenized ascomycetes

Abstract The Red Queen dynamic is often brought into play for antagonistic relationships. However, the coevolutionary effects of mutualistic interactions, which predict slower evolution for interacting organisms (Red King), have been investigated to a lesser extent. Lichens are a stable, mutualistic relationship of fungi and cyanobacteria and/or algae, which originated several times independently during the evolution of fungi. Therefore, they represent a suitable system to investigate the coevolutionary effect of mutualism on the fungal genome. We measured substitution rates and selective pressure of about 2000 protein‐coding genes (plus the rDNA region) in two different classes of Ascomycota, each consisting of closely related lineages of lichenized and non‐lichenized fungi. Our results show that independent lichenized clades are characterized by significantly slower rates for both synonymous and non‐synonymous substitutions. We hypothesize that this evolutionary pattern is connected to the lichen life cycle (longer generation time of lichenized fungi) rather than a result of different selection strengths, which is described as the main driver for the Red Kind dynamic. This first empirical evidence of slower evolution in lichens provides an important insight on how biotic cooperative interactions are able to shape the evolution of symbiotic organisms.


| INTRODUC TI ON
Coevolution can substantially shape the evolution of organisms involved in intimate ecological connections that range from antagonistic to mutualistic relationships. In its essence, coevolution is a reciprocal evolutionary change induced by interacting species (Thompson, 2014). Possibly, every biotic interaction within the food web involves a certain degree of interdependence resulting in coevolutionary patterns, as any change in a species will influence one or more connected species. If the relationship is tight enough, and the reciprocally induced evolutionary changes last long enough, coevolutionary effects can become apparent. In this context, Van Valen's (1973, 1977 Red Queen dynamic described how biotic interactions can influence evolution. This hypothesis was initially formulated to explain extinction patterns recurring in fossil records of major taxa, but was later extended by the author emphasizing the importance of competitive biotic interactions in a macro-evolutionary framework. In this hypothesis, coevolution was described as an evolutionary action-reaction cycle, which is characterized by the fluctuations of the relative fitness of two antagonist species. This cycle leads to an arms race regulated by natural selection that eventually accelerates evolutionary rates. Many authors further revised the Red Queen dynamic theory (Brockhurst et al., 2014;Strotz et al., 2018) broadening its original meaning (Morran et al., 2011), confirming (Kerfoot & Weider, 2004), or challenging it (Gokhale et al., 2013;Wei & Kennett, 1983); some of these studies used model simulations (Dercole et al., 2010;Rabajante et al., 2015), or experimental systems (Decaestecker et al., 2007;Paterson et al., 2010), at different organizational (e.g., community, population), temporal, and taxonomic scales (Finnegan et al., 2008;Liow et al., 2011).
The incredibly diversified literature inspired by Van Valen's original hypothesis resulted in a wide concept of the Red Queen dynamic that will be used in this study: coevolution as a driving force that can accelerate evolution (Delaye et al., 2018;Pal et al., 2007;Paterson et al., 2010) and/or modify the selective pressure acting on the coevolving species and their genes (Ejsmond & Radwan, 2015). Though abiotic interactions play a prevalent role as a selective constraint at the largest time and spatial scales (Benton, 2009(Benton, , 2010Venditti et al., 2010), biotic interactions can also have a relevant role in stable environments, for long-lasting, specific associations, such as symbioses. Evidence of biotic relationships as an important long-term selective force was found in host-parasite interactions, such as a New Zealand snail and its trematode parasites (Dybdahl & Lively, 1998), a plant-fungus association (Thrall et al., 2012), and a bacteria-ant association (Degnan et al., 2004(Degnan et al., , 2005. In contrast to the accelerated evolution in host-parasite interactions due to the Red Queen dynamic, the so-called Red King dynamic (Bergstrom & Lachmann, 2003) hypothesizes slower evolutionary rate as beneficial for mutualistic interactions in relevant classes of mutualistic interactions (Veller et al., 2017). Although empirical evidence for the Red King dynamic is still lacking, theoretical studies modeled Red Queen/King dynamics, evaluating parameters such as mutation rate, population size, selection strength, and generation time to understand what conditions can favor a slower evolving symbiont in mutualistic symbioses (Damore & Gore, 2011;Gao et al., 2015;Gokhale & Traulsen, 2012;Veller et al., 2017).
Molecular evolutionary rate measurements (e.g., nucleotide substitution rates) have been extensively used to test relevant evolutionary hypotheses involving lifestyles (Bromham et al., 2013), to compare large taxonomic groups (Buschiazzo et al., 2012;De la Torre et al., 2017;Wang et al., 2010), or to identify conditions likely responsible for rate shifts . Mutualistic symbioses have been investigated to test evolutionary hypotheses using substitution rates (Arab et al., 2020;Rubin & Moreau, 2016), but attention has been rarely focused on the lichen symbiosis (Lumbsch et al., 2008;Lutzoni & Pagel, 1997;Zoller & Lutzoni, 2003). These lichen studies used multi-gene datasets, but no study so far addressed differences in substitution rates in lichens using genome-scale data.
The lichen symbiosis is a stable, successful mutualistic association between at least one fungus (the mycobiont) and one or several photosynthetic partners (green algae and/or cyanobacteria: the photobionts). However, the definition of the lichen symbiosis ranged from a controlled parasitism (Ahmadjian, 1993) to mutualism, and it is still subjected to relevant extensions and revisions (Hawksworth & Grube, 2020). These symbioses developed multiple times independently along the evolutionary history of fungi (Schoch et al., 2009); moreover, the mycobiont is-with rare exceptions-an obligate symbiont, whereas the photobiont is usually not entirely dependent on the mycobiont for survival (Nash, 2008;Wedin et al., 2004). For these reasons, lichenized fungi are a suitable system to explore possible genomic consequences of a mutualistic lifestyle.
We are here using genome-scale data to test three specific hypotheses: (i) The evolutionary rate of lichen mycobionts differs from the rates of non-lichenized fungi; (ii) this change in evolutionary rates is due to a different selective pressure acting on mycobionts in comparison with non-lichenized fungi; and (iii) specific genes are under positive selection in a scenario of general slower or faster evolution.

| Taxon sampling
A total of eight lichen-forming fungal species and 11 nonlichenized fungal species were included in this study. Two datasets, corresponding to two independent lichenization events, which occurred in two Ascomycota classes, were prepared. In  Figure S1), and in the polytomy necessary to identify the constrained trees used for rate analyses as unrooted (Figure 1). In both datasets, an equal number of samples in lichenized and nonlichenized comparisons were used, in order to avoid a possible node-density bias (Hugall & Lee, 2007;Venditti et al., 2006). Illumina's NextSeq Platform. High-molecular-weight DNA isolation and long-read sequencing on a Nanopore GridIONx5 sequencer of Astrothelium subdiscretum were done as described before for the lichen fungal culture of Physcia stellaris (Wilken et al., 2020).
In addition to protein-coding genes retrieved by BUSCO, rDNA regions of each genome assembly were extracted. For each genome, F I G U R E 1 Constrained topologies used to run the PAML rate analyses. (a) Dataset A (Dothideomycetes); (b-e) the four topologies used to compare the two lichenized clades to the non-lichenized clade in dataset B (Eurotiomycetes). Lichen clade in green, non-lichenized in gray the assemblies were aligned with blastn to the 18S, ITSRefSeq, and 28S fungal databases (BLAST v2.11.0+; https://www.ncbi.nlm.nih. gov). All identified scaffolds with rDNA regions were then aligned to sequences from the NCBI nucleotide database to delimit the specific rDNA region. For assemblies that only contained a partial rDNA region or without BLAST hits, rDNA regions were reconstructed using the raw reads with GRAbB (Brankovics et al., 2016). Incomplete rDNA sequences from the previous BLAST step or ribosomal markers of the same species from NCBI were used as bait to assemble complete rDNA regions. A further scaffolding step (when needed) and trimming of the poorly aligned positions were then carried out manually. GC content was calculated after the filtering steps both for protein-coding and for rDNA alignments.

| Molecular phylogeny
Nucleotide alignments of genes that were longer than 900 bp (300 codons), after filtering, were concatenated using FASconCAT-G (Kück & Longo, 2014). The concatenated alignment was used to calculate a maximum-likelihood tree with IQ-TREE 2 (Minh et al., 2020) using the GTR+G substitution model. The fast-bootstrapping option implemented in IQ-TREE 2 was used to calculate 1000 bootstrap replicates. The phylogenetic relationship inferred ( Figure S1) was used for the unrooted (pruned, in dataset B; Figure 1) constrained topology in subsequent rate analyses.

| Rates of molecular evolution
Nucleotide substitution rates were measured in baseml or codeml (PAML v4.7e; Yang, 2007) using nucleotide and amino acid alignments longer than 900 bp/300 codons/300 amino acids. The topology was constrained on the base of the precalculated ML tree. In baseml, branch lengths were calculated for each nucleotide alignment separately with parameters: model = 7, Mgene = 0, clock = 0, fix_alpha = 0, Malpha = 0,ncatG = 10, and cleandata = 0. In order to partition the protein-coding genes by codon position, the same alignments used in the previous analysis were converted to the phylip format to exploit the options G (multiple partitions in the alignment) and C (partition by codon position), together with Mgene = 1, which calculates a separate set of branch lengths for each partition (codon position). Ribosomal DNA markers (18S, ITSs with 5.8S and 28S) were analyzed by baseml using nucleotide substitution settings (see above) and partitioning the alignments by locus.
Codon model analyses were performed using the extension for codeml ete3 evol (ETE3; Huerta-Cepas et al., 2016). ω is the ratio between non-synonymous (dN) and synonymous substitution rates (dS). Nested branch models M0 (same ω for the entire tree), b_free (two ω), b_free (three ω), and fb (one ω each tree branch; freeratio model) (Yang & Nielsen, 2002) were run on each alignment to evaluate what model better fits our data and to estimate dN, dS, and ω. In the b_free (two ω) model, the lichenized clade and nonlichenized clade have the same ω, but the outgroup has a different one. In the b_free (three ω) model, a different ω parameter was assigned to the outgroup and to the lichenized and non-lichenized clade, respectively. Highly divergent sequences can easily lead to a saturation of synonymous substitution estimation; therefore, the data from the free-ratio codon model were strictly filtered; only the genes with dS values lower or equal to three (dS ≤ 3) (Yang, 2014) were retained. In addition, all genes with ω > 10 were discarded (filtered dataset), since large ω values are very likely due to assembly or annotation errors that caused dS values to tend toward zero (Rubin & Moreau, 2016).
Branch-site models bsA and bsA1 (null model) were applied (Zhang et al., 2005)  From all PAML output files, long-term evolutionary rates were calculated for each clade (lichenized, non-lichenized) in the trees by averaging branch lengths from the tips to the common ancestor node in a tip-to-root fashion (Barraclough & Savolainen, 2001;Lanfear et al., 2010Lanfear et al., , 2013. This procedure was applied to both nucleotide (using codon or not) and amino acid rate estimations. This method of rate calculation was adopted to avoid the bias introduced when non-independent samples are used in comparative analyses (Felsenstein, 1985). Since the sum of all branches (from the tips, i.e., the present species, to the common ancestor node) represents the same timespan for the considered clades, we did not calibrate the tree to obtain absolute substitution rates.

| Statistical analysis
The distributions of gene rates and ω were compared in Prism 8.3.0 to assess global differences in genome evolutionary rates.
Nonparametric test was selected after testing the normality of the distributions of gene rate with D'Agostino-Pearson and Shapiro-Wilk tests. Therefore, we applied the Wilcoxon matched-pairs signed rank rest (nonparametric equivalent of the paired t test) to compare the distributions of averaged tip-to-root values of lichenized and non-lichenized clades. The distributions of rates were paired by gene.
The four nested branch models and the two nested branch-site models, used to test the presence of different selective pressure and positively selected genes, were compared in pairs by the likelihoodratio test (LRT) using a chi-square distribution.

| RE SULTS
Genome assemblies were generated from Illumina short reads or Nanopore long reads for Astrothelium subdiscretum (Table 1, boldfaced). The total length of the assemblies was between 30 and 40 Mb and in line with the expected genome sizes for filamentous ascomycetes. Genome assembly statistics represented by the contig number, total assembly size, and the N50 value highlight good contiguity and a completeness of 89%-94% evaluated with BUSCO at the class level. The phylum-level universal ortholog percentage was in the range of 94%-98% (Table 2). Using long-read sequencing for A. subdiscretum resulted in a similar genome assembly as the genomes assembled from short reads. Particularly, the short read assembly of P. massariospora outperformed every other assembly (including the long-read assembly of A. subdiscretum) and resulted in 41 contigs and an N50 value of 1,416,161 bp. Since most genomes assemblies were already in sufficient quality with only Illumina sequencing, we refrained from additional Nanopore sequencing for other fungal genomes than A. subdiscretum.
Gene models were extracted from these genome assemblies and were fully supported ( Figure S1), as expected for such a large supermatrix and a low number of tips.
The strict filtering was applied to exclude potentially saturated markers (dS ≤ 3, ω < 10) for the codon free-ratio model. In data-  Figure 2). The complete rate distributions in Figure 2 show higher density for lichenized clades at the median, as samples in these clades are more closely related than the samples in non-lichenized clades; however, the range of the rate distributions is similar. While the majority of genes were slower evolving in the lichenized clades, 12.9%-26.3% of the analyzed genes showed a faster substitution rate (Table 3). In addition to the nucleotide substitution rates, we measured median values of the amino acid replacement rate, which were also significantly lower in the lichenized clades than in the non-lichenized clades (Wilcoxon's test, p < .0001; Table 3, Figure S2). Furthermore, there were significantly slower substitution rates of the lichenized lineages in each codon position (Table 3), with p < .0001 for all the comparisons except one (Wilcoxon's test, p < .01). We also measured codon position rates in the strictly filtered dataset. When rate differences between lichen and non-lichen genes occurred, the strict filtering of genes (about 5%-10% genes survived) determined a generalized decrease in the substitution values (of about 20%-50%), which was expected, as fast-evolving genes (prone to saturation) were excluded. However, the filtering also determined a biased rate proportion between lichenized and non-lichenized

TA B L E 3 (Continued)
rates on each codon position (data not shown), but more apparent on third codon positions, which contributes the most to synonymous substitutions. Therefore, differences for third codon positions in filtered datasets were not always significant (Table 3).
We also measured the nucleotide substitution rate of ribosomal markers (18S, ITSs with 5.8S and 28S), of which most were faster evolving in lichens contrary to our findings in most protein-coding genes. Only the slower evolving 18S gene (0.019 subs/site) and ITS region (0.140 subs/site) rates in Verrucariales were slower when compared to the non-lichenized clade rate (Table 3).
Branch codon models were tested pairwise by LRT (p < .01) (Table S1) to assess which model was able to fit best our data. The most parameter-rich model fb (free-ratio) was the one passing the LRT for the largest fraction of genes, when tested against the 2ω or M0 null models (M0-fb and 2ω-fb in Table S1). However, it provided a better fit only for a smaller fraction of genes when the null model already accounts for different ω parameters between lichenized and non-lichenized lineages (3ω-fb in Table S1). A comparison between the 2ω and 3ω models allowed to reject the null hypothesis for the majority of genes, except for the comparisons performed on dataset B (Pyrenulales) (2ω-3ω in Table S1). For these genes in dataset B, a simpler model using less ω parameters fitted better.
Based on the results of the LRT, we chose the free-ratio branch codon model to calculate dS, dN, and ω and to compare their distributions using the complete and strictly filtered datasets (Figure 3).
dS and dN were significantly lower (Wilcoxon's test, p < .0001) in lichenized clades than in non-lichenized clades (Table 3). In the strictly filtered datasets, the removal of most of faster evolving genes made the difference for synonymous substitutions not significant (Wilcoxon's test, p > .05) except for one of the comparisons in dataset B (p < .01). The removal of extreme dS values strongly influenced the estimation of ω, which is significantly higher for lichens in the complete datasets, and has instead lower median value when the strict filtering is applied (Table 3, Figure 3b,d,f). The filtering approach was applied to the free-ratio codon model, as it is known as sensitive to substitution saturation on the third codon position  (Yang, 2014); it was also used in nucleotide substitution model analysis, when dataset was partitioned by codon position.
The GC content of protein-coding genes and rDNA was lower in lichenized clades than in non-lichenized clades (Table 3)

| DISCUSS ION
In this study, we detected overall slower evolutionary rates in a representative part of the protein-coding genes from the genomes of two distantly related lineages of lichenized fungi, when compared to their non-lichenized sister clades. Since these two lichenized clades evolved from independent lichenization events, these findings provide strong preliminary evidence of a convergence toward slower F I G U R E 3 dN (a, c, e) (substitutions/ site) and ω (b, d, f) distributions before ("complete") and after ("filtered") the strict filtering. Violin plots in a and b correspond to dataset A, c and d correspond to dataset B (Pyrenulales), and e and f correspond to dataset B (Verrucariales). Only one of the two comparisons performed for dataset B is reported (K. petricola, C. Psammophila), the E. sideris, and C. epimyces comparison is reported in Figure S3. Green plots represent lichenized lineages, and gray plots represent non-lichenized lineages; median value is represented by the white dot, the black bar shows the interquartile range, black line shows lower/upper adjacent value, and violin shows the probability density of the distribution evolution, possibly under the influence of the symbiotic lifestyle and its ecological implications.
Fungi are well known for their ability to form diverse associations with photosynthetic organisms, to the point that a generalized, latent capacity of symbiosis between fungal and algal partners was verified for non-symbiotic species (Hom & Murray, 2014). However, it is less clear how this peculiar lifestyle, once in place, can influence the evolution of fungi involved in mutualistic symbioses. Few studies directly investigated the possible consequences of a mutualistic lifestyle on evolutionary rates (Lutzoni & Pagel, 1997;Rubin & Moreau, 2016), highlighting relevant differences in substitution rates for mutualistic lineages. Only one of these investigations has been conducted to verify the possible connection between the switch to a lichenized lifestyle and an evolutionary rate change (Lutzoni & Pagel, 1997), although lichen symbiosis is a successful association, with almost 20% of currently known fungi adopting this lifestyle (Lumbsch & Rikkinen, 2017). Lutzoni and Pagel (1997) detected an increase in evolutionary rates in ribosomal markers for independent events of lichenization and concluded that lichenized fungi could have these elevated evolutionary rates due to higher UV exposure than in non-lichenized relatives with subterraneous vegetative hyphae.
We measured a lower GC content across protein-coding genes and ribosomal regions in lichens, which may indicate a C-to-T mutation bias that could be caused by increased UV radiation (Ikehata & Ono, 2011) due to their exposed lifestyle. We also detected higher substitution rates in ribosomal regions, but as an exception to overall reduced rates in protein-coding genes. This pattern of the ribosomal region not following the trend of the genome was also found for other organisms (Mitterboeck et al., 2016;Su & Hu, 2012), but it remains unclear why the ribosomal regions evolved differently than many protein-coding genes (also when only considering the third codon position).
To determine whether our results based on protein-coding genes can be classified under a broad definition of the Red King dynamic is not straightforward, given the theory is not completely settled on this topic, and also because our results on evolutionary rates do not contain any information about the relative benefits the bionts receive from being in a symbiosis (Bergstrom & Lachmann, 2003).
In addition, Veller et al. (2017) identified several symbiosis classes in which slower evolution could be beneficial (i.e., the Red King).
They assessed the impact of biological parameters such as generation time, mutation rate, selection strength, and population size in population models. The model indicated that mutation rate has a relevant role only for antagonistic symbiosis (Red Queen effect; i.e., faster evolution is more successful), but not for mutualist symbionts.
However, it was also shown that depending if the mutualism has a small or large benefit for the bionts, evolutionary rate parameters such as longer generation time, lower selection strength, and smaller population size can have a short-term and/or long-term advantages for mutualistic symbioses. Some of these described evolutionary rate parameters leading to a Red King effect (i.e., slower evolution is more successful) may be also applicable to the discussion of the result we found for lichenized fungi.
The "universal" protein-coding marker genes (BUSCO genes) used in this study are predictably under strict purifying selection (ω → 0) for both lichenized and non-lichenized lineages. However, the complete dataset identified lichens as having a slightly less strict purifying selection (higher ω) acting on the genes we tested, which can be beneficial in some mutualistic symbiosis (Red Queen). The opposite trend was identified for the filtered dataset, which produced biased ω values, due to the drastically diminished sample size and the exclusion of extreme dS values (mostly present in non-lichenized samples [data not shown]). However, completely different selection strengths (e.g., positive selection), acting on genomic regions other than the one studied, and involved in the establishment or functioning of the lichen symbiosis, cannot be excluded.
Our measurements are limited to two clades of lichenized fungi.
In a general scenario of reduced evolutionary rates, these two clades of lichenized lineages had more genes with (few) sites subjected to positive selection, or neutrally evolving (ω = 1). However, the detection of such sites was consistent only for a small fraction of genes when lichenized clades were compared to a different sister clade. Therefore, the changed sites cannot be attributed with confidence to positive selection or neutral evolution acting on lichens.
Moreover, these models are thought to lack detection power under synonymous substitutions saturation (Gharib & Robinson-Rechavi, 2013), which was the case for the divergent sequences we used.
Lower evolutionary rates in lichenized fungi could be a consequence of the lichen biology and ecology. Lichens are thought to have long generation times, as indirectly confirmed by their generally low growth rates (Armstrong, 1983;Fortuna & Tretiach, 2018) and by direct estimations (Høistad & Gjerde, 2011). Lichen growth can be constrained by the carbon production of a relatively small population of algae (Scheidegger & Goward, 2002). But even in axenic cultures, where nutrient-rich culture media are used, lichenized fungi often exhibit slow growth rates in comparison with many other filamentous ascomycetes with different lifestyles. Slow growth rates and longer generation times can provide a possible explanation for the lower substitution rates that we detected in lichenized clades and could have been contributed to the success and stability of the lichen symbiosis (Red King). Such an association between long generation time and slow evolutionary rate was also highlighted in other organisms (Welch et al., 2008) and in Ascomycota at a subphylum level (Shen et al., 2020). An attempt to assess a relationship between these two characteristics was made by Lanfear et al. (2013) who detected traits in plants that can influence their evolutionary rate.
However, this study used measurements (e.g., plant height), which are unavailable to us for the lichens in this study. Although there are no such data for the species we used in this study, it is reasonable to think that most lichens have long generation times by low growth rates.
Another evolutionary rate parameter that can slow evolution and benefit the Red King effect is the size of a population.
Population size can be estimated from genomic data. However, these population size estimations rely on multiple genomic samples belonging to the same species (or closely related species). The analyses are usually conducted on neutrally evolving sequences, which exclude coding regions (Gronau et al., 2011). Unfortunately, the number of samples included in this study did not allow reliable estimation of the population sizes. Moreover, the non-lichenized groups available for comparisons were rich in lifestyles, such as pathogenic or parasitic lifestyles that could have a strong effect on evolutionary rates. In particular, pathogens are often subject to accelerated evolutionary rates as a result of the Red Queen dynamic (Papkou et al., 2019;Paterson et al., 2010). This limited our selection, and we had, for example, to exclude the recently described order Phaeomoniellales (Chen et al., 2015) from dataset B as it is mostly composed of phytopathogenic and endophytic species. An important aspect of studies on the Red King dynamic is the rate relationship between two bionts in the same symbiosis.
For lichens, we currently lack information about the rates of the corresponding photosynthetic partners. The only experimental data about relative rates in lichens were provided by Zoller and Lutzoni (2003) who verified higher rDNA substitution rates in mycobiont Omphalina, a basidiolichen, when compared to its photobiont Coccomyxa. Since we focused on the genomic evolution of lichenized fungi in this study, we only sequenced mycobiont cultures. Future studies should include the photobionts to allow an investigation of relative rates in the lichen symbiosis.
Despite some limitations, our analyses provided the first evidence of slower evolutionary rates of lichen mycobiont genomes.
This shift in evolutionary rates was often hypothesized for lichens, but never tested. Given the limited sampling this study allowed, further research involving other lichenized lineages, and other symbiotic systems (e.g., mycorrhizae) will be necessary to generalize this possible convergence toward slower evolution. This empirical evidence provides nevertheless important initial insights on how biotic cooperative interactions can shape the evolution of symbiotic organisms.

ACK N OWLED G M ENTS
This work was supported by the Science Innovation Award at the Field Museum. We would like to thank Sabine Huhndorf for the help of subculturing lichen mycobionts and providing suitable biomass for molecular analyses.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.