Genetic variation and factors affecting the genetic structure of the lichenicolous fungus Heterocephalacria bachmannii (Filobasidiales, Basidiomycota)

Heterocephalacria bachmannii is a lichenicolous fungus that takes as hosts numerous lichen species of the genus Cladonia. In the present study we analyze whether the geographical distance, the host species or the host secondary metabolites determine the genetic structure of this parasite. To address the question, populations mainly from the Southern Europe, Southern Finland and the Azores were sampled. The specimens were collected from 20 different host species representing ten chemotypes. Three loci, ITS rDNA, LSU rDNA and mtSSU, were sequenced. The genetic structure was assessed by AMOVA, redundance analyses and Bayesian clustering methods. The results indicated that the host species and the host secondary metabolites are the most influential factors over the genetic structure of this lichenicolous fungus. In addition, the genetic structure of H. bachmannii was compared with that of one of its hosts, Cladonia rangiformis. The population structure of parasite and host were discordant. The contents in phenolic compounds and fatty acids of C. rangiformis were quantified in order to test whether it had some influence on the genetic structure of the species. But no correlation was found with the genetic clusters of H. bachmannii.


Introduction
The term "lichenicolous fungi" refers to a polyphyletic group of fungi specialized in living on lichens, whether as parasites, as commensals, or as saprophytes [1][2]. It is widely accepted that most of the species of lichenicolous fungi are highly specific, restricting themselves to a few host species that generally belong to a single genus [2][3]. Nevertheless, little is known about the factors that determine the host specificity and the genetic structure of the lichenicolous fungi. Very few researches using molecular tools have addressed these questions thus far. Molina

Sampling and loci selection
The specimens were collected between July 2014 and July 2015 in areas belonging to three different biogeographical regions: the Azores Islands (Macaronesian region), Southern Europe (Mediterranean region) and Southern Finland (Hemiboreal region). On every locality, one to five specimens of H. bachmannii per host species were collected. The collected specimens were spaced at least 5 m from each other. To complete the sampling, specimens deposited in the herbaria H, LE, and MACB were selected, including eight specimens from America, and two from Asia. The new collections were deposited in H. In total, DNA sequences from 123 specimens of H. bachmannii were obtained ( Table 1). The Cladonia host species were identified according to Ahti & Stenroos [20], by morphological and chemical study of the specimens. Three loci were selected for the population-based study of H. bachmannii: ITS rDNA, LSU rDNA, and mtSSU. Of these, ITS rDNA and LSU rDNA were selected according to the results of Millanes et al. [6], who recommend the use of these loci as barcodes in Tremellomycetes. The following loci were tested to determine which of them was the most informative at population level: l41 (with the primers L41F/L41R), IGS rDNA (with LR12R/5SRNA), atp6 (with ATP6-1/ATP6-2), rpb2 (with bRPB2-6F/bRPB2-11R), rpb1 (with RPB1-Af/RPB1-CR), ef1α (with 983F/1567). So far no sequences of H. bachmannii were obtained for any of these loci. Either the amplification failed or the obtained sequences corresponded to the lichen mycobiont. Therefore we decided to use mtSSU, the only additional locus that worked.
In addition, DNA was extracted from one of the host species, Cladonia rangiformis Hoffm., trying to amplify the same loci used for H. bachmannii, such as recomended by Vienne et al. [21]. However, mtSSU was uninformative for C. rangiformis, whereby we amplified IGS rDNA, a highly variable locus that we had previously tested in several species of the genus Cladonia [22].

DNA extraction, PCRs and sequencing
We cannot discard the coinfection of a single thallus of Cladonia by several strains of H. bachmannii. For this reason, we selected a single gall of each sample of H. bachmannii for the DNA extraction. The total DNA was extracted using E.Z.N.A. forensic DNA kit (Omega Bio-tek, Georgia, U.S.A.), according to the manufacturer's instructions. The DNA was eluted in the final step in 100 μl of elution buffer included in the kit. PCRs were carried out with Ready-to-Go-PCR Beads (GEHealthcare Life Sciences, UK). The volume of the reaction was 25 μl with 3 μl of extracted DNA. The primers used were: ITS1F/BasLSU3-3 [23][24] for ITS rDNA, BasLSU3-5/ LR5 [24][25] for LSU rDNA, and MS1/MS2 [26] for mtSSU (primer sequences in S4 Table). The PCR programs for ITS rDNA and LSU rDNA were the same used by Millanes et al. [24], and for mtSSU: 95˚C for 5 min; 35 cycles of 95˚C for 30 s, 50˚C for 30 s and 72˚C for 1 min; with a final extension at 72˚C for 10 min. The primers used to amplify C. rangiformis were: ITS1F/ITS4 [23,26] for ITS rDNA, LROR/LR5 [25] for LSU rDNA, and IGSf/IGSr [27] for IGS rDNA. The PCR programs for ITS rDNA and IGS rDNA were described in Pino-Bodas et al. [28], and for LSU rDNA 95˚C for 5 min; 30 cycles of 95˚C for 30 s, 55˚C for 30 s and 72˚C for 1 min; with a final extension at 72˚C for 10 min. PCR products were purified with E.Z.N.A. Ultra-Sep Gel Extraction Kit or ExoSAP-IT (USB Corporation, OH, USA). The sequencing was performed at Macrogen Europe service (www.macrogen.com), with the same primers as used for the PCR.

Secondary metabolites of the hosts
All the host specimens were studied by thin layer chromatography (TLC) according to the standardized procedures [29][30], with the solvent systems A and B. The secondary metabolites of C. rangiformis were additionally studied by high-performance liquid chromatography (HPLC) and ultra-performace liquid chromatography-mass spectrometry (UPLC-MS) according to [9,[31][32]. The fatty acids of C. rangiformis were studied by gas chromatography-mass spectrometry (GC-MS). Secondary metabolites extractions, HPLC, GC-MS and UPLC-MS protocols are included as Supporting Information (S1 File).

Phylogenetic analyses and haplotypes networks
The sequences were assembled in Sequencher 4.1.4 (GeneCodes, Ann Arbor, MI). The alignments were implemented in MAFFT [33] and BIOEDIT [34]. BLAST searches were done to verify the identity of the sequences. A few sequences corresponded to other Tremellomycetes genera and were deleted. A phylogenetic analysis based on ITS rDNA and LSU rDNA was carried out to test whether H. bachmannii is monophyletic. To this end the following species were selected: two species of Hetereocephalacria, two species of Syzygospora, one species of Piskurozyma, three species of Filobasidium, three of Goffeauzyma (voucher specimens in Supporting Information, S1 Table). As outgroup Cystofilobasidium bisporidii, C. capitatum and C. ferigula were selected according to the results of Millanes et al. [26] and Weiss et al. [35]. The ambiguous regions were removed using Gblock [36] with the less stringent options. Each region and the combined dataset were analyzed by maximum likelihood (ML). The ML analyses were implemented using RAxML 7.0.3 [37] assuming the GTRGAMMA model. The node support was estimated with rapid bootstrap algorithm, using 1000 pseudoreplicates. Congruence between the loci was tested following Lutzoni et al. [38], manually checking the clades with at least 70% bootstrap support. No incongruity was detected.
Haplotype networks for each locus under statistical parsimony were constructed in TCS 1.21 [39], considering the gaps as a 5th character. Genetic structure analyses Summary statisticals including haplotype diversity and nucleotide diversity were calculated in DnaSP v.5 [40]. In order to assess the contribution of the different potential factors (host species, geographical origin and host secondary metabolites) to the overall genetic variation of H. bachmannii, analyses of molecular variance (AMOVA) were conducted in ARLEQUIN V3.5 [41]. In these analyses only the groups containing more than five specimens were considered.
To examine the relative contribution of these three factors, redundancy analyses (RDA) and partial redundancy analyses (pRDA) were run in R (R Development core Team 2017), using the Vegan package [42]. For these analyses a binary matrix with the haplotypes of H. bachmannii was used as dependent matrix. Other three binary matrices containing data of host species, host secondary metabolites (classified in chemotypes), and the geographical origin (The Azores, Southern Europe, Southern Finland, America, Asia) were used as explanatory matrices. The variation explained by each variable group was estimated using adjusted R2 because the number of variables in each matrix was not the same. The statistical significance was assessed using a permutation-based ANOVA test with 2000 permutations. A bayesian clustering algorithm to estimate genetically homogeneous groups was implemented in STRUCTURE 2.3.4 [43]. Although this method was implemented to infer the population structure using unlinked markers, several studies have proved that the SNPS from sequences data as independent loci are also suited [44][45]. The analyses were run assuming an admixture model, without consideration of locality or host species of the specimens and allele frequencies independently modelled. The analysis was performed with five runs per K value (number of clusters), the range of K was from 1 to 20 with 100.000 iterations discarded as burnin and 2000.000 iterations kept for each replicate. The optimum value of K was calculated with the Evanno method in STRUCTURE HARVESTER [46][47]. The results of different runs were combined using CLUMPP 1.1.2 [48] and the barplots were generated with DISTRUCT 1.1 [49]. To assess whether the specimens of H. bachmannii on different species hosts, on different chemotypes or from different geographical regions were randomly distributed across the clusters, Chi-square tests were performed (http://www.physics.csbsju.edu/stats/). Only the groups (host species, secondary metabolites and geographical regions) with a frequency ! 5 were included.

Host-parasite comparisons
To compare the genetic structure of the parasite and the host, we selected Cladonia rangiformis as host species for the following reasons: 1) the largest number of specimens of H. bachmannii was found on this host; 2) it has two different chemotypes, one with atranorin, rangiformic, and nor-rangiformic acids; another that, in addition to the above compounds, contains fumarprotocetraric and protocetraric acids, 3) the largest number of haplotypes of H. bachmannii found growing on this host.
The genetic variation for host and parasite was calculated in DnaSP v.5. Three approaches were used in order to compare the genetic variation of H. bachmannii and C. rangiformis. Firstly, we performed partial Mantel test between pairwise [Fst/(1-Fst)] matrices of H. bachmannii and C. rangiformis correcting with geographic distances with 2000 random permutations using the VEGAN package in R. The pairwise geographical distance matrix was calculated with euclidean distances from the geographic coordinates using dist() function. Three different comparisons were tested: ITS-H. bachmannii/ITS-C. rangiformis, LSU-H. bachmannii/LSU-C. rangiformis and the combined dataset-H. bachmannii/combined dataset-C. rangiformis. Secondly, an analysis was run in STRUCTURE 2.3.4 with the same conditions specified above for the specimens of H. bachmannii and C. rangiformis. The congruence between the clusters of host and parasite defined by STRUCTURE was assessed by computing the proportion of individuals assigned together in the same cluster in both analyses. Contingency table analyses (http://www.physics.csbsju.edu/stats/) were used to test the association between the host chemotype and clusters generated by STRUCTURE. Variance analyses (ANOVA) in STATGRAPHICS 5.1 were carried out to study whether the clusters of H. bachmannii differed in the amount of phenolic compounds or fatty acids. Several ANOVAs were carried out, considering the total amount of substances (total amount of phenolic compounds, and total amount of fatty acids), the total amount of different subsets of fatty acids (saturated fatty acids, unsaturated acids and free fatty acids) and the amount of every phenolic substance separately (atranorin, fumarprotocetraric acid and protocetraric acid) and the amount of most abundant fatty acids separately (linoleic, oleic, palmitic and stearic acids). Bartlett test was used to check the variance homogeneity while Kolmogorov test checked the normality of the variables. All the variables had a homogeneous variance. Three variables were not normal (phenolic compounds total contents, fumarprotocetraric acid contents and protocetraric acid contents) and were analyzed by Kruskal-Wallis.
Finally, simple AMOVA analyses were made to study the genetic differentiation of H. bachmannii and C. rangiformis among geographic regions and host chemotypes.

Results
The 123 specimens of H. bachmannii were found on 20 host species (Table 1) with ten different chemotypes ( Table 1). The chemotype containing fumarprotocetraric and protocetraric acids was the most abundant, found on 47 host specimens from ten species. HPLC and UPLC-MS analyses revealed 14 compounds for C. rangiformis (S1 File). The total contents of phenolic compounds varied from 19.1 to 61.5 mg/g of dry weight (d.w.). Atranorin was the major compound in all the populations except in one of them, where the contents of fumarprotocetraric acid exceeded that of atranorin. By means of GC-MS a total of 16 fatty acids methyl esters (FAME) derived from triglycerides and additionally two free fatty acids (FFA) were identified from C. rangiformis (S1 File).
In total 304, new sequences of H. bachmannii and 90 of C. rangiformis were generated (S5 Table). In the phylogenetic analysis, all the specimens of H. bachmanii formed one strongly supported clade (S1 Fig), into which another lichenicolous species, H. physciacearum was grouped. The haplotype and nucleotide diversity values are shown in Table 2.
For each of the three loci, the haplotypes were connected in a unique network with 95% of confidence (Fig 1). A total of 28 ITS rDNA haplotypes, 24 of LSU, and six of mtSSU were found. The distribution of haplotypes among host species revealed that six ITS rDNA haplotypes, nine LSU rDNA haplotypes, and two mtSSU haplotypes were present on multiple host species, while 22 ITS rDNA haplotypes, 19 LSU rDNA haplotypes and four mtSSU haplotypes restricted themselves to a unique host species. Five haplotypes of ITS rDNA, seven of LSU rDNA, and two of mtSSU were present in several geographical regions. Southern Europe showed the highest number of haplotypes: 19 in ITS rDNA, 14 in LSU rDNA, and three in mtSSU ( Table 2). In the ITS rDNA haplotype network, the most abundant haplotype represented 32.1% of all the specimens, and it was found on several host species (Fig 1). A group of haplotypes restricted to C. rangiformis was separated by four mutational steps (Fig 1A). Another haplotype found on C. granulosa also separated from the rest by five mutational steps. In LSU rDNA network, the most abundant three haplotypes represented 44% of all the specimens, and they were found on several host species (Fig 1B). One haplotype restricted only to C. granulosa and another one only to C. macilenta were the most distant. In mtSSU haplotype network, the most abundant haplotype represented 58% of individuals, and it was present on 11 host species (Fig 1C).

Genetic structure of Heterocephalacria bachmannii
The results of AMOVA analyses are shown in Table 3. Most of the varition was found among the chemotypes (33.3-73.1%) and among host species (29.0-68.5%) while the genetic differentiation among geographical regions explained less percentage of variation (7.0-40.3%). Fig 2 shows the results of the redundancy analyses. The pRDAs show that the proportion of the genetic variation explained by the chemotypes, the host species or the geographical origin, when the other factors are kept under control, was very small. The greatest proportion of genetic variation was explained by the host species in conjunction with the chemotype (0.13 in ITS rDNA, 0.13 in LSU rDNA and 0.53 in mtSSU). The geographic origin on its own explained only a little genetic variation in ITS rDNA and LSU rDNA, while in mtSSU it explained a greater proportion than each of the other two factors. All the three variables, jointly taken, explained only a small proportion of the genetic variation in the three loci.
The clustering algorithm implemented in STRUCTURE revealed that the optimal number of clusters was two (highest values of ΔK were obtained for K = 2). The 90.9% of individuals were assigned to one cluster with membership coefficients > 0.7. Ten speciemens could not be assigned to any of the clusters (membership coefficients < 0.7). The inviduals were nonrandomly distributed in the two genetic clusters (S2 Table). Different association with respect to the host species (Chi-square = 68.4, d.f. = 7, P-value 0.001), chemotypes (Chi-square = 78.5, d.f. = 6, P-value 0.001) and the geographical region (Chi-square = 14.4, d.f. = 4, Pvalue 0.01) was found. One cluster contained specimens on Cladonia rangiformis (n = 26), C. furcata (n = 1) and C. cervicornis (n = 1). All the specimens of this cluster were collected in Southern Europe and in the Azores. The other cluster was associated with the others host species, although five specimes were also associated with C. rangiformis as the host species (Fig 3). The specimens assembled in this cluster were collected in all geographic regions and chemotypes.

Genetic variation of Heterocephalacria bachmannii
Comparison of Heterocephalacria bachmannii with Cladonia rangiformis Table 4 shows the genetic variation of H. bachmannii (collected on C. rangiformis), as well as that of the host species itself. The results of partial Mantel test are shown in Table 5. Cladonia rangiformis and H. bachmannii pairwaise Fst values were not significantly correlated.
The STRUCTURE analysis carried out with H. bachmannii samples collected on C. rangiformis showed a ΔK peak at a K = 2, as did the STRUCTURE analysis of C. rangiformis (Fig 3). In both analyses all the specimens were assigned to one of the two clusters with values > 0.8, but the cluster composition was different. The clusters defined by STRUCTURE for H. bachmannii and C. rangiformis were less congruents, only 43.8% of specimens were assigned to the same cluster. The contingency table analyses indicated that there is not a significant association between the host chemotypes and the H. bachmannii clusters (S3 Table) generated in STRUCTURE (Chi-square = 0.406, d.f. = 1, P-value = 0.524). ANOVA and Kruskal-Wallis analyses did not find significant differences in the contents of the host secondary metabolites between the two clusters of H. bachmannii ( Table 6).
The AMOVA analyses are shown in the Table 7. Significant genetic differenciation was found between H. bachmannii on different host chemotypes and geographical regions.   However, the variance explained by the geographical origin (31.7% in ITS rDNA, 15.9% in LSU rDNA and 58.40% in mtSSU) was higher than the explained by the chemotypes (16.8% in ITS rDNA, 13.6 in LSU rDNA and 15.7% in mtSSU). The genetic variance of C. rangiformis from different geographical origins was only significantly different in LSU rDNA analysis, while no significant differences were found in any loci between specimens with different chemotypes.

Discussion
Owing to the high number of host specific lichenicolous fungi, it has been assumed that these fungi and their lichen hosts have coevolved [2]. Coevolution depends on the genetic variation and on the genetic structure of interacting species [50]. In a coevolution process, the parasite populations would be expected to have a structure determined by the hosts, and with a low gene flow among the host species [51]. Thus far, however, the investigations on the relations between lichenicolous fungi and lichens are still too poor. The genetic variation of H. bachmannii has been studied here for the first time in a broad geographic range. The data indicate that the genetic differentiation of H. bachmannii is related to the host secondary metabolites and the host species, while this species presents few geographic isolation. This study supports the presence of two genetic clusters of H. bachmannii specialized to different hosts.

Genetic variation and genetic structure of Heterocephalacria bachmannii
The number of haplotypes in ITS rDNA and LSU were greater in southern Europe than in the other regions, an expected result considering that the sampling and the number of host species were higher in this region than in the other ones (Table 1). Our results also show mtSSU was less variable than the nuclear loci, as to the number of haplotypes, haplotype diversity, and nucleotide diversity. This result is concordant with prior studies [52][53][54][55][56] conducted in other fungi, the authors considered that the mutational rate in the mitocondrial genome is lower than in the nuclear genome. Nevertheless, this result could be also attributed to the fact that the number of mtSSU sequences obtained was smaller than the number of sequences of the other two loci. To check whether this explanation was plausible, we calculated the number of haplotypes for the nuclear loci taking into account only the specimens for which a sequence of mtSSU was obtained (data not shown). However, the number of haplotypes (20 in ITS rDNA and 15 in LSU rDNA) continued being smaller in mtSSU. According to our results, the host species and the host secondary metabolites are the most relevant factors in the genetic structure of H. bachmannii, while the populations of H. bachmannii from different geographical regions show slight genetic differentiation. It supports the presence of gene flow among the populations separate by long distances. Similar results have been found for other lichenicolous fungi of the genus Tremella [5,57]; while the opposite occurs in Marchandiomyces corallinus [4]. According to Werth et al. [5] the generalist lichenicolous fungi might be geographically structured, while the more specialist species might be structured by the host species. Therefore, H. bachmannii is expected to be structured by its hosts, all of which belong to the genus Cladonia. In other basidiomyceteous non-lichenicolous fungi, populations not geographically structured have been frecuently found [58][59][60][61][62]. It is assumed that fungal spores have an ability for long distance dispersal. However in some cases, no correlation has been found between the genetic and geographic distances [59] and other hypotheses have been proposed. Kretzen et al. [59] proposed two hypotheses to explain the absence of geographic structure in ectomycorrhizal fungi of the genus Rhizopogon. The first holds that competence can prevent offspring from establishing close to parents. The second hypothesis asserts that the mating systems could strongly support the outcrossing, whereby the spores coming from other populations would have a higher probability of survival success. Owing to our limited knowledge about the reproductive biology of H. bachmanii, we can not be sure that either hypothesis explains our results. Therefore, further studies are required to determine dispersal mechanisms of this species. The effect of the lichen secondary metabolites on the lichenicolous fungi had been previously studied. Lawrey [63] found that certain phenolic compounds inhibit the growth of these fungi. Other lichenicolous fungi can only colonize the host species if another fungus has previously degraded those compounds [13]. It is also proved that the thalli of parasitized lichens have a lower concentration of phenolic substances than the non-parasitized ones [64]. However, ours is the first study in which genetic divergence among the populations of a lichenicolous fungus has been proved to be associated with the host secondary metabolites. Werth et al. [5] demonstrated that the host species was the factor that explained most of the data variance in T. lobariacearum, suggesting that the host species could create a selective environment that only certain strains would have the ability to infect. This particular environment might be created by the secondary metabolites of the host species [5]. Therefore it was to be highly expected that the genetic structure of the parasite should be strongly influenced by the secondary Genetic variation of Heterocephalacria bachmannii metabolites of the host. Even though the AMOVA results indicated that secondary metabolites explained slightly more genetic variation of H. bachmannii than the host species, in the redundance analyses the conjuntion of both factors explained most of the genetic variation. The lack of resolution in our analyses may be due to the low intraspecific chemical variation of the hosts infected by H. bachmannii. Although several chemotypes are known for many of the host species studied [22], only one parasitized host species had several chemotypes (C. rangiformis), what makes difficult to separate both effects, "host species" and "host secondary metabolites". We do not know if our sampling was biased towards one chemotype, or if H. bachmanii in some Cladonia species only parasitizes one chemotype. Despite this parasite is not a rare species [17][18], there are no data about its distribution on the chemotypes within the host species and the total range of secondary metabolites it tolerates. It is also remarkable that the redundacy analyses showed high unexplained genetic variation in H. bachmannii. Therefore, other envarionmental factors could also be important in the genetic structure of this lichenicolous fungus.

Discordant genetic structure between Heterocephalacria bachmannii and Cladonia rangiformis
The results also point out that the genetic structures of the host and the parasite are different. The populations of H. bachmannii were genetically more variable in the Azores than in Southern Europe, while the populations of C. rangiformis were more variable in Southern Europe. In addition, the differentiation among populations from different geographical regions was higher for H. bachmannii than for C. rangiformis (Table 7), which reveals a lower gene flow among the populations of the parasite than among those of the host. The differences between the genetic structure of the parasite and their hosts are frequent [65][66][67], which could be due to differences in evolutionary rates, differences in dispersion rates or movement of H. bachmannii among the different potential hosts. Mantel test did not reveal any significant correlation between the genetic distances of the host and the parasite, which indicate the absence of codispersion. Therefore, the hypothesis that the galls of H. bahchmannii were dispersed together with thallus fragments of C. rangiformis is not supported by our data. The potential movement among different host species could explain why the geographical structure of H. bachmannii is more accentuated when only the specimens on C. rangiformis are considered. Other plausible explanation would be that the lineage of H. bachmannii (specialized on C. rangiformis) colonized once the Azores and subsequently has extended across the potential hosts in the islands.
The host chemotypes have less influence on the genetic structure of H. bachmannii on C. rangiformis (Table 7) than expected. These chemotypes differ in the presence or absence of fumarprotocetraric acid, a compound synthesized by a number of Cladonia species [19][20] that act as hosts for H. bachmannii. Most likely a large number of genotypes of H. bachmannii are tolerant to this substance. The genetic structure of the lichenicolous fungi is neither shaped by the amount of phenolic compounds nor by the concentration of fatty acids (Table 6), indicating that the two genetic clusters tolerate a wide range of these compounds.
No evidence of geographical isolation was found between the populations of the Southern Europe and those of the Azores in C. rangiformis. Cladonia species present different genetic structure patterns. While Printzen & Ekman [68] found that populations of C. subcervicornis separated by a few kilometers were genetically isolated, in other species the populations were geographically weakly structured [69][70][71]. The most extreme case was that of C. arbuscula and C. mitis, two species with a bipolar distribution, the specimens of which appear intermixed in the northern and the southern hemispheres [72]. It is assumed that the predominant reproduction in C. rangiformis is asexual by dispersion of thallus fragments, because the apothecia are rare. The propagule size can directly influence the population structure of the species. Since several studies show that thallus fragments are dispersed at close range by wind [73][74], one would expect that the species with this type of reproduction should have a marked geographic structure in its populations. However, this is not always the case; in species of Lobaria with a varying propagule size, no correlation has been found between the propagule size and the geographic structure [75]. In C. rangiformis, gene flow between the populations of the two geographic regions could have two explanations. First, long distance dispersion of thallus fragments by the wind or other vectors. In fact, some studies have proved that birds can act as dispersion vectors for lichen fragments [76][77]. Second, that although the sexual reproduction is less frequent, the ascospores (which easily disperse long distance by the wind) could have a higher establishment success than the thallus fragments. The success of the propagules in developing new thalli is highly depending on the substrate on which they fall, the substrate formed by bryophytes being the most suitable in C. mitis [78].
Supporting information S1 File. Secondary metabolites of Cladonia rangiformis. The analytical protocols, results, tables and references are decribed. Phenolic compounds were identified using HPLC and UPLC-MS, and fatty acids using GC-MS.  Table. GenBank accession numbers of Heterocephalacria bachmannii and Cladonia rangiformis. (XLSX)