Fungal-Algal Association Drives Lichen Mutualistic Symbiosis: A Case Study With Trebouxia-Related Lichens

1. Introduction

Lichens are important components of Fungi, accounting for about 20% of known fungal species [1]. Lichenization events have happened independently several times in the evolutionary history of Ascomycota and Basidiomycota [2,3]. The peculiar mutualistic symbiosis of lichens, composed of fungi, green algae, or cyanobacteria [4,5], helps them predominantly exist in various extreme environments such as polar, plateau, and desert regions [6], playing a crucial role in the early formation of terrestrial ecosystems [7].

Up to now, there are nearly 20,000 accepted lichen-forming fungi (LFF) species [1], while lichen-forming algae (LFA) species are no more than 200 [8]. The composition of LFF and LFA in lichen symbiosis is not symmetrical [9], however, the association between LFF and LFA is far more complex than the ratio of species number. Both biotic and abiotic factors including cophylogeny of LFF and LFA, geography, and reproductive strategy of LFF, influence the formation of fungal-algal pairings in the lichen symbiosis.

Various degrees of availability, selectivity, and specificity between LFF and LFA are the most frequent descriptions for driving the lichen associations. Specificity is defined as the narrow taxonomic range of interacting partners, whereas selectivity refers to the frequency of association with distinct partners [10]. Generally, compared with LFA-LFF association with higher selectivity, higher specificity of LFF-LFA association has been discovered at any taxonomic level. For example, at the order level, the LFF Arthoniales, Ostropales, Pyrenulales, and Trypetheliales, mainly distributed in the tropics, are preferentially associated with green algae Trentepohliales, and Lecanorales and Teloschistales are mainly associated with green algae Trebouxiales [11]. Most LFF genera and species choose one green alga or cyanobacterial genus and one specific species as their LFA, respectively [12,13,14], except for some scarce examples such as species in Lobaria, Peltigera, and Pseudocyphellaria of Peltigerales, which can associate with two different LFA Domains like cyanobacteria and green algae [14,15]. Until now, almost all researches on symbiont association patterns were performed using a specific lichen taxon as a model, in which genus level and family level associations are most common [8,12,15,16,17]. It is poorly known the association pattern between LFF and LFA when distantly related fungi are considered.

Different lichen taxa usually show different responses to diverse geography on the LFF-LFA association. Some LFF are not affected by the geographical factors that show strong specificity to the LFA across remarkable global ranges, such as the associations between the genus Oropogon and limited Trebouxia OTUs, Omphalina and Coccomyxa, and Evernia mesomorpha and Trebouxia jamesii s.l. [12,18,19]. While in some lichens, geographic distributions lead to different patterns of symbioses, such as Stereocaulon and Sticta [16,17,20]. The ecological preference of LFA would indirectly affect lichen distribution, different LFF associate with the same LFA generally exhibiting similar ecological properties [21,22]. Asexual reproduction is proceeded by vegetative propagules, such as isidia, soredia, and lobules in lichens, by which lichens colonize new habitats by co-dispersing LFF and LFA (vertical transmission). Sexual reproduction depends on the spores of LFF, and new symbioses form until LFF meets the suitable LFA (horizontal transmission) [23]. Asexual reproduction would be apt to structure stable lichen mutualistic associations but possibly result in low LFA diversity. In contrast, the propagation distance would be farther, and high LFA diversity occurs through sexual reproduction [21]. Nevertheless, these studies are either still not considering distantly related fungi or only limited to a certain geographical zone, it is not clear the association pattern between LFF and LFA when multi types of geographical zones as well as different LFF orders are considered simultaneously.

Among the known lichen symbionts, 50 % –70% are symbiotic with green algae Trebouxiaceae, of which the LFA genus Trebouxia is the most common [8]. There is an ongoing interest in understanding the diversity and evolutionary relationships in Trebouxia-associated lichens [12]. Recently, a total of 109 to 113 candidate species lineages were determined, a four-fold increase from the 29 formally described Trebouxia species, thereby underestimating Trebouxia diversity [5]. Nonetheless, very few among the nearly 1,600 samples in these studies were collected from China, encompassing a cryptic or unknown diversity of Trebouxia considering its extensive ecogeographical range. Therefore, it provides an impetus for intensively studying the diversity of Trebouxia-associated lichens in China to understand their evolutionary pattern. Furthermore, Trebouxia-associated lichens occur in different fungal orders and across different climatic zones.

In this study, we chose Trebouxia-associated lichen symbionts as a model, considering different taxonomic levels (including order) and distinct biomes (mainly from China) at the same time. A total of 200 samples were included here, among which we found 41 species of lichenized fungi, associating with 16 species of trebouxioid green algae, of which 62% previously unreported. We demonstrated that species identity of the fungal and algal partner had the greatest effect on the outcome of the symbiosis. Partner specificity was higher in tropical regions. Co-phylogenetic analyses indicated congruent phylogenies of fungi and algae, in which loss’ and ‘failure to diverge’ of the algal partners were observed to contribute most to the coevolutionary mechanisms. Overall, we demonstrated a biotic factor (partner identity) is the more important driver of symbiosis of Trebouxia spp. and associated fungal partners over a broader taxonomic and ecological context, which broadens our knowledge of fungi-algae symbiotic patterns.

2. Results

2.1. Molecular sequence data and phylogenetic analysis

We generated 1,238 DNA sequences for this study, including 196 ITS, 188 nuLSU, 176 mtSSU, and 136 RPB1 sequences from LFF and 186 ITS, 170 nuSSU, and 186 rbcL sequences from LFA (Table S1). The BLASTn of ITS sequences against the GenBank database indicated the LFF was grouped into 41 species, 29 genera, 14 families, eight orders, and two classes in Ascomycota, of which 40 species belonging to Lecanoromycetes, and one species to Candelariomycetes, by the high identity of ITS sequences with the described species. The species delimitation was confirmed with their phylogenetic topology using the ITS dataset (Figure S1) and the four-locus concatenated phylogeny (Figure S2).

All the LFA belong to Trebouxia based on their ITS sequence analysis using the BLASTn tool. The phylogenetic topology of Trebouxia was constructed based on the combined ITS and rbcL sequences by incorporating sequences from the study by Muggia, Nelsen, Kirika, Barreno, Beck, Lindgren, Lumbsch and Leavitt [5] (Figure S3). Three DNA species delimitation analyses (GMYC, bPTP, and ABGD) were used to confirm the species boundaries (Figure S4). Sixteen species were recovered, and ten were undescribed (Figures S3 and S5).

2.2. Mutualistic association pattern between LFF and LFA

Fifty-one associations were observed between 41 LFF and 16 LFA (Figure 1; Table S1). At the class level, the LFF species in Lecanoromycetes were associated with all 16 LFA species. At the order level, the LFF orders, including samples of more than two species, were associated with multiple LFA species, that is, the species of Lecanorales, Caliciales, Teloschistales, and Pertusariales, were associated with ten, six, three, and two LFA species, respectively. Similarly, at the family level, the LFF species of Parmeliaceae, Megasporaceae, Teloschistaceae, and Physciaceae, including samples of more than two species, were associated with six, four, three, and two LFA species, respectively. At the genus level, Circinaria was the only genus that included samples of more than two species, which was associated with two LFA species. Each LFF species was associated with one to four LFA species at the species level. Diplotomma alboatrum had the most potent selectivity but the poorest specificity to LFA, with the association ratio being 1 vs. 4.

Each Trebouxia species was compatible with one to 11 LFF species, and 11 of 16 Trebouxia species were associated with more than one LFF (Figure 1). Five LFA species, viz. Trebouxia arboricola, T. cretacea, T. impressa, T. jamesii, and Trebouxia sp. 10 were associated with four to ten LFF species across multiple families, genera, and species, among which the first four LFA species even were associated with multiple species distributed into different orders of LFF. Trebouxia cretacea had the most robust selectivity but the poorest specificity to LFF, which accepted 11 LFF species across two classes, seven orders, seven families, and eight genera.

Four pairs of associations between LFF and LFA have unique specificity no matter the source of the samples, i.e., Hypogymnia hypotrypa vs. Trebouxia sp. 5, Menegazzia subsimilis vs. T. sp. 4, Parmotrema tinctorum vs. T. corticola, and Psoroma fruticulosum vs. T. sp. 6. Trebouxia sp. 8 was found to symbiosis with two LFF species, namely Heterodermia japonica and H. speciose; however, 95% samples (21 of 22) were H. japonica. Most of these specific associations happened in tropic/subtropic zones, indicating that the LFF-LFA association tends to be more specific in warmer climates. The mutualistic associations corresponding to the three climate zones are summarized in Table 1. Interaction network analyses further revealed that the tropic and subtropic networks had lower web connectance (0.15), linkage density (2.13), and nestedness (22.92) and were segregated into more compartments (6) compared with the cold and temperate networks.

2.3. Correlation among LFF, LFA, and ecological factors

We performed variation partitioning analysis to assess the relative contributions of climatic variables, geographical distance, LFF or LFA, distribution types, reproductive modes, and host specificity to shape the mutualistic association of lichens (Figure 2). For the LFA, climatic variables, geographical distance, LFF, and the other factors explained 36% of the variation (Figure 2a). The LFF explained the most significant proportion of the variation with 12% independent effect and 14% in combination with other variables. Distribution types, reproduction modes, and host specificity independently explained 3% of the variation, whereas 16% was shared with other variables. Geographical distance and climatic variables independently explained 1% and 2% of the variation, and 5% and 10% shared with other variables, respectively. On the other hand, climatic variables, geographical distance, LFA, and other factors explained 28% of the variation of the LFF (Figure 2b). The LFA explained the most significant proportion of variation with 11% independent effect and 12% in combination with other variables. The climatic variables independently explained 3% of the total variation and 7% shared with other variables, followed by 2% and 9% of variation explained by distribution types, reproduction modes, and host specificity independently and in combination with other variables, respectively. The geographical distance explained 3% when combined with other variables while zero independently.

2.4. Cophylogenetic analyses

The cophylogenetic analyses were performed using the well-supported LFF and LFA phylogenies. The global-fit tests (ParaFit and PACo) supported overall congruence between tree topologies of fungal hosts and Trebouxia species (ParaFit Global = 0.08, P = 0.001; PACo m² global value = 2.73, P = 0, Figure S6), and rejected the null hypothesis of random association. Sixteen out of the 51 mutualistic associations (31.4%) contributed significantly to Parafit Global based on ParaFit1 values and 18 associations (35.3%) based on ParaFit2 values (P ≤ 0.05).

Most cost scenarios tested in Jane 4.0 supported significant congruence between LFF and LFA phylogenies. All random solutions were worse than the solution reconstructed by the program, except for costs H and Q (Table 2). Among the significant reconstructions, cost regime B yielded the lowest overall cost that suggested one cospeciation, eight duplications, six duplications with host-switches, 90 losses, and 35 failures to diverge events leading to the congruent phylogenies between LFF and LFA (Figure S7). The roles of losses and failure to diverge were vital and all the cost regimes showed the same solution regardless of the penalized values.

3. Discussion

Accurate species delimitation of symbiotic partners is a crucial prerequisite. Otherwise, a biased identity will result in masked interactions [15]. We performed phylogenetic analysis using the dataset in Muggia, Nelsen, Kirika, Barreno, Beck, Lindgren, Lumbsch and Leavitt [5] as the backbone to construct the phylogeny of Trebouxia species. Besides, we used three species delimitation approaches to define the algal species boundary. The results corresponded well with those of Muggia, Nelsen, Kirika, Barreno, Beck, Lindgren, Lumbsch and Leavitt [5] (Figure S3). Sixteen Trebouxia species were delimitated, of which ten (62.5%) are potentially new species to science, which supports the cryptic Trebouxia species diversity has been dramatically underestimated. Some inconsistent delimited results occurred that would make the potential Trebouxia species diversity higher. For example, within the Heterodermia japonica/H. speciosa clade, four or six Trebouxia species were recognized by GMYC (ITS), bPTP (ITS), ABGD (ITS + rbcL), and bPTP (ITS + rbcL), whereas only one species was adopted by GMYC (ITS + rbcL), ABGD (ITS), and Muggia, Nelsen, Kirika, Barreno, Beck, Lindgren, Lumbsch and Leavitt [5] (Figures S2, S4 and S5). We chose the last delimitation considering some samples in this clade growing in the same locality, even on the same bark, such as QJ2-1, QJ2-2, and QJ2-8. In a clade composed of 18 samples representing eight lichen species (Bryoplaca tetraspora, Ochrolechia tartarea, Gondwania regalis, Megaspora verrucosa, Lecanora fuscobrunnea, Physcia caesia, Physcia dubia, and Rhizoplaca chrysoleuca), all the delimitation methods accepted only one Trebouxia species except for GMYC (ITS) supporting two (Figures S4,A5). Hence, only one Trebouxia species was accepted according to the rule of the majority.

Our data showed that 32 (78%) of the 41 LFF species were associated with only one LFA species. In comparison, 11 (69%) of 16 LFA species accepted more than one LFF species (Figure 1). These asymmetric mutualistic association patterns indicated higher selectivity and weaker specificity from the perspective of LFA choosing LFF over the reverse. Previous interaction network studies have hypothesized that the asymmetry may result from the uneven distribution of fungal hosts and algal species, called an abundance-asymmetry hypothesis [24,25]. This hypothesis was under further development that species abundance would determine interspecific interaction frequency and strength and the resulting asymmetry structure if individuals in a community interact randomly [26]. The known species abundance differs about one hundred times between LFF and LFA (Kroken & Taylor, 2000; Lücking et al., 2017), which inevitably impacts the mutualistic association patterns. The reproduction modes also affect the asymmetric mutualistic association patterns. LFF have both sexual and asexual reproduction modes [27]. Macrolichens, like most species included in our studies, generally develop asexual structures like isidia and soredia (Bowler & Rundel, 1975), by which LFF disperse together with their LFA. Asexual reproduction was thought to be the only propagation strategy in Trebouxia-containing lichens [28]. Therefore, LFF tends to act as a specialist in mutualistic association. However, for LFA, Trebouxia species are seldom found in a free-living state and even cannot be alive when lichens die; therefore, they tend to act as generalists in the mutualistic association chosen by LFF, especially in a harsh environment [29,30].

Furthermore, in our study, different symbiotic partners showed different association patterns across different climatic zones and ecoregions. The network analyses supported that boreal/arctic and temperate zones had higher web connectance, linkage density, and nestedness than the tropic/subtropic zones, indicating higher reciprocal specificity between the symbiotic partners in the tropics, which is consistent with the finding of Singh, et al. [31]. More than 56% of mutualistic associations are specifically one-to-one in the tropics compared to 10%–12% in boreal/arctic and temperate zones. There is a higher selectivity of LFA in arctic/temperate regions than tropic/subtropic regions as it associates with eight to ten LFF in the former regions but at most five LFF in the latter regions, suggesting frequent mycobiont switches happened especially in harsh environments such as Antarctic [10,32]. On the contrary, the patterns of LFF specificities toward LFA remain unchanged across different ecoregions and climatic zones, in which most LFF maintain relationships with one LFA except for a very few cases, such as Circinaria tortuosa, Diplotomma alboatrum, Gondwania regalis, Lobothallia semisterilis, Ochrolechia tartarea, Parmotrema clavuliferum, and Ramalina calicaris associating two or three LFA species. Among all the mutualistic associations, the most apparent modularity can be well shown by Hypogymnia hypotrypa because both LFF and LFA acted as one-to-one specialists. This high reciprocal specificity is most possibly related to their vertical transmissions via soredia and the limited niche of LFA, which was also found in Evernia mesomorpha, belonging to the same family as Hypogymnia hypotrypa [13,19].

The symbiotic partners have the most significant impact that drives mutualism, with 11%–12% independent effects (12%–14% in combination with other variables). The contribution of LFF and LFA to their corresponding partners is based on the genetic distance of their phylogenies. Therefore, the effects of the symbiotic partners could be understood as their cophylogenetic trend; that is, the effect value here is roughly equivalent to the degree of cophylogenetic evolution. Variation partitioning analyses showed that the geographical location (GPS) and 19 climate factors related to the temperature and precipitation contributed no more than 3%. Similar results were found in many other symbiosis analyses, in which the symbiotic partners mainly drove the mutualism, whereas ecological mechanisms and geographical distributions such as extreme environments or cooler climate zones partly explained the variation [12,16,23].

There are still 64% and 72% unknown interpretations that need further exploration. One of these unknown factors is historical climatological and geological processes that have turned out to affect the current species distribution patterns [33]. Whether or to what extent this factor brings the effect is not known in this study. Besides, the research model in this study is Trebouxia-associated lichens that distribute across different genera, families, and orders, which behave in much more complex interactions than the studies focusing on species or genus level. A much higher explanation of the variation within one fungal genus could be easily obtained [16,31].

Cophylogeny evaluates the dependency of two associated groups of organisms to illuminate how ecological and evolutionary processes affect species diversification [34]. Our study rejected the null hypothesis of a random association between LFF and LFA, as previous studies reported [5,17,35,36], indicating the significant congruence between Trebouxia algae and its associated LFF. Loss and failure of LFA to diverge played a much more critical role than cospeciation, duplication, and host switch driving the coevolution in lichen mutualism (Figure S7). It has been verified that cospeciation between symbionts seldom happened except for those with strictly vertical transmission instead of the reciprocal selection by mutualistic partners such as lichens [37]. Host-switches, failure of the LFA to diverge, and losses are commonly found to be the prevalent events shaping lichen symbionts, such as Protoparmelia-Trebouxia, Sticta-Symbiochloris, and Cladonia-Trebouxia [17,31,38]. Failure of the LFA to diverge with its fungal hosts is known to be responsible for the occurrence of the generalist algal species [17,39], as shown by one Trebouxia species associated with multiple fungal species across the different taxonomic levels in our study. Losses of LFA are a consequence of extinction or incomplete lineage sorting that has been reported as another significant event during the coevolutionary process of lichens [31]. It is usually caused by the inability of the algae to parasite the fungal hosts when the host speciates incipiently with small population size [37,39]. Furthermore, this study's evolutionary events leading to Trebouxia-associated symbionts do not correspond to the climate zones. In contrast, macroclimate may influence unique lichen association patterns, such as Protoparmelia-Trebouxia symbionts driven by different evolutionary mechanisms in different climatic regions [31].

Our study highlights the LFF-LFA association in Trebouxia-related lichens. The sampling scope is not restricted to species, genus, family or order but across a broader taxonomic scale covering a comprehensive framework of lichen phylogeny. In addition, our study involves different climatic zones and ecological types, including local, regional to global scales other than a specific geographical scale. Although the sampling size of this study is not very large that cannot present the whole picture of the Trebouxia-related lichens, this study gave us an overview of the LFF-LFA mutualistic associations on larger taxonomic and ecological scales. More climate zones, ecological types, taxa, and sample sizes should be considered in the future exploration of mechanisms underlying the coevolution of mutualistic associations.

4. Materials and Methods

4.1. Taxon sampling

Two hundred Trebouxia lichen samples were collected from a wide range of geographical scales from the mainland (forest and sand desert) and Taiwan (forest) of China, Japan (forest), Israel (sand desert), and Antarctic (cold desert) across three climate zones (temperate, tropic/subtropic, and cold) in 2014–2018 (Figure 3, Table S1). The samples represented a wide range of taxonomical scales across different orders, families, genera, and species, mostly belonging to the class Lecanoromycetes.

4.2. DNA extraction, amplification, and sequencing

Total genomic DNA was extracted from lichen thalli using a modified CTAB method [40]. For LFA, we amplified the nuclear internal transcribed spacer region (ITS), the nuclear small subunit ribosomal RNA gene (SSU), and the chloroplast ribulose-bisphosphate carboxylase-RuBisCO gene (rbcL). For LFF, four nuclear markers were amplified, i.e., ITS, the large nuclear subunit ribosomal RNA gene (nuLSU), the largest subunit of RNA polymerase II (RPB1), and the small mitochondrial subunit ribosomal RNA gene (mtSSU). PCR amplification was performed in 25 μl reaction volume consisting of 12.5 μl of 2 × Taq Mix DNA polymerase (CWBIO, China), 9.5 μl of ddH₂O, 1 μl each (10 μM) of the primer set (forward and reverse), and 1 μl DNA template. Primers and PCR conditions are listed in Table S2. PCR products were examined on 0.8% agarose gels, purified using a gel purification kit (Shanghai Huashun Bioengineering Corporation, China) following the manufacturer’s instructions, and sequenced using ABI 3730 XL Sequencer by Shanghai BioSune Corporation of China.

4.3. DNA alignment and phylogenetic analysis

DNA sequence contigs were assembled and manually edited using the SeqMan module of the DNASTAR Lasergene Packages 7 (Madison, WI, USA). Sequences of each locus were aligned separately with MAFFT v 7.471 [41] using the G-INS-i algorithm. Gaps were treated as missing data, and ambiguously aligned regions were excluded using Gblocks v0.91b [42]. These sequence data have been submitted to the GenBank databases under accession number in Table S1. Sequence alignments were deposited at TreeBase (http://purl.org/phylo/treebase/phylows/study/TB2:S30069?x-access-code=9b532bf68bf8b1b507058d77e93697e4&format=html)

Individual gene topologies of LFA and LFF were reconstructed using the RAxML v.8 [43] program implemented on the CIPRES Science Gateway v.3.3 [44] with a GTRGAMMA model and 1,000 bootstrap (BS) replicates. Before concatenating the data sets, the individual gene topologies were checked for congruence. A conflict was significant when two data partitions supported conflicting monophyletic groups with ML bootstrap values of ≥70% in both trees [45]. Independent gene data sets were combined to achieve maximum phylogenetic resolution without evidence of conflicting nodes.

Bayesian inference was performed on the LFA and LFF concatenated data sets in MrBayes v.3.2.7. [46] on the CIPRES Science Gateway v.3.3. The best-fitting models of nucleotide substitutions were selected using the corrected Akaike information criterion (AICc) as suggested by jModelTest [47]. Two parallel Markov chain Monte Carlo (MCMC) runs were performed, each using four chains and 10,000,000 generations, sampling trees every 10,000^th generation. A 50% majority-rule consensus tree was generated from the combined sample of both runs after discarding the initial 25% as burn-in.

The phylogenetic trees were visualized using FigTree v.1.4.0 [48]. All clades with BS of ≥70% and posterior probability (PP) of ≥0.95 were considered significantly supported.

4.4. Species delimitation

The LFF identification was mainly based on DNA sequences combined with the phenotypic examination. The phylogenetic relationship of LFF species was first assessed with ITS trees using sequences generated newly for this study and downloaded from GenBank with ML and Bayesian methods and re-evaluated by the four-locus phylogeny. To determine the identity of LFA and their phylogenetic positions, we incorporated ITS and rbcL sequences from this study in the Trebouxia dataset from Muggia, Nelsen, Kirika, Barreno, Beck, Lindgren, Lumbsch and Leavitt [5] for phylogenetic analysis. The identity of Trebouxia species was then tested using three species delimitation approaches, i.e., Automatic Barcoding Gap Discovery (ABGD) [49], a Bayesian implementation of the PTP approach (bPTP) [50], and the coalescent-based General Mixed Yule Coalescent (GMYC) approach [51].

ABGD infers a model-based confidence limit for intraspecific divergence and detects the barcode gap as the first significant gap beyond this limit to partition the data. The primary data partitions are then recursively split to get finer partitions until no further partitions can be detected [49]. The ABGD analyses were performed at the web interface (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html) for each data set of ITS, SSU, and rbcL separately. Genetic distances were calculated using the JC69 model, with a prior P ranging from 0.001 to 0.1 and a relative gap width ranging from 1.0 to 1.5. The bPTP analyses were conducted on the bPTP Web Server (http://species.h-its.org/ptp/) based on unrooted phylogenetic trees inferred by Maximum Likelihood. The analyses were run for 200,000 generations, using 0.1 burn-in and 100 thinning. Both ML and Bayesian solutions were examined. For the GMYC analyses, ultrametric trees were estimated from the combined sequence dataset using the BEAST v2.6.3 program [52]. We ran MCMC chains for 15 million generations under the coalescent model with constant population size and a constant clock as the tree prior. Default values were used for the remaining priors. The chain mixture and convergence were evaluated in TRACER v.1.7.2 [53] after removing 10% of the samples as burn-in. The effective sample size (ESS) values greater than 200 were considered a good indicator. The consensus tree was generated using TreeAnnotator 1.8.2 [52] after discarding the first 3,000 trees. The GMYC analysis was executed under the single-threshold model using the SPLITS package [51] available for R 4.0.5 [54].

4.5. Interaction network analyses

The mutualistic association network of Trebouxia species and their associated fungal symbionts (Table S3) was represented in R using the function plotweb in the bipartite package based on the species trees [55]. Bipartite network analyses were conducted to measure the complexity of the interactions with five commonly used qualitative indexes, namely connectance, linkage density, links per species, number of compartments, and nestedness, using the functions networklevel in the R package bipartite [55,56].

4.6. Variation partitioning analysis

The relative effects of climatic variables, geographical distance, the symbiotic partner, and other factors, such as distribution types, reproduction modes, and host specificity (Table S4) on the variance in LFA and LFF diversity were analyzed by variation partitioning analysis based on redundancy analysis, using the varpart function in the vegan package [57]. The phylogenetic distances of LFA or LFF were calculated using the Patristic software with their phylogeny as the input [58], and the first ten PcoA axes were used as a response variable. The 19 climatic variables for the locations of 200 samples (Table S4) were obtained from the Worldclim database (www.worldclim.org) with a grid cell resolution of 2.5 min, using the software DIVA-GIS 7.5 [59]. These 19 environmental variables were implemented forward selection using the function RDA in vegan with the distribution of LFA and LFF, respectively. Geographical distance values (latitude and longitude) were transformed to the principal coordinates of neighbor matrices (PCNM) vectors representing the geographical distances at various spatial scales. PCNM vectors were calculated based on the pairwise geographical distances obtained by the function pcnm in the vegan package. All four PCNMs were used for the analysis. The presence/absence matrix of 30 variables, including distribution types, reproduction modes, and host specificity, was combined as a mixed factor. Similarly, forward selection with the distribution of LFA and LFF was applied, respectively.

4.7. Cophylogenetic analyses

Two distance-based methods, including PataFit (Legendre et al., 2002) and Procrustean Approach to Cophylogeny (PACo) [60], and an event-based method, i.e., Jane 4.0 [61], were used to assess the cophylogenetic patterns between LFF and LFA. One sample corresponding to each LFF-LFA association was chosen for the cophylogenetic analysis to assess potential cophylogenetic patterns.

Distance-based analysis was performed to assess the overall congruence between host and parasite phylogenies accompanying quantification of the relative contribution of individual host-parasite associations to the overall congruence. ParaFit was implemented in R package APE [62]. The ParaFitGlobal fit between parasite and host phylogenies was computed, and statistical significance was tested by 999 permutations to assess whether parasites are randomly associated with their hosts. If the observed ParaFitGlobal statistic is larger than the randomized statistic in more than 95% of cases, the null hypothesis of independent evolution can be rejected, and congruence can be assumed. The contribution of individual host-parasite associations is conducted using the ParaFitLink1 and ParaFitLink2 tests. A significant link suggests that a particular host-parasite association contributes to the global congruence between the host and parasite trees. PACo analysis was performed in R package APE and VEGAN [57], in which the dependence of algal phylogeny on the fungal host phylogeny was tested through 100,000 permutations. The contribution of each host-parasite association to the global fit is measured using jackknife estimation of their respective squared residuals and confidence intervals associated with each host-parasite association.

The event-based analyses were performed using Jane 4.0, which supports polytomies and accommodates parasites with multiple hosts. The topology-based program assigns costs to five coevolutionary events: cospeciation, duplication, host switch, loss/extinction, and failure to diverge [61]. It is difficult to estimate the relative cost of events; instead, a default event cost scheme (cospeciation = 0, duplication = 1, host switch = 2, loss/extinction = 1, failure to diverge = 1), as well as other 16 cost regimes derived from the default regime were selected to test the overall costs of reconstructions in Jane 4.0. All regimes were analyzed with a population size of 23, and the number of generations has been set to 45, as Conow et al. (2010) recommended. The statistical significance of reconstructions was evaluated with 100 random tips mapping permutations.