Phylogeographic analyses of an epiphytic foliose lichen show multiple dispersal events westward from the Hengduan Mountains of Yunnan into the Himalayas

Abstract Lobaria pindarensis is an endemic species of the Himalayas and the Hengduan Mountains. Little information is available on the phylogeography genetics and colonization history of this species or how its distribution patterns changed in response to the orographic history of the Himalayas and Hengduan Mountains. Based on samples covering a major part of the species' distribution range, we used 443 newly generated sequences of nine loci for molecular coalescent analyses in order to reconstruct the evolutionary history of L. pindarensis, and to reconstruct the species' ancestral phylogeographic distributions using Bayesian binary MCMC analyses. The results suggest that current populations originated from the Yunnan region of the Hengduan Mountains in the middle Pliocene, and that the Himalayas of Bhutan were colonized by a lineage that diverged from Yunnan ca. 2.72 Ma. The analysis additionally indicates that the Nepal and Xizang areas of the Himalayas were colonized from Yunnan as well, and that there was later a second dispersal event from Yunnan to Bhutan. We conclude that the change in climate and habitat related to the continuous uplift of the Himalayas and the Hengduan Mountains in the late Pliocene and middle Pleistocene influenced the geographic distribution pattern of L. pindarensis.


| INTRODUC TI ON
Many regional lichen floras include endemic taxa, even though most lichen species have broad, yet often fragmented, geographic distribution ranges (Galloway, 1996). It remains unclear whether current biogeographical patterns reflect long-distance dispersal or fragmentation of historically continuous ranges. Sexually and asexually produced lichen propagules lack special morphological adaptations for long-distance dispersal (Dettki et al., 2000;Heinken, 1999;Sillett et al., 2000). Therefore, it has been assumed that lichens are often dispersal-limited (Armstrong, 1990;Bailey, 1976). However, some studies indicate that lichens can disperse over very large distances, explaining the bipolar distribution patterns of some species (Fernández-Mendoza et al., 2011;Garrido-Benavent et al., 2018;Myllys et al., 2003;Wirtz et al., 2008). Molecular studies can help us to understand population history, determine the genetic differentiation among lichen populations, and trace past migration events (Sork & Werth, 2014;Werth et al., 2021).
Lobaria pindarensis Räsänen (Lobarioideae, Peltigeraceae) ( Figure 1) is an epiphytic foliose lichen associated with a green algal photobiont. This lichen species has developed a complex reproduction strategy including propagation by fungal ascospores sexually produced in disk-shaped fruiting bodies and by different types of asexual (symbiotic) dispersal units such as isidia, lobules, and thallus fragments (Büdel & Scheidegger, 1996;Scheidegger, 1995). The sexual cycle of the photobiont is suppressed while in symbiosis and only the mycobiont goes through sexual reproduction, resulting in the development of ascospores (Malachowski et al., 1980).

Studies of the geographic genetics of lichen-forming fungi
have led to different conclusions regarding species and geographic scale (Cassie & Piercey-Normore, 2008;Scheidegger et al., 2012;Werth, 2010Werth, , 2011. Although several large-scale studies have been conducted at an intercontinental scale or within European ecosystems (Buschbom, 2007;Dal Grande et al., 2010;Fernández-Mendoza et al., 2011;Geml et al., 2010;Otalora et al., 2015;Walser et al., 2003Walser et al., , 2005Werth et al., 2021), very few population genetic studies of Lobaria species have been carried out in the Himalayas and Hengduan Mts.  and Devkota, Dymytrova, et al. (2019) studied the genetic diversity of L. pindarensis populations throughout the species' distribution range in Nepal, applying 17 fungus-specific and nine alga-specific microsatellite loci, and showed that genetic diversity, allelic richness and gene pool composition were significantly influenced by elevation.
Here, we investigated the genetic variation at nine nuclear loci in the mycobiont of the epiphytic lichen Lobaria pindarensis from the Hengduan Mts. and the Himalayas. We aimed to identify the species' region of origin, evolutionary history, and range expansion by combining information related to geographic and climatic events in the Himalayas and the Hengduan Mts. We present data on the phylogeography of this endemic lichen species in the Himalayas and the Hengduan Mts., which can be used as a reference for future research on the phylogeography of lichens. In addition, three specimens from Nepal were collected by Devkota and Scheidegger and discussed in a previous study .

| Sampling of specimens for genetic analyses
Thalli of Lobaria pindarensis were sampled during 2017-2019. In total, 53 thallus fragments of L. pindarensis were sampled from different collection sites (Table S1)

| DNA extraction, PCR amplification, and sequencing
Genomic DNA was extracted from freshly collected and frozen herbarium specimens. About 15 mg of visually uncontaminated lichen thallus was used from each specimen for molecular analyses. To increase the phylogenetic information, we employed six loci (Lpi02, Lpi09, Lpi10, Lpi11, Lpi14, and Lpi19) that contain Simple Sequence Repeats (SSRs or microsatellites) Devkota et al., 2014). Additionally, two coding regions (partial sequences of two single-copy loci: the RNA polymerase II gene, RPB2, and the translation elongation factor-1α gene, EF-1α) and the noncoding nuclear ribosomal internal transcribed spacers, nrITS, enhanced our dataset to nine loci.
The DNA isolation of all specimens, the PCR and the cycle sequencing of the EF-1α, ITS and RPB2 loci were performed as described by Cornejo and Scheidegger (2015). The PCRs used for the amplification of the SSRs containing loci followed the conditions described by Devkota et al. (2014). These amplicons Nevertheless, all but Lpi02, Lpi11, and RPB2 datasets still contained ambiguously aligned regions, which were processed with the software Gblocks 0.91b on the freely available platform phylo geny.fr (Castresana, 2000;Dereeper et al., 2008Dereeper et al., , 2010Devier et al., 2010) enabling smaller final blocks, gaps on the final blocks and fewer strict flanking positions.
All sequences were aligned and edited using Geneious version 7.1.9 (https://www.genei ous.com). All newly produced sequences were checked using the BLASTN suite of the National Center for Biotechnology Information (NCBI) website (http://www.ncbi.nlm. nih.gov/BLAST/) to verify their close relatives and preclude potential contaminants (Ye et al., 2006). Matrices were aligned with MAFFT (version 7) web service (http://mafft.cbrc.jp/align ment/serve r/ index.html) (Katoh et al., 2019). Specimens used in this study, along with voucher information from the GenBank accession numbers, are listed in Table S2. Phylogenetic and molecular evolutionary analyses were conducted using MEGA version 6 (Tamura et al., 2013). The ambiguously aligned regions were arranged manually for the phylogenetic analyses. The gene fragments were combined using Geneious 7.1.9 for phylogenetic analysis, on the premise that no well-supported (BS > 70%, Nuhn et al., 2013) conflict was detected.
Maximum likelihood (ML) and Bayesian inference (BI) analyses were performed using RaxML v. 7.2.6 (Stamatakis, 2006) and MrBayes 3.1.2 (Ronquist & Huelsenbeck, 2003), respectively. We calculated the nine gene-trees under a maximum likelihood criterion. Models of the DNA sequence evolution for each locus were selected with the program jModelTest 3.7 (Posada, 2008) by using the Akaike information criterion (Akaike, 1973)

| Generating a time-calibrated tree of Lobaria pindarensis
We used a Bayesian method, implemented in the programs BEAUti and BEAST (both version 1.8.4; Drummond et al., 2012), to estimate a time-tree. Secondary calibrations (calibrations based on the results of previous molecular dating studies)  are applied in divergence time analyses. Bayesian inference (BI) based on Markov chain Monte Carlo methods (Yang & Rannala, 1997) was performed using MrBayes 3.1.2 (Ronquist & Huelsenbeck, 2003). The optimal model of molecular evolution and gamma rate heterogeneity was determined as implemented in jModelTest (Posada, 2008) by using Akaike

Information Criterion (AIC). The Markov chain Monte Carlo (MCMC)
algorithm was run for 20,000,000 generations with one cold and three heated chains, staring from random trees. Runs were repeated twice.
We selected a Yule speciation model for the species-tree prior, and the population size model was set to piecewise linear and constant root.

| Ancestral area reconstruction
The distribution range of Lobaria pindarensis was divided into four  Yu et al., 2015). In RASP, the BBM method inputs a posterior distribution of Bayesian inference, in this study the consensus tree from BEAST, to reconstruct the possible ancestral distributions of given nodes via a hierarchical Bayesian approach (Ronquist & Huelsenbeck, 2003). Using the R package BioGeoBEARS, we collapsed the phylogeny to a monophyletic population tree (Matzke, 2013) whose tips are monophyletic populations instead of samples. Outgroups were excluded, as the input phylogeny should include only monophyletic groups (Matzke, 2013) without an outgroup (Yu et al., 2015). The maximum number of areas occupied at each node was set to four. To account for phylogenetic uncertainty, 5 million generations of the MCMC chains were run, with sampling every 100 generations.
The phylogenetic analysis revealed a well-supported topology in which populations generally clustered according to their geographic proximity ( Figure 2). The phylogenetic tree indicated that the clade a from Yunnan, firstly branched off at node 2 to become a sister lineage for the rest of the lineages (Figure 2). All interior nodes of each region and subregion had bootstrap support/posterior probabilities higher than 70%/0.9. According to the BEAST analyses, the stem

| Ancestral area reconstruction
The ancestral area reconstruction analyses suggested that the first range expansion of Lobaria pindarensis occurred in the to Bhutan.
The Qinghai-Tibetan Plateau (QTP) has undergone uplift and peneplanation since the beginning of the Tertiary (Xu & Zhang, 1981).
The plateau was probably at about 1000 m a.s.l. during the Pliocene (Cao et al., 1981;Mao et al., 2021). At that time, the climate of southern Xizang was relatively warm and wet, but the northern regions were more arid and could not support glaciers (Spicer, 2017;Xu, 1992). Compared with other mountain ranges bordering the QTP, the Hengduan Mts. are younger, having uplifted over the last 8 Myr, and they showed a very high species diversification rate during this period (Xing & Ree, 2017). During the continuing uplift the basins within the QTP reached 3000 m a.s.l., and the mountains reached ∼4000-4500 m a.s.l. at the end of the Early Pleistocene (Zheng et al., 2002). These elevation ranges may have enabled species migration in what was likely an ecologically diverse region, F I G U R E 2 Time-calibrated phylogenetic tree derived from BEAST based on nine loci (ITS, EF-1α, RPB2, Lpi02, Lpi09, Lpi10, Lpi11, Lppi14, and Lpi19) for Lobaria pindarensis. Bayesian divergence time estimates of the main nodes are listed on the left. ML bootstrap support values/ Bayesian posterior probabilities greater than 70%/0.9 are indicated. Blue bars represent the 95% highest posterior density intervals for node (mean) ages. With different colors indicating the distribution of the major clades. The color bands on the right side: Purple (Hengduan mts.) and yellow (Himalayas). The outgroup-species L. devkotae was used to root the tree. The major geographic events related to the uplift of the Himalayas and the Hengduan mts. Are indicated at the bottom. The climate had become substantially colder (Zheng et al., 2002), which may have created additional forest habitat enabling increased connectivity between the Hengduan Mts. and the Himalayas.

F I G U R E 3 Ancestral area reconstruction for
Based on our field investigations, it is clear that infrastructure construction and forest management activities in these regions considerably reduce and threaten L. pindarensis habitats in areas that have recently been under intensive development. Further, the use of L. pindarensis in traditional medicine and food by cultures across the Himalayas and Hengduan Mts. (Yang et al., 2021) contributes to its sensitive/threatened status. Effective conservation strategies for the forest landscapes of the Himalayas and Hengduan Mts. must therefore ensure that these unique cradles of biodiversity do not turn into graves for biodiversity (Rangel et al., 2018), including for lichens like L. pindarensis.

ACK N OWLED G M ENTS
We are grateful to the Lichen Herbarium of the Kunming Institute of Botany, CAS for providing sample support and supporting the fieldwork collection. We thank Melissa Dawes for careful linguistic revision of the manuscript and Veronika Zengerer for

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequencing data: GenBank accession numbers and voucher information are listed in Table S1 and S2.