Stressed out: the effects of heat stress and parasitism on gene expression of the lichen-forming fungus Lobaria pulmonaria

Abstract

Gene expression variation can be partitioned into different components (regulatory, genetic and acclimatory effects) but for lichen-forming fungi, the relative importance of each of these effects is unclear. Here, we studied gene expression in the lichen-forming fungus Lobaria pulmonaria in response to thermal stress and parasitism by the lichenicolous fungus Plectocarpon lichenum. Our experimental procedure was to acclimate lichen thalli to 4 °C over three weeks and then expose them to 15 °C and 25 °C for 2 hours each, sampling infected and visually asymptomatic thalli at each temperature. Quantitative Real-Time PCR was utilized to quantify gene expression of six candidate genes, normalizing expression values with two reference genes. We found that two genes encoding heat shock proteins (hsp88 and hsp98), two polyketide synthase genes (rPKS1, nrPKS3) and elongation factor 1-1-α (efa) were upregulated at higher temperatures. Moreover, we observed higher expression of hsp98 at 25 °C in samples infected by P. lichenum than in uninfected samples. Finally, in partial redundancy analyses, most of the explained variation in gene expression was related to temperature treatment; genetic variation and long-term acclimatization to sites contributed far less. Hence, regulatory effects (i.e. direct adjustments of gene expression in response to the temperature change) dominated over genetic and acclimatory effects in the gene expression variability of L. pulmonaria. This study suggests that L. pulmonaria could become a valuable lichen model for studying heat shock protein responses in vivo.

Introduction

Throughout the history of life, organisms have been challenged to survive in habitats that are not stable but subjected to fluctuations in important abiotic conditions such as temperature, humidity, pH and UV-light (MacKenzie et al. 2001, 2004; Evans et al. 2013; Hamann et al. 2016). In order to deal with those changing conditions, the ability to regulate the expression of stress-related genes is vital (Evans et al. 2013). Investigations of both eukaryotes and prokaryotes have shown that gene expression plays a crucial role in the tolerance of extreme conditions such as drought (Wang et al. 2015; Carmo et al. 2019), temperature and salinity stress (Jamil et al. 2011; Che et al. 2013; Zhang et al. 2017), as well as exposure to toxins (Whitehead et al. 2010). Understanding the mechanisms and different pathways of this gene-expression response to stressful conditions is important for obtaining better insights into survival mechanisms and the interplay of organisms with their environment. Environmental stress response has been the subject of various studies in many different organisms (Mizoguchi et al. 1997; Gasch 2007; Dixon et al. 2020; Terhorst et al. 2020). In fungi, environmental stress response was first described in Saccharomyces cerevisiae (Gasch et al. 2000). Stress genes play an important role in carbohydrate metabolism, response to oxidative stress, intracellular signalling, DNA-damage repair and protein metabolism, especially protein folding (Gasch 2007; O'Meara et al. 2019).

One common environmental stressor, which organisms are confronted with, is thermal stress (Arshad et al. 2017). In most habitats, organisms have to deal with more or less rapidly changing temperatures. However, responses to thermal stress have also become an important issue due to the rapid increase in temperatures, and higher fluctuations and extremes, because of global climate change (IPCC 2021). Global mean surface temperatures will continue to increase in the first half of the 21st century, with the level of increase depending on the quantity of future man-made CO₂-emissions (IPCC 2021). Heat shock response represents one of the important mechanisms for organisms to adapt to stressful conditions at the cellular level (O'Meara et al. 2019).

In response to environmental stress, gene expression needs to be regulated to a new cellular equilibrium to ensure cell survival. We hereafter refer to the variation in gene expression that is involved in maintaining cellular equilibrium under temperature stress as ‘regulatory effects’. In order to respond to and survive thermal stress, organisms need to be able to sense heat; and as a response, they need to conduct an adequate regulation of genes that might prevent or reduce the damage caused by high temperatures. In general, heat stress can be sensed through two effects: first, the accumulation of denatured proteins which results in the activation of a heat shock factor (Franzmann et al. 2008); second, changes in thermosensitive structures such as DNA, RNA, proteins or lipids that serve as primary sensors. These can either have a direct effect or activate signal transduction pathways such as the very conserved mitogen-activated kinase (MAPK) pathways, important in the stress responses of filamentous ascomycetes (Hagiwara et al. 2016). The first reaction initiated by these signalling pathways can include fast responses such as the use of previously synthesized proteins or the regulation of channels and transporters. Whereas the main heat shock response is carried out through gene regulation leading to a major change in transcriptional patterns after a few minutes (Albrecht et al. 2010; Roncarati & Scarlato 2017). Many genes are simultaneously downregulated under stress conditions (e.g. those involved in cell-cycle, RNA metabolism and synthesis of proteins), some reaching several maxima in expression over a period of two hours or fluctuating over time (Albrecht et al. 2010; de Nadal et al. 2011; Takahashi et al. 2017).

An important reaction to thermal stress is the expression of genes encoding heat shock proteins (HSPs). HSPs are able to unfold and refold proteins which become misfolded because of heat exposure (Albrecht et al. 2010). The heat-induced upregulation of HSPs has been shown in many organisms including prokaryotes and eukaryotes, revealing many HSP families that interact and regulate each other in different pathways (Plesofsky-Vig & Brambl 1998; Miot et al. 2011; Smith et al. 2012; Li & Buchner 2013; Park et al. 2015). The heat shock protein gene hsp88 of an entomopathogenic fungus has been cloned and characterized by Park et al. (2014). Under thermal stress, hsp88 was 15–55-fold upregulated in the lichen-forming fungus Peltigera membranacea (Steinhäuser et al. 2016). An important heat shock protein gene in Aspergillus fumigatus is hsp98 (Do et al. 2009), and this gene was also upregulated under thermal stress in P. membranacea (Steinhäuser et al. 2016).

While the increased expression of heat shock protein genes is a universal and well-known reaction to environmental stressors, another reaction that could possibly be linked to stressful conditions is the production of polyketides in fungi (Timsina et al. 2013). Polyketides are secondary metabolites featuring antimicrobial, antitumour, immunosuppressive, antifungal and antiparasitic properties and they are therefore not only of great relevance for pharmaceutical purposes (Nivina et al. 2019), but also of interest for answering physiological and ecological questions. Polyketides have been suspected to protect organisms from environmental stresses such as high light levels and drought, or from herbivory and fungal parasites (Lawrey 1986, 1989; Torzilli et al. 1999; Gauslaa & McEvoy 2005; Timsina et al. 2013). The biosynthesis of polyketides out of 2-, 3- or 4-carbon compounds is catalyzed by polyketide synthases (PKSs), which are large multi-enzyme systems with a molecular weight of up to 10 000 kDa (Khosla et al. 1999). Type I PKSs are large proteins consisting of several functional domains and type III PKSs are simpler enzymes catalyzing the formation of a product within a single active site (Nivina et al. 2019). Non-reducing PKSs characteristically catalyze the synthesis of aromatic polyphenols but fungal reducing PKSs reduce beta-carbons with different domains to form reduced aromatic rings or aliphatic rings, for example macrolides (Bertrand & Sorensen 2018). Generally, there is a connection between polyketide production in lichens and abiotic conditions such as nutrient supply, substratum pH and light, with the production being higher under stressful conditions and negatively correlated with growth (Armaleo et al. 2008; Timsina et al. 2013). Thus, it is conceivable that heat stress would lead to an upregulation of polyketide synthase genes, causing a corresponding increase in polyketide production. In the lichen-forming fungus Lobaria pulmonaria (L.) Hoffm. (lichenized ascomycetes, Peltigerales), three major carbon-based secondary compounds are produced by PKS genes: stictic, constictic and norstictic acids, as well as some chemically related minor compounds (Bidussi et al. 2013; Gauslaa et al. 2013). The depsidones, norstictic and stictic acid, are produced via the acetate-polymalonate pathway (Ranković & Kosanić 2019). Some lichen secondary compounds, including those of L. pulmonaria, have anti-herbivore and antibiotic properties (Suleyman et al. 2003; Asplund & Gauslaa 2008; Nybakken et al. 2010). Some secondary compounds such as lecanoric acid may also have antifungal properties, preventing lichen colonization by certain lichenicolous fungi (Lawrey 1989, 2000; Lawrey & Diederich 2003), and some may be useful as anti-cancer drugs (Shrestha & St. Clair 2013; Dar et al. 2021; Yang et al. 2021).

The lichenicolous fungus Plectocarpon lichenum (Sommerf.) D. Hawksw. forms conspicuous darkish brown structures on thalli of Lobaria pulmonaria; these structures represent stromata made from a combination of hyphae of the lichenicolous fungus and of its lichen host (Bergmann & Werth 2017). A recent study based on qPCR found that the mycelium of this lichenicolous fungus is localized mainly in the stromata, with only a very low signal being detected directly adjacent to stromata (Bergmann & Werth 2017). Areas including stromata have on average twice the biomass when compared to adjacent asymptomatic thallus parts, and thalli infected by P. lichenum most often contain many stromata (Bergmann & Werth 2017). Thus, it is conceivable that P. lichenum taps substantially into the overall carbon pool of L. pulmonaria. Thalli of L. pulmonaria infected by P. lichenum were found to have a significantly reduced amount of carbon-based secondary compounds (Asplund et al. 2016). Similarly, in Lobarina scrobiculata (Scop.) Nyl. ex Cromb., polyketide concentration was reduced to less than half in thalli infected by the lichenicolous fungus Plectocarpon scrobiculatae Diederich & Etayo, when compared to uninfected thalli (Merinero et al. 2015). Either infections by Plectocarpon lead to an overall downregulation of PKS genes in the parasitized thalli, or the lichenicolous fungi might degrade the lichen's secondary compounds with extracellular enzymes (Lawrey 2000). The first hypothesis can be tested by an analysis of differential expression of PKS genes.

Abiotic conditions such as different habitats can also influence gene expression (e.g. MacFarlane & Kershaw 1980; Cheviron et al. 2008; Whitehead et al. 2010; Steinhäuser et al. 2016). Habitat-related differential gene expression could be composed of both genetic and acclimatory factors (Cheviron et al. 2008; Whitehead et al. 2010; Palumbi et al. 2014). If the differences in gene expression are caused by long-term physiological acclimatization effects, they should vanish after acclimation to common conditions in the laboratory, or in a common garden experiment. Lichen populations grown in the laboratory or a common garden can, however, retain site-specific differences in gene expression (Steinhäuser et al. 2016) or physiological state (MacFarlane & Kershaw 1980; Schipperges et al. 1995). These studies suggest that there might be a substantial genetic component to variation in gene expression. However, the relative importance of the genetic component has not yet been scrutinized.

The main aim of this study was to obtain a better understanding of gene expression variation in response to increased temperatures and its partitioning into different factors in the lichen-forming fungus L. pulmonaria. At the onset of our study, it was not known at which temperature heat shock is induced in L. pulmonaria. Therefore, we first investigated the expression patterns of L. pulmonaria heat shock protein and polyketide synthase genes exposed to different temperatures to quantify the regulatory component of gene expression variation. The specific question we asked was, does thermal stress caused by a temperature increase from 4 °C to 15 °C and then to 25 °C result in differential expression of heat shock protein and polyketide synthase genes?

Given that earlier studies indicated that the concentration of lichen secondary metabolites was reduced in Lobaria pulmonaria thalli parasitized by P. lichenum (Asplund et al. 2016, 2018), we hypothesized that the presence of the lichenicolous fungus P. lichenum would have an effect on the expression of polyketide synthase genes, leading to their down-regulation (biotic component of gene expression variation). However, since polyketide production may increase due to environmental stress, we expected higher gene expression in polyketide synthases under thermal stress conditions.

Furthermore, we examined whether physiological long-term acclimatization had a long-lasting effect on the physiological state of individuals, persisting as collecting site-related differences even after acclimation to common laboratory conditions (acclimatory component). To address this issue in our study, we compared thalli of L. pulmonaria from a population in Austria with one in Tenerife after acclimating them to common laboratory conditions. Finally, we related gene expression variation to genetic distance to quantify the genetic component of gene expression variation. To assess the relative roles of the regulatory, acclimatory, biotic, and genetic components of gene expression variation, a variance partitioning approach was used.

Materials and Methods

Collection of lichen samples

Samples were collected in February 2015 from a site in Austria (AU7) and a site in Tenerife (ST7). AU7 was chosen as one of four populations of Lobaria pulmonaria described in the literature, located in the Ennstaler Alps at Tamischbachgraben (47°38′N, 14°41′E) at c. 700 m above sea level (Scheidegger et al. 2012). Five thalli (AU7-01–AU7-05) of similar size were collected from trunks of sycamore maple (Acer pseudoplatanus L.). In order to collect different genotypes, the thalli were taken from trees at a distance of at least 20 m. Site ST7 was located in a pine (Pinus canariensis C. Sm. ex DC.) forest in Tenerife, the Canary Islands (28°24.51096′N, 16°25.06404′W, 1560 m a.s.l.); this site is frequently exposed to fog. From this site, ten thalli with Plectocarpon lichenum (ST7-11–ST7-20) and ten without (ST7-01–ST7-10) were gathered. Samples with Plectocarpon infection contained stromata visible to the naked eye. Samples were collected at a distance of at least 10 m from each other. All thalli were stored dry and in darkness at a temperature of c. 4 °C for 5 days until the beginning of the experiment.

Acclimation phase and temperature treatment

The thalli were placed in Petri dishes lined with filter paper, which was previously rinsed with distilled water to create a neutral substratum for the lichens. In order to allow them to acclimate, the lichens were grown in a styrofoam box for 3 weeks in a cold room at 4 °C under constant light conditions of 62.4 lx (in the middle of the box) to 38.4 lx (on the edge of the box). To achieve as equal conditions as possible, the samples in the middle and on the edge were swapped periodically. They were watered frequently with dH₂O, but allowed to dry out every few days in order to avoid mould and to simulate the natural change of metabolically active and inactive phases due to re- and dehydration. At the end of the acclimation period at 4 °C, tissue samples were taken for RNA extractions from fully hydrated lobes by cutting off 5 × 5 mm pieces from the edge of each thallus and placing them in ice-cooled RNA stabilization solution (3.53 M ammonium sulphate, 16.7 mM sodium citrate, 13.3 mM EDTA, pH = 5.2) (De Wit et al. 2012). From the samples ST7-11–ST7-20, only areas with no visible growth of P. lichenum were used for sampling since we were investigating if presence of the lichenicolous fungus influenced gene expression of the entire thallus, rather than just the symptomatic thallus parts. The light source and the lichen thalli were then moved to a plant growth chamber in which the thalli were exposed to higher temperatures, first to 15 °C and then to 25 °C. Each temperature treatment was kept for 3 h prior to tissue sampling as described above; this experimental set-up was similar to the one used by Steinhäuser et al. (2016).

RNA extraction and reverse transcription to cDNA

The samples taken after exposure to 4 °C and 15 °C were immediately used for RNA extraction, while the 25 °C samples were stored overnight at −80 °C. Successful RNA extractions for Lobaria pulmonaria have been reported using TRI Reagent (Doering et al. 2014), therefore this extraction chemistry was used. Samples were homogenized in TRI Reagent (Sigma Aldrich) using Tissue Lyser II (Qiagen) with a 3 mm stainless steel polishing bead (Kugel Pompel, Austria). RNA extraction was performed using 2 ml Heavy Gel Phase Lock gel tubes (5Prime) based on the manufacturer's instructions. RNA concentration was quantified using a P-Class NanoPhotometer (Implen). The RNA concentrations ranged from 120–620 ng/μl. To remove the remaining genomic DNA from the samples, a digest of genomic DNA was performed with the RNase-Free DNase Set (Qiagen). The RNA was pipetted to a mix of DNase I, RDD buffer and RNase free water and then the mix was incubated in a thermocycler (AlphaMetrix Biotech) at 37 °C for 15 min and at 75 °C for 5 min. After the DNA digest, the RNA concentration was quantified again and all samples were diluted to the same concentration (100 ng/μl) to enable quantitative comparisons.

For cDNA synthesis, 20 μl of digested RNA were pipetted to a mix of 4 μl 10× RT random primers, 1.6 μl dNTP mix (4 mM each), 4 μl 10× RT buffer, 2 μl MultiScribe Reverse Transcriptase (100 U) and 8.4 μl RNase-free water, using reagents and protocols provided with the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems). The assays were incubated at 25 °C for 10 min, at 37 °C for 120 min, at 85 °C for 5 min and then cooled to 4 °C. After the cDNA synthesis, the samples were diluted with 160 μl of RNase free H₂O to reach a final cDNA concentration of 10 ng/μl.

Selection of genes

As reference genes, we utilized two that play a vital role in metabolism which should have constant expression between different temperatures: beta-tubulin (bet) and glyceraldehyde 3-phosphate dehydrogenase (gpd). These genes had previously been validated in other studies of lichen-forming fungi (Joneson et al. 2011; Miao et al. 2012).

For the candidate genes, we focused on genes relevant under stressful conditions (especially genes encoding heat shock proteins that are likely to change in expression due to increasing temperatures) and we chose reducing and non-reducing types of polyketide synthase I. Table 1 shows the list of genes we considered as reference or candidate genes with putative functions. A BLASTX analysis was performed to verify the identity of loci (Altschul et al. 1997). Loci were selected based on PCR amplification (specific amplification, i.e. a single amplicon) and qPCR results. The sequences of all tested loci were deposited in GenBank (Accessions KX866397–KX866407).

Table 1. Reference and candidate genes used to study Lobaria pulmonaria gene expression variation in response to increased temperatures. Table headings are as follows: GenBank Accession (Accession); gene abbreviation (Gene); gene name (Name); alignment coordinates of blast hit on the L. pulmonaria genome (Coord. LPU); name of gene model from the L. pulmonaria Scotland JGI v1.0 reference genome (Gene model LPU); protein ID associated with L. pulmonaria gene model (ProteinID); KOG functional class assignment (KOG class); description of KOG function (KOG descr.); KOG ID; number of exons (No. exons); e-value from BLASTN analysis against the L. pulmonaria reference genome (LPU e-value); percent identity of blast hit to L. pulmonaria reference genome (Id LPU) (%). The loci nrPKS3 and nrPKS3' are exons of the same polyketide synthase gene.

As candidate and reference genes, we considered only conserved regions of the genes based on 454 genomic data of the mycobiont Lobaria pulmonaria (C. Scheidegger, unpublished data). Further information on the 454 data is provided in Werth et al. (2013); the multispore mycobiont culture F2, which was used to obtain the data, is described in Widmer et al. (2010) and Cornejo et al. (2015). Using partially sequenced mycobiont genomic data, we obtained genomic sequences of heat shock protein and PKS genes of L. pulmonaria, based on sequence similarity with protein sequences and a DNA sequence from GenBank. The following protein sequences were used to find L. pulmonaria sequences of PKS genes: ABV71377 (L. pulmonaria), BAN29051 (Lobaria orientalis (Asahina) Yoshim.), ABV71378 (L. scrobiculata), AEE87273, ADF28669, AEE87274, ADF28670, AEE65376, AEE65375, AEE65373, ADF28668, AEE65377, AEE65374, AEE65372 (Peltigera membranacea (Ach.) Nyl.), and a DNA sequence (EF363900, L. pulmonaria). The following GenBank Accessions were used to find stress genes: ACV03836 (Msn2, Aspergillus parasiticus Speare), EDN02919 (hsp88, Ajellomyces capsulatus (Kwon-Chung) McGinnis & Katz), EDP56763 (heat shock protein gene hsp98/hsp104/ClpA, Aspergillus fumigatus Fresen.), EYE93161 (putative signal peptide peptidase, a gene involved in signal transduction in Aspergillus ruber Thom & Church), AAR30137 (putative histidine kinase HHK2p, Fusarium verticillioides (Sacc.) Nirenberg), and elongation factor 1-α (AFQ55277), which has been shown to function as a molecular chaperone upregulated under heat conditions and salt stress in plants (Shin et al. 2009). Reference genes were obtained through sequences of β-tubulin (AFJ45056, P. membranacea) and glyceraldehyde 3-phosphate dehydrogenase (AFJ45057, P. membranacea). Only blast hits with an e-value < 10⁻⁴⁰ were retained. After inspecting alignments, we selected genes with a high similarity to hsp88, hsp98, putative signal peptide peptidase, putative histidine kinase HHK2p, reducing and non-reducing types of PKS I, actin, β-tubulin, glyceraldehyde 3-phosphate dehydrogenase, and elongation factor 1-α.

The Lobaria pulmonaria Scotland v.1.0 reference genome was released on JGI after we performed our experiment. To assess the correspondence of our gene set to NCBI gene models and annotations, we blasted each gene against the Lobaria pulmonaria genome on JGI MycoCosm (Table 1) (https://genome.jgi.doe.gov/pages/blast-query.jsf?db=Lobpul1; accessed 12 June 2018).

Primer design and efficiency

To design primers, we focused on an amplicon length of c. 100–200 base pairs and a primer length of 18–26 base pairs. The primers were designed using the NCBI Primer-BLAST software (Ye et al. 2012) and checked for melting point (optimum: 60–61 °C) and self-complementarity (< 5) with OligoAnalyzer 3.1 (Integrated DNA Technologies, Coralville, USA). Primers were obtained from Microsynth (Balgach, Switzerland), diluted to a concentration of 5 μM and first tested in a normal PCR. The PCR products were run on a 1% agarose gel stained with Midori Green at 80 V for 20 min in 1× TAE buffer. Since not all of the primers amplified sufficiently, we chose those producing the best results and showing specific amplification (a single amplicon) for the final qPCR experiments. In accordance with the MIQE guidelines (Bustin et al. 2009), all final primers were tested for their efficiency using LinReg 11.0. This software performs a linear regression analysis with the raw data for all replicate reactions of a primer (including the amplification data from all 40 cycles of a qPCR run) and calculates the primer efficiency.

RT-qPCR procedure

The qPCR was performed using 96-well optical PCR plates and seals (LabConsulting) and KAPA SYBR FAST qPCR Kit Universal (KAPA Biosystems). Each well contained a total reaction volume of 10 μl, consisting of 5 μl 2× KAPA SYBR FAST qPCR Master Mix Universal, 0.2 μl 50× ROX Low, 2.8 μl nuclease-free water, 250 nM of each forward and reverse primer, and 10 ng cDNA (1 μl). The qPCR was run on a 7500 Real-Time PCR System (Applied Biosystems). Cycling conditions were started with 3 min at 95 °C in order to activate the hot start polymerase, followed by 40 amplification cycles consisting of 15 s denaturation at 95 °C and 1 min annealing/extension at 60 °C.

The entire experiment was run once, and at the end material was harvested for RNA extractions. For each sample included in the qPCRs, we made a technical duplicate, which was preferably run on the same plate. These technical duplicates used the same cDNA and were performed to account for pipetting error in the qPCR. We also ran at least two non-template controls (NTC) per locus on each plate to detect potential contamination (NTCs with a cycle threshold (Ct)-value < 34). Technical duplicates varying by more than 1 cycle in their Ct-values were repeated, except for those with a Ct-value > 30, for which a difference of more than 1 cycle is not unusual due to the low RNA concentration.

Processing of qPCR data

Ct-values resulting from qPCR were standardized by the reference genes, and the resulting values (ΔCt) were used for data analysis. The cycle threshold Ct is defined as the number of PCR-cycles necessary for the fluorescent signal of a sample to exceed a predefined threshold (0.2), which allows a relative comparison of the original amount of cDNA copies of a gene. The earlier in a qPCR reaction the threshold cycle is reached, the higher the initial mRNA quantity. In order to minimize variation, we created the geometric mean of the Ct-values of each technical duplicate and used it for further calculations (to simplify, from here on referred to as the Ct-value). Then, for each candidate gene in each individual, a ΔCt-value was calculated according to the MIQE guidelines (Bustin et al. 2009). We subtracted the geometric mean of the reference genes from the Ct-value of the candidate gene: ΔCt = Ct_{candidate gene} − geomean (Ct_{reference gene 1}, Ct_{reference gene 2}).

In order to illustrate differential gene expression, we then used the ΔCt to create the relative expression (relative quantity = RQ) of each candidate gene. Here, we used the individual with the lowest expression as reference sample and calculated a ΔΔCt, from which the relative expression was calculated as follows: ΔΔCt = ΔCt − ΔCt_{reference sample}; RQ = 2^-ΔΔCt (Pfaffl 2001). Using RQ, we created charts to allow a visual inspection of gene expression as a function of temperature, site and presence of Plectocarpon.

Generation and processing of microsatellite data

To investigate the genetic component of gene expression, microsatellite data were generated so that genetic relationships among individuals could be inferred. For each individual used in the experiment, microsatellite data for the loci MS4, LPu37451, LPu28, LPu25, LPu09, LPu23, LPu17457, LPu39912, LPu13707, LPu15 and LPu04843 (Walser et al. 2003; Widmer et al. 2010; Werth et al. 2013) were generated by Cecilia Ronnås at the University of Graz and genotyped by SW using the microsatellite plugin implemented in Geneious v. 6.1.6. Individual genetic distance calculation followed the methods of Kosman & Leonard (2005) and the BIONJ algorithm, an improved version of the neighbour-joining algorithm, was used to generate an unrooted tree (Gascuel 1997); these algorithms were implemented in the R packages PopgenReport (Adamack & Gruber 2014) and ape (Paradis et al. 2004; Paradis 2006), and analyses were run in R v. 4.0.2 (R Core Team 2018).

Data analysis

For each putative reference gene, stability of expression was assessed over all studied samples and experimental conditions using boxplots. Additionally, NormFinder v. 0953 (Andersen et al. 2004) was used to quantify the stability of expression for the reference genes. The NormFinder program identifies genes with optimal normalization among a set of candidate genes. The lowest stability value indicates the most stable expression within the gene set examined, having the least variation within and among groups (Andersen et al. 2004).

Statistical analysis was performed in R v. 3.2.2 (R Core Team 2018). We tested for statistically significant differences in temperature and site using a multifactorial ANOVA of linear mixed effect models, with temperature and site as fixed factors and individual as random factor. If statistical significance was found, Tukey's post-hoc tests were used to calculate the P-values for comparisons between the three temperatures and/or between sites. In order to eliminate unintended factors, only individuals without P. lichenum from the ST7 population were used for comparisons between sites. To examine the difference between individuals with and without P. lichenum within the ST7 population, Student's t-tests were applied. We partitioned the variance in gene expression onto temperature, site, genetic factors and Plectocarpon infection in a partial redundancy analysis framework. First, a principal component analysis was performed on the microsatellite data to reduce their dimensionality. To do so, microsatellite alleles were coded as binary variables for each studied sample and a principal component analysis (PCA) was performed with the ‘princomp’ function in R v. 4.2.0. A total of 10 PCA axes were retained, explaining 80% of the variation in the microsatellite data, and these were included in (partial) redundancy analyses which were implemented in the package vegan (Oksanen et al. 2016). The aim of the redundancy analysis was to determine how much of the variance in gene expression of Lobaria pulmonaria was explained by genetic versus other factors (temperature, site, Plectocarpon lichenum infection).

Results

Verification of gene identities and expression stability

As expected, the 454 DNA sequences of Lobaria pulmonaria used to design primers matched with parts of the Lobaria pulmonaria Scotland v. 1.0 reference genome with identities of 99–100% (Table 1). Our gene names matched the KOG descriptions in the annotations of the L. pulmonaria genome for bet, efa, and gpd. Moreover, as expected, hsp88 and hsp98 were chaperones according to the Lobaria pulmonaria Scotland genome annotation. The PKS genes were annotated as ‘fatty acid synthase and related’ proteins in the Lobaria pulmonaria Scotland v.1.0 reference genome.

The efficiency of all primer pairs was ≥ 88% (Table 2). The stability values of bet and gpd were assessed with NormFinder software and found to be 0.014 and 0.015; hence these genes were stable in expression.

Table 2. Reference and candidate genes for Lobaria pulmonaria, used to study gene expression responses to increased temperature. Gene names are presented, with forward and reverse primer sequences and primer efficiency (Eff.).

Effects of Plectocarpon lichenum infection

Comparing the gene expression patterns of individuals with and without P. lichenum from site ST7, a significant difference was found in only one gene. While showing no difference in the 4 °C (Student's t-test: P = 0.7605, see Table 3) and 15 °C temperature treatments (Student's t-test: P = 0.4527, see Table 3), expression of the heat shock protein gene hsp98 was significantly higher at 25 °C for individuals infected with P. lichenum (Student's t-test: P = 0.0102). None of the other genes were differentially expressed between individuals with or without P. lichenum (Student's t-test: P > 0.1; Fig. 1C).

Fig. 1. Relative expression of mycobiont genes in thalli of the epiphytic lichen Lobaria pulmonaria from sampling sites AU7 (Austria) and ST7 (Spain, Tenerife) at 4 °C, 15 °C and 25 °C. For ST7, thalli with (ST7_Plect.) and without stromata of the lichenicolous fungus Plectocarpon lichenum were compared. The thallus with the lowest expression was used as a reference sample and set to one. The loci nrPKS3 and nrPKS3' represent two exons of the same gene. The letters ‘a’ and ‘b’ indicate a significant expression difference between samples infected with P. lichenum and those not infected. In colour online.

Table 3. P-values of Student's t-tests for the differences in gene expression between individuals of the ST7 (a site in Tenerife, Spain) population of Lobaria pulmonaria with and without Plectocarpon lichenum infection at 4 °C, 15 °C and 25 °C. Statistically significant values are given in bold.

Effects of temperature and collecting site

In all genes tested, a significant difference in gene expression due to increased temperatures was observed (ANOVA: P < 0.009; see Table 4, Fig. 1). There was a positive correlation of temperature and gene expression, except for efa in the AU7 population (Fig. 1A).

Table 4. P-values of ANOVA, using a linear mixed effects model with temperature and habitat as fixed factors, and site and lichen individual (Lobaria pulmonaria) as random factors, for differences in the expression of the heat shock protein genes (hsp88 and hsp98), elongation factor 1-α (efa) and the polyketide synthase genes (rPKS1, nrPKS3 and nrPKS3'). Statistically significant values are given in bold.

Since in all genes significant differences in gene expression due to increased temperature were found, Tukey's honest significance test was performed to find out at which temperatures exactly differential expression took place. There was a significant difference in gene expression of both heat shock protein genes hsp88 and hsp98 (Fig. 1B & C) with every temperature increase (Tukey's test: P < 0.002), being highly significant (Tukey's test: P < 0.0001) between the 4 °C and 25 °C temperature treatments (Table 5).

Table 5. P-values of Tukey's honest significance test for differences in the expression of the heat shock protein genes (hsp88 and hsp98), the elongation factor 1-α (efa) and the polyketide synthase genes (rPKS1, nrPKS3 and nrPKS3') of Lobaria pulmonaria, due to temperature treatments at 4 °C, 15 °C and 25 °C. Statistically significant values are given in bold.

The polyketide synthase genes rPKS1, nrPKS3 and nrPKS3' (Fig. 1D–F) were upregulated at the temperature increase from 4 °C to 15 °C, as well as at 4 °C to 25 °C (Tukey's test: P < 0.008), but did not show a significant difference at 25 °C compared to 15 °C (Tukey's test: P > 0.05; Table 5). In efa, significant upregulation was found only at 25 °C compared to 4 °C (P < 0.03; see Table 5). For efa, hsp98 and rPKS1, there was differential expression not only between temperatures but also between sites (ANOVA: P < 0.02; Table 4). For nrPKS3', a significant interaction between temperature and site was observed (ANOVA: P = 0.0115; Table 4). In AU7, an upregulation of nrPKS3' took place at 15 °C compared to 4 °C (Tukey's test: P = 0.0050) and at 25 °C compared to 4 °C (Tukey's test: P = 0.0007; Fig. 1, Table 6). In ST7, however, there was already a high expression of nrPKS3' at 4 °C, which did not increase in the 15 °C temperature treatment (Tukey's test: P = 1); while there was an upregulation at 25 °C, this was only near significant in Tukey's test (P < 0.1; Table 6).

Table 6. P-values of Tukey's honest significance test for differences in the expression of the polyketide synthase gene nrPKS3', due to the temperature treatments at 4 °C, 15 °C and 25 °C in Lobaria pulmonaria individuals from the sites AU7 (Austria) and ST7 (Spain, Tenerife) and both sites combined. Statistically significant values are given in bold.

Genetic distance among samples of Lobaria pulmonaria

Analysis of microsatellites indicated that both the Austrian and the Spanish population of L. pulmonaria were genetically diverse, with Austrian samples clustering together in the unrooted BIONJ tree (Fig. 2).

Fig. 2. Unrooted BIONJ neighbour-joining tree (see Gascuel 1997) for 11 microsatellite loci of the 25 Lobaria pulmonaria samples from Austria (AU7) and Spain (ST7) included in the gene expression experiment. Branches containing Austrian samples are shown in grey.

Partitioning of variance in gene expression data

Using redundancy analysis, 59.7% of the variance in gene expression was explained by regulatory (temperature), acclimatory (site), genetic and biotic (Plectocarpon-infection) effects. A total of 40.3% of the total variance was unexplained. Regulatory effects were the most important, with variation in gene expression due to temperature increase accounting for 81.4% of the explained variance (site = 2.9%, Plectocarpon-infection = 0.5%; Fig. 3). A total of 11.8% of the explained variance was attributed to genetic factors. Covariance among variable sets amounted to 3.4% of the explained variance. In other words, temperature treatment explained seven times more variance than genetic distance, 28 times more variance than acclimation to collecting site, and 156 times more variance than Plectocarpon-infection.

Fig. 3. Partitioning variance in gene expression of the lichen-forming fungus Lobaria pulmonaria onto four variable sets in partial redundancy analyses: regulatory (three different temperatures), genetic (10 principal components based on 11 microsatellite loci), acclimatory (site of origin, Austria or Spain) and biotic components (presence or absence of Plectocarpon lichenum infection). Covariance refers to variance shared among the variable sets. The percentage of explained variance is illustrated.

Discussion

Expression stability of reference genes

Our study provides two new reference genes for qPCR studies of Lobaria pulmonaria. The genes bet and gpd were stable in their expression and did not vary with temperature, therefore fulfilling the criteria for use as reference genes (Bustin et al. 2009).

Effects of Plectocarpon lichenum infection

The overall effect of Plectocarpon lichenum infection on variance in gene expression was low. However, the heat shock protein gene hsp98 showed significant infection-related differential expression in L. pulmonaria. Pathogen attack is known to induce upregulation of heat shock responses in plants (Aranda et al. 1996; Havelda & Maule 2000; Chivasa et al. 2005; Andrási et al. 2021). There is a lack of knowledge of how fungi, including lichenized species, react to pathogen attack but they seem to possess the genetic mechanisms required to detect and respond to pathogens (Uehling et al. 2017).

Effects of temperature and collecting site

The main hypothesis in our study was confirmed, that thermal stress influences the expression of candidate genes for stress response. Playing an important role in refolding of denatured proteins (Miot et al. 2011; Li & Buchner 2013), most heat shock protein genes are upregulated at least in the first response to thermal stress (Plesofsky-Vig & Brambl 1998; Che et al. 2013; Park et al. 2015; Steinhäuser et al. 2016). The heat shock protein genes of the lichen-forming fungus Lobaria pulmonaria were indeed significantly upregulated after the temperature increases: a heat shock response took place. Simultaneously with the heat shock response, the PKS genes showed a significant upregulation with every temperature increase. Since stress-induced polyketide production has been observed in bacteria (Auckloo et al. 2017) and in lichen-forming fungi (Armaleo et al. 2008; Timsina et al. 2013), an upregulation of PKS genes was anticipated. Little is known about the conditions under which fungal PKS genes are upregulated or by which biosynthetic genes fungal metabolites are produced (Kim et al. 2021), but the importance of these compounds for lichen tolerance of stressful biotic or abiotic conditions has previously been emphasized (Huneck 1999).

Interestingly, elongation factor 1-α (efa) showed upregulation with each temperature increase in L. pulmonaria. This gene is involved in protein biosynthesis and specifically in chain elongation by recruiting t-RNAs to ribosomes (Anand et al. 2003). While this gene has been used as a reference gene for qPCR because of its stable expression, for example in potato (Nicot et al. 2005) and cod (Aursnes et al. 2011), there is evidence that it is heat-induced in plants (Nikolaou et al. 2009; Momčilović et al. 2016; Sun et al. 2020), where it may also function as a molecular chaperone involved in protein degradation (Talapatra et al. 2002; Shin et al. 2009). Under higher temperatures, this gene may therefore be upregulated in lichenized fungi, presumably to also function as a molecular chaperone.

We found a heat shock response in L. pulmonaria even at moderate temperatures (i.e. 15 °C and 25 °C); there was an upregulation of both hsp88 and hsp98 with every temperature increase. In its natural growth habitat, L. pulmonaria is wet and physiologically active mostly at temperatures up to 15 °C (Pannewitz et al. 2003). Apparently, even moderate temperatures can provoke heat shock reactions in cold-adapted L. pulmonaria, although the effect was much less pronounced at 15 °C than at 25 °C. Others have found a temperature of 25 °C to be sufficient to induce severe stress conditions in Peltigera scabrosa (MacFarlane & Kershaw 1980). The fungal gene hsp88, encoding a heat shock protein similar to the hsp110 family (Plesofsky-Vig & Brambl 1998), was strongly induced at 25 °C in AU7. Although the expression was distinctly higher and there was no overlap among standard errors, the difference between the sites was not statistically significant. This might be caused by the high variance due to the small sample size of AU7. The gene hsp98, which encodes a prominent heat shock protein (Vassilev et al. 1992), showed less upregulation, although there was a significant difference between sites, mainly with the 15 °C treatment in AU7 showing higher gene expression. This might indicate that individuals from Austria are more sensitive to heat stress than those from Tenerife.

Response to high temperature may potentially affect many physiological processes, including growth and resistance to pathogens. For example, in plants, increased temperatures lead to suppressed immunity to pathogens, since higher temperatures can shift the allocation of heat shock proteins from defense responses to heat stress responses (Lee et al. 2012; Dangi et al. 2018; Janda et al. 2019). It is conceivable that heat-stressed lichens possess a lower ability to defend themselves against pathogens for the same reason. A temperature-dependent reduced defense could potentially modify interactions with lichenicolous fungi, making them increasingly more antagonistic. Furthermore, defense mechanisms against herbivores could also become weakened, which could lead to decreased survival rates.

Timsina et al. (2013) reported an increase of lichen polyketide content in Ramalina dilacerata under stressful conditions and, in general, polyketide content of lichens is thought to confer increased tolerance to biotic and abiotic stressors (Huneck 1999). In the PKS genes included in this study, expression increased significantly with the temperature rise from 4 °C to 15 °C, as well as highly significantly from 4 °C to 25 °C. While these results are promising, more work is needed to characterize the functions of PKS genes in lichens and the pathways producing specific secondary compounds (Kim et al. 2021).

Our data exhibited a small effect of collecting site, which represents the remaining effect of physiological long-term acclimatization to sites after laboratory acclimation. This finding is consistent with the results of Steinhäuser et al. (2016), who also found collecting site-related differential expression in Peltigera membranacea after three weeks of acclimation to cold in the laboratory. Collecting site-related different physiological responses to heat stress were also found in Peltigera canina (MacFarlane & Kershaw 1980). Our two collecting sites are situated in different climatic zones where the local environmental conditions should be rather different (Pannewitz et al. 2003).

We found a significantly stronger induction of rPKS1 in individuals from Austria compared to those from Tenerife which, together with the stronger induced heat shock protein gene expression in Austria, indicates that the gene response can vary in magnitude between populations. Profound gene expression differences between populations were also reported for Peltigera membranacea exposed to increases in temperature (Steinhäuser et al. 2016). In our study of L. pulmonaria, the residual acclimatory effects were nevertheless small, representing only 2.9% of the explained variance. This is not surprising as the thalli were acclimated to cold for three weeks, and lichens can acclimate their photosynthesis to changed conditions within a few days (Kershaw 1977; MacKenzie et al. 2004).

As expected, the variance in gene expression of L. pulmonaria in response to thermal stress appeared to be mainly due to the manipulated variable in our laboratory experiment, temperature; thus, the response reflects mostly an adjustment to thermal stress to maintain cellular functions. That this regulatory component of variation dominates in gene expression variation is perhaps not overly surprising in a mutualistic lichen symbiosis, where a fine-tuned physiological equilibrium between mycobiont and photobiont must be maintained to ensure the long-term persistence of the association. Our finding that genetic differences represent, with a total of 11.8% of the explained variance, the second largest component of gene expression variation in response to thermal stress in L. pulmonaria, and that acclimation explained only 2.8% is remarkable because it implies that the three week acclimation treatment to 4 °C removed most differences in gene expression due to long-term physiological acclimatization to the sites of origin in Austria and Tenerife, if any larger acclimatory differences existed in the first place. In our study, we did not quantify the maximum (initial) acclimation effect, since our first sample was taken after several weeks of acclimation to cold conditions in the laboratory. Other studies have found seasonal light acclimation of photosynthesis in L. pulmonaria (Schofield et al. 2003) which occurs via macromolecular allocation to chlorophyll and RuBisCo protein (MacKenzie et al. 2004). Such acclimation to changes in ambient light and temperature can occur immediately in lichens, over as little time as two days (Kershaw 1977, 1985; MacKenzie et al. 2004). Within the three week laboratory acclimation period, the samples should therefore have become completely acclimated to cold.

As much as 40.3% of the total variance in gene expression data was not explained by the factors covered in our study. This finding is not surprising, given that gene expression data tend to have a large stochastic component, even for populations of clonal cells under standardized conditions (McAdams & Arkin 1997; Elowitz et al. 2002; Blake et al. 2003; Kærn et al. 2005). Much greater variance would be expected in data gathered from natural populations where individuals may differ in genomic background, physiological acclimatization, phenotype, age, reproductive state, and other factors. Differences among individuals might contribute to some of the unexplained variation in gene expression. Substantial inter-individual variation in gene expression has also been reported for another Peltigeralean lichen, Peltigera membranacea (Steinhäuser et al. 2016).

Conclusions

The lichen-forming fungus Lobaria pulmonaria may provide an interesting model for in vivo studies of heat shock responses. Overall, our results show clearly that gene expression variation in L. pulmonaria under thermal stress is substantially influenced by the abiotic environment (temperature), with regulatory effects predominating (i.e. direct responses to elevated temperature). Lichen-forming fungi have evolved powerful molecular pathways to withstand environmental fluctuations and stress, and heat shock responses are a critical component conveying stress tolerance. Our results suggest that the colonization of thalli by lichenicolous fungi might have an influence on the mycobiont's heat shock responses; abiotic and biotic factors appear to cause cumulative effects. While L. pulmonaria has the molecular machinery to counteract short-term thermal stress, its persistence in a given landscape depends on the overall long-term positive carbon balance, which can be compromised by warmer temperatures leading to increased respiration rates and by reduced precipitation during summer, both predicted for Central Europe in connection with global climate change (Middelkoop et al. 2001; Ahrens et al. 2014; IPCC 2021). These topics deserve more attention in future work.