Anthropogenic climate change has led to unprecedented shifts in temperature across many ecosystems. In a context of rapid environmental changes, acclimation is an important process as it may influence the capacity of organisms to survive under novel thermal conditions. Mechanisms of acclimation could involve upregulation of stress response genes involved in protein folding, DNA damage repair and the regulation of signal transduction genes, along with a simultaneous downregulation of genes involved in growth or cell cycle, in order to maintain cellular functions and equilibria. We transplanted Lobaria pulmonaria lichens originating from different forests to determine the relative effects of long-term acclimation and genetic factors on the variability in expression of mycobiont and photobiont genes. We found a strong response of mycobiont and photobiont to high temperatures, regardless of sample origin. The green-algal photobiont had an overall lower response than the mycobiont. The gene expression of both symbionts was also influenced by acclimation to transplantation sites and by genetic factors. Lobaria pulmonaria seems to have evolved powerful molecular pathways to deal with environmental fluctuations and stress and can acclimate to new habitats by transcriptomic convergence. Although L. pulmonaria has the molecular machinery to counteract short-term thermal stress, survival of lichens like L. pulmonaria depends mostly on their long-term positive carbon balance, which can be compromised by warmer temperatures and reduced precipitation, and both these outcomes have been predicted for Central Europe in connection with global climate change.

We established reciprocal transplants between four Austrian sites: Koralpe (K), Kleinsölk Valley (S) Tamischbach Canyon (T) and Raab Canyon (R), which received additional transplants from six Swiss sites (Table 1) (Fig 1B, Fig S1). In early summer 2015, Silke Werth (SW) set up reciprocal transplants between four ‘recipient’ sites in Austria (Lowland; Southern, Central and Northern Alps) differing in elevation and macroclimatic regime (K, S, R, T) (Fig 1B, Fig S1, Table S1). In 2015, lichen lobes originating from two sites in the Swiss Jura mountains were transplanted to the recipient sites in Austria, and the transplantation was repeated in 2016 with lobes from two additional Swiss Jura Mountain sites. Furthermore, in early summer 2016, additional reciprocal transplants were made between the Austrian sites. In each case, parts of large thalli were collected from Austria (A) or Switzerland (W), split into four lobes, and each lobe was transplanted to a different recipient site. After acclimatization to the recipient sites over at least three vegetation periods, the transplants were harvested on two consecutive days in August of 2019 during a time period with constant weather conditions. They were grown under controlled, cold (4 °C) conditions with medium light (150-200 PAR) for 24 hours to adjust to the lab conditions so that the results would not simply reflect the prevailing local weather or light conditions during sampling. At a temperature of 4 °C, we took three samples of marginal tissue (ca. 25 mm2) from a given lobe of each fully hydrated lichen to include lichen lobes originating from the two Origin Sites (W and A) that were transplanted to each of the four transplantation sites (K, R, S, T). The samples were immediately flash-frozen in liquid nitrogen and stored at -80 °C until further processing. We transferred the lichen lobes into a growth chamber at 15°C and we kept them there for two hours and took another three samples of the lichen lobes exposed to this temperature. Lastly, we increased the temperature to 25 °C and kept the lichen lobes for two more hours at this temperature before collecting three additional samples from each lobe (Fig 1B) RNA extraction was performed according to the manufacturer’s protocol using the RNeasy Plant Mini Kit (Qiagen, Stockach, Germany). After elution of extracted RNA in RNase-free H2O, RNA concentration and quality were checked using first a NanoDrop ND-1000 UV/Vis-Spectrophotometer (Thermo Scientific, Carlsbad CA, USA) and then a Qubit 2.0 fluorometer with the DNA HS assay kit (Qiagen, Stockach, Germany). First, RNA concentration and quality were checked using a high-sensitivity RNA chip (Agilent Bioanalyzer). RNA sequencing was performed using an adapted version of the prime-seq protocol (step-by-step protocol: dx.doi.org/10.17504/protocols.io.s9veh66), but with adjustments as follows. Raw fastq files were processed using zUMIs software (Parekh et al. 2018).The normalized counts (vst) were used to perform principal components analysis (PCA) (R version 3.6.3). The Adonis test (a randomization/Monte Carlo permutation test was utilized to test for significant relationships between gene expression variation in samples exposed to different temperatures and from different transplantation sites. Differential expression analysis was carried out by DESeq2 in R. We used the Wald test for hypothesis testing when comparing differentially express transcripts between two groups, and only transcripts with a log2-fold change >1.5 and adjusted Benjamini-Hochberg p-value <0.0001 were considered differentially expressed. We used the adjusted Benjamini-Hochberg p-value and such low p-value levels to reduce the false positive rate. The study design allowed us to identify transcripts that were up and downregulated in the lichen thalli, among the three temperature treatments (acclimation effects), among the two origins (genetic effects), or among the four transplantation sites (acclimatization effects). We made three Wald test comparations for the different temperatures: 1) 4 °C vs 25 °C 2) 4 °C vs 15 °C 3) 15 °C vs 25 °C. Six different Wald test comparations between the transplantation sites: 1) K vs S 2) K vs R 3) K vs T 4) S vs R 5) S vs T 6) T vs R. In addition, we made a Wald test comparison between the two L. pulmonaria samples origin: A vs W. These Wald test comparations were made for the mycobiont and for the photobiont. While DESeq2 is useful for conducting Wald pairwise comparisons, the Weighted Gene Co-Expression Network Analysis (WGCNA) (Langfelder & Horvath 2008) package in R allows to examine patterns of gene expression across all treatments simultaneously. WGCNA was used to generate a correlation network to identify clusters of similarly expressed transcripts that clustered into gene modules, using a cutoff of minimum 30 transcripts per module. A module’s so-called ‘eigengene’ summarizes the overall module expression profile. Modules with highly correlated eigengenes were merged using a threshold value by the ‘dynamicMergeCut’ WGCNA function in R and the clusters were summarized by the eigengene. Then we related the modules to the explanatory variables as temperature treatment (acclimation effects), transplantation site (acclimatization effects), and origin site (genetic effects) using eigengene network methodology. A heatmap was created to visualize all significant correlations between modules and treatments in WGCNA. Then we identified transcripts with high significance for each of the explanatory variables for the mycobiont and the photobiont separately. Gene ontologies of the Lobaria pulmonaria reference genome hosted at DOE-JGI, https://mycocosm.jgi.doe.gov/LobpulSc1/LobpulSc1.home.html (date accessed: 17/07/2020) were used to match the differentially expressed (DE) transcripts with protein and GO-term information. Lists of gene ontology (GO) terms were generated for the differentially expressed transcripts found with DESeq2. When there was no GO term information available, we annotated the respective transcript manually by homology searches via the translated nucleotide BLAST algorithm (blastx) (Altschul et al. 1997) against the Uniprot database (https://www.uniprot.org/blast/). We searched for homology of L. pulmonaria mycobiont transcript sequences to protein sequences available for several whole-genome sequenced species of fungi. Then we searched for homology of the S. reticulata photobiont’s transcript sequences to protein sequences available for algae. We used the target databases‘ fungi’ and ‘algae’ with default algorithm parameters (expect threshold=10, matrix=auto, filtering=none, gapped=yes, hits=250) for the Blast searches. We considered a locus homologous if the e-value returned was smaller than 1.0 * 10-5. We only used this method for the differentially expressed transcripts found with DESeq2 analysis, if we did not find any GO term information available for those transcripts. The GO term gene information for the WGCNA is presented in the Supporting Information (Tables S1- S2).

We upload all the gene data counts needed for our analysis

We included the R script with our analysis

We included all the DE genes list from the Wald test comparing variables. We also included the list of the total genes found in the modules of the WGCNA.

The data with M is from the mycobiont and the P is from the photobiont