In total, 335 herbarium specimens deposited in B, H, HBG, LD, S, UGDA and UPS were used for morphological, chemical and ecological niche modelling (ENM) study: 61 of P. discordans, 113 of P. pinnatifida and 161 of P. omphalodes. A total of 34 specimens were selected for molecular study using the nuclear internal transcribed spacer region (ITS rDNA). Thirty four ITS rDNA sequences of the mycobionts and 17 ITS rDNA sequences of their photobionts were newly generated (Table 1). Additionally, 22 sequences from 10 Parmelia taxa and 67 representative sequences of Trebouxia OTUs, as proposed by Leavitt et al. (2015), were downloaded from GenBank. The specimens deposited in MAF herbarium, which sequences were also used here, have been morphologically and chemically analysed. Newly obtained ITS rDNA sequences were subjected to BLAST search (Altschul et al. 1997) in order to check their identity. All sequences have been deposited in GenBank (see Table 1).

The upper surfaces of all specimens were examined to determine the type of pseudocyphellae orientation such as: only marginal, marginal with few laminal in older parts of thalli and marginal and laminal in young and older parts of thalli. Pseudocyphellae were analysed on the whole thalli surfaces. Moreover, the length (distance between points of lobe branching) and width (distance between two adjacent lobe edges at the point of their branching) of lobes were also measured. Based on morphology and chemistry (see below), the studied specimens were divided into groups, which are characterised in Table 2. From each group (see Table 2) the samples were selected for DNA analysis.

Secondary lichen compounds were identified using thin-layer chromatography (TLC) in solvents A and C (Orange et al. 2001). The presence or absence of fatty acids was checked on two types of TLC plates: glass and aluminium. In order to check the differences in the concentration of lobaric acid in different parts of thalli, samples from marginal and central parts of thalli were analysed using TLC.

Total genomic DNA was extracted using the Sherlock AX Kit (A&A Biotechnology, Poland) in accordance with the manufacturer’s protocol, with slight modifications described by Ossowska et al. (2018).

Fungal ITS rDNA was amplified using the primers ITS1F and ITS4A (White et al. 1990; Gardes and Bruns 1993), while algal ITS rDNA was amplified using the following primers: LR3, ITS4M, ITS1T, ITS4T and AL1500bf (Friedl and Rokitta 1997; Kroken and Taylor 2000; Helms et al. 2001; Guzow-Krzemińska 2006). Amplification was performed in a total volume of 25 μl containing 1.0 μl of 10 μM of each primer, 12.5 μl of Start-Warm HS-PCR Mix Polymerase (A&A Biotechnology, Poland), 1.0 μl of dimethyl sulphoxide (DMSO), 3.0 μl of template DNA (~10–100 ng) and water.

The amplifications were performed in an Eppendorf thermocycler and carried out using the following programme: for fungal ITS rDNA marker: initial denaturation at 94 °C for 3 min and 33 cycles of: 94 °C for 30 sec; annealing at 52 °C for 45 sec; extension at 72 °C for 1 min and final extension at 72 °C for 10 min. For green-algal ITS: initial denaturation at 94 °C for 3 min and 35 cycles of: 94 °C for 45 sec; annealing at 55 °C for 45 sec; extension at 72 °C for 90 sec and final extension at 72 °C for 7 min.

The PCR products were purified using Wizard SV Gel and PCR Clean Up System (Promega, US), according to the manufacturer’s instruction. The cleaned DNA was sequenced using Macrogen sequencing service (http://www.macrogen.com).

The newly generated mycobiont sequences, together with selected representatives of Parmelia spp., were automatically aligned in Seaview (Galtier et al. 1996; Gouy et al. 2010) using the algorithm MUSCLE (Edgar 2004), followed by manual correction and elimination of terminal ends. Then, selection of unambiguously aligned positions was performed using Gblocks 0.91b (Castresana 2000) employing less stringent conditions. The final alignment of mycobionts consisted of 58 ITS rDNA sequences and 444 characters. A sequence of P. sulcata (JN118597) was used as an outgroup.

The newly generated photobiont sequences, together with representative Trebouxia OTUs, downloaded from Dryad database (Dryad Digital Repository) (Leavitt et al. 2015) and described in Leavitt et al. (2015), were automatically aligned using MAFFT – Multiple Alignment using Fast Fourier Transform (Katoh et al. 2002), as implemented in UGENE (Okonechnikov et al. 2012). It was followed with a selection of unambiguously aligned positions using Gblocks 0.91b (Castresana 2000) with less stringent settings (i.e. allowing smaller final blocks, gap positions within the final blocks and less strict flanking positions).

The final alignment of photobionts consisted of 84 ITS rDNA sequences and 580 characters. The names of operational taxonomic units (OTU) for Trebouxia ITS rDNA sequences were given according to Leavitt et al. (2015).

The GTR+I+G best-fit evolutionary model was selected for the mycobiont dataset, based on Akaike Information Criterion (AIC) (Akaike 1973) as implemented in MrModelTest2 (Nylander 2004). For photobionts, we used Partition Finder 2 (Lanfear et al. 2016), implemented at CIPRES Science Gateway (Miller et al. 2010) to determine the best substitution model for each partition under Akaike Information Criterion (AIC) and greedy search algorithm (Lanfear et al. 2012). Two different models were found for partitions, i.e. TRNEF+I+G for 5.8S and GTR+I+G+X for both ITS regions.

Bayesian analysis was carried out using the Metropolis-coupled Markov chain Monte Carlo (MCMCMC) method by using the Markov chain Monte Carlo (MCMC) method, in MrBayes v. 3.2.6 (Huelsenbeck and Ronquist 2001; Ronquist and Huelsenbeck 2003) on the CIPRES Web Portal (Miller et al. 2010) using best models. Two parallel MCMCMC runs were performed, each using four independent chains and 2 million generations for the mycobiont tree and 10 million generations for the photobiont tree, sampling every 1000^th tree. Tracer v. 1.6 (Rambaut and Drummond 2007) was used by plotting the log-likelihood values of the sample points against generation time. Convergence between runs was also verified using the Potential Scale Reduction Factor (PSRF) with all values equal or close to 1.000. Posterior Probabilities (PP) were determined by calculating a majority-rule consensus tree after discarding the initial 25% trees of each chain as the burn-in.

A Maximum Likelihood (ML) analysis was performed using RAxML-HPC2 v.8.2.10 (Stamatakis 2014) with 1000 ML bootstrap iterations (BS) and the GTRGAMMAI model for both analyses.

Phylogenetic trees were visualised using FigTree v. 1.4.2 (Rambaut 2012). Since the RAxML tree did not contradict the Bayesian tree topology for the strongly supported branches, only the latter was shown with the bootstrap support values, together with posterior probabilities of the Bayesian analysis (Figures 1, 2). BS ≥ 70 and PP ≥ 0.95 were considered to be significant and are shown near these branches.

Figure 1. 

Phylogenetic relationships of Parmelia discordans, P. omphalodes and P. pinnatifida, based on Bayesian analysis of the ITS rDNA dataset. Posterior probabilities and maximum likelihood bootstrap values are shown near the internal branches. Newly generated sequences are described with herbarium numbers following the species names. GenBank Accession numbers of sequences downloaded from GenBank follow the species names. Clades with Parmelia discordans, P. omphalodes and P. pinnatifida are highlighted.

Figure 2. 

Phylogenetic placement of Trebouxia photobionts from selected Parmelia spp., based on Bayesian analysis of the ITS rDNA dataset. Posterior probabilities and maximum likelihood bootstrap values are shown near the internal branches. Newly generated sequences are in bold, with collecting numbers preceding the species names. Representative Trebouxia OTUs, as described in Leavitt et al. (2015), were downloaded from Dryad database (Dryad Digital Repository, Leavitt et al. 2015). Clades with photobionts from Parmelia discordans, P. omphalodes and P. pinnatifida are highlighted.

Sequences of ITS rDNA from specimens belonging to P. discordans and P. omphalodes were aligned using Seaview software (Galtier et al. 1996; Gouy et al. 2010) and the terminal ends were deleted. The alignment consisted of 13 sequences and 463 sites. The TCS network (Clement et al. 2002) was created using PopART software (http://popart.otago.ac.nz) (Figure 3).

Figure 3. 

Haplotype network showing relationships between ITS rDNA sequences from Parmelia discordans and P. omphalodes. The names of species are followed with herbarium numbers of specimens or GenBank Accession Numbers. Mutational changes are presented as numbers in brackets near lines between haplotypes.

To evaluate the similarity of niches occupied by all studied taxa, ecological niche modelling (ENM) was applied.

The database of localities of P. discordans, P. omphalodes and P. pinnatifida was compiled, based on information provided on labels of herbarium specimens. The geographic coordinates provided on the herbarium sheet labels were verified. If there were no information about the latitude and longitude on the herbarium sheet label, we followed the description of the collection site and assigned coordinates as precisely as possible to this location. Google Earth (Google Inc.) was used to validate all gathered information. In total, 61 records of P. discordans, 161 of P. omphalodes and 113 of P. pinnatifida were used to perform ENM analysis (Figure 4 and Suppl. material 1: Table S1).

Figure 4. 

Localities of Parmelia discordans (red), P. omphalodes (blue) and P. pinnatifida (green) used in ENM analysis.

The maximum entropy method, as implemented in Maxent version 3.3.2 software, was used to create models of the suitable niche distribution (Phillips et al. 2004, 2006). This application has been proved to provide the most robust response across the number of environmental variables tested (Duque-Lazo et al. 2016) and it has been shown to work better with a small number of samples than with other approaches (Hernandez et al. 2006). MaxEnt settings previously used in research where limited samples were available (e.g. Pietras and Kolanowska 2019) were used in our computations. To assess the high level of specificity of the analysis, the maximum iterations of the optimisation algorithm were established as 10000 and the convergence threshold as 0.00001. The neutral (= 1) regularisation multiplier value and auto features were used. The “random seed” option was used for selecting training points. The run was performed with 1000 bootstrap replications and the default logistic model was used. The Area Under the Receiver Operating Characteristic (AUC) was used to evaluate the reliability of analyses. This is a commonly used threshold independent metric for evaluation of species distribution models (Hosmer and Lemeshow 2000; Elith et al. 2006; Evangelista et al. 2008) which was also used in studies involving a small number of samples (Pietras and Kolanowska 2019). Using more specific metrics, which could evaluate the possible overfitting of the model, would require implementing absence points and, in the case of our study object, such a dataset could not be prepared due to the lack of comprehensive studies on the distribution of genus representatives.

Twelve bioclimatic variables in 2.5 minutes developed by Hijmans et al. (2005; http://www.worldclim.org) were used as input data (Table 3). The study area which was used to evaluate the global identity of niches occupied by P. discordans, P. omphalodes and P. pinnatifida extended from 86.583°N to 17.83°N. As some previous studies (Barve et al. 2011) indicated that usage of a restricted area in ENM analysis is more reliable than calculating habitat suitability on the global scale, the similarity of niches occupied in America was calculated for an area that extended from 180°W to 31.749°W and from 85.292°N to 17.833°N and the study area of all three species occurring in Eurasia was reduced to 84.83–17.83°N and 17.833°W-180°E.

The differences amongst the niches occupied by the populations of three studied lichens were evaluated using the niche identity indices: Schoener’s D (D) and I statistic (I) as available in ENMTools v1.3 (Schoener 1968; Warren et al. 2008, 2010). Additionally, the predicted niche occupancy (PNO) profiles were plotted to visualise differences in the preferred climatic factors amongst all taxa. PNO integrates species probability (suitability) distributions derived with MaxEnt with respect to a single climatic variable (Heibl and Calenge 2015).

Principal components analysis (PCA) was performed to explain the general variation pattern amongst the studied species, based on 12 bioclimatic factors used in ENM analysis. Statistical computations were performed with the programme PAST v. 3.0 (Hammer et al. 2001).

Trees of similar topologies were generated using the maximum likelihood method (RaxML; best tree likelihood LnL = −1512.540166) and the Bayesian approach (BA; harmonic mean was −1667.09). The Bayesian tree is presented in Figure 1 with added bootstrap supports from the RaxML analysis and posterior probabilities from the BA. The phylogenetic analyses showed that, despite morphological similarities of species, the P. omphalodes group is not monophyletic. Specimens are separated into three distinct clades. One clade (0.99 PP) is related to P. imbricaria Goward et al. (Figure 1). In this clade, specimens containing salazinic acid, but variable in fatty and lobaric acids content (Table 2), are grouped with sequences labelled as P. pinnatifida, downloaded from GenBank. Analysis of morphological features revealed that all specimens in this clade have predominantly marginal pseudocyphellae. Specimens with similar chemical variation (Table 2), but having both marginal and laminal pseudocyphellae and, thus, referable to P. omphalodes, form two distinct clades (Figure 1), one containing the majority of the studied specimens and also the sequences downloaded from GenBank (1 PP and 79 BS) and the second (1 PP and 95 BS) grouping only two samples (specimens S-F238139 and UGDA L-23632). The latter clade consists of specimens indistinguishable in all morphological and chemical features from other specimens of P. omphalodes used in this study. This lineage may represent a cryptic species, but more specimens and additional molecular markers are necessary to be analysed before it is described.

Within the larger clade of P. omphalodes, four sequences obtained from specimens containing protocetraric acid and determined as P. discordans are nested. Three of those specimens form a highly supported lineage (1 PP and 93 BS), while the fourth sample of P. discordans is placed outside this subclade (Figure 1). Moreover, to better understand the phylogenetic position and genetic variation of the ITS rDNA marker within P. omphalodes s.l., we generated a haplotype network for specimens of both P. discordans and P. omphalodes (Figure 3). There is no significant difference between specimens of those two taxa, except two samples of P. omphalodes (specimens S-F238139 and UGDA L-23632) representing the second lineage found in our study (see above), that differ from other representatives of this species in at least 10 sites. One specimen of P. discordans (S-F252494) shares the same haplotype with P. omphalodes (AY036998), which differs from other haplotypes of the former taxon in 5 sites. Moreover, three other specimens of P. discordans share the same haplotype, which differs from haplotypes of P. omphalodes in at least 3 positions.

So far, the taxonomy of P. omphalodes group was unclear. Kurokawa (1976) recognised three species within this group: P. discordans, P. omphalodes and P. pinnatifida, whereas Skult (1984) classified P. discordans and P. pinnatifida as subspecies within P. omphalodes. On the other hand, Hale (1987) recognised two species, P. discordans and P. omphalodes. However, our results agree to a certain point with those presented by Molina et al. (2004) and Thell et al. (2008), who showed that P. pinnatifida is a taxon well-separated from P. omphalodes. In the case of P. discordans, Thell et al. (2008) used only a single sequence of this species (AY583212), which was nested within the P. omphalodes clade. In the discussion, those authors concluded that the status of P. discordans as a separate taxon required further molecular analyses (Thell et al. 2008). In our study, sequences of P. discordans are also nested in the clade of P. omphalodes. Perhaps the former should be synonymised with P. omphalodes, as some specimens of both taxa share the same ITS rDNA haplotypes (Figure 3). However, the final conclusions should await more data from other molecular markers as the use of a single genetic marker to delimit species might be inappropriate (e.g. Leavitt et al. 2011, 2013a; Pino-Bodas et al. 2013). However, in the case of many taxonomic groups, ITS rDNA helps to discriminate species, for example, in Parmeliaceae, including Parmelia, and has been shown to be effective and proposed to be used as a primary fungal barcode (e.g. Crespo and Lumbsch 2010; Leavitt et al. 2014; Divakar et al. 2016; Corsie et al. 2019).

The distinguishing character between P. omphalodes and P. pinnatifida is the development of pseudocyphellae; however, the determination of the type and orientation of pseudocyphellae requires checking of the entire thallus surface, not only marginal or central parts of the thalli. We concluded that P. pinnatifida has mostly marginal pseudocyphellae forming white rims around lobes margins (Figure 5C), in some samples with few laminal ones in older parts of thalli. Laminal pseudocyphellae, in this species, predominantly start at the edge of lobes and are connected to the marginal pseudocyphellae and only very few are separated from the marginal ones (Figures 5C, D). Thalli of P. omphalodes always have marginal and laminal pseudocyphellae and, in the case of the latter, many are not connected to the margins of lobes (Figure 5B). We also checked the orientation of pseudocyphellae in P. discordans. In young thalli, they may be exclusively marginal, but in most cases laminal ones are also developed (Figure 5A), as in the case of P. omphalodes.

Figure 5. 

A Parmelia discordans, with marginal and laminal pseudocyphellae, laminal pseudocyphellae mostly not connected with marginal ones (S F-252494) B P. omphalodes, with marginal and laminal pseudocyphellae, laminal pseudocyphellae mostly not connected with marginal ones (S F-252845) C P. pinnatifida, with marginal pseudocyphellae (UGDA L-24298) D P. pinnatifida, with marginal and laminal pseudocyphellae, laminal pseudocyphellae starting predominantly from pseudocyphellae formed at the edge of lobes (S F-239397). Scale bars: 200 μm (A, B, D), 150 μm (C).

The presence of lobaric and fatty acids cannot be treated as diagnostic for the separation of P. omphalodes and P. pinnatifida, as it does not correspond with molecular data. Until now, P. pinnatifida was characterised as a species lacking lobaric acid (Kurokawa 1976; Skult 1984; Molina et al. 2004; Ossowska and Kukwa 2016). In this study, the specimens with morphology of pseudocyphellae typical for this species and with or without lobaric acid are grouped in one clade. The same variation in the presence of lobaric acid was noted in P. omphalodes, which was reported as constantly containing this substance (Kurokawa 1976; Skult 1984; Ossowska and Kukwa 2016). A similar issue was noted in the P. saxatilis group. The presence or absence of lobaric acid was treated as a diagnostic character to differentiate species (e.g. Feuerer and Thell 2002; Molina et al. 2004; Thell et al. 2011; Ossowska et al. 2014), but the recent results obtained by Thell et al. (2017), Ossowska et al. (2018), Corsie et al. (2019) and Haugan and Timdal (2019), revealed that the production of this substance is variable, for example, P. serrana A. Crespo et al., typically lacking lobaric acid, may also produce this substance (Ossowska et al. 2018; Corsie et al. 2019; Haugan and Timdal 2019). Similar variation in lobaric acid production was also observed in Stereocaulon condensatum Hoffm. (Oset 2014). Moreover, lobaric acid was detectable in P. omphalodes and P. pinnatifida only when lobes from the central parts of the thalli were taken for TLC.

Kurokawa (1976) reported that P. omphalodes and P. pinnatifida also differ in the production of fatty acids (absent in P. omphalodes, present in P. pinnatifida), but both species also showed intraspecific variation in this character (Table 2). Moreover, the detection of fatty acids may differ due to the type of TLC plates used. The glass TLC plates are better suited for the detection of these substances than aluminium plates (Orange et al. 2001) and, for example, protolichesterinic acid was undetectable on aluminium plates, but visible on glass plates.

Morphological and chemical characteristics of all taxa of the group are summarised in Table 4 and the determination key is presented below (see also Table 2).

Trees of similar topologies were generated using maximum likelihood (RaxML; best tree likelihood LnL = -7013.073328) and Bayesian analysis (BA; harmonic mean was -6996.31). The Bayesian tree is presented in Figure 2 with added bootstrap supports from RaxML and posterior probabilities from BA. The phylogenetic analyses showed that photobionts of P. discordans, P. omphalodes and P. pinnatifida belong to the Trebouxia S clade (T. simplex/letharii/jamesii group) sensu Leavitt et al. (2015) and represent at least five different lineages (Figure 2). The most common photobiont in the species analysed in this work is Trebouxia OTU S02, which was found in one specimen of P. discordans and most specimens of P. pinnatifida (Figure 6). Additionally, we detected Trebouxia OTU S04 in a single specimen of P. pinnatifida (UGDA L-24293) and one specimen of this species (S-F252763) has an unnamed Trebouxia species (SUn2). Therefore, P. pinnatifida associates with at least three different photobiont taxa of which, based on the BLAST search, OTU S04 seems to be very rare. We also found some variation in photobionts of P. omphalodes which associates with two lineages of Trebouxia, i.e. OTU S05 (two specimens) and an unnamed Trebouxia lineage (three specimens) (SUn1), closely related to the photobiont present in one sample of P. pinnatifida (S-F252763). Moreover, Trebouxia OTU S05 was also detected in P. discordans. In Leavitt et al. (2015), it was reported that, based on 98% sequence similarity, Parmelia species form associations with Trebouxia OTU I02, belonging to the T. impressa/galapagensis group, but this group of photobionts might only be characteristic for P. saxatilis and P. sulcata groups, as we have not found this lineage in the studied specimens.

Figure 6. 

Association network between lichen mycobionts of P. omphalodes group (i.e. Parmelia discordans, P. omphalodes and P. pinnatifida) and photobiont OTUs. The line width is proportional to the number of specimens forming the association with the particular OTU. SUn1 and SUn2 represent unnamed lineages of Trebouxia belonging to clade S.

According to Beck et al. (2002), ‘selectivity’ refers to the taxonomic range of partners that are selected by one of the bionts, while ‘specificity’ should be used for the symbiotic association and depends on the range and taxonomic relatedness of acceptable partners. Lichens with high selectivity may associate with a limited number of photobionts. Numerous mycobionts, belonging to Parmeliaceae, have been shown to associate with identical species of Trebouxia, while others exhibited higher photobiont flexibility (e.g. Kroken and Taylor 2000; Ohmura et al. 2006, 2018; Doering and Piercey-Normore 2009; Leavitt et al. 2013b, 2015; Lindgren et al. 2014). Our results indicate that taxa from P. omphalodes group are moderately selective in their photobionts choice, as these taxa associate with at least two or three Trebouxia lineages (Figure 6).

Lichens that reproduce sexually via independent dispersal of fungal spores, undergo a process of re-lichenisation. This means that the germinating spore of the mycobiont can easily exchange its autotrophic partner, in contrast to asexually reproducing lichens distributing both partners together, which allows continuation of the symbiosis without the need to re-associate with another biont (Beck et al. 1998, 2002; Romeike et al. 2002; Sanders and Lücking 2002). However, even asexually reproducing lichens, such as the Lepraria species, have been shown to switch their algal partners (Nelsen and Gargas 2008). Moreover, in populations of Physconia grisea (Lam.) Poelt with a vegetative propagation strategy, mycobionts associate with more than one photobiont genotype (Wornik and Grube 2010). It was also reported that both sexual and vegetative reproduction allows lichens to generate almost the same amount of diversity to adapt to their environments (Cao et al. 2015). Moreover, Protoparmeliopsis muralis (Schreb.) M.Choisy, which does not produce vegetative propagules, exhibited a low selectivity level (Guzow-Krzemińska 2006; Muggia et al. 2013); however, P. muralis has wider geographical distribution and occurs on a wider range of substrata and ecological conditions than taxa from the analysed group.

The ecological ’lichen guilds‘ hypothesis, i.e. communities of lichens growing on the same type of habitat and forming associations with the same photobiont species, have been proposed for cyanobacterial lichens (Rikkinen et al. 2002). This hypothesis was tested by Peksa and Škaloud (2011) for the eukaryotic genus Asterochloris Tschermak-Woess. These authors showed that ecological niches available to lichens may be limited by algal preferences for environmental factors and thus can lead to the existence of specific lichen guilds, but their results were based only on selected species of Lepraria Ach. and Stereocaulon Hoffm. On the other hand, results obtained by Leavitt et al. (2015) indicated that ecologically specialised lichens from different genera form associations with different Trebouxia OTUs in the same habitat. Moreover, observations made by Deduke and Piercey-Normore (2015) for species of Xanthoparmelia (Vain.) Hale, growing on different rock types, did not support the photobiont guild hypothesis. However, they suggested that the range of rock substrata type in their study may have been too narrow to differentiate algal preference. On the other hand, they indicated that Peksa and Škaloud (2011) compared broadly defined types of substrata (defined as a ‘bark of tree’ and ‘rock’).

In this study, we found that the most common photobiont in P. pinnatifida was Trebouxia OTU S02. All samples of P. pinnatifida were collected from rocks; however, some authors previously reported the same Treboxia OTU S02 from terricolous, saxicolous and corticolous Parmeliaceae (i.e. genera Cetraria Ach., Melanohalea O.Blanco et al., Montanelia Divakar et al., Protoparmelia M.Choisy and Rhizoplaca Zopf and species Xanthoparmelia coloradoensis Hale and Vulpicida juniperinus (L.) J.-E.Mattsson & M.J.Lai) (Lindgren et al. 2014; Leavitt et al. 2015; Singh et al. 2017), but it may also occur in lichen genera representing other families, for example, Chaenotheca (Th.Fr.) Th.Fr., Circinaria Link and Umbilicaria Hoffm. (Beck 2002; Romeike et al. 2002; Molins et al. 2018). On the other hand, Trebouxia OTU S04, which corresponds to T. jamesii (UBT-86.156C3), was identified in a single specimen of P. pinnatifida (UGDA L-24293). It was previously reported exclusively from corticolous Melanohalea and Bryoria species (Lindgren et al. 2014; Leavitt et al. 2015) and seems to be very rare or at least rarely sampled, as it is poorly represented in GenBank. Moreover, the unnamed lineage of Trebouxia (SUn2) was detected in a single specimen of P. pinnatifida and, based on 99% identity, we found that it may also associate with, for example, Bryoria simplicior (Vain.) Brodo & D.Hawksw., Cetraria aculeata (Schreber) Fr., Evernia divaricata L. (Ach.) (Piercey-Normore 2009; Domaschke et al. 2012; Lindgren et al. 2014). Some variation in photobionts was also found in specimens of P. omphalodes which associate with Trebouxia OTU S05 and an unnamed lineage (SUn1). Leavitt et al. (2015) reported Trebouxia OTU S05, which corresponds to Trebouxia suecica (SAG2207), from terricolous and corticolous Parmeliaceae (i.e. Cetraria aculeata (Schreber) Fr., Letharia vulpina (L.) Hue and Melanohalea spp.). Photobionts, very similar to Trebouxia OTU S05 (100% identity), were additionally found in, for example, Bryoria fremontii (Tuck.) Brodo & D.Hawksw., Lasallia hispanica (Frey) Sancho & Crespo, Lecanora rupicola (L.) Zahlbr. and Tephromela atra (Huds.) Hafellner (Blaha et al. 2006; Lindgren et al. 2014; Muggia et al 2014; Paul et al. 2018). Moreover, the unnamed lineage of Treboxia (SUn1) was detected in three specimens of P. omphalodes and, based on 99% identity, we found that it may also associate with, for example, Bryoria spp., Cetraria spp., Evernia mesomorpha Nyl. Flavocetraria nivalis (L.) Kärnefelt & A.Thell and Vulpicida pinastrii (Scop.) J.-E.Mattsson & M.J.Lai (Opanowicz and Grube 2004; Piercey-Normore 2009; Lindgren et al. 2014; Onuț-Brännström et al. 2018). Therefore, the results obtained, based on our dataset, do not support the ecological guild hypothesis; however, our sampling was rather limited and we did not analyse co-occurring species. Although the type of substrata seems not to correspond to any of Trebouxia OTUs, bioclimatic factors, such as annual mean temperature, maximum temperature of warmest month or precipitation, may influence the patterns of photobionts distribution. However, to perform such an analysis, a larger set of specimens should be examined.

Interestingly, although P. omphalodes was found to associate with two lineages of Trebouxia photobionts (i.e. OTU S05 and an unidentified lineage SUn1), it does not associate with Trebouxia OTU S02, which, on the other hand, was found to associate with P. discordans (two samples). However, P. discordans also associates with Trebouxia OTU S05. As those species differ in morphology and chemistry, we suggest that those differences might be related to the photobiont type. Although some researchers did not find any correlation between different chemotypes and the associated photobionts (e.g. Blaha et al. 2006; Lindgren et al. 2014), recent studies suggested that the production of certain secondary metabolites might be triggered by the environment, for example, climate, edaphic factors or associated symbionts (e.g. Spribille et al. 2016; Lutsak et al. 2017). However, due to limited sampling, we cannot confirm this hypothesis for Parmelia spp. analysed in this study.

The created models, derived from MaxEnt, received high AUC scores, indicating high reliability of analyses (Table 5). Generated maps of distribution of suitable niches of the three lichen species were wider than the known geographical range of these lichens (Figures 7–9).

Figure 7. 

Distribution of suitable niches of P. discordans (A), P. omphalodes (B) and P. pinnatifida (C) in the Northern Hemisphere.

Figure 8. 

Distribution of suitable niches of P. omphalodes (A) and P. pinnatifida (B) in America.

Figure 9. 

Distribution of suitable niches of P. discordans (A), P. omphalodes (B) and P. pinnatifida (C) in Eurasia.

The distribution of P. discordans is limited mainly by precipitation of the driest month (bio14), but two other factors that can influence the occurrence of this taxon, varied in analyses conducted for the Northern Hemisphere and Eurasia separately. While in the former analysis, annual mean temperature (bio1) and mean diurnal range (bio2) gave important contributions to the model, the latter analysis indicated maximum temperature of the warmest month (bio5) and temperature seasonality (bio4) as significant limiting factors. Additionally, in cases of P. omphalodes and P. pinnatifida, different variables gave various contributions to the models created for different study areas. Mean diurnal range (bio2) was the crucial limiting factor for Eurasian populations of P. omphalodes, while within the American range of this species, its occurrence depends on precipitation of the driest month (bio14). For the American distribution of P. pinnatifida, the annual mean temperature (bio1) significantly influenced the model and the distribution of Eurasian populations appears limited by the maximum temperature of the warmest month (bio5) (Table 6).

The PCA diagram (Figure 10) showed that the highest bioclimatic variation is observed in P. omphalodes and that niches of P. discordans and P. pinnatifida are embedded in this highly flexible bioclimatic tolerance of P. omphalodes. The overall high similarity in bioclimatic preferences of all three studied taxa is presented in PNO profiles created for various geographic areas (Suppl. material 2: Figure S2, Suppl. material 3: Figure S3, Suppl. material 4: Figure S4). On a global scale, P. pinnatifida and P. omphalodes occupy similar niches (D = 0.581, I = 0.840), while bioclimatic preferences of P. discordans are more similar to P. omphalodes than to P. pinnatifida (Table 7). In the American range, P. omphalodes and P. pinnatifida occupy very similar habitats (D = 0.821, I = 0.968; Table 8). Within Eurasian populations, the highest similarity is observed for P. omphalodes and P. discordans (D = 0.587, I = 0.828); however, P. pinnatifida and P. omphalodes also occupy similar niches (D = 0.564, I = 0.820; Table 9).

Figure 10. 

Principal components analysis (PCA) of P. discordans (red), P. omphalodes (blue) and P. pinnatifida (green), based on the bioclimatic factors from individuals.

According to published data (Sanders and Lücking 2002; Büdel and Scheidegger 2008), lichens without vegetative propagules, dispersing both bionts independently, require the contact of the mycobiont with a compatible photobiont species in suitable environmental conditions to establish new thalli. Results of ecological niche modelling, presented here, confirmed that species from the analysed group occupy similar niches. In Figure 2, one sequence of photobionts, associating with P. discordans, belong to Trebouxia OTU S05 and the second to Trebouxia OTU S02. The latter is the most common photobiont of P. pinnatifida which, on the other hand, was also found to associate with Trebouxia OTU S04 and an unnamed Trebouxia lineage SUn2. However, none of photobionts from P. omphalodes belongs to Trebouxia OTU S02 and OTU S04, but this taxon associates with two lineages of Trebouxia photobionts (i.e. OTU S05 and an unnamed lineage SUn1). These results show that, despite the species from P. omphalodes group differing in associated photobiont species, they exhibit similar niche preferece.

PCA (Figure 10) results showed that P. omphalodes is characterised by the highest bioclimatic variation in comparison with other species from the P. omphalodes group. On the other hand, the ENM method has shown that the potential distribution of P. omphalodes is wider than its known current occurrence range (Figures 4, 6–8). The absence of this taxon in the potential niches may be caused by the lack of suitable photobiont species in those areas or that the model did not capture the relevant variation and so overestimates the niche. Two Trebouxia lineages are found in this species, i.e. OTU S05 and an unnamed lineage. Such flexibility in the photobiont choice may facilitate the mycobiont colonisation of new niches; however, some of those photobionts may be relatively rare. Trebouxia OTU S05, which corresponds to the generalist Trebouxia suecica, was previously reported from numerous terricolous and corticolous species in temperate, boreal and alpine climates, while the unnamed lineage of Trebouxia (SUn1, Table 10), present in three specimens, probably also occurs in selected terricolous and corticolous species (Table 10). Probably the latter is characterised by narrower ecological amplitude, but it needs further studies. On the other hand, P. pinnatifida forms associations with three Trebouxia lineages, i.e. OTUs S02 and S04 and an unnamed lineage (SUn2, Table 10). Most photobiont sequences from P. pinnatifida were grouped in OTU S02 clade. They were collected from different localities in Poland (Beskidy Mts, Sudety Mts, Stołowe Mts), Norway and Sweden. Moreover, the same Trebouxia OTU S02 was found in terricolous, saxicolous and corticolous lichens (e.g. Leavitt et al. 2015). It suggests that Trebouxia OTU S02 has a broad ecological amplitude and worldwide distribution. Therefore P. pinnatifida may also have wider geographical distribution than current data suggest. The absence of those species in some localities may be caused by the lack of unambiguous morphological and chemical features necessary for their identification. For this reason, herbarium material from the group P. omphalodes requires re-determination. On the other hand, the possible overestimation of the MaxEnt models may be due to additional, ecological factors (e.g. interaction with other organisms) which were not included in our analyses, but limit the distribution of the studied lichens.