Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Adaptive differentiation coincides with local bioclimatic conditions along an elevational cline in populations of a lichen-forming fungus

BMC Evolutionary Biology volume 17, Article number: 93 (2017) Cite this article

Abstract

Background

Many fungal species occur across a variety of habitats. Particularly lichens, fungi forming symbioses with photosynthetic partners, have evolved remarkable tolerances for environmental extremes. Despite their ecological importance and ubiquity, little is known about the genetic basis of adaption in lichen populations. Here we studied patterns of genome-wide differentiation in the lichen-forming fungus Lasallia pustulata along an altitudinal gradient in the Mediterranean region. We resequenced six populations as pools and identified highly differentiated genomic regions. We then detected gene-environment correlations while controlling for shared population history and pooled sequencing bias, and performed ecophysiological experiments to assess fitness differences of individuals from different environments.

Results

We detected two strongly differentiated genetic clusters linked to Mediterranean and temperate-oceanic climate, and an admixture zone, which coincided with the transition between the two bioclimates. High altitude individuals showed ecophysiological adaptations to wetter and more shaded conditions. Highly differentiated genome regions contained a number of genes associated with stress response, local environmental adaptation, and sexual reproduction.

Conclusions

Taken together our results provide evidence for a complex interplay between demographic history and spatially varying selection acting on a number of key biological processes, suggesting a scenario of ecological speciation.

Background

Fungi are diverse and ubiquitous, having evolved over time to occupy a wide range of ecological niches. Some fungal species are exceptionally proficient at surviving a broad range of environmental conditions. In nature, these species can inhabit latitudinal and altitudinal clines that span considerable temperature ranges [1]. Individuals from species with broad ecological amplitudes exhibit local adaptation when divergent selection is strong relative to the rate of gene flow [2, 3]. Locally adapted individuals show higher fitness in their focal environment relative to immigrants. Despite the wealth of studies investigating genetic structure and dispersal of fungi, the processes that shape adaptive genetic polymorphism in wild populations are not well understood.

Environmental factors have been shown to be drivers of local adaptation in diverse fungi. For example, temperature differences are responsible for the maintenance of differentially adapted populations of pathogens [4]. Temperature changes were also shown to be drivers of adaptation in natural populations of the saprotrophic fungus Neurospora crassa, resulting in genomic islands of differentiation involved in cold-response and circadian rhythm [5]. Differences in climate and soil salinity correlated with regions of extreme genomic divergence between coastal and montane populations of an ectomycorrhizal basidiomycete [6]. Furthermore, experimental evidence suggests that salinity and temperature are drivers of ecological isolation in experimentally derived lineages of baker’s yeast Saccharomyces cerevisiae [7]. In some cases, ecological divergence can promote reproductive barriers between populations of differentially adapted ecotypes [8], as is the case in the filamentous fungus Neurospora crassa [9]. Overall however, adaptive genomic polymorphisms have been investigated only in a small set of fungal species and life styles.

Lichen-forming fungi, which constitute about half of the described ascomycetes, are a nutritionally specialized group of fungi that form obligate symbiotic associations with green algae and bacteria [10]. Lichens are suitable to study local adaptation because they can tolerate extreme environmental conditions and sustain growth despite frequent cycles of desiccation and rehydration, low nutrient availability, and large fluctuations in temperature (e.g., [11, 12]). The distributional ranges of many lichens span broad climatic ranges [13–16]. Furthermore, long-lived, sessile organisms such as lichens experience strong selection pressures [17]. This may lead to reduced survival of maladapted individuals, and create steeper genetic gradients between differentially selected populations [18]. Environmental stressors, such as drought and high-light conditions, have been shown to trigger physiological adjustments in different lichens [19–21]. The genetic bases of these adaptive responses are currently poorly understood. Many species of lichenized fungi show genetic differentiation among populations despite ongoing gene flow, even across thousands of kilometers. This suggests a role of spatially varying selection in maintaining biogeographic structure (reviewed in [22]). Interestingly, one study reported gene pool associations with altitude and interpreted this as evidence for climate-driven local adaptation [23]. However, no study so far has specifically addressed adaptive diversity in geographically close (<20 km), but ecologically distant lichen populations.

Altitudinal gradients are suited to study local adaptation because ecological transitions are typically steep and occur at relatively short distances, thus limiting the confounding effect of distinct regional evolutionary histories [24]. Moreover, altitudinal gradients are also climatic gradients, characterized by decreasing temperature and atmospheric pressure, increasing relative air humidity, rainfall, and solar radiation with increasing altitude [25]. Thus, adaptation along altitudinal gradients can be explored as a proxy model for genomic responses to climate change [26–29].

Here we report on the population genomics of a lichen-forming ascomycete along an altitudinal gradient in the Mediterranean region. As model, we chose Lasallia pustulata (Umbilicariaceae), a species with a distribution from southern Europe to northern Scandinavia, which forms dense populations on exposed, siliceous rocks [30]. Using genomic data from geographically close populations along a steep altitudinal gradient in northern Sardinia (Italy), we analyzed whether genetic clusters were present, and whether relatedness between clusters was correlated with signatures of local adaptation. Heat, drought, and radiation stress constitute determining factors in the composition of biological communities inhabiting rocky outcrops and boulders in Mediterranean mountains [31]. Therefore we tested the hypothesis that environmental factors shape genome-wide population differentiation in lichenized fungi which occur across different bioclimatic regions. Specifically, we addressed the following questions: i) what is the genome-wide population structure and connectivity between geographically close populations along an elevation gradient?, ii) what are putative functions of highly differentiated genes between the genetic clusters?, iii) what are putative functions of the genes showing strong correlation with local climatic factors?, and iv) do individuals belonging to different genetic clusters (and environments) display fitness differences?

Methods

Study organism, study site, and sampling

Lasallia pustulata is a foliose, rock-inhabiting, haploid lichen-forming ascomycete. Individuals are attached to the substrate with a central holdfast. L. pustulata has a mixed strategy of asexual and sexual cycles. Asexual reproduction via isidia – macroscopic dispersal units containing both symbionts that break off the mother thallus – is the typical (and in most populations only) way of reproduction. Isidia are considered to be detached from the thallus mainly by raindrops, and dispersed over short distances [30]. We collected samples from six populations along an altitudinal gradient in the Limbara massif (Sardinia, Italy). The transect extended from Lake Coghinas (population 1: 176 m a.s.l.) to Punta Balestrieri (population 6: 1303 m a.s.l.) covering a linear distance of ~13.5 km. Intermediate populations were located at 297 m a.s.l. (population 2), 588 m a.s.l. (population 3), 842 m a.s.l. (population 4), and 1125 m a.s.l. (population 5). The maximal linear distance between populations was ~9 km (Fig. 1). Populations were located on horizontal or gently sloping, fully sun-exposed rock faces in scattered Paleozoic granitic outcrops, and covered an area of ~50 m2. For each population, we collected 100 thallus pieces of ~8 mm in diameter. Our sampling design aimed at capturing the maximal diversity present at the sites. The minimum distance between sampled individuals was 50 cm to maximize the inclusion of different genets. Samples were collected with sterile tools and transferred into sterile 2 ml tubes.

Fig. 1
figure 1

Location of the study site in north-eastern Sardinia (Italy) and location of six populations of L. pustulata along the elevation cline. The numbers in the color legend refer to the bioclimatic zones by [32]: 17 - lower mesomediterranean-upper dry-weak euoceanic, 20 - lower mesomediterranean-lower sub-humid-weak euoceanic, 26 - upper mesomediterranean-lower sub-humid-weak euoceanic, 28 - upper mesomediterranean-upper sub-humid-weak euoceanic, 30 - upper mesomediterranean-lower humid-weak euoceanic, 36 - upper mesotemperate (sub-Mediterranean)-lower humid-weak euoceanic, 37 - upper mesotemperate (sub-Mediterranean)-lower humid-weak semi-continental, 40 - lower supratemperate-lower humid-weak semi-continental, 41 - lower supratemperate-lower hyperhumid-weak semi-continental. The color gradient represents the bioclimatic profile of the cline, ranging from Mediterranean pluviseasonal oceanic (M) (populations 1–4) to temperate oceanic, submediterranean variant (TOm) (population 5), to temperate oceanic (TO) (population 6) climate

Full size image

The study area encompasses three distinct bioclimates: i) Mediterranean pluviseasonal oceanic (populations 1–4), ii) temperate oceanic submediterranean (population 5), and iii) temperate oceanic (population 6) (Fig. 1) [32]. Temperature data collected between May 28th 2014 and June 2nd 2015 from loggers positioned at the level of L. pustulata thalli (2 loggers per population) indicate that localities 1 to 4 have higher summer temperatures and are less prone to freezing during winter than populations 5 and 6. Logger data also showed that temperatures of the rock surfaces to which thalli are attached frequently exceeded 50 °C in summer. The populations therefore experienced seasonal temperature fluctuations on the order of >40 °C (see Additional file 1).

DNA extraction and genome resequencing

For each population, we extracted genomic DNA separately from each individual using a CTAB-based method [33]. DNA concentration was measured with a Qubit fluorometer (dsDNA BR, Invitrogen). A pooled sample was created for each population containing equal amounts of DNA from each sample (Pool-seq). Library preparation (200–300 bp insert size), sequencing on an Illumina HiSeq2000 with 100 bp paired-end chemistry at ~90x coverage per population, as well as tags and adaptor removal were performed by GenXPro GmbH (Frankfurt am Main, Germany).

Genome annotation

As reference genome, we used the draft assembly of L. pustulata available at the European Nucleotide Archive (http://www.ebi.ac.uk/ena/data/view/GCA_000938525.1). The draft genome is composed of 3891 scaffolds (average length of 10Kbp) for a total length of 39.2 Mbp and an N50 scaffold of size 104.3Kbp. The genome is ~92% complete according to an assessment with the software BUSCO 2.0 [34] and a lineage-specific set of Ascomycota single-copy orthologs.

For gene model prediction we used both ab-initio based methods and RNA-Seq derived transcript mapping onto the assembled genome following the method described in [35]. For this purpose, total RNA was isolated from a thallus of L. pustulata collected near Orscholz (Saarland, Germany; N49.5012, E6.5440) in July 2013 using the method by [36], and purified using the RNeasy MinElute Clean-up Kit (Qiagen). Paired-end sequencing was performed using Illumina MiSeq (2x250 bp) by StarSEQ (Mainz, Germany). RNA-Seq data was quality-filtered using Trimmomatic [37], with a length cutoff of 200 and a quality cutoff of 20 in a window of 5 bp.

We used Blast2GO [38] to annotate the predicted protein sequences with gene ontology (GO) terms and protein names using NCBI's nr database at an E-value cut-off of 1x10−3, and default weighting parameters. We also annotated each protein with InterPro domains using InterProScan [39].

SNP analysis

We filtered out reads shorter than 80 bp, reads with N's, and reads with average base quality scores less than 26 along with their pairs using FastQFS [40]. Trimmed paired-end reads of each pool were mapped to the L. pustulata genome using BWA-MEM [41] and default parameters. Unambiguously aligned reads with a minimum mapping quality of 20 were extracted with SAMtools v1.18 [42]. Reads were sorted and duplicates were marked with Picard using the tools SortSam.jar and MarkDuplicates.jar. Single nucleotide polymorphisms were called with SAMtools (mpileup, [43]). Indels were detected and masked with PoPoolation [44] using the scripts identify-genomic-indel-regions.pl (−−min-count 2 −−indel-window 5) and filter-pileup-by-gtf.pl. The synchronized file was converted into a gene-based synchronized file using the script create-genewise-sync.pl in PoPoolation2 [45]. The coverage for each population was reduced to a uniform coverage of 30 with PoPoolation2 using the sync-file and the script subsample-synchronized.pl (−−without-replacement).

Population genetics analyses

To characterize genome-wide patterns of variation, we estimated three population genetic parameters by accounting for pooling: i) π, a measure of polymorphism in a sample of sequences scaled to their length, ii) Watterson's θ (θW), a measure of the number of segregating sites, and iii) Tajima's D, a measure of the skew of allele frequency distribution. All estimates were calculated in non-overlapping 10-kb windows across the genome using PoPoolation [44], assuming a minimum count of two. Differences in genetic diversity among populations were tested using linear mixed effect models in R 3.2.2 [46, 47]. For each diversity measure, models included population as a fixed effect predictor and incorporated scaffolds as a random effect across populations. Pairwise population comparisons were then obtained from post hoc Tukey contrasts of the respective model predictors [48].

To identify strongly differentiated alleles, we adopted an empirical outlier approach. Genetic differentiation (FST) was calculated with fst-sliding.pl in PoPoolation2. We only considered SNPs with a minimum count of 4, a minimum quality of 20, and falling into the upper 0.5% tail of the FST distribution, corresponding to an FST threshold of 1.0. Highly differentiated SNPs were further inspected with Fisher’s exact test [49] to identify significant allele frequencies differences between population pairs using the script fisher-test.pl in PoPoolation2 and a Bonferroni-corrected p-value of 0.003. In addition, we estimated average FST across all polymorphic SNPs for each gene and only considered those falling into the upper 5% tail of distribution to be truly differentiated. Based on this analysis, we calculated the percentage of overlap between SNP-based and gene-based lists. SNPs were classified as genic and non-genic loci. Genic SNPs were further classified as exonic, coding, and intronic. Non-genic SNPs located in the 600 bp 5'-flanking sequence of each gene were considered putative promoter SNPs.

To visualize groups of populations with varying degrees of similarity to one another, we first obtained a reduced set of pairwise FST distance matrices based on sample quantiles (0.975, 0.75, 0.5, 0.25, 0.025) of the full set of distances across all polymorphic SNPs. The resulting set of pairwise (quantile) distances was then jointly analyzed using a three-way generalization of classical multidimensional scaling (DISTATIS, [50]). DISTATIS calculates a compromise distance space from the weighted average of all cross-product matrices derived from the set of quantile distance matrices. This compromise can then be used to visualize positional relations among populations 1 – 6 based on their genetic distances. Moreover, we obtained 95% confidence intervals around each population's compromise position using bootstrap resampling [51].

To reconstruct the historical relationships among populations using their current genome-wide allele frequencies, we used TreeMix v1.12 [52]. We created a maximum likelihood phylogeny of the populations based on all polymorphic SNPs, using blocks of 500 SNPs to account for linkage disequilibrium. To test for the presence of admixture and migration among populations, we calculated f3 and f4 statistics. These statistics are formal tests for admixture as they detect correlations in allele frequencies that are not compatible with population evolution following a bifurcating tree. To calculate f3 and f4 statistics we used the threepop and fourpop functions in TreeMix on all possible triplets and tetraplets population groups. To further explore population relationships, we calculated neighbor joining trees for the correlation matrix obtained with Bayenv2.0 (see below) and for the matrix of pairwise FST values calculated across all polymorphic sites using the package ape in R [46].

To corroborate the hypothesis of admixture, we estimated gene flow among the three major groups (A: populations 1 to 4, B: population 5, C: population 6) across multiple intergenic polymorphic loci using the coalescent-based method MIGRATE-N 3.2.6 [53]. We estimated relative effective population size (Ne) according to the relation θ = Neμ, assuming identical but unknown mutation rates (μ) in all populations. Haplotypes (20 per population) were obtained by parsing the output files of PoPoolation2 for genomic regions shorter than 90 bp, and containing three or more SNPs. To minimize the chance that the loci are linked, we retained only one such region per scaffold. This information was then used to extract the individual Illumina reads covering this region from the mapping file for each population. The reads were aligned and trimmed to the informative region. The alignments were filtered for a minimum coverage of 15x, maximum coverage of 100x and minimum length of 20 bp in all populations. This resulted in a data set of 5880 sequences from 49 loci, covering in total 1763 bp. All previously described steps were performed with custom Python scripts (see [54]). Bayesian estimates of number of migrants (Nm) and θ were obtained under an unconstrained migration model with variable θ using MIGRATE-N 3.2.6 [53] for each pair of genetic clusters separately. We used a uniform prior on both θ (0.0-0.40) and Nm (0.0-600). A Metropolis-coupled Monte-Carlo chain with static heating (1.0, 1.5, 3, 1 × 106) was run for 1.8 × 106 generations, recording every 600th step after a burn-in period of 6 × 104 generations. Convergence was monitored with Tracer (http://beast.bio.ed.ac.uk/). All effective sample sizes of the MCMC chain were larger than 104.

Environmental association analysis

To identify candidate loci for altitude specific adaptation, we correlated allele frequencies of populations with the environment using Bayenv2.0 [55]. Bayenv2.0 models the sampling error of pooled sequencing, and accounts for the confounding effect of neutral, demographic signals. To summarize the climate along the cline, we used elevation and 19 bioclimatic variables from the WorldClim database [56]. Correlation among these was checked using the function rcorr in R [46]. Elevation was strongly correlated with all bioclimatic variables. Variation in the bioclimatic variables and elevation was thus summarized using PCA, resulting in one composite climate variable explaining 97.6% of the variance, hereafter referred to as env1 (see Additional file 2). To build the reference covariance matrix of population allele frequencies, we used a subsample of 10,000 polymorphic SNPs based on 1,000,000 MCMC iterations. To ensure convergence we estimated a second matrix from another subsample of 10,000 SNPs. We explored gene-environment correlation by estimating the statistic Z between allele frequencies of all SNPs and env1 per population for 200,000 iterations. SNPs with the highest score possible for Z (i.e., Z = 0.5) were considered as showing strong support for a non-zero correlation.

Gene ontology enrichment

Gene Ontology (GO) term enrichment is a technique for interpreting the functions of a set of genes making use of the GO system of classification (http://www.geneontology.org/). In this system genes are assigned to predefined bins depending on their functional characteristics in a species-independent manner. An enrichment analysis will find which GO terms are over-represented in a given data set using the annotations for that gene set. We used the R package topGO [57] to search for an enrichment of different GO categories. The analysis was performed for i) the set of genes containing SNPs falling into the upper 0.5% tail of the SNP-derived FST distribution, ii) the subset of these that are differentially fixed between the low altitude (population 1 to 4) and high altitude (population 6) genetic clusters, and iii) the set of genes containing the Baynev2.0 top 1% SNPs. All genes with a GO annotation were used as background. Significance for each GO-identifier was computed with Fisher’s exact test at α = 5%. We used the ‘elim’ method in topGO to iteratively remove genes mapped to significant GO terms from more general terms, thus reducing the rate of false positives. Only GOs with more than three associated genes were considered. We used the REVIGO tool [58] to produce summaries of non-redundant GO terms grouped into functional categories.

Ecophysiology

We performed ecophysiological experiments to assess differences in physiological traits in the populations. Measurements were performed on three samples per population. For this analysis, we randomly selected specimens with a minimum diameter of 6 cm to have enough material to perform replicate measurements. To explore the genetic relatedness of the individuals we genotyped each specimen at six loci, covering a total of approximately 4.1 Kbp. We selected three of the loci from genes containing top 0.5% differentiated SNPs, and three from those with top Bayenv2.0 SNPs. For primers and genetic characteristics of the loci see Additional file 3.

We investigated the thalli for differences in biomass and chlorophyll content per surface area. To calculate the specific thallus area (mm2/mg), we first determined the thallus size by photographing wetted thalli on scale paper using a binocular microscope and the AxioVision software (Carl Zeiss, Jena, Germany). We determined the dry weights (DW) of these thalli by weighing after 3 days of oven drying at 60 °C. We also measured thallus chlorophyll content according to [59]. Statistical significance of differences in biomass and chlorophyll content between groups was determined using a Mann–Whitney test.

To characterize the physiological response to different light conditions, and different thallus water contents, we conducted CO2 gas exchange measurements using a portable mini cuvette system (GFS 3000, Walz Company, Effeltrich, Germany). We measured the response of net photosynthesis (NP) and dark respiration (DR) to thallus water content (WC) for a subsample of six thalli, representing the two main genetic groups present along the gradient. We measured complete desiccation cycles (from water saturated to air dry thalli) at saturating light (750 μmol photons m−2s−1), ambient CO2, at 17 °C (within the optimal temperature range for CO2-gas exchange of this species). We weighed the samples between each measurement and later extrapolated WC as a percentage of DW. We determined DW after 5 days in a desiccator over silica gel. We considered ninety percent of maximum NP to be a reasonable estimate for optimal water saturation. We measured the samples at 3–15 h intervals for NP and DR. Immediately after each measurement we removed the samples from the cuvettes and determined their weight to calculate the decrease in WC. This process was continued for 97 h with each sample until CO2 exchange ceased due to complete drying. Statistical significance of differences in TWC and maximum NP between groups was determined using a t-test.

Results

Reference gene set

We identified a total of 8268 genes in the L. pustulata genome. In total 5747 genes were assigned a GO term.

Genome-wide variation

After adapter and quality trimming, we obtained 179,809,145 paired-end reads (26.6–32.1 million per pool, total of 16.2 GB, average coverage per pool: 89.91x).

To examine sequence variation, we used two estimates of nucleotide diversity, π and θW. When averaged for all 10-kb windows across the genome, estimates of π were highest in population 5 (π = 0.006 ± 0.004), and lowest in population 4 (π = 0.004 ± 0.003). Estimates of θW were instead highest in population 2 (θW = 0.005 ± 0.003), and lowest in population 4 (θW = 0.003 ± 0.003). To examine deviation from neutrality, we calculated Tajima's D in 10-kb windows across the genome. Average D deviated from neutrality and differed significantly among populations. D was negative in populations 2 (D = −0.839 ± 0.713) and 6 (D = −0.260 ± 0.322), and positive in all other populations, being highest in population 5 (D = 1.112 ± 0.54) (Table 1, Additional file 4, Additional file 5).

Table 1 Sampling locations of Lasallia pustulata populations and the mean and standard deviation for three standard population genetic parameters, π, Watterson's θ (θW), and Tajima's D in nonoverlapping 10-kb windows across the genome of L. pustulata
Full size table

Patterns of genetic differentiation

Among the six populations we identified 722,401 polymorphic SNPs (Table 2). Mean pairwise FST based on all polymorphic sites was moderate with an average of 0.124, and ranging from 0.044 (Pool1 vs. Pool4) to 0.236 (Pool2 vs. Pool6) (Fig. 2).

Table 2 Number of variants for all annotated features of the L. pustulata genome
Full size table
Fig. 2
figure 2

Pairwise FST comparisons of populations 1 to 6 across a total of 722,401 SNPs. Boxes extend from the first to the third quartiles, with a horizontal line indicating the median. The horizontal line across the graph indicates the top 0.5 quantile (FST = 1)

Full size image

A total of 30,571 SNPs located in 2944 genes fell into the upper 0.5% tail of the distribution. Of these, 4170 SNPs in 595 genes were differentially fixed between the low (1 to 4) and high-altitude (6) population clusters. When calculating average FST across all polymorphic sites within a given gene, 2413 genes fell into the top 5% tail. Our SNP-based approach detected 72.9% of the highly differentiated genes.

We found strong genetic structure separating lower altitude populations (populations 1 to 4) from the other populations. The multidimensional scaling of the FST quantile distance SNP matrix illustrated the close genetic affinities of populations 1 to 4, with population 5 occupying an intermediate position between these and population 6 (Fig. 3a). The tree-based analyses showed a similar, well-resolved structure (Fig. 3b), with the tree derived from the Bayenv2.0 correlation matrix showing longer internal branches (see Additional file 6). Bayenv2.0- and FST-matrices were highly correlated (Mantel-test, r = 0.999, P = 0.001; see Additional file 7).

Fig. 3
figure 3

a Compromise configuration of categories between populations based on the FST quantile distance matrix for 722,401 polymorphic SNPs, with 95% tolerance ellipses. b Tree inferred with TreeMix for the six population of L. pustulata based on 722,401 polymorphic SNPs. Numbers at the branching points are support values from bootstrapping based on 1000 runs

Full size image

Using the threepop test we found clear evidence of admixture at the level of population 5. All population triplets having population 5 as the admixed group displayed significantly negative f3 values (see Additional file 8). This is in accord with the higher nucleotide diversity and positive Tajima's D values for population 5 (see Additional file 4). We also found support for migration among populations as 35 out of 45 four-population tests rejected all possible tree topologies without migration (|z| > 3, i.e., p < 0.001). We inferred higher significance for pairs grouping together population 5 and 6 with one of the lower altitude populations, respectively (see Additional file 8), which is also in accordance with the proposed admixture scenario.

To further describe the migration pattern, we performed estimations of migration rates among the major genetic groups with MIGRATE-N 3.2.6. Mutation-scaled effective population sizes varied between groups, ranging from θ = ~0.054 in group A (pop. 1–4) to θ = ~0.065 in group B (pop. 5) (see Additional file 9). Migration rates varied by several orders of magnitude. Results supported the hypothesis that the genetic diversity of population 5 is the result of admixture from the other two genetic groups, while gene flow rates between the other groups are negligible in comparison. This is in line with the FST- and tree-based analyses.

Candidate loci for local adaptation

To identify candidate loci for altitude specific adaptation, we correlated allele frequencies with an environmental variable summarizing altitude and climate using Bayenv2.0. A total of 2978 SNPs showed the highest score possible for Z and were located in 616 genes.

At the SNP level, the overlap between the Bayenv2.0- and FST-based approaches was low (3.02%, 90 SNPs). Of these, 42 SNPs were located in 39 genes (see Additional file 10). Genes containing top Bayenv2.0 SNPs matched 216 of the FST-based top 5% differentiated genes.

Functional inference of candidate genes

Gene set enrichment analysis of the 2944 genes containing top 0.5% SNPs indicated the presence of 62 enriched biological processes (see Additional file 11a). These involve many pathways, some of which are centrally important for stress response, cell growth, carbohydrate transport, and both asexual and sexual reproduction. For example among the significantly enriched categories we found biological processed like response to abiotic stimulus, growth, gene expression, RNA processing, translation, fungal-type cell wall polysaccharide biosynthetic process, catabolic processes, protein N-linked glycosylation, trehalose biosynthetic process, and developmental process involved in reproduction. GO enrichment of the 595 genes containing SNPs differentially fixed between the low (1 to 4) and high-altitude (6) genetic groups resulted in 38 enriched biological processes, including sexual and asexual reproduction, trehalose biosynthesis, growth, response to oxidative stresses, cell wall and ribosome biogenesis, and gene expression (see Additional file 11b).

GO enrichment of the top 1% Bayenv2.0 SNPs indicated that 23 biological processes were enriched (Table 3, see Additional file 12). Among the biological processes likely involved in adaptation to altitude we found localization, signal transduction, DNA repair, lipid modification, histone methylation, catabolic processes, and cell-wall biogenesis.

Table 3 GO enriched categories for top 1% Z Bayenv2.0 environmentally associated SNPs
Full size table

Ecophysiology

Multi-locus genotyping grouped the samples into two genetic groups, one composed of all thalli from population 6 (G1), and one composed of thalli from all remaining sites (G2). Differentiated SNPs in all six markers were monomorphic between the two groups (see Additional file 3, Additional file 13). The genetic separation coincided with differences in anatomy and physiological responses to thallus water content (WC; Fig. 4) and high light conditions (see Additional file 14). First, G1 thalli had higher biomass (p = 0.002), and higher chlorophyll a + b content per surface area unit (p = 0.002) (Table 4). Second, G1 thalli needed higher WC for reaching maximal NP (0.5498 ± 0.104 mm H2O) compared to G2 thalli (0.175 ± 0.0 mm H2O; Fig. 4). In addition, G1 thalli reached their maximum NP rates (>90% of max value) at lower light intensity (see Additional file 14) (p = 0.003). In relation to surface area unit, G1 thalli fixed almost three times as much CO2. In relation to thallus dry weight G1 and G2 specimens did not significantly differ in their CO2 fixation rates (p = 0.65) (see Additional file 14, Additional file 15).

Fig. 4
figure 4

Photosynthetic CO2 gas exchange of L. pustulata highland (population 6; blue) and lowland population (populations 1 to 5; red) related to thallus water content (TWC). TWC is expressed as mm “precipitation”. Polynomic regressions lines are indicated with their r2 value. Circles = net CO2 uptake, triangles = dark respiration

Full size image
Table 4 Dry weight and chlorophyll a + b content of the different populations (N = 3 for each of the populations)
Full size table

Discussion

Temperature and precipitation drive large-scale distribution patterns of lichens. It is thus expected that much of the signal of adaptation among lichen populations should occur along these gradients [60]. Here we presented the first genome-based analysis of population differentiation associated with an environmental gradient for a lichen-forming fungus. The studied populations underwent several generations of asexual reproduction, as only rarely sexual structures (apothecia) were observed in any of the sites. Thus, by applying Pool-seq resequencing to large population samples, we were able to track the frequencies of diverged long-lived lineages in each of the sites, and describe genome-wide population divergence in relation to changes in altitude.

Our study revealed significant differentiation and structure of L. pustulata populations. We found two genetic clusters along the gradient. One cluster is predominant at low elevations (up to ~800 m a.s.l.), while the other is predominant at high elevations (~1300 m a.s.l.). Our data suggest extensive admixture of these clusters at ~1100 m a.s.l. Given the high number of reciprocally fixed SNPs between the clusters and the high levels of clonal propagation in this fungal species, possible explanations for the observed pattern include ancient divergence, and a combination of limited gene flow, long generation times, and strong environmental filtering. Unfortunately there are no estimates for sexual or asexual generation times in L. pustulata from which we could calculate the age of the split between the clusters. Thus it is currently impossible to formally distinguish between the above scenarios. Ancient population splits and high genomic divergence have been reported in non-lichenized fungi. For example, strong genomic divergence, evidence for ancient population splits and introgression between subpopulations were inferred for the human pathogens Coccidioides immitis and C. posadasii [61]. Strong genetic structure was also found in locally adapted subpopulations of Neurospora crassa [5]. Numerous theoretical and empirical studies suggest that strong population divergence among continuously distributed populations may be caused by selection along environmental gradients promoting adaptation to different environmental conditions and ultimately impeding gene flow [62–64]. Interestingly, the genetic clusters of L. pustulata correspond to the major bioclimatic zones covered by our transect, with their admixture zone coinciding with the transition between the Mediterranean and the temperate-oceanic climate. Therefore, environmental filters likely contribute to the observed genetic structure.

To search for loci putatively involved in environmental adaptation, we detected allele-climate associations. Heat, drought, and intense light belong to the selective forces that cause differentiation among and within plant species in Mediterranean ecosystems [65]. We thus expected to find diversifying selection in those loci associated with pathways of the fungal environmental stress response (ESR). The ESR is in fact a common feature in the response of fungi to different environments, and it is responsible for initiating gene expression that protects the cell against stress [66, 67]. In yeasts, the ESR includes ~900 genes and requires a coordinated effort from multiple pathways, including signal transduction molecules, enzymes involved in cellwall biogenesis and maintenance, genes responsible for regulation of transcription, post-translational modification, and enzymes with proteolytic or antioxidant activities [68]. We found representative genes of each of the above pathways in our set of candidates for altitudinal adaptation. One example is alpha-ketoglutarate-dependent dioxygenase, a gene that is involved in the catalysis of taurine. Taurine is a solute required in osmoregulation, and has been linked to the survival of the fungus Ochroconis mirabilis in different habitats [69]. In the same functional category, we found candidate SNPs in genes such as a flavin-binding monooxygenase, and in the putative essential subunit of U3-containing 90S preribosome (NOP9), which has been reported among the 71 essential genes required for oxidative stress tolerance in Saccharomyces cerevisiae [70]. Additionally we found the calcium channel subunit Cch1, which has been reported to be involved in the Ca2+ release in response to exogenous oxidative stress in yeast [71], and a putative flap endonuclease, which is part of the base-excision repair pathway that removes lesions resulting from exposure to reactive oxygen species in yeast [72, 73]. Among the candidates involved in thermal stress response, we found the small heat shock protein Hsp20, glutamate carboxypeptidase, and a putative thermotolerance protein. Other candidates for diversifying selection include genes involved in the regulation of gene transcription, in particular Isw1, a gene that functions in parallel with the NuA4 and Swr1 complexes to regulate stress-induced gene transcription via chromatin remodeling in yeast [74]; interestingly, also Swr1 was detected among the candidates. Furthermore we found genes putatively involved in cellwall integrity and filamentous growth pathways (septin [75]), and in UV-damage response (UV-damage endonuclease [76, 77]).

Several of the top differentiated genes containing environmentally associated SNPs were also involved in the ESR [78]. In particular we found genes linked to signal transduction and cell-wall integrity pathways, such as the transmembrane cellwall sensor Wsc4, a MAPKKK-cascade protein kinase, a two-component osmosensing histidine kinase, and a calcium/calmodulin-dependent kinase (CAMK). CAMK proteins were reported to be involved in thermotolerance and oxidative stress survival in Neurospora crassa [79]. The presence of a putative trehalose-phosphate synthase, and of a number of ribosomal proteins genes, RNA helicases, and thioredoxin among the candidates suggests that many genes involved in redox homeostasis and also important for cold-shock response are putative targets of selection [80]. It is therefore tempting to speculate that adaptation to different temperatures, in particular cold-shock, is a driver of population differentiation in L. pustulata, given a difference of ~6 °C in the mean annual temperature between top and bottom of the cline, and frequent winter frost above 1100 m a.s.l.

The overlap between environmentally associated SNPs and FST-based outliers was low in our study. As in other studies of local adaptation [81–83], different approaches to identify candidate loci yielded different sets of candidates. This is probably due to the effects of population structure, and the parameters used for the environmental correlation. The limited overlap between sets of outliers indicates that selection along the gradient occurs mainly at the scale of local populations, and only partially at the evolutionary scale of the ancestral genetic groups [81]. Another reason for this difference may be attributed to the presence of environmental drivers that do not covary with chosen environmental factors, either biotic (e.g., interactions with photosynthetic partners, bacteria, or pathogens, and intra- and inter-specific competition), or abiotic (e.g., cloud cover, wind speed). In addition, covariance of population structure with the environment has been shown to make the method correcting for neutral population structure over-conservative [81]. Thus, candidates identified via methods that do not adjust for population structure should not be ignored, just treated carefully as their interpretation is necessarily post hoc. Future studies of L. pustulata will have to include more populations and replicate independent clines to fully disentangle demography from local selection.

Many candidate genes have known roles in stress response and growth regulation, and it is thus tempting to hypothesize that variation at these loci might affect fitness-related traits. Our ecophysiological experiments showed the presence of genetic lineages with differential fitness under different environmental conditions along the cline. The high altitude group seems to be better adapted to more humid conditions and to lower light intensities than the lowland group. In particular, samples from this population had thicker thalli and thus more biomass per area unit. In poikilohydric organisms such as lichens, an increased fungal biomass may be beneficial at wetter high-elevation sites and where winds speed up drying, especially on the exposed rocky outcrops where the species lives. This is because a higher fungal biomass may lead to i) higher mechanical stability against mechanical damage, ii) prolongation of the wet phase for an increase of the active period [84], and iii) acceptance of higher thallus water content [85]. The latter point is supported by our finding that high-altitude individuals need more water for maximal net photosynthesis than low-altitude ones. Moreover, under similar light and humidity conditions, low-altitude individuals would eventually die because of respiration rates exceeding carbon fixation rates. Structural changes towards improved thallus hydration in relation to improved photosynthetic exploitation were shown to drive acclimation in populations from different slopes in Ramalina capitata [86], between vagrant and attached morphs of Cetraria aculeata [21], between shaded and exposed populations in the Antarctic endemic Catillaria corymbosa [87], and between populations from different biomes in Psora decipiens [88]. Thallus thickness is a property controlled by the mycobiont and thus is related to the mycobiont’s response capability to environmental conditions [19]. Future genome-wide studies on the photobiont are required in order to elucidate the relative role of each symbiont in shaping the lichen’s response to the environment. Studies on the ecological benefits of phenotypic plasticity in lichens are still in their infancy, and the genomic basis of physiologically-relevant traits is far from being understood. Our results suggest that a surprisingly high, possibly adaptive, genetic diversity is responsible for anatomical and physiological differences between ecomorphs of a morphologically homogeneous lichen species. Further genomic research and physiological experiments based on replicate populations from different geographic areas, and comparisons with other lichen species will enable us to test the aforementioned hypothesis.

Adaptation to divergent environments promotes environmental specialization and reproductive isolation among fungal populations [5, 9]. In fact, the exclusion of immigrating individuals due to higher fitness of local genotypes (i.e., isolation by adaptation) can also lead to reproductive isolation resulting in little or no effective gene flow between geographically close populations [89]. This process, known as ecological speciation, interprets reproductive isolation as a by-product of adaptation to divergent environments. Our analysis of fixed variation between the highly differentiated genetic clusters of L. pustulata showed that several genes involved in sexual reproduction were significantly enriched, supporting the hypothesis of reproductive isolation. This is interesting, because we did not observe morphological evidence for sexual reproduction, such as presence of apothecia. Additional indication for reproductive isolation stems from the observation that genomic divergence between clusters is not limited to a few genomic areas, but rather widely dispersed across the genome. It involves many pathways centrally important for the response to environmental signals and stress, gene expression, growth, and metabolism. Overall, the high level of genomic divergence, and the presence of physiological differences between genetic clusters suggest the existence of L. pustulata ecotypes adapted to the Mediterranean and temperate-oceanic bioclimatic zone.

Conclusions

Pool-seq genome resequencing is a cost-effective and powerful approach to assess allelic diversity in populations, and to identify genes that are potentially under selection [90–92]. However, this method also has limitations, mainly associated with the increased impact of sequencing errors and resampling of alleles. To circumvent these issues, we created equimolar genomic libraries based on high sampling density, sequenced each population at high coverage, and used strict thresholds for sequence quality filtering. Furthermore, we based the calculations of genetic diversity measures on subsampled data to avoid coverage bias, and used analysis tools specifically adapted to Pool-seq data. An intrinsic limitation of pooled sequencing is the loss of linkage disequilibrium information. The proposed candidates may not be under selection but may have been detected because of being linked to the actual targets of selection, or to other loci that diverged because of different evolutionary processes (e.g., genetic drift). This limitation can only be overcome by using different sequencing strategies, such as genome resequencing of individuals.

Our exploratory work into the genomics of adaptation of a lichen-forming fungus reveals numerous loci and pathways putatively involved in environmental adaptation, including many loci shown in other fungi to be linked to temperature and UV-radiation stress response. Such genes provide excellent targets for further investigations. Future studies based on individual genotyping, possibly including replicate populations from different regions, additional physiological analyses including more samples, and quantitative trait locus mapping experiments of the candidate genes in controlled and field settings will help to elucidate the drivers of local adaptation in this and other fungal species.

References

  1. Gladieux P, Ropars J, Badouin H, Branca A, Aguileta G, de Vienne DM, et al. Fungal evolutionary genomics provides insight into the mechanisms of adaptive divergence in eukaryotes. Mol Ecol. 2014;23:753–73.

    Article  PubMed  Google Scholar 

  2. Savolainen O. The genomic basis of local climatic adaptation. Science. 2011;334:49–50.

    Article  CAS  PubMed  Google Scholar 

  3. Pauls SU, Nowak C, Bálint M, Pfenninger M. The impact of global climate change on genetic diversity within populations and species. Mol Ecol. 2013;22:925–46.

    Article  PubMed  Google Scholar 

  4. Mboup M, Bahri B, Leconte M, De Vallavieille-Pope C, Kaltz O, Enjalbert J. Genetic structure and local adaptation of European wheat yellow rust populations: The role of temperature-specific adaptation. Evol Appl. 2012;5:341–52.

    Article  PubMed  Google Scholar 

  5. Ellison C, Hall C, Kowbel D. Population genomics and local adaptation in wild isolates of a model microbial eukaryote. Proc Natl Acad Sci U S A. 2011;108:2831–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Branco S, Gladieux P, Ellison CE, Kuo A, LaButti K, Lipzen A, et al. Genetic isolation between two recently diverged populations of a symbiotic fungus. Mol Ecol. 2015;24:2747–58.

    Article  CAS  PubMed  Google Scholar 

  7. Dettman JR, Sirjusingh C, Kohn LM, Anderson JB. Incipient speciation by divergent adaptation and antagonistic epistasis in yeast. Nature. 2007;447:585–8.

    Article  CAS  PubMed  Google Scholar 

  8. Nosil P. Divergent host plant adaptation and reproductive isolation between ecotypes of Timema cristinae walking sticks. Amer Nat. 2007;169:151–62.

    Article  Google Scholar 

  9. Dettman JR, Anderson JB, Kohn LM. Divergent adaptation promotes reproductive isolation among experimental populations of the filamentous fungus Neurospora. BMC Evol Biol. 2008;8:35.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Honegger R. The Lichen Symbiosis—What is so Spectacular about it? Lichenol. 1998;30:193–212.

    Article  Google Scholar 

  11. Honegger R. Functional aspects of the lichen symbiosis. Annu Rev Plant Physiol Plant Mol Biol. 1991;42:553–78.

    Article  CAS  Google Scholar 

  12. Seymour FA, Crittenden PD, Dyer PS. Sex in the extremes: lichen-forming fungi. Mycologist. 2005;19:51–8.

    Article  Google Scholar 

  13. Printzen C, Ekman S, Tønsberg T. Phylogeography of Cavernularia hultenii: evidence of slow genetic drift in a widely disjunct lichen. Mol Ecol. 2003;12:1473–86.

    Article  CAS  PubMed  Google Scholar 

  14. Printzen C, Domaschke S, Fernández-Mendoza F, Pérez-Ortega S. Biogeography and ecology of Cetraria aculeata, a widely distributed lichen with a bipolar distribution. MycoKeys. 2013;6:33–53.

    Article  Google Scholar 

  15. Widmer I, Dal Grande F, Excoffier L, Holderegger R, Keller C, Mikryukov VS, et al. European phylogeography of the epiphytic lichen fungus Lobaria pulmonaria and its green algal symbiont. Mol Ecol. 2012;21:5827–44.

    Article  PubMed  Google Scholar 

  16. Sork VL, Werth S. Phylogeography of Ramalina menziesii, a widely distributed lichen-forming fungus in western North America. Mol Ecol. 2014;23:2326–39.

    Article  PubMed  Google Scholar 

  17. Sork VL, Gugger PF, Chen JM, Werth S. Evolutionary lessons from California plant phylogeography. Proc Natl Acad Sci U S A. 2016;113:8064–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Jaramillo-Correa J-P, Rodriguez-Quilon I, Grivet D, Lepoittevin C, Sebastiani F, Heuertz M, et al. Molecular proxies for climate maladaptation in a long-lived tree (Pinus pinaster Aiton, Pinaceae). Genetics. 2014;199:793–807.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Büdel B, Scheidegger C. Thallus morphology and anatomy. In: Nash T, editor. Lichen Biology. Cambridge: Cambridge University Press; 1996. p. 40–69.

    Google Scholar 

  20. Kappen L, Schroeter B, Green TGA, Seppelt RD. Chlorophylla fluorescence and CO2 exchange of Umbilicaria aprina under extreme light stress in the cold. Oecologia. 1998;113:325–31.

    Article  CAS  PubMed  Google Scholar 

  21. Pérez-Ortega S, Fernández-Mendoza F, Raggio J, Vivas M, Ascaso C, Sancho LG, et al. Extreme phenotypic variation in Cetraria aculeata (lichenized Ascomycota) adaptation or incidental modification? Ann Bot. 2012;109:1133–48.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Werth S. Population genetics of lichen-forming fungi – a review. Lichenol. 2010;42:499–519.

    Article  Google Scholar 

  23. Scheidegger C, Bilovitz PO, Werth S, Widmer I, Mayrhofer H. Hitchhiking with forests: Population genetics of the epiphytic lichen Lobaria pulmonaria in primeval and managed forests in southeastern Europe. Ecol Evol. 2012;2:2223–40.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Keller I, Alexander JM, Holderegger R, Edwards PJ. Widespread phenotypic and genetic divergence along altitudinal gradients in animals. J Evol Biol. 2013;26:2527–43.

    Article  CAS  PubMed  Google Scholar 

  25. Körner C. The use of ‘altitude’ in ecological research. Trends Ecol Evol. 2007;22:569–74.

    Article  PubMed  Google Scholar 

  26. Berner D, Körner C, Blanckenhorn WU. Grasshopper populations across 2000 m of altitude: Is there life history adaptation? Ecography. 2004;27:733–40.

    Article  Google Scholar 

  27. Bonin A, Taberlet P, Miaud C, Pompanon F. Explorative genome scan to detect candidate loci for adaptation along a gradient of altitude in the common frog (Rana temporaria). Mol Biol Evol. 2006;23:773–83.

    Article  CAS  PubMed  Google Scholar 

  28. Stapley J, Reger J, Feulner PGD, Smadja C, Galindo J, Ekblom R, et al. Adaptation genomics: The next generation. Trends Ecol Evol. 2010;25:705–12.

    Article  PubMed  Google Scholar 

  29. Hübner S, Rashkovetsky E, Kim YB, Oh JH, Michalak K, Weiner D, et al. Genome differentiation of Drosophila melanogaster from a microclimate contrast in Evolution Canyon, Israel. Proc Natl Acad Sci U S A. 2013;110:21059–64.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Hestmark G. Sex, size, competition and escape - strategies of reproduction and dispersal in Lasallia pustulata (Umbilicariaceae, Ascomycetes). Oecologia. 1992;92:305–12.

    Article  PubMed  Google Scholar 

  31. Giordani P, Incerti G, Rizzi G, Rellini I, Nimis PL, Modenesi P. Functional traits of cryptogams in Mediterranean ecosystems are driven by water, light and substrate interactions. J Veg Sci. 2013;25:778–92.

    Article  Google Scholar 

  32. Canu S, Rosati L, Fiori M, Motroni A, Filigheddu R, Farris E. Bioclimate map of Sardinia (Italy). J Maps. 2014;11:711–8.

    Article  Google Scholar 

  33. Cubero OF, Crespo A. Isolation of nucleic acids from lichens. In: Kranner I, Beckett R, Varma A, editors. Protocols in Lichenology. Berlin, Heidelberg: Springer Lab Manuals, Springer; 2002. p. 381–91.

    Chapter  Google Scholar 

  34. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015. doi:10.1093/bioinformatics/btv351.

    PubMed  Google Scholar 

  35. Sharma R, Xia X, Cano LM, Evangelisti E, Kemen E, Judelson H, et al. Genome analyses of the sunflower pathogen Plasmopara halstedii provide insights into effector evolution in downy mildews and Phytophthora. BMC Genomics. 2015;16:741.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Rubio-Piña JA, Zapata-Pérez O. Isolation of total RNA from tissues rich in polyphenols and polysaccharides of mangrove plants. Electron J Biotechnol. 2011;14:1–8.

    Google Scholar 

  37. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114-20.

  38. Conesa A, Götz S. Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int J Plant Genomics. 2008;619832. doi: 10.1155/2008/619832.

  39. Jones P, Binns D, Chang H-Y, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30:1236–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Sharma R, Thines M. FastQFS – a tool for evaluating and filtering paired-end sequencing data generated from high throughput sequencing. Mycol Prog. 2015;14:60.

    Article  Google Scholar 

  41. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Li R, Li Y, Fang X, Yang H, Wang J, Kristiansen K, et al. SNP detection for massively parallel whole-genome resequencing. Genome Res. 2009;19:1124–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Kofler R, Orozco-terWengel P, De Maio N, Pandey RV, Nolte V, Futschik A, et al. PoPoolation: a toolbox for population genetic analysis of next generation sequencing data from pooled individuals. PLoS One. 2011;6, e15925.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Kofler R, Pandey RV, Schlötterer C. PoPoolation2: identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq). Bioinformatics. 2011;27:3435–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. R Core Team. R: A Language and Environment for Statistical Computing. 2015.

    Google Scholar 

  47. Bates D, Maechler M, Bolker B, Walker S. Linear Mixed-Effects Models using ‘Eigen’ and S4. 2015.

    Google Scholar 

  48. Bretz F, Hothorn T, Westfall P, Hothorn T, Westfall P. Multiple comparisons using R. 2011.

    Google Scholar 

  49. Fisher RA. On the interpretation of χ2 from contingency tables, and the calculation of P. J R Stat Soc. 1922;85:87–94.

    Article  Google Scholar 

  50. Abdi H, O’Toole AJ, Valentin D, Edelman B. (2005). DISTATIS: the analysis of multiple distance matrices. Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05). 2005;1063-6919/05.

  51. Abdi H, Dunlop JP, Williams LJ. How to compute reliability estimates and display confidence and tolerance intervals for pattern classifiers using the Bootstrap and 3-way multidimensional scaling (DISTATIS). Neuroimage. 2009;45:89–95.

    Article  PubMed  Google Scholar 

  52. Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 2012;doi:10.1371/journal.pgen.1002967.

  53. Beerli P, Palczewski M. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics. 2010;185:313–26.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Pfenninger M, Patel S, Arias-Rodriguez L, Feldmeyer B, Riesch R, Plath M. Unique evolutionary trajectories in repeated adaptation to hydrogen sulphide-toxic habitats of a neotropical fish (Poecilia mexicana). Mol Ecol. 2015;24:5446–59.

    Article  PubMed  Google Scholar 

  55. Günther T, Coop G. Robust identification of local adaptation from allele frequencies. Genetics. 2013;195:205–20.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.

    Article  Google Scholar 

  57. Alexa A, Rahnenfuhrer J. topGO: Enrichment analysis for Gene Ontology. R package version 2.20.0. 2010.

    Google Scholar 

  58. Supek F, Bošnjak M, Škunca N, Šmuc T. Revigo summarizes and visualizes long lists of gene ontology terms. PLoS One. 2011;6, e21800.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Ronen R, Galun M. Pigment extraction from lichens with dimethyl sulfoxide (DMSO) and estimation of chlorophyll degradation. Environ Exp Bot. 1984;24:239–45.

    Article  CAS  Google Scholar 

  60. Françoise LD, Holger T, Marie-Laurence A, David D, Joël B. Oxidative stress regulation in lichens and its relevance for survival in coastal habitats. Adv Bot Res. 2014;71:467–503.

    Article  Google Scholar 

  61. Engelthaler DM, Roe CC, Hepp CM, Teixeira M, Driebe EM, Schupp JM, et al. Local population structure and patterns of Western Hemisphere dispersal for Coccidioides spp., the fungal cause of Valley Fever. Mbio. 2016;7:pii: e00550-16.

  62. Nosil P, Vines TH, Funk DJ. Perspective: Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution. 2005;59:705–19.

    PubMed  Google Scholar 

  63. Schluter D. Evidence for ecological speciation and its alternative. Science. 2009;323:737–41.

    Article  CAS  PubMed  Google Scholar 

  64. Keller I, Seehausen O. Thermal adaptation and ecological speciation. Mol Ecol. 2012;21:782–99.

    Article  CAS  PubMed  Google Scholar 

  65. Thompson JD. Plant Evolution in the Mediterranean. Oxford: Oxford Univ. Press; 2005.

    Book  Google Scholar 

  66. Robbins NM. The role of fungal stress responses in regulation of azole resistance. PhD Thesis. Toronto: University of Toronto; 2013.

    Google Scholar 

  67. Lasky JR, Des Marais DL, Lowry DB, Povolotskaya I, McKay JK, Richards JH, et al. Natural variation in abiotic stress responsive gene expression and local adaptation to climate in Arabidopsis thaliana. Mol Biol Evol. 2014;31:2283–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Gasch AP. The environmental stress response : a common yeast response to diverse environmental stresses. Top Curr Genet. 2002;1:11–70.

    Article  Google Scholar 

  69. Yew SM, Chan CL, Toh YF, Toh YF, Ngeow YF, Na SL, et al. The genome of newly classified Ochroconis mirabilis: Insights into fungal adaptation to different living conditions. BMC Genomics. 2016;17:1471–2164.

    Article  Google Scholar 

  70. Okada N, Ogawa J, Shima J. Comprehensive analysis of genes involved in the oxidative stress tolerance using yeast heterozygous deletion collection. FEMS Yeast Res. 2014;14:425–34.

    Article  CAS  PubMed  Google Scholar 

  71. Popa CV, Dumitru I, Ruta LL, Danet AF, Farcasanu IC. Exogenous oxidative stress induces Ca2+ release in the yeast Saccharomyces cerevisiae. FEBS J. 2010;277:4027–38.

    Article  CAS  PubMed  Google Scholar 

  72. Fridovich I. Fundamental aspects of reactive oxygen species, or what’s the matter with oxygen? Ann NY Acad Sci. 1999;893:13–8.

    Article  CAS  PubMed  Google Scholar 

  73. Llopis S, Querol A, Heyken A, Hube B, Jespersen L, Fernández-Espinar MT, et al. Transcriptomics in human blood incubation reveals the importance of oxidative stress response in Saccharomyces cerevisiae clinical strains. BMC Genomics. 2012;13:1471–2164.

    Article  Google Scholar 

  74. Vander Heiden MG, Choy JS, VanderWeele DJ, Brace JL, Harris MH, Bauer DE, et al. Bcl-x(L) complements Saccharomyces cerevisiae genes that facilitate the switch from glycolytic to oxidative metabolism. J Biol Chem. 2002;277:44870–6.

    Article  CAS  PubMed  Google Scholar 

  75. Vargas-Muñiz JM, Renshaw H, Richards AD, Lamoth F, Soderblom EJ, Moseley MA, et al. The Aspergillus fumigatus septins play pleiotropic roles in septation, conidiation, and cell wall stress, but are dispensable for virulence. Fungal Genet Biol. 2015;81:41–51.

    Article  PubMed  PubMed Central  Google Scholar 

  76. Sinha RP, Häder DP. UV-induced DNA damage and repair: a review. Photochem Photobiol Sci. 2002;1:225–36.

    Article  CAS  PubMed  Google Scholar 

  77. Solhaug KA, Gauslaa Y, Nybakken L, Bilger W. UV-induction of sun-screening pigments in lichens. New Phytol. 2003;158:91–100.

    Article  CAS  Google Scholar 

  78. Fuchs BB, Mylonakis E. Our paths might cross: The role of the fungal cell wall integrity pathway in stress response and cross talk with other stress response pathways. Eukaryot Cell. 2009;8:1616–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Kraus PR, Heitman J. Coping with stress: Calmodulin and calcineurin in model and pathogenic fungi. Biochem Biophys Res Commun. 2003;311:1151–7.

    Article  CAS  PubMed  Google Scholar 

  80. Aguilera J, Randez-Gil F, Prieto JA. Cold response in Saccharomyces cerevisiae: New functions for old mechanisms. FEMS Microbiol Rev. 2007;31:327–41.

    Article  CAS  PubMed  Google Scholar 

  81. Eckert AJ, Bower AD, González-Martínez SC, Wegrzyn JL, Coop G, Neale DB. Back to nature: ecological genomics of loblolly pine (Pinus taeda, Pinaceae). Mol Ecol. 2010;19:3789–805.

    Article  CAS  PubMed  Google Scholar 

  82. Keller SR, Levsen N, Olson MS, Tiffin PL. Local adaptation in the flowering-time gene network of balsam poplar, Populus balsamifera L. Mol Biol Evol. 2012;29:3143–52.

    Article  CAS  PubMed  Google Scholar 

  83. Pyhäjärvi T, Hufford MB, Mezmouk S, Ross-Ibarra J. Complex patterns of local adaptation in teosinte. Genome Biol Evol. 2013;5:1594–609.

    Article  PubMed  PubMed Central  Google Scholar 

  84. Cowan IR, Lange OL, Green TGA. Carbon-dioxide exchange in lichens: determination of transport and carboxylation characteristic. Planta. 1992;187:282–94.

    Article  CAS  PubMed  Google Scholar 

  85. Colesie C, Scheu S, Green TGA, Weber B, Wirth R, Büdel B. The advantage of growing on moss: facilitative effects on photosynthetic performance and growth in the cyanobacterial lichen Peltigera rufescens. Oecologia. 2012;169:599–607.

    Article  PubMed  Google Scholar 

  86. Pintado A, Valladares F, Sancho LG. Exploring phenotypic plasticity in the lichen Ramalina capitata: morphology, water relations and chlorophyll content in north- and south-facing populations. Ann Bot. 1997;80:345–53.

    Article  Google Scholar 

  87. Sojo F, Valladares F, Sancho LG. Structural and physiological plasticity of the lichen Catillaria corymbosa in different microhabitats of the maritime antarctica. Bryologist. 1997;100:171–9.

    Article  Google Scholar 

  88. Colesie C, Williams L, Büdel B. Internal thallus water status in the soil crust lichen Psora decipiens is optimised via a high phenotypic plasticity. Lichenol. 2017;In press.

  89. Hoffmann A, Griffin P, Dillon S, Catullo R, Rane R, Byrne M, et al. A framework for incorporating evolutionary genomics into biodiversity conservation and management. Clim Chang Responses. 2015;2:doi: 10.1186/s40665-014-0009-x.

  90. Futschik A, Schlötterer C. The next generation of molecular markers from massively parallel sequencing of pooled DNA samples. Genetics. 2010;186:207–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Fabian DK, Kapun M, Nolte V, Kofler R, Schmidt PS, Schlötterer C, et al. Genome-wide patterns of latitudinal differentiation among populations of Drosophila melanogaster from North America. Mol Ecol. 2012;21:4748–69.

    Article  PubMed  PubMed Central  Google Scholar 

  92. Gautier M, Foucaud J, Gharbi K, Cézard T, Galan M, Loiseau A, et al. Estimation of population allele frequencies from next-generation sequencing data: pool-versus individual-based genotyping. Mol Ecol. 2013;22:3766–79.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We are grateful to Ingo Ebersberger and Bastian Greshake (Frankfurt) for support with genome assembly. We thank Bob O'Hara and Garima Singh (Frankfurt) for helpful discussions. Eva-Maria Gerstner assisted with ArcGIS, Claus Weiland with cluster access, and Deepak Kumar Gupta (all Frankfurt) with Perl scripts. Anna Sadowska-Deś (Frankfurt) provided valuable information on field sites. We thank the following for granting access to the bioclimatic data of the sites: Simona Canu, Michele Fiori, Andrea Motroni (Dipartimento Meteoclimatico of ARPA Sardegna), Emmanuele Farris (University of Sassari), and Leonardo Rosati (University of Basilicata). Environmental association analysis was performed on the LOEWE cluster of the Center for Scientific Computing (CSC) of the Goethe University in Frankfurt am Main.

Funding

The work was funded by the research funding program Landes-Offensive zur Entwicklung Wissenschaftlich-Oekonomischer Exzellenz (LOEWE) of Hesse’s Ministry of Higher Education, Research, and the Arts through the Senckenberg Biodiversity and Climate Research Centre (BiK-F) and the Integrative Fungal Research Program (IPF).

Availability of data and materials

The annotated genome sequence is available at the European Nucleotide Archive under the accession numbers FWEW01000001 to FWEW01003891.Illumina sequence data are available at the European Nucleotide Archive under the accession numbers ERR1110280 to ERR1110285.

Sync files from Popoolation2, DNA sequences linked to Additional file 3, and the complete lists of FST and Bayenv2.0 outliers are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.255hj.

Authors’ contributions

FDG and IS designed research; FDG, JO, and BB collected the data; RS, BM, AM contributed scripts; FDG, RS, AM, GR and BB analyzed data; and FDG, MP, and IS wrote the manuscript with contributions from the other authors. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not required.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations

  1. Senckenberg Biodiversity and Climate Research Centre (BiK-F), Senckenberganlage 25, 60325, Frankfurt am Main, Germany

    Francesco Dal Grande, Anjuli Meiser, Gregor Rolshausen, Bagdevi Mishra, Marco Thines, Jürgen Otte, Markus Pfenninger & Imke Schmitt

  2. Institut für Ökologie, Evolution und Diversität, Goethe-Universität Frankfurt, Max-von-Laue-Str. 9, 60438, Frankfurt am Main, Germany

    Rahul Sharma, Anjuli Meiser, Markus Pfenninger & Imke Schmitt

  3. Plant Ecology and Systematics, Biology Department, University of Kaiserslautern, 67653, Kaiserslautern, Germany

    Burkhard Büdel

Authors
  1. Francesco Dal Grande
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Rahul Sharma
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Anjuli Meiser
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Gregor Rolshausen
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Burkhard Büdel
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Bagdevi Mishra
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Marco Thines
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jürgen Otte
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Markus Pfenninger
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Imke Schmitt
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding authors

Correspondence to Francesco Dal Grande or Imke Schmitt.

Additional files

Additional file 1:

Temperature data (in °C) collected between May 28th 2014 and June 2nd 2015 from loggers positioned at the level of L. pustulata thalli (two loggers per population). (PDF 113 kb)

Additional file 2:

Plot of the first two axes of a PCA showing high correlation between 19 bioclimatic variables and altitude for the six populations along the cline. (PDF 524 kb)

Additional file 3:

Genetic characteristics of the six loci used to genotype 18 L. pustulata samples for anatomical and physiological measurements. (PDF 113 kb)

Additional file 4:

Distribution of nucleotide diversity, measured as Tajima's D, across 309 non-overlapping sliding windows of 10-kb. For the six scatterplots on the left, the data are plotted in ascending order of mean diversity across the six populations (solid line) together with the respective 95% confidence interval (grey envelope). Black crosses depict population-specific outlier SNPs. For the boxplot on the right, the data are summarized for each population across all windows. (PDF 517 kb)

Additional file 5:

a) The mean and standard deviation for three standard population genetic parameters, π, Watterson's θ (θW), and Tajima's D in non-overlapping 10-kb windows across the genome of L. pustulata. b) Pairwise Tukey contrasts based on linear mixed effect models for each measure and population. (PDF 446 kb)

Additional file 6:

NJ trees of the correlation matrix obtained with Bayenv2.0 and the pairwise FST values. (PDF 165 kb)

Additional file 7:

Average pairwise FST (lower triangle) and pairwise correlations of Bayenv2.0 correlation matrix among allele frequencies (upper triangle). (PDF 81 kb)

Additional file 8:

Three-(a) and four-(b) population test of the L. pustulata populations. The test was performed on 722,401 high confidence SNPs. Bold values indicate significant f3 (if z < 0) or f4 statistics (if |z| < 3). (PDF 236 kb)

Additional file 9:

Parameter estimates from MIGRATE-N based on 49 loci for pairs of three genetic groups of L. pustulata found along the altitudinal gradient (A: populations 1 to 4, B: population 5; C: population 6). (PDF 110 kb)

Additional file 10:

Allele frequency and annotation for the 42 environmentally associated (top Z Bayenv2.0 score), highly differentiated (top 0.5% FST distribution) SNPs. (PDF 126 kb)

Additional file 11:

GO enriched categories for a) top 0.5% highly differentiated SNPs, b) top 0.5% highly differentiated SNPs differentially fixed between populations 1–4 and population 6. (PDF 228 kb)

Additional file 12:

List of significantly enriched GO terms for environmentally associated genes. Gene names in bold indicate genes containing top Z Bayenv2.0 SNPs. (PDF 103 kb)

Additional file 13:

Neighbor-joining tree of the aligned concatenated gene fragments of the six loci at which we genotyped 18 thalli L. pustulata used for the anatomical and ecophysiological measurements. Samples in bold and underlined (N = 6) were used for the gas exchange experiments. (PDF 165 kb)

Additional file 14:

Light (a), temperature (b) curves for 6 thalli representing two genetic groups (1: pop. 6, high altitude; 2: populations 1 to 5, low altitude). (PDF 556 kb)

Additional file 15:

Net photosynthesis at optimal thallus water content and optimal temperature for genetic group A (pop. 6, high altitude) and B (populations 1 to 5, low altitude). (PDF 94 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Dal Grande, F., Sharma, R., Meiser, A. et al. Adaptive differentiation coincides with local bioclimatic conditions along an elevational cline in populations of a lichen-forming fungus. BMC Evol Biol 17, 93 (2017). https://doi.org/10.1186/s12862-017-0929-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-017-0929-8

Keywords

BMC Ecology and Evolution

ISSN: 2730-7182

Contact us