Limitations of Species Delimitation Based on Phylogenetic Analyses: A Case Study in the Hypogymnia hypotrypa Group (Parmeliaceae, Ascomycota)

Delimiting species boundaries among closely related lineages often requires a range of independent data sets and analytical approaches. Similar to other organismal groups, robust species circumscriptions in fungi are increasingly investigated within an empirical framework. Here we attempt to delimit species boundaries in a closely related clade of lichen-forming fungi endemic to Asia, the Hypogymnia hypotrypa group (Parmeliaceae). In the current classification, the Hypogymnia hypotrypa group includes two species: H. hypotrypa and H. flavida, which are separated based on distinctive reproductive modes, the former producing soredia but absent in the latter. We reexamined the relationship between these two species using phenotypic characters and molecular sequence data (ITS, GPD, and MCM7 sequences) to address species boundaries in this group. In addition to morphological investigations, we used Bayesian clustering to identify potential genetic groups in the H. hypotrypa/H. flavida clade. We also used a variety of empirical, sequence-based species delimitation approaches, including: the “Automatic Barcode Gap Discovery” (ABGD), the Poisson tree process model (PTP), the General Mixed Yule Coalescent (GMYC), and the multispecies coalescent approach BPP. Different species delimitation scenarios were compared using Bayes factors delimitation analysis, in addition to comparisons of pairwise genetic distances, pairwise fixation indices (FST). The majority of the species delimitation analyses implemented in this study failed to support H. hypotrypa and H. flavida as distinct lineages, as did the Bayesian clustering analysis. However, strong support for the evolutionary independence of H. hypotrypa and H. flavida was inferred using BPP and further supported by Bayes factor delimitation. In spite of rigorous morphological comparisons and a wide range of sequence-based approaches to delimit species, species boundaries in the H. hypotrypa group remain uncertain. This study reveals the potential limitations of relying on distinct reproductive strategies as diagnostic taxonomic characters for Hypogymnia and also the challenges of using popular sequence-based species delimitation methods in groups with recent diversification histories.


Introduction
Molecular sequence data have had a pronounced effect on our understanding of species boundaries, especially in organisms with relatively simple morphologies and considerable variability of phenotypic characters, such as lichen-forming fungi. Similar to most major biological groups, identifying the appropriate character sets is one of the greatest challenges with empirical species delimitation in lichen-forming fungi [1][2][3][4][5][6]. However, fungi generally have few taxonomically informative traits, in comparison to other major clades of life, and intraspecific variation makes accurate taxonomic circumscriptions more difficult. Hence, molecular data now play a prominent role in circumscribing fungal species. Cryptic species are often identified using molecular data, and in some cases cryptic species are corroborated by formerly overlooked phenotypic characters [7][8][9][10][11]. In other cases, some species-level lineages were shown to consist of chemically or morphologically polymorphic individuals that were previously regarded as separate taxa [12][13][14].
Differences in reproductive strategies have traditionally played an important role in circumscribing lichen-forming fungal species, with populations forming asexual diaspores (such as powdery soralia or corticated isidia) being separated at the species level from populations lacking those and exhibiting ascomata [15,16]. However, this classification has been challenged [17,18] and several molecular studies have shown that the taxonomic importance of reproductive mode has probably been over-emphasized in a number of fungal groups [3,[19][20][21][22][23]. This is partly due to a correlation of reproductive mode and environmental modulation. The development of reproductive structures is often correlated with the ontogeny of lichen thalli, since it has been found that some lichen species use a mixed strategy of early asexual and late sexual reproduction [24]. Besides, macro-and microclimatic variables are also reported to affect reproductive capacity, for example isidia (one type of asexual reproductive structure in lichens) are produced in higher frequency under greater microclimate stress (higher radiation and temperatures, lower humidity) [24]. Some lichenologists found a positive correlation between production of apothecia and microclimatic conditions [25,26], and Seymour et al. [27] showed that lichens more frequently produce sexual structures in hostile environments.
During our studies on the genus Hypogymnia (Parmeliaceae) in China, the Hypogymnia hypotrypa species group drew our attention as an important case study for assessing the importance of reproductive strategies to delimit species of lichen-forming fungi. The currently recognized species pair includes the sorediate H. hypotrypa and esorediate H. flavida. Both taxa are characterized within the genus by a large thallus and yellowish color of the upper surface.
Because species delimitation based on presence or absence of soredia has been shown to be incongruent with phylogenetic relationships in some lichenized fungi, and a correlation between reproductive mode and environmental conditions has also been observed, we reexamined the relationship between H. hypotrypa and H. flavida. We reassessed phenotypic characters and generated molecular data to delimit species boundaries in this group. The phenotypic analysis was based on morphological, anatomical, and chemical characters, whereas the molecular data included sequences from the nuclear ribosomal internal transcribed spacer region (ITS) and two protein-coding nuclear markers, GPD and MCM7. Specifically, our study attempts to assess: (1) whether the presence vs. absence of sorediais diagnostic of two separate lineages in the group, (2) whether presence vs. absence of soredia is related to geography or elevation, and (3) whether other phenotypic characters can be associated with lineages recovered in the phylogenetic analyses.

Phenotypic study
Over 500 specimens of H. hypotrypa and H. flavida collected throughout the species distributions, China (including Taiwan), Japan and Russia, were examined for this study, also including the lectotype of H. hypotrypa (BM) and holotype of H. flavida (OSC). No specific permissions were required for these locations/activities. The field studies didn't involve endangered or protected species. A Geographic Information System (GIS) was used to analyze the geographic distribution of H. hypotrypa and H. flavida, based on the locality information of examined specimens.

DNA extraction, PCR amplification, and sequencing
Seventy-four lichen specimens of seven Hypogymnia species were selected for DNA extraction based on availability of fresh materials suitable for DNA extraction. The collection information of these specimens is listed in S1 Table, including the latitudes and longitudes of all localities. A total of 62 specimens represented the Hypogymnia hypotropa group were collected from a broad geographic range-China (including Taiwan), Japan and Russia-to ensure the range of phenotypic variation. All sequences used in the analyses were newly generated for this study, except sequences of Arctoparmelia centrifuga, Brodoa intestiniformis, Letharia spp. and Pseudevernia spp. that were chosen as outgroup and downloaded from GenBank.
The extraction procedure followed a modified CTAB method [36]. The internal transcribed spacer of nuclear ribosomal DNA (nrDNA ITS), and fragments of protein-coding genes GPD (glyceraldehyde 3-phosphate dehydrogenase) and MCM7 (mini-chromosome maintenance proteins) were chosen as the genetic markers. Primers used for the PCR amplifications were listed in Table 1.
Reactions were carried out in 50 μl reaction volume and the components used were 3 μl total DNA, 1 μl each primer (10 μM), 25 μl 2×Taq MasterMix (CWBIO, China) and 20 μl ddH 2 O. PCR amplifications were carried out in a Biometra T-Gradient thermal cycler, following conditions: initial heating step for 5 min at 95°C, followed by 35 cycles of 30 s at 94°C, 30 s at 56°C for amplifying ITS and GPD sequences or 54°C for amplifying MCM7 sequences, and 1 min 30 s at 72°C, a final extension step of 8 min at 72°C was added, after which the samples were kept at 4°C. Negative controls were prepared for each amplification series. PCR products were purified using gel purification kit (Shanghai Huashun Bioengineering Corporation, China) following the manufacturer's instructions. PCR products were sequenced using ABI 3730 XL Sequencer by Shanghai BioSune Corporation of China.

Multiple sequence alignments and data analysis
Sequences were aligned using ClustalW Multiple Alignment [42] in BioEdit 7.2.5 [43]. The alignment files were transformed into phylip format in SeaView 4 [44,45]. Pairwise genetic distances were separated into intraspecific and interspecific parameters and calculated to characterize both intra-and interspecific variation within and between H. hypotrypa and H. flavida. Pairwise distances can be viewed as a rough measure for the overall sequence divergence, and an intra-interspecific threshold of ca. 0.015-0.017 substitutions per site has been proposed for species in Parmeliaceae [46]. Pairwise genetic distances were computed for the ITS locus using the general time-reversible model in PAUP Ã [47] for each nominal taxon individually-H. flavida and H. hypotrypa-and all pairwise interspecific comparisons. Genetic distance were then exported from PAUP Ã and the distribution and mean of pairwise distance were calculated. Pairwise distances between different haplotypes were reported as the number of nucleotide substitutions per site (s/s).
Phylogeny of Hypogymnia hypotrypa group. Conflicts were not detected in the three single-gene trees, and the three data sets were concatenated. Phylogenetic analyses were performed using RAxML-HPC BlackBox 8.2.6 [50] and MrBayes 3.2.6 [52,53] on the Cipres Science Gateway (http://www.phylo.org; [51]). In the ML analysis, the default GTR + G model was used as the substitution model with 1000 pseudoreplicates. The data was partitioned according to the different genes. For gpd and MCM7 data were also partitioned by codon position. In the Bayesian analysis, the best model for the three single genes had been found in advance with PartitionFinder v1.1.1 [54]. The ITS and the two protein-coding genes data sets were partitioned by the length of sequences and codon position, respectively. Two parallel

Population genetic analyses and Bayesian clustering
The program SITES [55] was used to assess genetic differentiation and polymorphisms within and between the two traditionally circumscribed taxa in the H. hypotrypa group, the number of fixed differences, shared polymorphisms and pairwise fixation indices (F ST ) [56]. To measure genetic differentiation, we used the program DnaSP V5.10.1 [57]. Furthermore, an intra-interspecific threshold of ca. 0.015-0.017 substitutions per site has been proposed for species in Parmeliaceae [46], and comparisons of pairwise genetic distances were made within and between H. flavida and H. hypotrypa. Bayesian clustering implemented in the program STRUCTURE v.2.3.2 [58,59] was used to assign specimens to genetic clusters. All constant nucleotide position characters in the concatenated multi-locus sequence alignment were excluded, and the data matrix for STRUC-TURE was comprised of only variable nucleotide position characters (SNPs). Indels and 'N's were ignored for the purpose of SNP identification. Individual population assignments were inferred for K values ranging from 1-5; with 10 replicate runs for each K value. Each run consisted of 50,000 burn-in generations, followed by 50,000 iterations using the admixture options.
For ABGD we used default parameters except for using a Pmax at 0.01 and a relative gap width of 1.5, with the model Jukes-Cantor (JC69). The PTP model is intended for delimiting species in single-locus molecular phylogenies, and provides an objective approach for delimiting putative species boundaries that are consistent with the phylogenetic species criteria. We used the bPTP web server (http://species.h-its.org, [67]) to delimit putative species groups using the ITS topology as the input tree and implementing default settings.
We employed the GMYC approach [62,63] to test whether the data support a scenario supporting all samples in the H. hypotrypa/flavida group as belonging to a single species or not. The GMYC method aims at detecting shifts in branching rates between intra-and interspecific relationships. Within a likelihood framework it uses chronograms to compare a null model under which the whole sample belongs to a single species and hence follows a coalescent process and an alternative general mixed Yule coalescent (GMYC) model. The latter combines equations describing branching patterns within and among lineages. A likelihood ratio test (LRT) is used to evaluate whether the null model can be significantly rejected. If the GMYC model fits the data significantly better than the null model, the threshold T allows estimating the number of species present in the data set. The outgroup samples were excluded from the data set. The GMYC analysis based on the ITS sequences was then run online (http://species.hits.org/gmyc/), employing a single and multiple threshold methods.
The multispecies coalescent model implemented in the program BPP v3.2 [64][65][66] was used to assess support for the separation of the sampled Hypogymnia species. BPP incorporates coalescent theory and phylogenetic uncertainty into parameter estimation; and the posterior distribution for species delimitation models is sampled using a reversible-jump Markov Chain Monte Carlo (rjMCMC) method. We used the unguided species delimitation algorithm (' A11'; [68]). This algorithm explores different species delimitation models and different species phylogenies, with fixed specimen assignments to populations. Specimens were assigned to either H. hypotrypa or H. flavida based on the conventional phenotype-based descriptions (sorediate vs. esorediate). The program attempts to merge populations into one species, and uses the nearest neighbor interchange (NNI) or subtree pruning and regrafting (SPR) algorithms to change the species tree topology [69]. We used Prior 0, equal probabilities for the labeled histories, to assign probabilities to the models. Rates were allowed to vary among loci (locus rate = 1), and the analyses were set for automatic fine-tune adjustments. Multiple analyses using different combinations of the theta (θ) and tau (τ) priors spanning a range of possible population sizes and divergence times were performed for each genus. The rjMCMC analysis was run for 200,000 generations, sampling every 2 generations discarding the first 10% as burn-in. Each analysis was run twice using a different search algorithm (algorithm 0 or 1) to confirm consistency between runs. Speciation probabilities greater than 0.95 were considered supported species delimitations.
Given that different species delimitation analyses supported different species scenarios for the H. hypotrya /flavida group (see Results), the most likely hypothesis of species boundaries was inferred by comparing marginal likelihoods using Bayes factor delimitation (BFD) test [70]. The optimal species delimitation scenario was evaluated by comparing marginal likelihoods using the BFD framework described previously [70]. We calculated marginal likelihood estimates (MLEs) for three species delimitation scenarios: (i) assigning specimens within the H. hypotrypa/flavida group to two separate species based on traditional, phenotype-based identifications; (ii) lumping all H. hypotrypa/flavida specimens into a single putative lineage; and (iii) assigning specimens within the H. hypotrypa/flavida group to two separate candidate species-level lineages based on the results of the PTP analysis (see Results). All other sampled Hypogymnia species were assigned species-level membership based on morphological identifications.
For each of the three species delimitation scenarios we reconstructed a species tree using Ã BEAST v1.8.3 [71]. Substitution models for each of the three loci were chosen using Partition-Finder [54], as described above. We selected a birth-death model for the species tree prior; the population size model was set to piecewise linear and constant root. Ã BEAST analyses were performed using 20,000,000 generations, sampling every 1000 generations, and the first 25% of each run was discarded as burn-in. MLEs were estimated using the stepping-stone method [72], with 100 path steps, a chain length of 100,000 generations and likelihoods saved every 100 generations. Bayes factors were then calculated as described by Grummer et al. [70], with 2lnBf >10 being considered as 'decisive' support for a hypothesis.

Phenotypic studies
All sampled specimens from the H. hypotrypa group (H. hypotrypa and H. flavida) were identical in the anatomical structure and chemical substances, both of which developed internally heteromerous thalli: prosoplectenchymatous upper cortex, photobiont layer, medulla and prosoplectenchymatous lower cortex with similar thicknesses. However, in rare instances some lobes tip of H. flavida lacked obvious dorsoventrality (Fig 1D), resulting in two upper cortex layers and two algal layers. In chemistry, the H. hypotrypa group constantly contained usnic acid, physodalic acid, and protocetraric acid; some also contained 3-hydroxyphysodic acid. The only apparent phenotypic differences between H. flavida and H. hypotrypa were in regards to lobe morphology and presence of soredia. Although soredia were present in all H. hypotrypa specimens, in many cases, the soredia were distributed along the cracks of the upper surface and hence could easily be overlooked (Fig 1C). Compared with H. hypotrypa, H. flavida had a broader range of variation in lobe morphology. In addition to the broad and richly branched lobes typical of H. hypotrypa (Fig 1A-1E), the lobes of H. flavida were occasionally found to be fingerlike and sparsely branched (Fig 1D). Production of apothecia was observed in both H. hypotrypa and H. flavida (Fig 1E and 1F).

Geographic distribution
Both H. hypotrypa and H. flavida usually grow in alpine to montane habitats in eastern Asia, although each species is known to occur across a broader altitudinal range. Based on the analysis of over 500 specimens in our study and in agreement with previous results [33], H. hypotrypa has a broader geographic distribution and wider altitudinal range. Hypogymnia flavida can be found between 2150 m to 4300 m, while H. hypotrypa is found at an altitude between 65 m to 4300 m. We confirm the occurrence of H. flavida in China (including Taiwan), and H. hypotrypa in China, Japan and Russia (Fig 2). Hypogymnia hypotrypa has also been reported from Taiwan and North Korea [33,[73][74][75][76][77], but we have not seen this material and thus cannot confirm the identity of these collections.

Genetic differentiation and Bayesian clustering
No fixed differences in nucleotide positions were observed between H. hypotrypa and H. flavida in any of the three sampled loci (Table 2). F ST indices were calculated to assess the degree of genetic isolation within H. hypotrypa group, which can vary from 0 (complete panmixis) to 1 (complete isolation between populations). In our study, the values of F ST were relatively low, ranging from 0.035 to 0.276. The shared polymorphisms revealed 5-10 nucleotide shared by Results from the Bayesian clustering analysis of the H. hypotrypa/flavida group performed under the assumption of two distinct populations are shown in Fig 4. Specimens representing each traditionally circumscribed species were recovered in two distinct genetic clusters-'cluster 1' and 'cluster 2', with approximately 10% of samples specimens with evidence of admixed genomes. However, the majority of H. flavida specimens were assigned membership to 'cluster 1', while those identified as H. hypotrypa were assigned membership to 'cluster 2' (Fig 4). The information of samples from different localities assigning to 'cluster 1' and 'cluster 2' is seen in S1 Table.

Phylogenetic analyses
Single-locus maximum likelihood (ML) topologies and ML and Bayesian trees inferred from the concatenated, 3-locus data set (1859bp) are shown in S1-S5 Figs. In order to clearly depict relationships among H. hypotrypa and H. flavida specimens, cartoon topologies of the ITS and concatenated matrix are reported (Figs 5 and 6). Because the topology of ML and Bayesian trees are highly similar, the posterior probability values above 0.5 are noted directly after the bootstrap values at the nodes of the ML tree (Fig 6). The H. hypotrypa group formed a wellsupported clade (BS = 90, 100 and PP = 1) and was comprised of closely related specimens distinct from other sampled species of Hypogymnia species (Figs 5 and 6).
Within the H. hypotrypa/flavida clade, several samples of H. hypotrypa (blue branches) or H. flavida (Figs 5 and 6; yellow branches) clustered into small sub-branches, then intermixed each other. In some cases, samples of H. hypotrypa and H. flavida clustered together, forming separated sub-clades (orange branches). No formation of two well defined separated bigger clades corresponds to H. hypotrypa and H. flavida. No obvious relationship between clades andlarge-scale geographic distribution were seen, although we found small-scale geographic differentiation. For example, some samples from Shaanxi Province of China (CH-Sx, highlighted in red) often formed separated sub-clades. Samples from Japan (JA, highlighted in pale blue) and Russia (RU, highlighted in purple) had a closer relationships to each other than

Species delimitation analyses
The ABGD analysis based on nrDNA ITS provided evidence supporting one species delimitation scenario, e.g. all specimens within the H. hypotrypa/flavida group are inferred as conspecific (P = 0.0021-0.01). The ABGD circumscription of all specimens within the H. hypotrypa/ flavida group as a conspecific OTU was concordant to the lack of deep, well-supported phylogenetic substructure within this clade.
The In the GMYC analyses using the single and multiple threshold methods, the GMYC model was not significantly better than the null model of uniform (coalescent) branching rates. The likelihood ratio for the single threshold method analyzing on the ingroup (H. flavida and H. hypotrypa) was 1.3, and three or four clusters within H. hypotrypa group were included (S7 Fig). More than 10 clusters were shown when the multiple threshold method was performed, which seems rather unreasonable because most of time one cluster was only composed of two or three samples. It has previously been suggested that the single-threshold model outperforms the multiple-threshold version [78], and results from the multiple threshold GMYC analysis were not considered further.
In contrast to the ITS-based species delimitation analyses, the multispecies coalescent species delimitation method BPP provided strong support (posterior probability = 1.0) for the Individual population assignments were inferred using a STRUCTURE analysis of single nucleotide polymorphisms from multi-locus sequence data from specimens identified as H. flavida and H. hypotrypa under a model assuming two genetic groups. The horizontal axis gives specimen numbers. The vertical axis represents the inferred proportion of each individual's genome assigned to a genetic cluster, with assignment probability into the two different genetic clusters depicted with distinct colors-'cluster 1' shown in yellow and 'cluster 2' in blue. Specimens within each taxon are clustered by geographic region (SAA = Shaanxi; YN = Yunnan; XZG = Tibet; SC = Sichuan; TWN = Taiwan; JP = Japan; and RU = Russia). Population assignments for each specimen are reported in S1

Discussion
In this study we used an integrative approach to investigate species boundaries between two closely related lichen-forming fungal species in the genus Hypogymnia, H. flavida and H. hypotrypa. The production of vegetative reproductive propagules in H. hypotrypa, differences in lobe morphology, and variation in geographic distributions have traditionally separated species Species Delimitation in Hypogymnia hypotrypa Group in the H. hypotrypa group. Our morphological analyses of over 500 specimens supported the previous observations that H. hypotrypa and H. flavida differ in the presence or absence of soredia, wide or narrow lobes, and the former had a broader geographic distribution. However, in this study, analyses based on the DNA sequences data failed to provide a consensus on species boundaries in H. hypotrypa/flavida group.
The three species delimitation analyses based on ITS sequence data alone-ABGD, bPTP, and GMYC-indicated multiple scenarios of species boundaries in the H. hypotrypa/flavida group: one being that all members of this group are conspecific (ABGD), while bPTP and GMYC support multiple species-level lineages within this group. However, candidate species circumscribed using bPTP and GMYC did not correspond with the traditional diagnostic character of the presence or absence of soredia.
Although both bPTP and GMYC delimited multiple candidate species within the H. hypotrypa/flavida clade, the supporting evidence was not particularly robust. The Bayesian implementation of PTP provided only weak statistical support for the two species delimited in this group, with posterior probabilities << 0.5 (see Results), and the GMYC model was not significantly better than the null model of uniform branching rates. Similar to the results of the ABGD analysis which suggested H. hypotrypa and H. flavida to be conspecific, F ST values between H. hypotrypa and H. flavida were relatively low, e.g. from 0.035 to 0.276, suggesting little isolation between the two nominal species. Additionally, there were 10 shared polymorphisms at most, supporting the hypothesis that the nominal taxa H. hypotrypa and H. flavida do not form two distinct evolutionarily independent lineages.
In contrast, the coalescent-based species validation method BPP (see Results) and BFD tests (Table 3) provided decisive evidence supported the recognition of H. flavida and H. hypotrypa as distinct separate species. If the independence of these nominal taxa is legitimate, it is not tracked by the ITS marker, the formal barcoding marker for fungi [79], suggesting a recent diversification history for this clade. This result highlights a potential limitation of using single-locus datasets and phylogenetic species recognition criteria for groups with recent diversification histories and incomplete sorting among lineages [80]. However, the relatively high intraspecific genetic distances in both H. flavida and H. hypotrypa, with some pairwise comparisons exceeding the proposed threshold for species in Parmeliaceae [46], suggest the potential for more complex species delimitation scenarios. Recently, phylogenomic data has shown promise in resolving relationships in closely related lichen-forming fungal species groups with recent divergence histories [81], and we anticipate that genome-scale data will provide important insight and resolution into relationships in the H. hypotrypa/flavida group.
Species in the H. hypotrypa group were not recovered as monophyletic in phylogenetic analyses of multilocus DNA sequence data. Additional species delimitation analyses, genetic clustering, and comparisons of genetic differentiation indicated multiple possible scenarios of species boundaries in the H. hypotrypa group, e.g. conspecificity or multiple independent species. This raises the question, what are the existing limitations in delimitating species boundaries using molecular sequence data and phylogenetic analyses and what are the limitations of traditionally diagnostic phenotypic characters? In regards to our initial question about the taxonomic utility of the presence or absence of soredia, our data suggest that differences in reproductive strategies may not correspond to species boundaries with high fidelity (Fig 4; STRUCTURE). In some groups, molecular data suggested that the sorediate and non-sorediate taxa were conspecific, and the sorediate populations usually have a larger range (e.g., Usnea antarctica morph of U. aurantiaco-atra) [23]. Phenotypically, H. hypotrypa and H. flavida differ in the presence or absence of soredia, and they vary in production of wide or narrow lobes. Furthermore, H. hypotrypa has a broader geographic distribution. The geographical ranges of Hypogymnia hypotrypa and H. flavida, however, are more complex with esorediate morphs absent from Russia and Japan and sorediate morphs absent from Taiwan of China. While the former agrees with other studies, the absence of sorediate morphs from Taiwan is difficult to interpret and may be due to the fact that populations belonged to different species. Although H. hypotrypa had not been confirmed in Taiwan, our data indicate that H. flavida fromTaiwan has a closer relationship to H. hypotrypa than to specimens identifiable as H. flavida from other localities (Figs 5 and 6, S1-S5 Figs). This can be interpreted in several ways: (1) H. flavida of Taiwan represents intermediates by introgression between H. hypotrypa and H. flavida, (2) H. flavida from Taiwan is close to the ancestral state at the time of divergence of sorediate and esorediate lineages, but is currently reproductively isolated from both H. hypotrypa and continental H. flavida, or (3) the pattern represents a random deviation in an otherwise panmictic species complex. For any of these three scenarios, one could conclude that H. flavida is conspecific with H. hypotrypa. But both scenarios 1 or 2 are also compatible with a taxonomy that accepts two or more species, using a phylogenetic species concept.
However, this scenario contradicts the results of the BPP species validation analysis and BFD test, which support the traditional recognition of species based on the presence or absence of soredia to delimit the H. hypotrypa group. The presence or absence of soredia may generally correspond to distinct evolutionary lineages, e.g., H. hypotrypa and H. flavida, but may not be a consistent diagnostic feature (see Fig 4). The misspecification of individuals in coalescentbased species delimitation analyses, such as BPP, is not well understood. The strong support in BPP and BFD tests may reflect the general pattern of the presence or absence of soredia in each lineage, rather than an exclusive pattern in each group.
The influence of reproductive mode on distributional ranges of lichens is currently poorly understood [82]. Hypogymnia species with soredia tend to have broader transcontinental ranges than esorediate species [83]. Poelt [84] showed that sorediate populations are generally expected to have higher potential for long-distance dispersal and hence often have larger distributional ranges. The elevation range of the esorediate taxon H. flavida (2150-4300 m) is about half that of the sorediate form H. hypotrypa (65-4300 m), which occurs in high montane to subalpine or alpine habitats. Note that this vertical difference is exactly analogous to the broader distribution observed for sorediate counterparts to fertile species. Higher altitude habitats may be correlated with harsher environmental conditions. Ecological stress, including biotic and abiotic factors, as important correlation factors contributing to genomic and phenomic diversity in nature, and has been shown to bepositively correlated with increased sexuality (by means of meiospores) in lichens and soil microfungi [27,85,86]. Because H. flavida grows exclusively at higher elevations, it would be not surprising having some adaptive phenotype under this ecological stress, such as narrower finger-like lobes, differing from the most common wide lobes of both species, and depending on sexual reproduction but not on vegetative reproduction.
Previous studies of other lichen genera have suggested that some sorediate and esorediate populations likely belong to a single species [14,20,22,23,[87][88][89]. Our data from the H. hypotrypa group suggest a more complex relationship between esorediate and sorediate populations, including the presence of reproductively uniform species being closely related to lineages exhibiting different reproductive modes [90] or the presence of several sorediate populations each representing distinct lineages [91,92]. Despite the fact that in the majority of cases studied using molecular data sorediate and esorediate populations were found to represent variations within one species, no conclusions can be drawn a priori. The lack of a generalizable pattern in the taxonomic utility of differences in reproductive strategies demonstrates that each case requires careful consideration. The genus Hypogymnia is a prime example since it includes distinct lineages characterized by the morphology of soralia [93][94][95].
Supporting Information S1 Fig. The RAxML tree based on nrDNA ITS sequences. The numbers in each node represents bootstrap support value, and the numbers lower than 50 were not shown. The samples marked with '' were downloaded from GenBank, and others were newly generated for this analysis. The number of each sample is listed in S1 Table. (TIF) S2 Fig. The RAxML tree based on GPD sequences. The numbers in each node represents bootstrap support value, and the numbers lower than 50 were not shown. The samples marked with '' were downloaded from GenBank, and others were newly generated for this analysis. The number of each sample is listed in S1 Table. (TIF) S3 Fig. The RAxML tree based on MCM7 sequences. The numbers in each node represents bootstrap support value, and the numbers lower than 50 were not shown. The samples marked with '' were downloaded from GenBank, and other were newly generated for this analysis. The number of each sample is listed in S1 Table. (TIF)

S4 Fig. The
RAxML tree based on 3-loci concatenated sequences. The numbers in each node represents bootstrap support value, and the numbers lower than 50 were not shown. The samples marked with '' were downloaded from GenBank, and others were newly generated for this analysis. The number of each sample is listed in S1 Table. (TIF) S5 Fig. The Bayesian tree based on a concatenated 3-locus data matrix. The numbers in each node represents posterior probability value, and the numbers lower than 0.5 were not shown. The samples marked with '' were downloaded from GenBank, and others were newly generated for this analysis. The number of each sample is listed in S1 Table. (TIF)  Table. Specimens used for DNA extraction and sequences used in this study. (DOC)