Phylogenetic diversity of two geographically overlapping lichens: isolation by distance, environment, or fragmentation?

Phylogenetic diversification is a precursor to speciation, but the underlying patterns and processes are not well‐studied in lichens. Here we investigate what factors drive diversification in two tropical, morphologically similar macrolichens that occupy a similar range but differ in altitudinal and habitat preferences, testing for isolation by distance (IBD), environment (IBE), and fragmentation (IBF).


| INTRODUC TI ON
This past year, the world celebrated the 250th birthday of Alexander von Humboldt , whose work laid the foundation for the field of biogeography. Humboldt's travel to South America (Venezuela, Colombia and Ecuador), with Aimé Bonpland , between 1799 and 1803, and particularly the exploration of the northern Andes, resulted in a pioneering model of altitudinal zonation in the tropics in the Essai sur la Géographie des Plantes (von Humboldt & Bonpland, 1805). Lichens played an integral role in this concept: 'These stony and icy tops that the eye distinguishes barely above the clouds, are covered only by mosses and lichenous plants '. (von Humboldt & Bonpland, 1805:14). The authors also recognized lichen communities as one of their 15 plant physiognomies in the tropics (von Humboldt & Bonpland, 1805:31). During their exploration of the Colombian and Ecuadorian Andes, Humboldt and Bonpland must have come across the two most common species of the genus Sticta (Peltigeraceae), although their actual names, S. andina and S. scabrosa, were only recently established (Moncada et al., 2020). As we demonstrate in this paper, these two taxa exemplify the critical role of autecology in driving phylogenetic diversity, specifically the connection between altitudinal zonation and biogeography, first established by Humboldt. They thus represent an excellent case study to explore links between geo-and biodiversity in lichen fungi.
Phylogenetic diversification and speciation are often related to key innovations that allow a lineage to conquer a new area or niche space (Heard & Hauser, 1995;Schluter, 2000;Hagen & Kadereit 2003;Coyne & Orr, 2004;Kraichak et al., 2015). However, while a key innovation may explain the success of a lineage in terms of abundance and occupied range, it does not directly cause diversification or radiation into separate lineages. The latter also requires mechanisms of phylogenetic diversification, population fragmentation and reproductive isolation (Kisel & Barraclough, 2010).
Populations are generally expected to follow the classic model of isolation by distance (IBD; Wright, 1943), comparable to the concept of regional equilibrium (Hutchinson & Templeton, 1999). In this model, gene flow is assumed to be limited by distance and hence geographically more distant individuals are expected to be also genetically more different. In contrast, isolation by environment (IBE) denotes a correlation between genotypes and environmental parameters, with the consequence of gene flow being limited between ecological niches, even in close proximity. IBD appears to be more frequent in plants and IBE more common in animals (Sexton et  The survey by Sexton et al. (2014) did not consider fungi (including lichens), because of the few studies available and because mechanisms potentially facilitating IBE are demonstrably rare in fungi. The detection of isolation patterns in fungi is also challenging because the largely stochastic dispersal of typically small diaspores, including vegetative propagules in lichens (Tripp et al., 2016), causes a considerably amount of noise in population structure (Buschbom, 2007).
Since fungi behave more like plants in terms of reproduction and usually lack specific fertilization mechanisms such as vectorized pollinization, one would assume IBD to be the predominant mechanism driving phylogenetic diversification in these organisms. Indeed, the few available studies suggest that fungi mostly diversify through IBD (Allen et al., 2018;Branco et al., 2015;Liti et al., 2009).
While there are various ways for populations to become effectively isolated, either geographically or ecologically (Coyne & Orr, 2004;Losos & Ricklefs, 2009;Ryan et al., 2007;Wiens & Graham, 2005), geographic isolation is generally based on plate tectonics: fragmentation of land masses subsequently separated by water bodies, creation of oceanic islands through volcanic hot spots, and the uplift of mountain ranges along collision zones that generate a topographically disrupted terrain (Mila et al., 2010). Mountain uplift will act selectively upon lineages in dependence of their autecology, particularly their altitudinal preferences, leading to a third mechanism of isolation besides IBD and IBE, namely isolation by fragmentation (IBF; Hutchinson & Templeton, 1999). Lineages preferring high altitudes and occurring in disrupted terrain possibly diversify through IBF, whereas lineages occupying lowland regions with plane to undulate terrain likely diversify through IBD (Sexton et al., 2014).
Also, lineages adapted to more illuminated conditions and tolerating disturbances may be less susceptible to IBF than those restricted to well-preserved, closed vegetation. Thus, autecology will variously affect IBF and IBD and lead to different population structures, even if other factors such as lineage age and geographic range are comparable (Cannon et al., 2009;Malhi et al., 2013;Morley, 2011;Shimizu-Kimura et al., 2017).
IBF has been much less frequently studied than IBD and IBE (Sexton et al., 2014), and mostly in a context of anthropic land use changes and conservation (Ghazoul & McLeish, 2001;Saeki et al., 2018). Mountain ranges provide an environment highly susceptible to IBF, and effects of within-lineage phylogenetic diversification and speciation at high altitudes have been shown for a broad range of organisms, mostly plants and birds (Boucher et al., 2016;Givnish, 2010;Givnish et al., 2014Givnish et al., , 2015Hughes & Atkinson, 2015;Küper et al., 2004;Losos & Ricklefs, 2009;Madriñan et al., 2013;Ryan et al., 2007). Demonstrating effects of IBF in mountain ranges is challenging, as first a spatial model for patterns of fragmentation has to be established.
Here we use the two aforementioned species in the genus Sticta, S. andina and S. scabrosa, to test for evidence of IBD, IBE and IBF, in correlation with altitudinal zonation and habitat preferences. These lichens largely disperse by vegetative propagules and, although the symbiotic nature of such propagules leads to larger distribution ranges due to quick establishment and subsequent re-propagation (Tripp et al., 2016), the relatively large size of the propagules causes dispersal limitations, making these taxa excellent study objects to analyse scale-dependent phylogenetic diversification patterns. Both species were previously identified with the broadly defined name S. weigelii (Galloway, 1994(Galloway, , 1998a(Galloway, ,b, 2007, but are only distantly related to each other and not closely related to S. weigelii s. str. (Moncada et al., 2020). The two species have largely congruent geographic ranges, being broadly distributed in the Neotropics and Hawaii, and S. andina is also known from the Azores; both are comparatively abundant lichens. Their main differences are found in their autecology: S. andina is a tropical montane species found in well-preserved, upper montane cloud forest and paramo (usually above 2,500 m), whereas S. scabrosa occurs in tropical lowland to lower montane forest (usually below 1,500 m), often extending into open and anthropic vegetation, locally with a 'weedy' character. We therefore hypothesized that S. scabrosa follows IBD, whereas S. andina exhibits IBF. Both taxa are well-sampled at global level, representing the two largest species-level clades in the currently available global Sticta ITS phylogeny, with 164 and 180 OTUs, respectively, and therefore constitute an excellent model to approach this question.

| ITS barcoding
In this study, we focus on the fungal ITS barcoding marker (Schoch et al., 2012), which is highly effective in delimiting species in the genus Sticta and produces phylogenies congruent with other ribosomal and protein markers (Magain & Sérusiaux, 2015;Simon et al., 2018;Widhelm et al. 2018). Congruence with other markers and with phenotype characters of the ITS-based lineages also preclude the possibility of ITS paralogy; in general, reported cases of presumed ITS paralogy in fungi mostly do not reflect gene duplication but hybridization (Lücking & Hawksworth, 2018). The exclusive use of ITS sequences allowed to place the samples in a broad context including numerous species . The S. andina complex initially comprised 20 OTUs, all sampled in Colombia, believed to comprise three species informally labelled 'andina', 'colombiana' and 'paramuna' ; for this study, 144 ITS sequences were newly generated, resulting in a total of 164 OTUs from Mexico, Costa Rica, Colombia, Ecuador, Brazil, the Azores and Hawaii. Sticta scabrosa was initially based on 10 OTUs from Colombia and the Dominican Republic , while the current set included 165 new accessions, for a total of 175 OTUs from Mexico, Costa Rica, the Dominican Republic, Puerto Rico, Colombia, Brazil, Argentina, Galapagos and Hawaii (Appendix S1). DNA extraction, sequencing and sequence assembly was performed at the Field Museum, the University of Liège and the Botanischer Garten und Botanisches Museum, Freie Universität Berlin, with methodological details described elsewhere (Lücking et al., 2017;Magain & Sérusiaux, 2015;Simon et al., 2018).

| Phylogenetic analysis
To delimit the two species and to verify their position in the global Sticta phylogeny based on previous analyses Widhelm et al. 2018), the sequences of S. andina and S. scabrosa were aligned with 338 ITS sequences of Sticta from GenBank, using Lobaria pulmonaria as outgroup (Appendix S2). The alignment was assembled in BIOEDIT 7.2.5 (Hall, 1999) and sequences were in alignment with MAFFT 7.244 (Katoh et al., 2002(Katoh et al., , 2009) using the [--auto] option. The final alignment included 677 ingroup OTUs and was 626 bases long. Phylogenetic analysis was performed using maximum likelihood in RAxML 8.2.8 on the CIPRES Science Gateway (Miller et al., 2010;Stamatakis, 2015;Stamatakis et al., 2008), applying the GTR-Gamma model and 1,000 bootstrap replicates. The resulting tree was visualized in FIGTREE 1.4.2 (Drummond & Rambaut, 2007) and was used to delimit the two lineages. Subsequently, two subtrees were generated for each lineage, with the most closely related species as outgroup: for S. andina, we used 164 ingroup and six outgroup sequences (S. phyllidiata, S. squamifera), with a total alignment length of 543 sites (Appendix S3), whereas for S. scabrosa, we used

| ITS haplotype network and expansion analysis
ITS haplotype networks were elaborated for both S. andina and S. scabrosa using POPART 1.7 (Leigh & Bryant, 2015) on the separate species alignments (Appendix S3, S4) but excluding the outgroup sequences; we thereby employed the Templeton-Crandall-Sing method or TCS (Clement et al., 2002) and coded the main geographic origins of the specimens, and in the case of S. andina for Colombia also the departments.
To test for haplotype structure in terms of a potential recent expansion (selective sweep), we employed Tajima's D on subsets of the ITS alignments for both species with only specimens from continental Central and South America included. The test was performed using DnaSP 6 (Rozas et al., 2017). To evaluate the effect of terminal and internal gaps, we analysed both the original (raw alignments) and subalignments where sites with internal gaps (indels) were removed, sequences with long terminal gaps were removed and the remaining alignments were trimmed or short remaining terminal gaps in constant sites were filled with the corresponding bases. This resulted in two sets of alignments: andina continental original (158 terminals, 543 sites), andina continental edited (147 terminals, 511 sites), scabrosa continental original (78 terminals, 548 sites) and scabrosa continental edited (62 terminals, 503 sites).

| Environmental data and phenotype analysis
The 339  Universität Berlin, using standard microscopical techniques described in Moncada (2012) and Ranft et al. (2018). For all specimens, we also recorded altitude (Appendix S5). Using the georeference data, we assigned values for 19 climate variables to 337 specimens, using the bioclimatic variables BIOclim1 to BIOclim19 and the altitude layer (Appendix S6) from Worldclim in 2.5 arc minutes (http:// www.world clim.org; Fick & Hijmans, 2017;Hijmans et al., 2005), as well as the Landsat tree cover layer (http://glcf.umd.edu/data/ lands atTre ecover). Two specimens had to be excluded: for one (from Colombia), no precise georeference data could be obtained, and the other (from the Azores) did not have a close enough grid match among the Worldclim grid cells. In the case of Hawaii, we also analysed 56 additional, older collections housed at HAW for altitudinal and ecological differences. Differences in altitude, climate data and tree cover between the two taxa were statistically compared using the nonparametric Mann-Whitney U and the more sensitive Kolmogorov-Smirnov tests in STATISTICA TM 6.0.

| Isolation by distance, environment and fragmentation
In order to test for isolation by distance (IBD), environment (IBE) and/or fragmentation (IBF), we elaborated distance matrices for each of the parameters for both species under different scenarios, assessing the data globally, on a continental scale (to filter effects of rare long-distance dispersal e.g. to island biota), and for the northern Andeas, focusing on a region with disruptive topology.
Based on georeferenced data for all samples, pairwise geographic distances for samples within each species were calculated using Batch Distance Computations (Morse, 2011). Pairwise genetic distances based on the ITS barcoding marker were computed from the above species alignments using Kimura's two-parameter distance computed in DNADist 3.69 (Felsenstein, 1993(Felsenstein, , 2008. Correlations between pairwise geographic and genetic distances were tested for each species using a simple Mantel test as employed in zt (Bonnet & Van de Peer, 2002), with 10,000 permutations. The tests were employed both at global level (all known samples) and for continental samples in Central and South America only.
In order to evaluate isolation by environment, we used the com- developed a spatially constrained permutation procedure to replace the partial Mantel test. However, this promising approach can currently only be applied to site data (e.g. Cañedo-Argüelles et al., 2020, in which each site is unique and therefore the spatial vector contains no duplicates. For specimen data in which more than one specimen share the same site, this method is currently not applicable. To assess a potential effect of isolation by fragmentation in Sticta andina in the northern Andes, we developed a migration model for that species following reconstruction of the final uplift of the northern Andes (Hoorn et al., 2010). In Colombia, the Andes divide into three mountain ranges (Cordillera Occidental, Central and Oriental), which allows to hypothesize migration of altoandine species from south to north and subsequent fragmentation of populations along each of these three ranges. We consequently arranged the Ecuadorian and Colombian samples of the species into eight spatial clusters (A-H), allowing for dispersal parallel to each mountain range. 'Fragmentation distances' between each cluster were then calculated as the total of dispersal lines connection two clusters (Figure 1), resulting in distance scores between 0 (samples within the same cluster) and 4 (between clusters F and D, G or H). In this case, a simple Mantel test was applied to the distance matrices.
All distance matrices are available in a zipped file (Appendix S7). Appendix S8-S10). Due to lack of terminal resolution in the global ITS data set (due to the higher frequency of indels), the species most closely related to S. scabrosa, S. laselvae and S. pseudobeauvoisii, appear nested within S. scabrosa in the global tree, but in smaller data sets and multilocus analyses form a sister clade (Moncada, Reidy, et al., 2014;Widhelm et al. 2018).

| Phylogeny and phenotype
Sticta andina and S. scabrosa can be considered geographically congruent, as they exhibited rather similar distribution ranges, predominantly across the Neotropics and Hawaii, with S. scabrosa also including the Caribbean (Dominican Republic and Puerto Rico) and the Galapagos Islands, whereas S. andina was further present in the Azores (Figure 4). Sticta andina included both apotheciate and isidiate to phyllidiate specimens, as well as sun and shade forms with higher or lower concentration of brown melanine pigments (Figures 3a, 5a-f). Sticta scabrosa was uniformly phyllidiate but, besides most specimens having an uneven lobe surface, featured two additional, unique surface morphodemes in the Hawaiian subspecies, either faveolate-pitted or papillose (Figures 3b, 5g-l).
The internal topology of Sticta andina lacked distinct correlations between clades and morphodemes or geography, but exhibited some partial structure (Figure 3a, Appendix S9). Thus, clades I, IV and V were predominantly phyllidiate, with some apotheciate morphs nested within, whereas clades II and III mostly included isidiate specimens, with some apotheciate and phyllidiate specimens intermingled.
Some small subclades with peculiar propagule types (verruciform or dark isidia or small squamiform phyllidia) were found in clades II, III and V. In S. scabrosa, subspecies hawaiiensis formed a terminal clade (clade IV), but within that clade, the three surface morphodemes were not phylogenetically clustered (Figure 3b, Appendix S10).

| ITS haplotype network and expansion analysis
The ITS-based TCS haplotype networks for the two species showed substantial differences. Sticta andina exhibited a highly reticulate network with numerous interconnected haplotypes; while some haplotypes were more frequent than others, and no haplotype was particularly dominant (Figure 6; see also Appendix S3). There was no immediately obvious correlation between haplotypes and distribution ranges, although some patterns were discernible. Thus, specimens from the northern Andes in Colombia (Cundinamarca to Risaralda)

| Ecological differentiation
Occurrences of Sticta andina and S. scabrosa differed significantly in 13 bioclimatic variables for both statistical tests, including all temperature variables, and an additional five precipitation variables for the more sensitive Kolmogorov-Smirnov test (Table 1). Mean annual temperature in the grids with occurrence of S. andina was five degrees lower than in those corresponding to S. scabrosa. Temperature seasonality as well as mean temperature in the warmest and coldest month and the wettest, driest, warmest and coolest quarter was also significantly higher in grids with occurrence of S. scabrosa. Significant differences were further found in precipitation seasonality and precipitation in the driest month and driest quarter (Table 1). Mean tree cover in the grids with Sticta scabrosa was not significantly different with either test (Table 1).
A marked difference was found in altitudinal preferences, both for the altitudinal grid layer and the actually measured altitude at the

F I G U R E 3 Individual circle cladograms of Sticta andina (A) and
S. scabrosa (B) based on the ITS fungal barcoding locus, depicting reproductive morphodemes in S. andina and main distribution and lobe surface morphodemes in S. scabrosa. See Appendices S9 and S10 for fully labeled phylograms [Colour figure can be viewed at wileyonlinelibrary.com] sample localities (Table 1) (Figure 4). In the Azores, S. andina was extremely rare, being confined to little disturbed cloud forest on the island of Pico; it was not found on the nearby island of Faial, which also has remnants of these natural forests (Purvis & James, 1993).

| Isolation by distance, environment and fragmentation
At global level, S. scabrosa showed a highly significant correlation between genetic and geographic distances (IBD), whereas such a correlation was absent in S. andina (Table 2) (Table 2).
Both taxa exhibited a similar range of Kimura distances, mostly between 0 and 0.015, with S. scabrosa including some instances of slightly higher values, suggesting that both lineages were of similar evolutionary age.

| D ISCUSS I ON
The two studied species, Sticta andina and S. scabrosa, differed strongly in their phylogenetic structure. Sticta scabrosa had a predominant haplotype throughout most of its range except Hawaii. The high number of haplotypes in Colombia (9) Lücking 34607;h,Brazil,Lücking et al. 40073;i,Brazil,Lücking et al. 40042;j,Hawaii,Moncada et al. 6917). k. Damp, olive-green thallus with marginal and laminal phyllidia and minute surface papillae (Hawaii, Moncada et al. 7017). l. Brown thallus with distinctly foveolate lobe surface (Hawaii, Moncada et al. 7009) [Colour figure can be viewed at wileyonlinelibrary.com] northern Argentina (4), compared to the low number in Mexico (1), Costa Rica (1), Dominican Republic (1), Puerto Rico (2), Galapagos (1) and Hawaii (1), together with the more terminal position of specimens from Central America, the Caribbean, Galapagos and Hawaii, were consistent with the interpretation of continental South America being the evolutionary centre of this species, although this hypothesis needs to be tested with additional markers and specimens. The genetic structure of this species within its continental range, with a single dominant haplotype and several infrequent satellite haplotypes, indicates a recent, rather a fast expansion consistent with its ecological interpretation as a somewhat 'weedy' species; however, our analysis using Tajima's D also showed that this test is sensitive to how missing data and indels are treated by editing the underlying alignment prior to analysis. For instance, a single gap within a polymorphic site will cause exclusion of that site from the computation of Tajima's D.
Besides the main haplotype in Sticta scabrosa, the only other distinctive and predominant haplotype was the one in Hawaii, representing subspecies hawaiiensis (Moncada et al., 2020). In contrast to the rather uniform morphology of the neotropical populations, the Hawaiian material featured three distinct morphodemes: the common morphodeme with uneven lobe surface, a second morphodeme characterized by faveolate to pitted lobes especially towards the tips and a third morphodeme with surface papillae. The common morphodeme occurred throughout Hawaii, whereas the faveolate-pitted morphodeme was found on Kauai and rarely on Maui, and the papillose morphodeme was exclusive to Maui.
The phenomenon of morphological disparity combined with a rather uniform genotype is proper to island biotas and has been well-documented for plant radiations, such as the Hawaiian lobeliads in the family Campanulaceae and the Silversword alliance in the family Asteraceae (Baldwin et al., 1991;Baldwin & Sanderson, 1998;Barrier et al., 2001;Carlquist et al., 2003;Givnish et al., 2009). For lichens, it was recently shown for the genera Lobariella and Pseudocyphellaria in Hawaii, with closely related species developing distinctive morphologies (Lücking et al., 2017;Moncada, Reidy, et al., 2014).
In contrast to Sticta scabrosa, S. andina had a much larger number of frequent haplotypes, interconnected in reticulate fashion.
The highest number of haplotypes was found in the northern Andes (Colombia and Ecuador), decreasing towards North America in the recent past, leading to a phylogenetically isolated subspecies (Moncada et al., 2020 Of particular interest for the accidental dispersal of associated lichens is Cupressus lusitanica Mill., native to Mexico and northern Central America but widely naturalized in Central and South America including at high altitudes in Colombia (Little & Skolmen, 1989).
In Sticta andina, there was no significant correlation between genetic and geographic distances, neither at global nor at continental level, although at continental level, the correlation was marginally significant, indicating a weak effect of IBD. However, maximum genetic distances were found already at small geographic distances, corresponding to the scale of topographical isolation in the highly distrupted terrain of the northern Andes. This pattern was closest to the case III model defined by Hutchinson and Templeton (1999),  and Lobaria pulmonaria (Walser et al., 2005;Fernández-Mendoza et al., 2011;Allen et al., 2018), IBD was the predominant pattern detected. In plants and animals, where IBE appears to be the main pattern besides IBD (Sexton et al., 2014), IBF has also been demonstrated for high altitude lineages, such as andine orchids and paramo plants (Givnish, 2010;Givnish et al., 2014Givnish et al., , 2015Küper et al., 2004;Madriñan et al., 2013). A recent study on plants of the Cape Floristic Region (Verboom et al., 2015) found effects of IBF to be predominant among higher altitude lineages, as opposed to IBD combined with IBE in lower-altitude species.
Our findings also lend support to the notion that key innovations alone do not cause phylogenetic diversification but require additional steps of isolation. Given that Sticta andina and S. scabrosa have comparable, rather wide geographic ranges, both with similar abundances, one would consider them as equally 'successful' in evolutionary terms. Thus, both likely acquired key innovations leading to their range expansions compared to other species, although the nature of these is not known (though likely physiological). However, only in S. andina we found a level of genetic, phenotypic and ecological diversification potentially leading to subsequent speciation across its continental range in Central and South America, whereas S. scabrosa is genetically, morphologically and ecologically more homogeneous within the same area and is not likely to further speciate, except for long-distance, long-term isolation, as is the case with the Hawaiian subspecies.
In conclusion, Sticta andina and S. scabrosa constitute an excellent case study to highlight the link between altitudinal and habitat preferences, that is, autecological features and geographical metapopulation structure, which relates to species biogeography. The early work of von Humboldt and Bonpland (1805), largely performed in the region where the two lichens have their centres of distribution, already hinted at such a connection, although at the time the underlying evolutionary mechanisms were far from being recognized.

ACK N OWLED G EM ENTS
Funding for field and laboratory work for this study was provided by a grant from the National Science Foundation (NSF) to The