Circadian clock‐ and temperature‐associated genes contribute to overall genomic differentiation along elevation in lichenized fungi

Circadian regulation is linked to local environmental adaptation, and many species with broad climatic niches display variation in circadian genes. Here, we hypothesize that lichenizing fungi occupying different climate zones tune their metabolism to local environmental conditions with the help of their circadian systems. We study two species of the genus Umbilicaria occupying similar climatic niches (Mediterranean and the cold temperate) in different continents. Using homology to Neurospora crassa genes, we identify gene sets associated with circadian rhythms (11 core, 39 peripheral genes) as well as temperature response (37 genes). Nucleotide diversity of these genes is significantly correlated with mean annual temperature, minimum temperature of the coldest month and mean temperature of the coldest quarter. Furthermore, we identify altitudinal clines in allele frequencies in several non‐synonymous substitutions in core clock components, for example, white collar‐like, frh‐like and various ccg‐like genes. A dN/dS approach revealed a few significant peripheral clock‐ and temperature‐associated genes (e.g. ras‐1‐like, gna‐1‐like) that may play a role in fine‐tuning the circadian clock and temperature‐response machinery. An analysis of allele frequency changes demonstrated the strongest evidence for differentiation above the genomic background in the clock‐associated genes in U. pustulata. These results highlight the likely relevance of the circadian clock in environmental adaptation, particularly frost tolerance, of lichens. Whether or not the fungal clock modulates the symbiotic interaction within the lichen consortium remains to be investigated. We corroborate the finding of genetic variation in clock components along altitude—not only latitude—as has been reported in other species.


| INTRODUC TI ON 1.| Variation in circadian rhythms along latitudinal and elevational gradients
The circadian clock is a central molecular timekeeping mechanism that regulates biochemical processes and helps organisms to perceive and respond to abiotic and biotic environmental cues.Given this fundamental role of the circadian clock in environmental adaptation, it is not surprising that biological rhythms vary along latitudinal as well as elevation gradients (Hut et al., 2013).Signals of adaptation along environmental gradients, in particular latitudinal gradients, have been observed in the circadian clock components of a variety of lineages, from insects (Helfrich-Förster et al., 2020;Lindestad et al., 2022) and other animals (Moreno et al., 2021) to plants (Yerushalmi et al., 2011) and fungi (Ellison et al., 2011).

| Role of the circadian clock in temperature-dependent adaptation
Individual circadian clock components are known to be directly temperature-responsive.The central circadian oscillator in fungi is composed of the negative element frequency (frq) and its interactions with the White Collar Complex (WCC), composed of proteins encoded by the white collar-1 (wc-1) and white collar-2 (wc-2) genes (Collett et al., 2002;Froehlich et al., 2002), and the Frequency Interacting RNA Helicase (frh) gene, which helps to modulate FRQ degradation rates (Hurley et al., 2013;Shi et al., 2010).Expression levels of long and short isoforms of FRQ are known to be regulated by the thermosensitive splicing of intron 6, and overall translation of frq RNA is thermoregulated by the presence of six upstream ORFs that trap scanning ribosomes, regulating the expression of the main functional ORF (Colot et al., 2005;Diernfellner, 2005).
Aside from the temperature regulation of specific clock components, the circadian clock machinery generally provides fitness advantages.Strains of Neurospora discreta with habitat-specific circadian rhythms maintain higher fitness under their respective habitats (Koritala et al., 2020).Circadian rhythms have also been observed to tightly track with latitude and elevation in the model plant Arabidopsis thaliana (Michael et al., 2003;Salmela & Weinig, 2019;Vidigal et al., 2016), as well as in natural populations of Mimulus guttatus and domesticated populations of Glycine max (Greenham et al., 2017).

| Variation in circadian clock and temperature-response components
Temperature compensation is a core feature of circadian clock systems (Edwards et al., 2005;Gil & Park, 2019), and in many circadian systems, including those of fungi, the temperature-and clock-responsive mechanisms share many components (Hunt et al., 2012).In plants, temperature compensation mechanisms are known to modulate the function of the circadian clock under both stressful and non-stressful temperature deviations from the optimum growing temperature (Gil & Park, 2019).Furthermore, natural variation in a circadian clockassociated locus, FLOWERING LOCUS C (FLC), has been shown to mediate high temperature-specific clock responses (Edwards et al., 2006).

| Variation along environmental gradients in the lichen symbiosis
Many species of lichenized fungi have broad geographical and climatic ranges, and exhibit significant population structure (e.g.Allen et al., 2021;Fernández-Mendoza et al., 2011;Muggia et al., 2014).Species of the genus Umbilicaria have been previously extensively studied for their population turnovers along elevation gradients, along which species such as U. pustulata and U. phaea grow in a continuum from Mediterranean to cold temperate climate zones (Rolshausen et al., 2022).Part of the observed population differentiation in these species can be linked to environmental factors, evidenced, for example, by abrupt, parallel genomic breaks along elevation gradients, separating populations into high-and low-elevation genetic clusters, which correspond to different climate zones (Dal Grande et al., 2017;Rolshausen et al., 2020Rolshausen et al., , 2022)).U. pustulata and U. phaea exhibit different reproductive strategies; while both taxa have similar growth forms, habitats and ecological ranges, U. phaea produces fruiting bodies and reproduces sexually via fertile ascospores (Hestmark, 2004), while U. pustulata typically does not exhibit sexual structures and reproduces vegetatively via isidia, thallus outgrowths acting as clonal dispersal units for fungal hyphae and algal cells simultaneously (Hestmark, 1992).
Despite this difference in reproductive strategies, adaptation to different niches along elevation gradients may be linked to similar processes in both species, as evidenced by the presence or absence of particular enzymes in either climate zone (Merges et al., 2023).Adaptive responses of the circadian system to climate gradients are presently unknown in lichens.
In the case of evolutionarily more recent divergence, loci driving local adaptation responses and loci in close linkage to those loci under adaptation can be detected against a largely undifferentiated genomic background (see e.g., Nosil & Feder, 2012;Ravinet et al., 2017).In order to gain a clearer picture of what adaptive mechanisms drive differentiation in the case of highly divergent populations, like those of many lichen species, it may be advantageous to target known sources of climate adaptation.In many species, including fungi (Ellison et al., 2011), temperature changes have been shown to drive adaptation in cold-response and circadian clock-associated loci, and these loci may, therefore, provide a sensible starting point for investigating climate adaptation in lichens.

| Hypotheses
Previous work has elucidated the existing variation within species of Umbilicaria and shown that different ecotypes are associated with distinct climate zones along elevation gradients (Dal Grande et al., 2017;Rolshausen et al., 2022).The putative core clock components have been identified in several lichen species, including the Umbilicaria species above, and the functional conservation of the FRQ-WCC core circadian clock loop has been implicated by experimental evidence (Valim et al., 2022).To examine whether the circadian clock-and light-and temperature-mediated responses more generally-are associated with the distinct ecotypes observed in U. pustulata, and U. phaea populations, we asked the following questions: Is there variation in the core circadian clock genes along elevation gradients?If so, is this variation linked to functional loci (i.e. SNPs that lead to radical non-synonymous substitutions)?

| ME THODS
The lichen species used in this study, Umbilicaria pustulata and U. phaea, were sampled along five different gradients in either Southern Europe or California.In Europe, these sites were Mount Limbara (Sardinia, Italy; 6 U. pustulata populations, IT), described in Dal Grande After sampling, genomic DNA was extracted using a CTAB-based method and individuals pooled at equimolar concentrations were sequenced with Illumina paired-end chemistry at ~90× coverage per population (100 bp for IT and ES2 populations, 150 bp for ES1 and CA populations) before trimming, mapping reads to reference genomes, SNP calling and analysis of circadian clock-and temperatureassociated genes.All analyses were conducted either in a Unix environment or in R and using RStudio (RStudio Team, 2020; R Core Team, 2018).Figures were created with ggplot2 (Wickham, 2016).
Scripts of all analyses are available on GitHub at https:// github.com/ hvalim/ lichen_ circa dian_ gradi ents, and detailed information about each analysis in this study is provided in the Data S1.

| Core circadian clock genes display non-synonymous substitutions along climatic gradients
We previously identified putative homologues of the core circadian clock genes frequency (FRQ), white collar-1 (WC-1), white collar-2 (WC-2) and frequency-interacting RNA helicase (FRH) in 31 lichen-forming fungal species.Using a similar reciprocal BLAST-based approach, we F I G U R E 1 Variation of clock genes is lower within climate zones than across climate zones.Schematic representation (a) of pool-seq sampling scheme across three gradients for Umbilicaria pustulata and two gradients for U. phaea.Populations falling into a higher and a lower elevation climatic zones, either Mediterranean/cool temperate or cool temperate/alpine, were analysed for allele frequency differentiation, represented by the negative logarithms of p-values for pairwise Fisher's exact tests (FET) scores at each SNP for populations either across (red) or between (blue) climate zones.(b) U. pustulata core circadian clock gene loci with significant FET scores are almost exclusively found across climate zones, including those that are inside of the same elevation gradient (red).(c) U. phaea core circadian core clock gene loci within climate zones (blue) are more significant relative to U. pustulata, although the most significant SNPs are reliably found across climate zones (red).BIO1: mean annual temperature.ESi: Talavera-Puerto de Pico (Sistema Central, Spain; 3 populations).ESii: Sierra de Gredos (Sistema Central, Spain; 6 populations).IT: Mount Limbara (Sardinia, Italy; 6 populations).MJ: Mt.Jacinto (Mount San Jacinto State Park, CA, USA, 7 populations).SN: Sierra Nevada (CA, USA, 4 populations).Dashed lines denote significant values above FET scores of 2.99, equivalent to −log(0.05).
The populations of each gradient fall along a climate gradient ranging from a low elevation, Mediterranean climate zone to a high elevation, cold temperate climate zone (Figure 1a).To determine whether SNPs in the 11 putative core circadian clock genes varied along this climatic gradient, we compared pairwise FET (Fisher's exact test) scores across and within climate zones (Figure 1a), allowing us to identify significant peaks across climate zones that were independent of any variation within climate zones.For both U. pustulata (Figure 1b) and U. phaea (Figure 1c), pairwise FET results demonstrated that the core circadian clock genes have many loci that are more strongly differentiated across climate zones at the elevation extremes of all gradients ("Across," red lines) than between different gradients but within the same climate zones ("Within," blue lines).
Following this initial analysis, Spearman's correlation coefficients were calculated for each SNP in order to identify those loci with allele frequencies that correlated significantly with BIO1 (mean annual temperature) of the populations (Figure 2a), yielding 105 SNPs across 10 out of 13 homologues in U. pustulata and 61 SNPs in 10 out of 11 homologues in U. phaea that were significantly correlated to BIO1 in at least one gradient per species (Table 1).Of these SNPs, many were not correlated in all gradients per species (e.g. Figure 2c,e) or caused synonymous or non-radical substitutions (e.g. Figure 2d,g).23 SNPs in 8 homologues (U. pustulata) and 8 in 7 homologues (U. phaea) were found to cause non-synonymous substitutions (e.g. Figure 2b,f, Table 1, see Tables S1 and S2 for a full list of analysed SNPs).

| Nucleotide diversity (π) in clock-and temperature-associated genes significantly varies along elevation gradients
In previous analyses of pool-seq data from U. pustulata, ca.30% of genes (2521 out of 8268) could not be assigned a Gene Ontology (GO) term.In order to identify as many putatively circadian clockassociated and temperature-associated genes as possible, we made use of a reciprocal BLAST approach to extract homologues from the well-annotated Neurospora crassa GO term database (https:// cyc. pnnl.gov/ ) that were associated with the circadian clock or temperature responses (Figure 3a).This approach yielded 51 U. pustulata and Identifying substitutions linked to climate in circadian clock genes.This is a schematic representation of the pipeline used to discover putatively functional variation in clock genes along elevation (temperature gradients).Allele frequencies of core circadian clock gene loci were processed using Spearman's rank correlation coefficient to investigate the relationship between allele frequency and mean annual temperature (BIO1).(a) Loci that were significantly correlated (black line) in at least one gradient were extracted (arrows) and further analysed.Three loci each from U. pustulata FRH-1 (b-d) and from U. phaea SKI2-2 (e-g) are shown to illustrate results.Most loci did not vary in all gradients (b, e) and were removed from further analysis.Some loci varied in all gradients, but yielded synonymous substitutions, and thus were not considered further (c, f).Some loci yielded non-synonymous substitutions, sometimes accompanied by changes in the chemical properties of the substituted amino acid (i.e.radical substitutions) (d, g).These substitutions have the highest likelihood of causing functional change.A count of loci with significant correlations in all gradients that yielded radical non-synonymous substitutions, can be found in Table 1, and an expanded table with all significant loci and the chemical properties of their radical substitutions can be found in Table S1.S3) and 37 U. pustulata and U. phaea temperature-associated homologues each (Table S4).
In order to analyse sequence variation of these homologues along the climate gradient, we used Tajima's π as a measure of nucleotide diversity at the individual SNP level.We used the envfit function in the vegan package to fit the available climatic variables and determine pvalues of the correlation of each climatic variable with the nucleotide diversity of the circadian-and temperature-associated loci (Table S5, Figure S1).Of these, BIO1 (mean annual temperature), BIO6 (minimum temperature of coldest month) and BIO11 (mean temperature of coldest quarter) were significantly correlated with nucleotide diversity for both U. pustulata and U. phaea.We further performed Mantel tests in order to further test the correlation between nucleotide diversity with BIO1 at each site (Figure 3), revealing strong correlations for circadian clock-associated genes in U. pustulata (r: .3947,p-value: .0018, Figure 3c) and in U. phaea (r: .4772,p-value: .004, Figure 3d), similarly to temperature-associated genes in U. pustulata (r: .4711,p-value: 8e-04, Figure 3e) and in U. phaea (r: .4427,p-value: .0046, Figure 3f).

| Peripheral genes in the clock-and temperature-associated gene clusters display signatures of selection along elevation gradients
We next examined signatures of selection in the previously identified clock-and temperature-associated loci based on dN/dS ratios (ω).We We observed fewer significant results for U. pustulata genes than for U. phaea genes across all dN/dS analyses in HyPhy (Table 2).The two methods used for testing individual sites (MEME and FUBAR) demonstrated a low number of significant results for U. pustulata circadian (6% and 4%) and for temperature-associated genes (5.41%), while a majority of both U. phaea circadian (25.5% and 51%) and temperature-associated genes (56.8% and 81.1%) displayed significant results.Similarly, the gene-wide test (BUSTED) yielded a low number of significant results for U. pustulata circadian (4%) and temperature-associated genes (8.1%), while U. phaea circadian (11.8%) and temperature-associated TA B L E 1 Radical non-synonymous substitutions along climatic gradients of U. pustulata and U. phaea core circadian clock genes.
Note: PROT_ID: UniProt protein ID of the N. crassa proteins used for TBLASTN (protein → nucleotide).Gene Annot: gene ID from NCBI-submitted annotation.Spearman: number of SNPs with significant Spearman rank correlation results between allele frequency and mean annual temperature (BIO1) along the climatic gradient.Radical?: number of Spearman correlation-identified SNPs causing radical, non-synonymous substitutions at the extremes of the climatic gradient.Shaded cells represent genes that returned no BLAST result or contained no non-synonymous substitutions for that species.
genes (43.2%) had a higher number of significant results.The differences between U. pustulata and U. phaea were most extreme for the branch-site test (aBSREL), where no circadian or temperatureassociated genes had significant results for U. pustulata, while for U. phaea displayed significant branches in both circadian (7.8%, Figure S2) and temperature-associated genes (32.4%, Figure S3).
When looking at site-level results (FEL), nearly all circadian (98%) and temperature-associated genes (100%) in U. phaea yielded significant results for purifying selection.A smaller proportion of U. pustulata circadian (60%) and temperature-associated genes (46%) yielded significant results.Few genes displayed significant results for positive or diversifying selection in either U. pustulata circadian (0%) or temperature-associated genes (5.4%), while the number of significant results for U. phaea circadian (5.9%) and temperature-associated genes (32.4%) approximated those of other analyses like aBSREL.

| Allele frequency changes are significantly different between randomly selected genes and U. pustulata circadian-associated genes
To compare the level of differentiation observed for the circadian-  steeper increase in pairwise mean F ST .An analysis of gene-wise nucleotide diversity across the genome yielded mean π for U. pustulata ranging from 0.000649 to 0.00293, while that of U. phaea ranged from 0.00650 to 0.0227 (Figure S4a,b).Tajima's D for U. pustulata ranged from −0.869 to 1.10, and from −0.905 to 0.676 for U. phaea (Figure S4c,d).
We then performed a Cochran-Mantel-Haenszel (CMH) test, using the populations located at the extremes of each gradient for each species.We compared the −log 10 p-values of the CMH test for the circadian-and temperature-associated gene SNPs to the SNPs in 50 randomly selected genes throughout the genome: our analysis TA B L E 2 dN/dS analysis of circadian-and temperature-associated genes in U. pustulata and U. phaea.

| Core circadian clock genes display non-synonymous substitutions along climatic gradients
For both U. pustulata and U. phaea, we observed higher FET scores for loci in the "Across category" relative to the "Within category," demonstrating that there was greater variation between alleles at different extremes of the climatic gradients than between independent gradients at the same elevation and climate zone.Given the geographical distance between the gradients used in this study, these results imply a strong selective effect of climate on these loci: the U. pustulata gradients at Sierra de Gredos and Mount Limbara are 1200 km apart, while the U. phaea gradients at Mt Jacinto and Sierra Nevada are roughly 700 km apart.
We detected fixed alleles for each SNP (cold_fixed and warm_fixed columns, Tables S1 and S2) to assess whether the warm and cold alleles are fixed at any point in all of the gradients.For U. pustulata, cold and warm alleles were fixed more reliably at the two temperature extremes (e.g. Figure 2b).For U. phaea, allele distributions followed a bell-shaped distribution (e.g. Figure 2f).This outcome may be related to the amount of ecological variation present along the studied gradients: with an altitudinal range of 1900-2400 m, the U. phaea gradients are almost twice as long as the U. pustulata gradients (1000-1200 m), and may thus encompass more climatic diversity and a higher number of ecological zones.Stark differences in mean annual precipitation between the U. phaea gradients, with much lower mean annual precipitation in the Mt.Jacinto gradient than in the Sierra Nevada gradient (Figure S5b) than observed between the U. pustulata gradients (Figure S5a), may also contribute to the greater variation in U. phaea loci.
Another potential source of variation between the U. phaea and U. pustulata loci may be related to the age of the WorldClim 2.1 climate dataset used in this study .A previous study comparing the influence of WorldClim 1.4 (1960WorldClim 1.4 ( -1990) ) and WorldClim 2.1 (1970-2000) climate data on habitat suitability models found a noticeable level of mismatch between these two WorldClim versions (Cerasoli et al., 2022), which may also similarly affect the climate boundary zone predictions of previous studies on U. phaea and U. pustulata (e.g.Rolshausen et al., 2022).However, focusing on climatic zones that co-occur with elevation changes, as well as focusing on changes at the extremes of the gradients (as we do here) may help to insulate the dataset from some of the mismatch in the climate data.Nonetheless, repeating the analyses performed here with more current data, or with future projection data (such as CMIP Phase 62021-2040 20-year average predictions), may paint a different and more accurate picture of where the climate zone boundaries lie for each gradient.

| Signatures of selection in lichen circadian and temperature-associated genes
In our work, we analysed signatures of selection using various dN/ dS-based hypothesis tests using phylogenies.These tests, performed at any site on a particular branch or set of branches?(BUSTED).We find that a subset of circadian and temperature-associated genes in both U. phaea and U. pustulata display significant results (Table 2), although more genes are significant for U. phaea than for U. pustulata, with the exception of the site-level results (FEL).Of the genes that were found to be significant in at least one of the above categories, several are homologous to the same Neurospora gene (wc-1, ras-1, ccg-9).Most significant genes in U. phaea were not found to be significant for U. pustulata, while two of the genes found to be significant in U. pustulata had no discernible homologue in U. phaea (frq, vvd).
dN/dS approaches are used to analyse evolutionary patterns across species (Moreno et al., 2021) as well as within species (Kryazhimskiy & Plotkin, 2008;Peterson & Masel, 2009;Pond et al., 2006;Wilson et al., 2020).We were able to observe some evidence of selection across circadian-associated genes (Table S6), similar to that of the temperature-associated genes we analysed (Table S7).However, several circadian clock genes for which non-synonymous substitutions between warm-and cold-adapted were observed (Table 1) did not yield significant results in our dN/dS analyses.Some limitations have been elucidated for dN/dS approaches, particularly if the methods are applied to intra-specific or intra-population data (Kryazhimskiy & Plotkin, 2008).Despite this, dN/dS approaches continue to be wellcited in the literature, and new tools for the use of dN/dS metrics in intra-specific contexts continue to be developed (Wilson et al., 2020).
Our results provide another example that the application of dN/dS metrics in within-species contexts must be carefully considered, as these metrics may fail to pick up on potentially adaptive variation in genes with well-described adaptive functions.
Another possible source of genetic adaptation along climate gradients in lichen symbioses are co-evolutionary processes between the fungal partners and environmentally adapted strains of lichen photobionts.Environmental factors, such as climate and substrate pH, have been shown to be important drivers of photobiont distribution (e.g.Rolshausen et al., 2017;Škvorová et al., 2022;Steinová et al., 2022;Vančurová et al., 2021).

| Differences between U. pustulata and U. phaea variation in circadian clock-and temperature-associated genes
We observed in the FET analysis of 11 core circadian clock genes that both U. phaea and U. pustulata genes follow similar patterns of within-and across-climate zone genetic differentiation.However, U. phaea has a higher number of variable loci than U. pustulata, and yields a higher number of significant results in terms gene-wide, branch-site and individual site selection tests (Table 2).Furthermore, U. phaea populations show more differentiation between gradients ("Within," blue lines, Figure 1c).These differences may be linked to different rates of recombination in the two species, due to their differing reproductive strategies (asexual for U. pustulata and sexual for U. phaea).Asexual species of lichenized fungi have been shown to have lower genomic diversity than closely related sexually reproducing species (e.g.Grewe et al., 2018).However, asexually reproducing lichenized fungi still show substantial intraspecific differentiation (Dal Grande et al., 2017;Grewe et al., 2018;Onuţ-Brännström et al., 2017;Otálora et al., 2013).This is similar to non-lichenized fungi without discernible sexual structures (e.g.Geiser et al., 1998;Varga & Tóth, 2003), and has led to the assumption that most asexual fungi are not exclusively clonal (Taylor et al., 2015).
Our genetic diversity results seem to corroborate field observations and morphological analyses that U. phaea reproduces sexually while U. pustulata reproduces primarily asexually.First, U. pustulata (0.022) has far fewer SNPs/base than U. phaea (0.079); this difference is yet more pronounced in genic SNPs, with an average of 81 SNPs/gene for U. pustulata and 360 SNPs/gene for U. phaea.U. pustulata shows a stronger pattern of isolation by distance, with a more drastic increase in pairwise F ST in relation to elevational differences (Figure 4c vs. d).Our analysis of nucleotide diversity for all genes across the genomes of U. pustulata and U. phaea further reinforces these results, with the asexual U. pustulata having a lower mean π, which may also be related to differences in reproductive strategies observed between the two species.Empirical evidence for reduced levels of standing genetic variation in asexual versus sexual species, due to the lack of meiotic recombination in asexuals, is also reported from of natural populations of stick insects (Bast et al., 2018).A potential explanation for the differences would be related to the quality or depth of the sequencing; however, despite some minor differences, this is unlikely in our estimation to cause such large differences in the overall patterns observed.These differences are also intermixed between the species, with the ESi population from U. pustulata and the two U. phaea populations being sequenced using 150 bp paired-end chemistry while the other two U. pustulata populations were sequenced using 100 bp paired-end chemistry.Moreover, all populations for both species were sequenced to the same sequencing depth, 90×, and both species have roughly comparable genome assembly sizes (35-38 Mb) and gene numbers (7.5-9.5 k genes) according to PacBio metagenomic sequencing (Singh et al., 2021).
A final considered source of variation between U. phaea and U. pustulata circadian and temperature-associated genes is the much greater difference in latitude between the two species' gradients: the Sierra Nevada and Mt.Jacinto U. phaea gradients are 4.65° of latitude apart, while the U. pustulata gradients are roughly at the same latitude, being 0.5° apart.There is a wide body of literature demonstrating that circadian clock genes vary functionally along latitudinal gradients in animals (Kyriacou et al., 2008;Moreno et al., 2021;Pegoraro et al., 2022), plants (Greenham et al., 2017;Rees et al., 2021) and fungi (Koritala et al., 2020;Koritala & Lee, 2017), making this a potential source of confounding variation between the two U. phaea gradients in our analysis.However, there are no striking differences between U. phaea and U. pustulata in the numbers of loci identified that cause radical, non-synonymous substitutions in the 11 core circadian clock genes (Table 1), and there is no evidence that nucleotide diversity is strongly segregated between the two U. phaea gradients in the NMDS analysis in either the circadian-or temperature-associated loci (Figure 3).

| Circadian clock variation and cold tolerance in lichenized fungi
In this work, we identify temperature as a significant explanatory variable for both circadian and temperature-associated gene variation in two lichen-forming fungi.Most of our analyses are based on mean annual temperature, for which a turnover zone in algal symbiont identity at around 12°C mean annual temperature for all the gradients studied here has been characterized (Rolshausen et al., 2020).
We also observe, however, that other temperature-associated bioclimatic variables are significantly correlated to nucleotide diversity of circadian and temperature-associated genes, specifically BIO6 (coldest temperature of the coldest month) and BIO11 (mean temperature of coldest quarter), which are both significantly correlated with nucleotide diversity across both sets of genes and in both species (Table S5).This suggests that winter temperatures and frost events are important selective forces for these genes.Indeed, the biome borders of the Mediterranean biome, based on vascular plant communities, are best predicted by winter temperatures (mean temperature of coldest month <5°C) (Prentice et al., 1992).
Some species of lichenized fungi are known to be extremely tolerant of freezing stress (Kappen et al., 1996), but whether populations of the same species -which are subject to different climatic selection pressures-uniformly tolerate cold temperatures is less understood.
In controlled environments, lichen individuals survive below-freezing temperatures when they are dry, but not when they are hydrated (Honegger, 2007).However, lichens in the cold temperate biome found at the top of all five gradients studied here are often simultaneously subjected to freezing and wet conditions, for example, during spring.Given that freezing tolerance is a well-known adaptive trait in other sessile organisms such as plants, there may be intra-specific variation in the ability of lichenized fungi to survive chilling and freezing stress in the absence of dry conditions.Furthermore, comparisons between species of the same genus have shown that taxa with oceanic distributions (e.g.Lobaria virens) are less frost tolerant than those inhabiting more frost prone environments (e.g.Lobaria pulmonaria) (Solhaug et al., 2018).Our findings suggest that genes associated with the circadian clock and temperature response provide useful targets to better understand the mechanisms conferring cold tolerance in lichenized fungi, similarly to their role in plants (Dong et al., 2011).
phaea circadian clock-associated homologues (Table did this using HyPhy (Hypothesis Testing using Phylogenies), a tool for positive and diversifying selection analysis.Our HyPhy analysis consisted of running the following methods: BUSTED (Branch-site Unrestricted Statistical Test for Episodic Diversification), which provides a gene-wide test for positive selection; MEME (Mixed Effects Model of Evolution), which provides a site-specific test for positive selection; FUBAR (Fast, Unconstrained Baysian AppRoximation), which infers per-site probabilities of positive or negative selection; aBSREL (adaptive Branch-Site Random Effects Likelihood), which provides a test of whether a proportion of sites have evolved under positive selection for each branch in a phylogeny; and FEL (Fixed Effects Likelihood), which provides a site-specific test for selection pressure along an entire phylogeny.All methods were run using standard settings and using gene trees reconstructed based on chimeric nucleotide sequences taken from pool-seq data using a Maximum Likelihood approach via RAxML.
and temperature-associated genes with the genomic background in U. phaea and U. pustulata, we first investigated SNP density, isolation by distance and nucleotide diversity across the genome in both species.After filtering the dataset, U. pustulata yielded 798,029 SNPs in 35.7 Mb (0.022 SNPs/base), while U. phaea F I G U R E 3 Nucleotide diversity in clock-and temperature-associated genes varies significantly along climatic gradients.Diagram of Neurospora crassa GO terms used for BLAST-based extraction of putative circadian-(a) and temperature-associated (b) gene homologuesin Umbilicaria spp.NMDS plot of nucleotide diversity (π) for circadian-associated genes in populations of U. pustulata for three elevation gradients (c) and in populations of U. phaea for two elevation gradients (d) varies along BIO1 (mean annual temperature, degrees Celsius) gradients.NMDS plot of nucleotide diversity (π) for temperature-associated genes in populations of U. pustulata for three elevation gradients (e) and in populations of U. phaea for two elevation gradients (f) varies along BIO1 gradients (red: lower climate zone, Mediterranean; blue: higher climate zone, cold temperate; purple: intermediate climate zone).Mantel statistic performed based on Spearman's rank correlation and 9999 permutations.
770,905 SNPs in 35.1 Mb (0.079 SNPs/base).Both U. pustulata and U. phaea showed significant isolation by distance based on pairwise F ST (Figure 4c,d), with U. pustulata showing a circadian-associated SNPs were significantly more differentiated than the randomly-selected genes for U. pustulata (p = 1.253 × 10 −6 , Figure4a), while the circadian-associated SNPs in U. phaea and the temperature-associated SNPs in both species were not significantly differentiated from those randomly selected genes (Figure4a,b).
using the command line tool HyPhy, are broadly designed to ask the following sets of questions: (a) Are individual sites subject to pervasive positive or purifying selection?(FEL, FUBAR), (b) Are individual sites subject to episodic positive or purifying selection?(MEME), (c) Are individual branches subject to episodic positive or purifying selection?(aBSREL), and (d) Has a gene experienced positive selection Turnovers of algae along climatic gradients show that some lineages of the common lichen photobiont Trebouxia prefer specific temperature and humidity conditions (e.g., Vargas Castillo & Beck, 2012; Werth & Sork, 2014), including algal lineages associated with both fungal species studied here (Dal Grande et al., 2018; Rolshausen et al., 2020, 2022).It is thus possible that climate either influences the fungal genome directly, or the fungal genome changes in response to the presence of a different photobiont lineage.In both cases, climate would be the ultimate driver of variation in circadian clock-or temperature-associated loci.The broader question of whether fungal and algal partners in the lichen symbiosis have coordinated circadian systems, which contribute to local adaptation of the holobiont, remains to be investigated.
With respect to environmental selection of clock-and temperature-associated genes, we observe that different clock components and different temperature-associated genes are affected in two lichen-forming fungi, and conclude that the evolutionary route to climate-specific nucleotide configurations in these genes is species-specific.Overall, we observe that differentiation patterns for both circadian and temperature-associated genes in U. phaea and U. pustulata follow broad genomic patterns of differentiation.This may indicate that, for both of these species, allele frequency changes across the genome are sufficiently high enough to provide evidence of subspecies differentiation for populations at the gradient extremes.However, more work will be necessary to determine if rangewide high-and low-elevation individuals form genetically distinct groups in these species of lichen-forming fungi.Previous work on birds (Zosterops lateralis) has indicated that, as populations become more differentiated across the genome, localized genomic islands of differentiation(Turner et al., 2005) may become masked by overall divergence along the entire genome(Sendell-Price et al., 2020).

No. genes Whole-gene Per-site Branch-site Site-level
BUSTED: No. of genes with significant (<.1) p-values for BUSTED analysis.MEME/FUBAR: number of genes with sites considered significant (p < .1)by MEME/FUBAR analysis.aBSREL: number of genes with at least one branch considered significant (p < .1)by aBSREL analysis.FEL pos./ divers.: number of genes with sites/cluster of sites found to have significant evidence of positive or diversifying selection by FEL analysis.FEL neg.: number of genes with sites/cluster of sites found to have significant evidence of negative or purifying selection by FEL analysis.