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. Reference MacKenzie, MacDonald, Dubois and Campbell2001, Reference MacKenzie, Johnson and Campbell2004; Evans et al. Reference Evans, Chan, Menge and Hofmann2013; Hamann et al. Reference Hamann, Kesselring, Armbruster, Scheepens and Stöcklin2016). In order to deal with those changing conditions, the ability to regulate the expression of stress-related genes is vital (Evans et al. Reference Evans, Chan, Menge and Hofmann2013). 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. Reference Wang, Zhang, Zhou, Zhang and Wei2015; Carmo et al. Reference Carmo, Martins, Martins, Passos, Silva, Araujo, Brasileiro, Miller, Guimarães and Mehta2019), temperature and salinity stress (Jamil et al. Reference Jamil, Riaz, Ashraf and Foolad2011; Che et al. Reference Che, Song and Lin2013; Zhang et al. Reference Zhang, Zhang and Fan2017), as well as exposure to toxins (Whitehead et al. Reference Whitehead, Triant, Champlin and Nacci2010). 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. Reference Mizoguchi, Ichimura and Shinozaki1997; Gasch Reference Gasch2007; Dixon et al. Reference Dixon, Abbott and Matz2020; Terhorst et al. Reference Terhorst, Sandikci, Keller, Whittaker, Dunham and Amon2020). In fungi, environmental stress response was first described in Saccharomyces cerevisiae (Gasch et al. Reference Gasch, Spellman, Kao, Carmel-Harel, Eisen, Storz, Botstein and Brown2000). 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 Reference Gasch2007; O'Meara et al. Reference O'Meara, O'Meara, Polvi, Pourhaghighi, Liston, Lin, Veri, Emili, Gingras and Cowen2019).
One common environmental stressor, which organisms are confronted with, is thermal stress (Arshad et al. Reference Arshad, Farooq, Asch, Krishna, Prasad and Siddique2017). 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 Reference Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger, Caud, Chen, Goldfarb and Gomis2021). 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 CO2-emissions (IPCC Reference Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger, Caud, Chen, Goldfarb and Gomis2021). Heat shock response represents one of the important mechanisms for organisms to adapt to stressful conditions at the cellular level (O'Meara et al. Reference O'Meara, O'Meara, Polvi, Pourhaghighi, Liston, Lin, Veri, Emili, Gingras and Cowen2019).
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. Reference Franzmann, Menhorn, Walter and Buchner2008); 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. Reference Hagiwara, Sakamoto, Abe and Gomi2016). 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. Reference Albrecht, Guthke, Brakhage and Kniemeyer2010; Roncarati & Scarlato Reference Roncarati and Scarlato2017). 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. Reference Albrecht, Guthke, Brakhage and Kniemeyer2010; de Nadal et al. Reference de Nadal, Ammerer and Posas2011; Takahashi et al. Reference Takahashi, Kusuya, Hagiwara, Takahashi-Nakaguchi, Sakai and Gonoi2017).
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. Reference Albrecht, Guthke, Brakhage and Kniemeyer2010). 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 Reference Plesofsky-Vig and Brambl1998; Miot et al. Reference Miot, Reidy, Doyle, Hoskins, Johnston, Genest, Vitery, Masison and Wickner2011; Smith et al. Reference Smith, Burns, Shearer and Snell2012; Li & Buchner Reference Li and Buchner2013; Park et al. Reference Park, Lee, Kang, Kim and Kwak2015). The heat shock protein gene hsp88 of an entomopathogenic fungus has been cloned and characterized by Park et al. (Reference Park, Kim, Kim, Park, Son, Hong and Lee2014). Under thermal stress, hsp88 was 15–55-fold upregulated in the lichen-forming fungus Peltigera membranacea (Steinhäuser et al. Reference Steinhäuser, Andrésson, Pálsson and Werth2016). An important heat shock protein gene in Aspergillus fumigatus is hsp98 (Do et al. Reference Do, Yamaguchi and Miyano2009), and this gene was also upregulated under thermal stress in P. membranacea (Steinhäuser et al. Reference Steinhäuser, Andrésson, Pálsson and Werth2016).
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. Reference Timsina, Sorensen, Weihrauch and Piercey-Normore2013). 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. Reference Nivina, Yuet, Hsu and Khosla2019), 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 Reference Lawrey1986, Reference Lawrey1989; Torzilli et al. Reference Torzilli, Mikelson and Lawrey1999; Gauslaa & McEvoy Reference Gauslaa and McEvoy2005; Timsina et al. Reference Timsina, Sorensen, Weihrauch and Piercey-Normore2013). 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. Reference Khosla, Gokhale, Jacobsen and Cane1999). 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. Reference Nivina, Yuet, Hsu and Khosla2019). 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 Reference Bertrand and Sorensen2018). 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. Reference Armaleo, Zhang and Cheung2008; Timsina et al. Reference Timsina, Sorensen, Weihrauch and Piercey-Normore2013). 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. Reference Bidussi, Goward and Gauslaa2013; Gauslaa et al. Reference Gauslaa, Bidussi, Solhaug, Asplund and Larsson2013). The depsidones, norstictic and stictic acid, are produced via the acetate-polymalonate pathway (Ranković & Kosanić Reference Ranković, Kosanić and Ranković2019). Some lichen secondary compounds, including those of L. pulmonaria, have anti-herbivore and antibiotic properties (Suleyman et al. Reference Suleyman, Odabasoglu, Aslan, Cakir, Karagoz, Gocer, Halıcı and Bayir2003; Asplund & Gauslaa Reference Asplund and Gauslaa2008; Nybakken et al. Reference Nybakken, Helmersen, Gauslaa and Selas2010). Some secondary compounds such as lecanoric acid may also have antifungal properties, preventing lichen colonization by certain lichenicolous fungi (Lawrey Reference Lawrey1989, Reference Lawrey2000; Lawrey & Diederich Reference Lawrey and Diederich2003), and some may be useful as anti-cancer drugs (Shrestha & St. Clair Reference Shrestha and St. Clair2013; Dar et al. Reference Dar, Dar, Islam, Mangral, Dar, Singh, Verma and Haque2021; Yang et al. Reference Yang, Devkota, Wang and Scheidegger2021).
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 Reference Bergmann and Werth2017). 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 Reference Bergmann and Werth2017). 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 Reference Bergmann and Werth2017). 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. Reference Asplund, Gauslaa and Merinero2016). 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. Reference Merinero, Bidussi and Gauslaa2015). 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 Reference Lawrey2000). 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 Reference MacFarlane and Kershaw1980; Cheviron et al. Reference Cheviron, Whitehead and Brumfield2008; Whitehead et al. Reference Whitehead, Triant, Champlin and Nacci2010; Steinhäuser et al. Reference Steinhäuser, Andrésson, Pálsson and Werth2016). Habitat-related differential gene expression could be composed of both genetic and acclimatory factors (Cheviron et al. Reference Cheviron, Whitehead and Brumfield2008; Whitehead et al. Reference Whitehead, Triant, Champlin and Nacci2010; Palumbi et al. Reference Palumbi, Barshis, Traylor-Knowles and Bay2014). 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. Reference Steinhäuser, Andrésson, Pálsson and Werth2016) or physiological state (MacFarlane & Kershaw Reference MacFarlane and Kershaw1980; Schipperges et al. Reference Schipperges, Kappen and Sonesson1995). 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. Reference Asplund, Gauslaa and Merinero2016, Reference Asplund, Gauslaa and Merinero2018), 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. Reference Scheidegger, Bilovitz, Werth, Widmer and Mayrhofer2012). 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 dH2O, 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. Reference De Wit, Pespeni, Ladner, Barshis, Seneca, Jaris, Therkildsen, Morikawa and Palumbi2012). 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. (Reference Steinhäuser, Andrésson, Pálsson and Werth2016).
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. Reference Doering, Miao and Piercey-Normore2014), 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 H2O 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. Reference Joneson, Armaleo and Lutzoni2011; Miao et al. Reference Miao, Manoharan, Snæbjarnarson and Andrésson2012).
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. Reference Altschul, Madden, Schäffer, Zhang, Zhang, Miller and Lipman1997). 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).
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. (Reference Werth, Cornejo and Scheidegger2013); the multispore mycobiont culture F2, which was used to obtain the data, is described in Widmer et al. (Reference Widmer, Dal Grande, Cornejo and Scheidegger2010) and Cornejo et al. (Reference Cornejo, Scheidegger and Honegger2015). 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. Reference Shin, Moon, Park, Kim and Byun2009). 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−40 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. Reference Ye, Coulouris, Zaretskaya, Cutcutache, Rozen and Madden2012) 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. Reference Bustin, Benes, Garson, Hellemans, Huggett, Kubista, Mueller, Nolan, Pfaffl and Shipley2009), 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.
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. Reference Bustin, Benes, Garson, Hellemans, Huggett, Kubista, Mueller, Nolan, Pfaffl and Shipley2009). We subtracted the geometric mean of the reference genes from the Ct-value of the candidate gene: ΔCt = Ctcandidate gene − geomean (Ctreference gene 1, Ctreference 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 − ΔCtreference sample; RQ = 2-ΔΔCt (Pfaffl Reference Pfaffl2001). 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. Reference Walser, Sperisen, Soliva and Scheidegger2003; Widmer et al. Reference Widmer, Dal Grande, Cornejo and Scheidegger2010; Werth et al. Reference Werth, Cornejo and Scheidegger2013) 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 (Reference Kosman and Leonard2005) and the BIONJ algorithm, an improved version of the neighbour-joining algorithm, was used to generate an unrooted tree (Gascuel Reference Gascuel1997); these algorithms were implemented in the R packages PopgenReport (Adamack & Gruber Reference Adamack and Gruber2014) and ape (Paradis et al. Reference Paradis, Claude and Strimmer2004; Paradis Reference Paradis2006), and analyses were run in R v. 4.0.2 (R Core Team 2018).
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. Reference Andersen, Ledet-Jensen and Ørntoft2004) 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. Reference Andersen, Ledet-Jensen and Ørntoft2004).
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. Reference Oksanen, Blanchet, Kindt, Legendre, Minchin, O'Hara, Simpson, Solymos, Stevens and Wagner2016). 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).
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.
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).
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).
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).
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).
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).
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.
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. Reference Bustin, Benes, Garson, Hellemans, Huggett, Kubista, Mueller, Nolan, Pfaffl and Shipley2009).
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. Reference Aranda, Escaler, Wang and Maule1996; Havelda & Maule Reference Havelda and Maule2000; Chivasa et al. Reference Chivasa, Simon, Yu, Yalpani and Slabas2005; Andrási et al. Reference Andrási, Pettkó-Szandtner and Szabados2021). 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. Reference Uehling, Deveau and Paoletti2017).
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. Reference Miot, Reidy, Doyle, Hoskins, Johnston, Genest, Vitery, Masison and Wickner2011; Li & Buchner Reference Li and Buchner2013), most heat shock protein genes are upregulated at least in the first response to thermal stress (Plesofsky-Vig & Brambl Reference Plesofsky-Vig and Brambl1998; Che et al. Reference Che, Song and Lin2013; Park et al. Reference Park, Lee, Kang, Kim and Kwak2015; Steinhäuser et al. Reference Steinhäuser, Andrésson, Pálsson and Werth2016). 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. Reference Auckloo, Pan, Akhter, Wu, Wu and He2017) and in lichen-forming fungi (Armaleo et al. Reference Armaleo, Zhang and Cheung2008; Timsina et al. Reference Timsina, Sorensen, Weihrauch and Piercey-Normore2013), 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. Reference Kim, Jeong, Yun and Hur2021), but the importance of these compounds for lichen tolerance of stressful biotic or abiotic conditions has previously been emphasized (Huneck Reference Huneck1999).
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. Reference Anand, Chakraburtty, Marton, Hinnebusch and Kinzy2003). While this gene has been used as a reference gene for qPCR because of its stable expression, for example in potato (Nicot et al. Reference Nicot, Hausman, Hoffmann and Evers2005) and cod (Aursnes et al. Reference Aursnes, Rishovd, Karlsen and Gjøen2011), there is evidence that it is heat-induced in plants (Nikolaou et al. Reference Nikolaou, Agrafioti, Stumpf, Quinn, Stansfield and Brown2009; Momčilović et al. Reference Momčilović, Pantelić, Zdravković-Korać, Oljača, Rudić and Fu2016; Sun et al. Reference Sun, Ji, Jia, Huo, Si, Zeng, Zhang and Niu2020), where it may also function as a molecular chaperone involved in protein degradation (Talapatra et al. Reference Talapatra, Wagner and Thompson2002; Shin et al. Reference Shin, Moon, Park, Kim and Byun2009). 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. Reference Pannewitz, Schroeter, Scheidegger and Kappen2003). 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 Reference MacFarlane and Kershaw1980). The fungal gene hsp88, encoding a heat shock protein similar to the hsp110 family (Plesofsky-Vig & Brambl Reference Plesofsky-Vig and Brambl1998), 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. Reference Vassilev, Plesofsky-Vig and Brambl1992), 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. Reference Lee, Yun and Kwon2012; Dangi et al. Reference Dangi, Sharma, Khangwal and Shukla2018; Janda et al. Reference Janda, Lamparová, Zubíková, Burketová, Martinec and Krčková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. (Reference Timsina, Sorensen, Weihrauch and Piercey-Normore2013) 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 Reference Huneck1999). 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. Reference Kim, Jeong, Yun and Hur2021).
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. (Reference Steinhäuser, Andrésson, Pálsson and Werth2016), 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 Reference MacFarlane and Kershaw1980). Our two collecting sites are situated in different climatic zones where the local environmental conditions should be rather different (Pannewitz et al. Reference Pannewitz, Schroeter, Scheidegger and Kappen2003).
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. Reference Steinhäuser, Andrésson, Pálsson and Werth2016). 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 Reference Kershaw1977; MacKenzie et al. Reference MacKenzie, Johnson and Campbell2004).
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. Reference Schofield, Campbell, Funk and MacKenzie2003) which occurs via macromolecular allocation to chlorophyll and RuBisCo protein (MacKenzie et al. Reference MacKenzie, Johnson and Campbell2004). Such acclimation to changes in ambient light and temperature can occur immediately in lichens, over as little time as two days (Kershaw Reference Kershaw1977, Reference Kershaw1985; MacKenzie et al. Reference MacKenzie, Johnson and Campbell2004). 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 Reference McAdams and Arkin1997; Elowitz et al. Reference Elowitz, Levine, Siggia and Swain2002; Blake et al. Reference Blake, Kærn, Cantor and Collins2003; Kærn et al. Reference Kærn, Elston, Blake and Collins2005). 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. Reference Steinhäuser, Andrésson, Pálsson and Werth2016).
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. Reference Middelkoop, Daamen, Gellens, Grabs, Kwadijk, Lang, Parmet, Schadler, Schulla and Wilke2001; Ahrens et al. Reference Ahrens, Formayer, Gobiet, Heinrich, Hofstätter, Matulla, Prein, Truhetz, Anders, Haslinger, Kromp-Kolb, Nakicenovic, Steininger, Gobiet, Formayer, Köppl, Prettenthaler, Stötter and Schneider2014; IPCC Reference Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger, Caud, Chen, Goldfarb and Gomis2021). These topics deserve more attention in future work.
Sophie S. Steinhäuser provided advice on statistical analyses and Cecilia Ronnås generated the microsatellite data. MK thanks family and friends for their encouragement and support. This work was supported by Helmut Mayrhofer and the Institute of Plant Sciences, University of Graz, Austria, and by the Swiss National Science Foundation [grants 31003A-105830 and 31003A-127346 to CS]. We thank anonymous reviewers for very constructive comments on earlier drafts of this paper.
Silke Werth, 0000-0002-4981-7850.