Resolving the phylogenetic relationship between Parmotrema crinitum and Parmotrema perlatum populations

Abstract

The widespread species Parmotrema crinitum (Ach.) M. Choisy and Parmotrema perlatum (Huds.) M. Choisy are mainly distinguished by their reproductive strategies. While P. crinitum propagates by isidia, P. perlatum produces soredia. In this study, we aim to evaluate the phylogenetic relationship between both species and to critically examine their species boundaries. To this purpose, 46 samples belonging to P. crinitum and P. perlatum were used in our analysis, including 22 for which we studied the morphology and chemistry, before extracting their DNA. We used 35 sequences of the internal transcribed spacer region of nuclear ribosomal DNA (ITS) of Parmotrema perlatum from Europe and Africa (20 of which were newly generated), and 11 of Parmotrema crinitum from Europe, North America and North Africa (two newly generated). Additionally, 28 sequences of several species from Parmotrema were included in the ITS dataset. The ITS data matrix was analyzed using different approaches, such as traditional phylogeny (maximum likelihood and Bayesian analyses), genetic distances, automatic barcode gap discovery (ABGD) and the coalescent-based method poisson tree processes (PTP), in order to test congruence among results. Our results indicate that all samples referred to P. crinitum and P. perlatum nested in a well-supported monophyletic clade, but phylogenetic relationships among them remain unresolved. Delimitations inferred from PTP, ABGD and genetic distance analyses were comparable and suggested that P. crinitum and P. perlatum belong to the same lineage. Interestingly, two samples of P. perlatum separate in a different monophyletic clade, which is supported as a different lineage by all the analyses.

Introduction

Lichenized fungi form mutualistic relationships with photo-autotrophic organisms (photobionts), mainly green algae (Trebouxiophyceae and Ulvophyceae) and/or cyanobacteria. The lichen symbiosis has been highly successful within fungi, especially Ascomycota, with c. 19 400 currently accepted species (Lücking et al. 2017) and an estimated diversity of more than 28 000 species (Lücking et al. 2009; Leavitt et al. 2013). Additionally, lichens are commonly used to assess environmental disturbance, serving as bioindicators of air pollution, forest age and health, and climate change (Nimis et al. 2002; Crespo et al. 2004; Giordani & Brunialti 2015; Sujetovienė 2015; Sancho et al. 2019; Abas 2021).

Recognizing phylogenetic relationships and delimiting species in lichens are crucial for ecological and conservation studies, assessing biotic diversity, and identifying factors driving diversification. They are also important for future investigations because phylogenetic differences may not be fully reflected in the phenotype.

Traditionally, to infer taxonomic boundaries in lichen-forming fungi, thin-layer chromatography (Culberson 1972), morphology and the expression of signature secondary metabolites, and isolation and identification of lichen substances have been used (Huneck & Yoshimura 1996; Huneck 1999). However, these characters are highly variable and their homology has proved difficult to assess. For example, it has been shown that apothecial form and spore wall thickness have changed in parallel within Pertusaria. The Pertusaria-type ascus is plesiomorphic within the Pertusariaceae and thus cannot be used to circumscribe Pertusaria (Lumbsch & Schmitt 2001). However, with the use of molecular tools, it has been shown that phenotype-based taxonomy may not reflect natural groups, including cases in which morphologically distinct forms formerly recognized as distinct species are shown to represent a polymorphic species (see e.g. Boluda et al. 2019), and other cases where several morphologically cryptic species are masked within a single nominal taxon (Del-Prado et al. 2016).

Within lichen-forming Ascomycetes, Parmeliaceae (Lecanorales) constitutes one of the largest and best-studied families (Crespo et al. 2007, 2011; Thell et al. 2012; Divakar et al. 2015). This family is usually characterized morphologically by a certain type of ascoma ontogeny and the presence of an ascomatal structure known as a cupulate exciple (Henssen & Jahns 1974). Parmeliaceae also comprises species which are frequently used in biomonitoring studies of atmospheric pollution, such as Parmelia sulcata, Flavoparmelia caperata, Parmotrema perlatum and Punctelia subrudecta (Hawksworth & Rose 1970; Crespo et al. 2004; De La Cruz et al. 2018). The application of phylogenetic analysis based on molecular (DNA) characters to delimit species allows us to determine a posteriori which types of phenotypic characters are good predictors of phylogenetic species and demonstrate how these characters evolve in this family and in lichenized fungi in general. These molecular data have led to the recognition of morphologically cryptic species, such as Parmelia serrana (Molina et al. 2004), P. barrenoae (Divakar et al. 2005a), P. encryptata (Molina et al. 2011a) and Melanelixia californica (Divakar et al. 2010), and conversely also to the union of species traditionally regarded as morphologically distinct (see e.g. Boluda et al. 2019). Recently, coalescent-based species delimitation approaches have shown to be well suited to critically evaluate species boundaries in Parmeliaceae, as well as lichens in general (Leavitt et al. 2015). Furthermore, these methods can accurately display relationships when incomplete lineage sorting and gene tree heterogeneity hide phylogenetic relationships among species (Knowles & Carstens 2007; Camargo et al. 2012). Commonly used methods to critically evaluate species boundaries include the poisson tree processes (PTP) model (Zhang et al. 2013), the automatic barcode gap discovery (ABGD) (Puillandre et al. 2012), the general mixed Yule coalescent model (GMYC) (Pons et al. 2006; Monaghan et al. 2009), and SpedeSTEM (Ence & Carstens 2011).

Parmotrema (Massalongo 1860) is one of the largest genera in the parmelioid group of the family Parmeliaceae. It includes c. 300 described species with an apparent centre of speciation in the Pacific Islands, tropical and subtropical regions of South America (Spielmann & Marcelli 2020). The species of the genus are characterized by a pored epicortex, large thalli with broad lobes, a broad, naked marginal zone on the lower surface, and large, thick-walled, ellipsoid ascospores, sublageniform or filiform conidia (Elix 1993), and (commonly) marginal cilia (Hale 1974). Reproductive strategies vary among Parmotrema taxa. Sexual reproduction is restricted to characteristic fungal fruiting bodies (ascomata) producing ascospores. Ascospores are dispersed independently of the photosynthesizing partner (photobiont) and require reacquisition of the appropriate photobiont partner to re-establish the lichenized condition. Species within Parmotrema also commonly propagate asexually by means of vegetative diaspores, either isidia or soredia. These specialized vegetative reproductive propagules contain both fungal and algal symbionts, eliminating the need for independent acquisition of the appropriate photobiont.

In a molecular phylogeny of parmotremoid lichens (Blanco et al. 2005), a single sample of both Parmotrema perlatum (Huds.) M. Choisy and Parmotrema crinitum (Ach.) M. Choisy, originally from Portugal, were included and these formed a well-supported monophyletic group, indicating that these specimens could be conspecific. Therefore, a critical evaluation of species boundaries is necessary.

Parmotrema crinitum is characterized by the presence of coralloid branched, apically ciliate isidia or often eciliate isidia (Divakar & Upreti 2005), and the stictic acid complex in the medulla. According to Elix (1994), P. crinitum is a cosmopolitan species that is widespread in temperate, tropical regions and even in high humidity, sub-boreal forests (Elix 1994; Kurokawa & Lai 2001). Many European countries have reported the presence of P. crinitum (Jablońska et al. 2009) as have some Asian countries such as Japan (Yoshimura 1974), China (Wei 1991) and Taiwan (Wang-Yang & Lai 1976).

Parmotrema perlatum is a greenish grey foliose lichen, saxicolous or corticolous, loosely adnate, with rounded lobes. It can be recognized by its broad lobes, irregularly branched (7.0 cm), laterally overlapping, with frequent black cilia (0.20–)0.50–1.00 × 0.02(–0.10) mm, evenly distributed but less common in the lobe apices. Its soralia are marginal and linear, sometimes widely distributed or subcontinuous, concolorous with the thallus, and the medulla is white. The thallus contains atranorin, stictic acid, hypostictic acid, menegazziaic acid and norstictic acid (Spielmann & Marcelli 2009). Parmotrema perlatum is a widespread species in temperate regions of the Northern and Southern Hemispheres: Asia (Hale 1965; Kurokawa 1991; Kurokawa & Lai 2001), Europe (Hale 1965), Africa (Hale 1965; Swinscow & Krog 1988), North America (Hale 1965; Brodo et al. 2001; Nash & Elix 2002), Central America (Hale 1965), South America (Hale 1965), Australia (Kantvilas 2019) and New Zealand (Blanchon 2013).

Morphologically, P. crinitum and P. perlatum can easily be separated based on the characteristics associated with their different dispersal strategies. Parmotrema crinitum has apically ciliate isidia or often eciliate isidia, and P. perlatum has marginal soralia. Parmotrema perlatum is one of the most common and widely utilized lichens in the Ayurvedic system of medicine and has been overexploited for uses in traditional medicine in India (Kumar & Upreti 2001). In India, this species is currently considered threatened (Divakar & Upreti 2005).

The aim of this study was to evaluate the phylogenetic relationships between P. crinitum and P. perlatum, and elucidate the possible monophyly between both species. Additionally, PTP, ABGD and genetic distance analyses were included to critically assess species boundaries. The ITS data matrix consisted of a total of 74 sequences of Parmotrema species, including 11 of P. crinitum and 35 of P. perlatum, plus two outgroup sequences of Crespoa carneopruinata. Morphological and chemical features of all samples were critically examined.

Material and Methods

Taxon sampling

Sequence data of the ITS locus were analyzed from 76 specimens, of which 11 sequences were of Parmotrema crinitum (two individuals were newly sequenced and nine sequences downloaded from GenBank) and 35 of Parmotrema perlatum (20 individuals were sequenced in this study and 15 sequences downloaded from GenBank). The new specimens were collected from distant geographical regions throughout the species distributions (Table 1). Additionally, 30 ITS sequences from 30 specimens, belonging to 14 different species (6 were newly sequenced and 24 sequences downloaded from GenBank), are included in this study to test the monophyly of both species within the genus. This total included Crespoa carneopruinata which was selected as outgroup since it has previously been shown to be the sister group of Parmotrema (Divakar et al. 2015). Details of the material studied, including GenBank Accession numbers, are shown in Table 1. Species recognition was based mainly on morphological and chemical characters.

Table 1. Parmotrema species used in this study with GenBank Accession numbers of the ITS sequences and voucher information for the specimens. Crespoa carneopruinata was included in the phylogenetic analyses as outgroup. * = newly generated sequences from specimens for which morphological and chemical characters are given in Table S1 (available online).

Chemistry and morphology

Thallus morphology of all new specimens of Parmotrema crinitum and P. perlatum included in the molecular analyses was studied using a Nikon SMZ-1000 stereomicroscope to identify morphological characteristics of the thallus: lobe shape, isidia, soralia, cilia and rhizines. This is because P. crinitum and P. perlatum are traditionally differentiated based on these features. Photographs were taken with a Nikon 105 mm f/2.8D AF Micro-Nikkor lens coupled to a Nikon D90 camera.

Spot tests were carried out on the medulla with usual chemical reagents such as aqueous potassium hydroxide (K), Steiner's stable paraphenylenediamine (PD) and aqueous calcium hypochlorite (C).

Thin-layer chromatography (TLC) was carried out following Orange et al. (2010). We used TLC solvent system C (200 ml toluene/30 ml acetic acid) (Elix & Ernst-Russell 1993), with concentrated acetone extracts at 50 °C spotted onto silica gel 60 F254 aluminum sheets (Merck, Darmstadt). The aluminum sheets were dried for 10 min in an acetic acid atmosphere to maximize resolution.

DNA extraction, PCR and sequencing

Total DNA was extracted from a single, clean (under a dissecting microscope) lichen lobe using the DNeasy Plant Mini Kit (Qiagen, Barcelona) with a slight modification to the manufacturer's instruction (Crespo et al. 2001). Genomic DNA (5–25 ng, quantified using a quantitative PCR machine) was used for PCR amplifications of the ITS region. Standard PCR amplifications were conducted in 25 μl reaction volumes containing 12.5 μl of Master Mix (50 units/ml of Taq DNA polymerase supplied in a proprietary reaction buffer (pH 8.5), 400 μM dATP, 400 μM dGTP, 400 μM dCTP, 400 μM dTTP, 3 mM MgCl₂), and 1.5 μl of each primer at 10 μM, 4.5 μl of water (H₂O) and 5 μl of DNA template. Fungal nuclear internal transcribed spacer (ITS) rDNA was amplified using the primer pair ITS1F (5ʹ [CTT GGT CAT TTA GAG GAA GTA A] 3ʹ) (Gardes & Bruns 1993)/ LR1 (5ʹ [GGT TGG TTT CTT TTC CT] 3ʹ) (Vilgalys & Hester 1990). We also tried primer pair ITS1-LM (5ʹ [GAA CCT GCG GAA GGA TCA TT] 3ʹ) (Myllys et al. 2001) and ITS2-KL (5ʹ [ATG CTT AAG TTC AGC GGG TA] 3ʹ) (Lohtander et al. 1998), but we had better results with primer pair ITS1F/LR1. For any failed samples, we tried a nested PCR using the primer pair ITS1F/LR1 for the first PCR, and 5 μl of this PCR product as DNA template for the second PCR using ITS1-LM/ITS2-KL. With this nested PCR, we aimed to amplify a smaller region than with the previous primers. However, it did not yield any additional products.

Amplification was run in an automatic thermocycler (Techne Progene, Jepson Bolton & Co., Watford, Herts, UK) using the following parameters: initial denaturation at 94 °C for 5 min followed by 35 cycles at 94 °C for 1 min, 54 °C (ITS1F/LR1) for 1 min, and 72 °C for 1 min 30 s, with a final extension at 72 °C for 10 min. Amplification products were visualized on 1% agarose gel stained with SYBR green DNA (Life Technologies Corporation, Grand Island, New York, USA), and subsequently purified using ExoSAP-IT (GE Healthcare, Chalfont St Giles, UK) according to the manufacturer's instructions. Sequencing was performed using BigDye Terminator reaction kit (ABI PRISM, Applied Biosystems, Waltham, Massachusetts, USA). Cycle sequencing reactions were performed with the same sets of primers used in the amplification step. Sequencing reactions were electrophoresed on a 3730 DNA Analyzer (Applied Biosystems) at the Unidad de Genómica (Parque Científico de Madrid).

Sequence alignments and phylogenetic analyses

Sequence fragments generated for this study were assembled and edited using the program SeqMan v. 7 (Lasergene R, DNASTAR, Madison, Wisconsin, USA). Sequence identity was confirmed using the megaBLAST search function in GenBank. For the alignment we prepared a matrix of 480 base pairs (bp). Then we used the program MAFFT v. 7 (Katoh & Standley 2013) with the parameters set to default. If sequences had different lengths, only the part shared by all the sequences was used, and after manually removing ambiguously aligned nucleotide positions, we kept 451 unambiguously aligned base pairs for the final matrix that we used as input for the phylogenetic reconstruction. The alignment was analyzed using maximum likelihood (ML) and a Bayesian Markov chain Monte Carlo (B/MCMC) approach. The ML analysis was performed using an online version of the program RAxML-HPC BlackBox v. 8.2.8 (Stamatakis 2006; Stamatakis et al. 2008), implemented on the Cipres science gateway (https://www.phylo.org/portal2/home.action; Miller et al. 2010). We used a GTRGAMMA model, which includes a parameter (Γ) for rate heterogeneity among sites but chose not to include a parameter to estimate the proportion of invariable sites (Stamatakis 2006; Stamatakis et al. 2008). Support values were assessed using the ‘rapid bootstrapping’ option with 1000 replicates.

For the Bayesian reconstruction, we used MrBayes v. 3.2.1 (Ronquist & Huelsenbeck 2003). The analysis was performed assuming a discrete gamma distribution with six rate categories (GTR + G). The nucleotide-substitution model parameters were selected using the Akaike Information Criterion as implemented in jModelTest (Posada 2008), molecular clock not assumed. A run with 10 million generations, starting with a random tree and employing 12 simultaneous chains, was executed. Trees were saved to a file every 200th generation. The first 2 million generations (i.e. 20 000 trees) were deleted as the ‘burn-in’ of the chains. We plotted the log-likelihood scores of sample points against generations using Tracer v. 1.5 (Rambaut & Drummond 2007) and determined that stationarity had been achieved when the log-likelihood values of the sample points reached an equilibrium value (Huelsenbeck & Ronquist 2001). The trees obtained before stationarity were discarded. Posterior probabilities (PPs) were obtained from the 50% majority-rule consensus of sampled trees after excluding the initial 20% as burn-in. Only clades that received bootstrap support ≥ 70% in the ML analyses and PP ≥ 0.95 were considered strongly supported. The phylogenetic tree was drawn using FigTree v. 1.2.3 (Rambaut 2009) and modified with CorelDRAW v. 8.

Candidate species identification

In order to establish candidate species limits in the phylogenetic tree, three computational approaches were used:

1. 1) Poisson tree processes (PTP): this method does not require an ultra-metric tree, as the transition point between intra- and inter-specific branching rates is identified directly using the number of nucleotide substitutions (Zhang et al. 2013). PTP incorporates the number of substitutions in the model of speciation and assumes that the probability of a substitution giving rise to a speciation event follows a Poisson distribution. The branch lengths of the input tree are supposed to be generated by two independent classes of Poisson events, one corresponding to speciation and the other to coalescence. The ML phylogeny of the ITS dataset obtained with RAxML was used as the input trees to run PTP species delimitation analysis on the PTP web server (http://species.h-its.org/ptp/). We ran the PTP analysis for 100 000 MCMC generations, with a thinning value of 100 and a burn-in of 25%. Outgroup taxa were removed for species delimitation.

2. 2) Automatic barcode gap discovery (ABGD): this is an automatic procedure that sorts the sequences into hypothetical species based on the barcode gap. This method automatically finds the distance where the barcode gap is located (Puillandre et al. 2012). The ABGD method was carried out on the ITS dataset using the web interface at http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html. Default parameters were chosen using Kimura 2-parameter (K2P) distances that correct for transition rate bias (relative to transversions) in the substitution process (Kimura 1980). The default for the minimum relative gap width was set to different values between 0.1 and 0.15.

3. 3) Genetic distances: pairwise ML distances (given as the number of nucleotide substitutions per site) among the ITS rDNA sequences of Parmotrema crinitum and P. perlatum were calculated. Genetic distances were calculated with TREE-PUZZLE v. 5.2 (Schmidt et al. 2002) using the GTR model of nucleotide substitution, assuming a discrete gamma distribution with six rate categories. The program generates an output file which consists of a triangular matrix with all pairwise distances between all the samples included. This matrix was visualized with Microsoft Excel 2010 and genetic distances between different specimens of P. crinitum and P. perlatum were manually identified. Candidate species were proposed based on the threshold of 0.016 substitutions per site (s/s) which separates intra- and interspecific distances in parmelioid lichens (Del-Prado et al. 2010). Distance values in the matrix ≤ 0.016 s/s have been considered as the values between samples of a single species. By using the Excel filter, we separated values ≤ 0.016, providing for every specimen included in the analysis a group of specimens that share the values characterizing its species range.

Results

Morphological and chemical analyses

Morphological and chemical characters of the specimens (P. crinitum and P. perlatum) newly generated and included in the phylogenetic analysis are presented in Supplementary Material Table S1, available online. Parmotrema crinitum and P. perlatum could not be distinguished based on the colour of the thalli, presence of cilia (Fig. 1A & B) and colour of the lower surface (Fig. 1E & F). However, the reproductive structures allowed unequivocal separation into the two species; Parmotrema perlatum has soralia (Fig. 1C) and P. crinitum has isidia (Fig. 1D).

Fig. 1. Main morphological characters of P. perlatum (A, C, E) and P. crinitum. (B, D, F), as indicated with arrows. A & B, black, simple cilia. C, marginal soralia concolorous with the thallus. D, isidia with apical cilia. E & F, black lower surface with brown naked zone peripherally. In colour online.

The TLC test showed that both species contained the stictic acid complex (stictic acid, menegazzic acid and hypostictic acid) and atranorin. The spot tests of P. crinitum were similar to those of P. perlatum (Fig. 2, Supplementary Material Table S1). However, the samples Parmotrema perlatum 34 and P. perlatum 35 lacked menegazzic acid (nos 10 and 12 in Fig. 2). The spot test was congruent with other samples of P. perlatum, and the morphology of both samples showed the same characteristics of P. perlatum, although the cilia were not as abundant as is usual in this species.

Fig. 2. Thin-layer chromatography profile of Parmotrema species (included in the DNA analysis) in solvent system C. P is the control, Pleurosticta acetabulum. Lanes 1 & 7, P. dilatatum. Lane 2, P. flavotinctum. Lane 3, P. mellissii. Lane 4, P. paulense. Lanes 5, 6, 8, 11 & 19, P. crinitum. Lane 9, empty. Lanes 10 & 12, Parmotrema perlatum. Lanes 13, 14, 15, 16, 17, 18 & 20, P. perlatum; these have the same TLC result as P. crinitum but were distinguished from each other based on the morphology. a = norstictic acid; b = atranorin; c = protocetraric acid; d = unknown; e = hypostictic acid; f = stictic acid; g = menegazzic acid. In colour online.

Phylogenetic analysis

The ITS data matrix consisted of a total of 74 sequences of Parmotrema species, including 46 of P. crinitum and P. perlatum (11 sequences of P. crinitum and 35 sequences of P. perlatum), with two sequences of Crespoa carneopruinata were included as outgroup (Table 1). The final matrix used as input for the phylogenetic reconstruction contained 451 unambiguously aligned base pairs (bp). The tree reconstruction (Fig. 3) comprised the following species of Parmotrema, included to test the monophyly of P. crinitum and P. perlatum: Parmotrema clavuliferum, P. dilatatum, P. flavotinctum, P. grayanum, P. haitiense, P. hypoleucinum, P. mellissii, P. paulense, P. perforatum, P. pseudoreticulatum, P. reticulatum, P. robustum, P. subtinctorium and P. tinctorum.

Fig. 3. Phylogenetic relationships between Parmotrema crinitum and P. perlatum. The green (thickened) line indicates the clade supported by maximum likelihood (ML) and Bayesian analyses; * = clades supported only by ML. ** = clades supported only by Bayesian analysis. Alphanumeric labels indicate clades and clusters. Species delimitation scenarios obtained from PTP, ABGD and genetic distances are indicated in columns to the right. In colour online.

The phylogenetic trees estimated from ML and Bayesian methods did not show any well-supported conflict; ML topologies are presented with bootstrap and B/MCMC analysis with posterior probability (ML bootstrap ≥ 70%; PP ≥ 0.95 in B/MCMC analysis) (Fig. 3). The LnL value was −1938.068210 for ML, and −2184.428 for the Bayesian analysis. The phylogenetic analysis showed that all samples of Parmotrema included in this analysis formed a monophyletic group. Within this genus the samples of Parmotrema crinitum and P. perlatum were grouped in one monophyletic group (clade A), with the exception of two samples of P. perlatum that were grouped in a separate clade B. While clade A comprised specimens of P. crinitum and P. perlatum collected from various geographical regions, clade B included only two specimens of P. perlatum collected from Tenerife (Canary Islands) and Lisbon (Portugal).

Although four well-supported monophyletic clusters can be recognized in clade A, the phylogenetic relationships among them remained unresolved (Fig. 3).

Cluster A1 (Fig. 3) contained 15 specimens of Parmotrema perlatum, most of them collected from Morocco, except four from the Canary Islands and two from Portugal and Turkey. Cluster A2 included eight specimens of Parmotrema crinitum from different geographical regions (Table 1). Cluster A3 grouped 18 samples of Parmotrema perlatum from the Canary Islands (Tenerife, Palma, Gomera), and cluster A4 included three samples of P. crinitum from the Canary Islands.

Identifying candidate species

We used the same RaxML tree obtained from the phylogenetic analysis to illustrate the delimitation of putative species recognized by the different approaches conducted with the ITS dataset. ABGD, PTP and genetic distance analyses applied to the ITS dataset detected two candidate species that corresponded to the well-supported clades A (including P. crinitum + P. perlatum) and B (including P. perlatum) obtained in the phylogenetic analysis (Fig. 3).

Discussion

Previous molecular phylogenetic studies have proposed the monophyly of Parmotrema crinitum and Parmotrema perlatum (Blanco et al. 2006; Crespo et al. 2010). The primary goal of the present investigation was to evaluate the phylogenetic relationship between the two species using ML and Bayesian analyses. Additionally, we used different approaches to assess the species boundary, such as genetic distances based on the threshold of 0.015–0.017 s/s established by Del-Prado et al. (2010) to measure intra- and interspecific genetic distances in parmelioid lichens, and separate ranges of intra- and interspecific divergence. This threshold was established using both phylogenetically and morphologically well-delimited species (for details see Del-Prado et al. (2010)). We also used the poisson tree processes (PTP) model (Zhang et al. 2013) and the automatic barcode gap discovery (ABGD; Puillandre et al. 2012). Furthermore, we combined ecological, biogeographical, morphological and chemical data.

Most of the samples belonging to P. crinitum and P. perlatum were recovered in a well-supported monophyletic clade (clade A, Fig. 3), with the exception of two samples of P. perlatum, one from Tenerife (Canary Islands) and one from Lisbon (Portugal), that are separated in the monophyletic clade B, supported by ML and B/MCMC analysis. Morphologically, P. perlatum in clade A has more cilia compared to P. perlatum in clade B.

Thin-layer chromatography (TLC) revealed that the samples in clade B (nos 10 and 12 in Fig. 2) differ chemically from the samples of Parmotrema perlatum grouped in clade A. Several studies have shown that environmental factors, such as light, temperature, pH and culture media, can influence the secondary metabolism in lichens (BeGora & Fahselt 2001). However, we have samples from the same localities and same conditions as the samples in clade A. Furthermore, ABGD, PTP and genetic distance analyses supported clade B as a new sister group to clade A (Fig. 3), suggesting polyphyly, a common phenomenon in Parmeliaceae in general (Lumbsch & Leavitt 2011; Leavitt et al. 2016) and Parmotrema in particular (Divakar et al. 2005b; Del-Prado et al. 2016, 2019; Widhelm et al. 2016).

Clade A was split into 4 groups; however, the phylogenetic relationships among them remained unresolved (Fig. 3). Previous studies have shown that phylogenetic analyses alone are insufficient to explain phylogenetic relationships within Parmeliaceae. For example, based only on maximum parsimony and Bayesian analyses the sorediate P. sulcata was shown to belong to the same clade as the isidiate P. squarrosa (Molina et al. 2004), leaving the authors unable to reach any conclusion regarding species boundaries. Similarly, phylogenetic analyses were insufficient to resolve genetic variability among Parmelia saxatilis specimens; whereas samples from distant geographical regions formed a monophyletic group, samples from neighbouring localities were separated (Crespo et al. 2002; Molina et al. 2011b, 2017). A study on Usnea perpusilla demonstrated that it was necessary to use a combined approach with molecular and morphological data to assess species boundaries in closely related and morphologically variable species (Wirtz et al. 2008). For this reason, coalescent-based species delimitation analyses have been applied with the goal of explaining relationships among clades and delimiting species boundaries (Parnmen et al. 2012; Leavitt et al. 2013).

PTP, ABGD and genetic distance analyses supported clades A and B as putative distinct species. The values obtained were within the interspecific ranges (genetic distances as defined in Del-Prado et al. (2010)), which would support the two samples of P. perlatum in clade B as a separate species, different from clade A. Focusing on clade A, our species delimitation approaches (PTP, ABGD and genetic distance) supported P. crinitum and P. perlatum as one species. However, based only on one molecular marker (ITS) and the various phylogenetic analyses and species delimitation approaches used, can we consider Parmotrema perlatum conspecific with P. crinitum?

Similar cases have been reported in the Parmotrema perforatum species complex (Widhelm et al. 2016), and authors have suggested that the phylogenetic relationships between sexual and asexual populations of this species group could be more complex than previously assumed. They also suggest that traditional tools based on reproductive mode and secondary metabolites (Culberson & Culberson 1973) are no longer the key to identify species such as P. perforatum (Widhelm et al. 2016).

A previous study on Umbilicaria (Ott et al. 2004) focused on Umbilicaria kappenii and U. antarctica, which are distinguished only by their reproductive strategies. Umbilicaria antarctica propagates by thalloconidia and U. kappenii exhibits a variety of asexual propagules: soredia, adventive lobes and sorediate thallyles. To infer phylogenetic relationships between both species, the authors used molecular data from three loci. Results indicated that all samples morphologically referred to U. antarctica and U. kappenii form a monophyletic group and they proposed placing U. kappenii into synonymy with U. antarctica (Ott et al. 2004). In Parmelia, Molina et al. (2011b) rejected the previous hypothesis that P. sulcata and P. squarrosa form a monophyletic group (Molina et al. 2004) and based on phylogenetic analyses and species delimitation approaches confirmed that P. sulcata is not conspecific with P. squarrosa. In addition, P. squarrosa is a reproductively isolated lineage and genetic distances clearly separate this from other Parmelia species (Del-Prado et al. 2010). Within Usnea, studies have suggested that Usnea subfloridana was a secondary species derived from U. florida (Clerc 1984, 1987, 1997; Purvis et al. 1992). Usnea florida and U. subfloridana show many morphological similarities but they have different reproductive strategies. Usnea florida always displays many apothecia and produces no specialized asexual propagules. Usnea subfloridana has soralia, isidiomorphs and occasionally apothecia. Multilocus phylogenetic analyses based on sequences of the ITS, IGS and LSU regions of the nuclear ribosomal DNA and the gene coding for β-tubulin, Mcm7, RPB1 and RPB2 showed that specimens of the two morphospecies formed one monophyletic intermixed group, and not two groups corresponding to morphology (Articus et al. 2002; Mark et al. 2016). This topology was further confirmed with a coalescent-based species delimitation approach (Mark et al. 2016). Authors have suggested that they could be conspecific but taxonomic conclusions must await further study. Moreover, a recent study using RADseq data suggests that closely related lichen species may need genome-wide data to test their species boundaries (Grewe et al. 2018).

Traditionally it was thought that asexually reproducing species in lichens, and in filamentous fungi in general, was an evolutionary dead end (Normark et al. 2003). However, recent molecular studies have demonstrated that lineages with vegetative propagation can also present high genetic diversity (e.g. Parmelia sulcata (Molina et al. 2011b), Parmotrema reticulatum (Del-Prado et al. 2016)). However, even in the absence of sexual reproduction, lichens can exchange photobionts and this process could provide opportunities for gene transfer (Piercey-Normore 2006; Nelsen & Gargas 2008).

While our study confirmed the monophyly of an intermixed clade of Parmotrema crinitum and P. perlatum, the taxonomic conclusion must await additional studies including more markers. The phylogenetic tree lacked ML or B/MCMC support for other widely accepted Parmotrema species, such as P. pseudoreticulatum and P. reticulatum. A more comprehensive taxon sampling and additional molecular markers will therefore be needed before making a formal taxonomic decision on the status of clade B.