Genetic mechanisms and biological processes underlying host response to ophidiomycosis (snake fungal disease) inferred from tissue‐specific transcriptome analyses

Emerging infectious diseases in wildlife species caused by pathogenic fungi are of growing concern, yet crucial knowledge gaps remain for diseases with potentially large impacts. For example, there is detailed knowledge about host pathology and mechanisms underlying response for chytridiomycosis in amphibians and white‐nose syndrome in bats, but such information is lacking for other more recently described fungal infections. One such disease is ophidiomycosis, caused by the fungus Ophidiomyces ophidiicola, which has been identified in many species of snakes, yet the biological mechanisms and molecular changes occurring during infection are unknown. To gain this information, we performed a controlled experimental infection in captive Prairie rattlesnakes (Crotalus viridis) with O. ophidiicola at two different temperatures: 20 and 26°C. We then compared liver, kidney, and skin transcriptomes to assess tissue‐specific genetic responses to O. ophidiicola infection. Given previous histopathological studies and the fact that snakes are ectotherms, we expected highest fungal activity on skin and a significant impact of temperature on host response. Although we found fungal activity to be localized on skin, most of the differential gene expression occurred in internal tissues. Infected snakes at the lower temperature had the highest host mortality whereas two‐thirds of the infected snakes at the higher temperature survived. Our results suggest that ophidiomycosis is likely a systemic disease with long‐term effects on host response. Our analysis also identified candidate protein coding genes that are potentially involved in host response, providing genetic tools for studies of host response to ophidiomycosis in natural populations.

caused by the chytrid fungus Batrachochytrium dendrobatidis (Bd) (Berger et al., 1998), and white-nose syndrome (WNS) in hibernating bats caused by the fungus Pseudogymnoascus destructans (Frick et al., 2010).Bd infection has been documented in hundreds of amphibian species and is responsible for amphibian species declines worldwide (Berger et al., 1998;Fisher & Garner, 2020;Kilpatrick et al., 2010).Infection with Bd disrupts the integrity of amphibian skin, causing host mortality because skin is a critical organ for amphibians involved in physiological activities such as respiration, ion balance, hydration and defence against other pathogens (Rosenblum et al., 2013;Voyles et al., 2009).
Similarly, white-nose syndrome is one of the most damaging infectious disease epidemics in hibernating bats, having been identified in over 60 species in all temperate regions of northern hemisphere (Hoyt et al., 2021) and linked to extirpations of entire populations of many bat species (Frick et al., 2015;Hoyt et al., 2021).
In cases of WNS, the lesions caused by P. destructans affect the skin of the bats' ears, nose, wing and tail tissue, disrupting crucial regulatory functions like thermoregulation (Meteyer et al., 2009), which results in evaporative water loss and higher activity levels during hibernation (McGuire et al., 2017), leading to fat loss, starvation and eventual death (Hoyt et al., 2021;Reeder et al., 2012).
Fungal pathogens have been recognized in numerous additional taxa, including sea turtles (Sarmiento-Ramírez et al., 2010), corals (Kim & Harvell, 2004), lizards (Schumacher, 2003), honeybees (Vanengelsdorp et al., 2009) and various plants (Anderson et al., 2004;Brown & Hovmoller, 2002), but estimates of disease prevalence and understanding of host responses to such pathogens are lacking in most species (Fisher et al., 2012;Ghosh et al., 2018;Williams et al., 2002).The recent emergence and global spread of fungal pathogens have led to increased surveillance efforts for known pathogens and the need to understand how specific diseases impact infected individuals to understand the pathological mechanisms that underlie disease (Anderson et al., 2004;Daszak et al., 2000;Fisher et al., 2012;Williams et al., 2002).
One example of a potentially impactful yet poorly understood wildlife fungal disease is ophidiomycosis (snake fungal disease; SFD).Ophidiomycosis is a recently identified disease in snakes caused by the fungus Ophidiomyces ophidiicola that has been detected in many free-ranging and captive snake species and has the potential to cause widespread morbidity and mortality in snakes (Allender, Raudabaugh, et al., 2015;Lorch et al., 2015Lorch et al., , 2016;;Rajeev et al., 2009).Ophidiomyces ophidiicola is a generalist fungus and is known to infect a wide range of snake species with different life histories, irrespective of taxonomy or habitat (Burbrink et al., 2017).
Ophidiomycosis has been documented in more than 30 species of wild snakes in the United States and Europe (Lorch et al., 2016), with reported instances in Australia (Sigler et al., 2013) and Southeast Asia (Sun et al., 2022;Takami et al., 2021).Ophidiomycosis was first reported in the United States in the mid-2000s but was likely present as early as the 1940s (Allender et al., 2011;Clark et al., 2011;Lorch et al., 2021).A recent population genomic study of O. ophidiicola confirms that ophidiomycosis is a recently emerging infectious disease in North American snakes that has spread after multiple transcontinental introductions (Ladner et al., 2022).
Clinical signs of ophidiomycosis vary among individuals and range from lethargy, skin lesions and excessive shedding, to skin crusts, granulomas, corneal opacity and ulcers on the head and body in more severe cases (Lorch et al., 2016).However, the mode of infection, specific mechanisms of pathology and physiological responses by infected individuals are unclear despite the value of this information for understanding and mitigating the impact of this disease (Davy et al., 2021;Lorch et al., 2016).Specifically, the two key elements of ophidiomycosis pathology that are still unknown are: (1) how snakes respond to O. ophidiicola infection and (2) the influence of temperature on disease severity of infected individuals.
A previous histopathological study (Guthrie et al., 2016) identified O. ophidiicola infection as localized to skin (Lorch et al., 2015), so we expect the highest fungal activity on snake skin.Second, since snakes are ectotherms, meaning that their thermoregulation is dependent on external temperatures, we predict that temperature would have a crucial impact on overall host response to infection.
Field evidence suggests that many infected snakes prefer higher temperatures (McBride et al., 2015;McCoy et al., 2017), potentially enabling them to recover from ophidiomycosis by facilitating a more robust immune response (Lind et al., 2018).Transcriptomic studies of controlled exposure experiments offer an approach to uncover the impact of infection and the effects of external conditions on host response (see examples in Ghosh et al., 2020).By identifying differentially expressed genes between treated and control hosts, we can understand what biological processes are occurring differently in host tissues due to conditions such as infection status or temperature (Davy et al., 2020;Ju et al., 2002;McNew et al., 2022;Navajas et al., 2008;Reemers et al., 2009;Zhan et al., 2019).
Here, we studied the response of a snake host to ophidiomycosis by performing controlled O. ophidiicola exposure experiments in Prairie rattlesnakes (Crotalus viridis) at two different temperatures and then conducting transcriptomic analyses of gene expression in multiple tissues of infected and control individuals.Prairie rattlesnakes are a common and widely distributed snake species with known infections in the wild (Allender et al., 2020) and are also closely related to many other susceptible species (Burbrink et al., 2017), which makes them a good model to study ophidiomycosis (Dibadj et al., 2021).
We describe the changes in gene expression occurring in different host tissues and use this information to infer the possible biological mechanisms underlying O. ophidiicola infection.We also test whether temperature differences affect host response among infected hosts.
Finally, we describe a list of candidate genes putatively involved in host response to ophidiomycosis that could be helpful for genetic monitoring of vulnerable populations and species of snakes, as done for other wildlife diseases (Ghosh et al., 2020).

| Controlled experiment varying infection status and temperature
Twelve adult (six male, six female) Prairie rattlesnakes (C.viridis) were caught in the wild as adults or young adults from a single population in North Dakota at least 10 years prior to the study and were part of a long-term research collection at the University of Illinois since 2017.For the controlled O. ophidiicola exposure experiment, snakes were divided into two temperature groups, which were maintained at 20°C or 26°C following a 2-week acclimatization period (Table S1).
These temperatures were chosen as thermal thresholds for the snakes (Allender, Raudabaugh, et al., 2015).Four animals (two males, two females) from each temperature group were randomly selected for experimental infection, while two animals (one male, one female) were maintained as controls (Figure 1a).All snakes were free of clinical signs suggestive of ophidiomycosis based on criteria previously described (Baker et al., 2019) for over 24 months prior, at the start of the study and prior to inoculation.The O. ophidiicola-free status of all individuals was confirmed by testing skin swabs collected prior to inoculation using a specific qPCR assay (Allender, Bunick, et al., 2015); all tested negative for Ophidiomyces DNA.All snakes were offered pre-killed mice every 1-2 weeks, and water was provided ad libitum.
A single snake in the infected group at 26°C died in the first 2 weeks post infection due to bacterial sepsis (unrelated to ophidiomycosis) and was removed from the study.A pure O. ophidiicola isolate (UI-VDL # 12-34933) cultured from an eastern massasauga rattlesnake (Sistrurus catenatus) was injected intradermally (0.1 mL of a pure culture containing 109,000 CFUs) for snakes assigned to the infection group; control animals were inoculated with a similar volume of sterile saline.All injections were performed over the dorsal mid-body of each snake.Once weekly, physical examinations were performed, and skin swabs were collected for O. ophidiicola DNA quantification using the same qPCR assay described above.Additional details regarding animal procedures, fungal preparation, animal care during the study, qPCR assay methods, survival, clinical signs and histopathology are provided in Allender et al. (2020).
The study was conducted for 90 days post infection (dpi) and snakes were euthanized when they developed one or more of the clinical thresholds established by the IACUC protocol (10% or more weight loss in a 1-week period, lack of righting reflex and extreme lethargy and/or lack of eating, drinking or visual activity due to granuloma growth), or at the end of the 90-day study period in cases with no or less severe clinical signs.Humane euthanasia was performed by complete anaesthesia with ketamine injected intramuscularly, followed by sodium pentobarbital injected intravenously.All infected snakes in the 20°C group (N = 4) were euthanized prior to the end of the study due to the severity of clinical signs, while only one animal was euthanized in the 26°C group due to severity of clinical signs.
The remaining two inoculated animals in the 26°C group survived with mild-moderate clinical signs and were euthanized at the termination of the study (90 dpi).All control animals survived throughout the study period.All procedures were approved by the University of Illinois Institutional Animal Care and Use Committee (IACUC;Protocol: 19165).Following euthanasia, each animal received a necropsy examination separately to collect samples of liver, kidney and skin.We chose liver and kidney tissues for transcriptomic analysis as these tissues are involved in metabolic activities and previous studies in snakes have indicated disruption in normal metabolism as a sign of host response to ophidiomycosis (Agugliaro et al., 2019).
Additionally, skin has been shown histologically to be the site of fungal activity (Lorch et al., 2016), and thus, was also included for transcriptomic analysis.

| RNA extraction, construction, and sequencing of RNA-Seq libraries
RNA was extracted from snake tissue using a Qiagen RNEasy kit according to the manufacturer's recommendations (Qiagen RNEasy, Valencia, CA).RNA samples were analysed using spectrophotometry (NanoDrop 1000; Thermo Fisher) to quantify concentration and purity, then stored at −80°C.RNA-Seq libraries were prepared for each sample with Illumina's TruSeq Stranded mRNA-seq Sample Prep kit with Unique Dual Indexes to prevent index switching.Libraries were made by enriching mRNA by selecting polyA containing mRNA using poly-T oligo attached magnetic beads.Enriched mRNAs were then converted to cDNA and adaptors were ligated to prepare individual RNA-Seq libraries.The libraries were then pooled; quantitated by qPCR and (PE) reads (2 × 150 bp) were sequenced on one SP lane for 151 cycles from both ends of the fragments on an Illumina NovaSeq 6000 sequencer at the University of Illinois Roy J. Carver Biotechnology Center.Following sequencing, we removed one skin and one liver sample due to low number of raw sequences generated.Ultimately, we analysed sequence data from 10 liver, 11 kidney and 10 skin tissues from each snake with different temperature condition and infection status (Figure 1a; Table S1).

| Sequence processing and expression quantification
We removed adapter sequences and clipped poor-quality bases (quality score <20) from both ends of raw reads from the sequencer using Trimmomatic (Bolger et al., 2014), aligned filtered reads to the annotated C. viridis reference genome (UTA_CroVir_3.0;GenBank assembly accession: GCA_003400415.2) using hisat2 (Kim et al., 2015) with the -downstream-transcriptome-assembly option and reporting primary alignments.We next assembled transcripts for each sample using StringTie (Pertea et al., 2015) default parameters and C. viridis reference annotation file to guide assembly and merged sample transcripts using StringTie.Next, a transcript count matrix was next created with featureCounts (Liao et al., 2013) which excluded chimeric fragments and only retained mapped fragments in which both read pairs successfully aligned to the reference.We only counted reads that matched to the exons present in the merged annotation file.Ultimately, we generated a table of the count data which reports, for each sample, the number of sequence fragments that were assigned to each transcript.
To identify the number of total transcripts from snake tissues that had fungal origins, we mapped the aligned filtered reads to O. ophidiicola reference genome (GenBank assembly accession GCA_002167195.1) using hisat2.We expect more fungal transcripts to be present in transcriptome of tissues where fungus is active and thus, would have higher alignment rate to the O. ophidiicola reference genome.

| Differential gene expression analysis
We conducted differential gene expression (DGE) analyses from the counts of all assembled transcripts (N = 107,167) for each tissue separately using DESeq2 (Love et al., 2014) in R (Team, 2013) following the workflow as described by the authors (Love et al., 2014(Love et al., , 2016)).To determine how different samples clustered in each tissue due to overall gene expression, we used logtransformed count data of all transcripts and identified samples clusters using the find.clustersfunction in the R package adegenet (Jombart & Ahmed, 2011).Successive K-means were run with an increasing number of clusters (K), and for each model, a statistical measure of goodness of fit was estimated using Bayesian information criterion (BIC), and the optimum number of clusters was chosen for the model with lowest BIC.We also performed hierarchical clustering using dissimilarity matrix calculated from individual logtransformed count data for each tissue as a complementary approach to identify clusters.
To identify differentially expressed genes (DEGs) due to infection status and temperature treatments, we designed a F I G U R E 1 Methodological overview of the study and overall individual gene expression profiles.(a) Schematic showing the experimental design for the controlled Ophidiomyces ophidiicola infection trial to study the host response.Prairie rattlesnakes (Crotalus viridis) were exposed to the fungus O. ophiodiicola ('infected') or sham inoculations ('control') at either higher (26°C) or lower (20°C) temperatures.N = number of samples within each group.After the end of the study (90 dpi), RNA was extracted and sequenced from liver (N = 10), kidney (N = 11) and skin (N = 10) tissues of each individual.Differential gene expression (DGE) analysis was performed using infection status and temperature as fixed effects, and the interaction between infection and temperature.Each tissue (liver, kidney and skin) was analysed independently.Functional enrichment and pathway networks were then inferred from biological functions of differentially expressed genes (DEGs).(b-d) Tissue-specific expression profile and individual clustering based on normalized read count data of all expressed genes in (b) liver, (c) kidney, and (d) skin.O. ophidiicola infected snakes at the lower temperature (red circles) clustered independently and were separated on PC1 axis (x-axis) in all three cases except for one skin sample (shown by arrow in panel d).Infected snakes at the higher temperature (red triangles) had similar profiles as control snakes at the higher (blue triangles) or the lower (blue circles) temperature.The dotted circles represent the number of clusters identified using discriminant analysis (see Figure S1 and Section 2 for details).measured the effect between temperature while controlling for infection, that is, if temperature difference had an effect on DGE within the infected group.The two sexes were equally represented in different groups but due to small sample sizes within each group and the absence of sex tissues, we did not account for differential gene expression due to sex differences in this study.
To remove potential false positives due to low expression and high transcript count variability due to small number of biological replicates in the study, we applied log-fold-count (LFC) shrinkage (Zhu et al., 2019).We transformed LFC change data using 'adaptive shrinkage' from the ashr package (Stephens, 2017) in R to shrink the change in expression.We then removed transcripts with low transformed counts and only considered transcripts with a false discovery rate (FDR) adjusted p-value of <.05 to be significantly differentially expressed (Figures S13-S15).This approach minimized noise in count data due to high individual variation and is suitable for low sample sizes (Schurch et al., 2016;Zhu et al., 2019).Since shrinkage only captures DGE for genes with high effect sizes, only differentially expressed genes (DEGs) that are most biologically significant are retained (Conesa et al., 2016).
We used the contrast function in DESeq2 to identify genes that were significantly differentially expressed purely due to fixed effects (infection and temperature) or due to differences in temperatures between infected tissues (i.e.interaction; Figure S4 and Figure 5a).

| Weighted gene co-expression network analysis
Besides identifying individual genes that showed differential expression in each tissue due to differences in treatment and temperature conditions, we also identified networks of co-expressed genes (termed 'modules') from the normalized counts of all expressed genes (Zhang & Horvath, 2005).Estimating co-expression patterns can provide insights into the biological processes that underlie the complex cascade of events that lead to the phenotypic differences observed due to different treatment effects (Nacu et al., 2007).The nodes of these co-expression gene networks correspond to gene expression profiles, and edges between genes are determined by the pairwise correlations between the level of expression of each gene (Zhang & Horvath, 2005).
We used weighted gene co-expression network analysis (WGCNA) (Langfelder & Horvath, 2008)  N genes_skin = 44,807).We then merged modules that were correlated to each other with R 2 > .75 to get a final set of merged modules.Modules were identified by hierarchical clustering of signed Topological Overlap Matrix (TOM) (Yip & Horvath, 2007) using a soft threshold (β) to assign a connection weight to each gene pair and a minimum of 30 co-expressed genes to be assigned to a module (Harder et al., 2019).β was chosen as the lowest power for which the scale-free topology fit index of all genes reached 0.9 within each tissue type (β liver = 12; β kidney = 8; β skin = 12) as suggested by previous work (Godwin et al., 2021).
Next, we identified modules that were significantly associated (p < .05after applying false discovery rate method to correct for multiple testing) with the external temperature and infection treatments.We only retained modules that were significantly correlated to (a) either infection or (b) both infection and temperature.Within each significant module, we also identified the module membership and gene significance for each gene within the module.Module membership measures correlation between a gene's expression profile with the module eigengene of a given module.Highly connected genes within a module are more likely to have higher module membership values to the respective module.Gene significance is the measure of correlation of a gene with an external treatment factor (infection and temperature) and indicates the biological significance of a module gene with respect to the fixed effects of our experimental design.For many of our final modules, genes with higher module membership had a positive and significant correlation with gene significance meaning that highly connected genes within a module are more likely to be associated with experimental treatments.Finally, we matched the genes belonging to final set of modules that were also differentially expressed due to the infection treatment only.

| Functional enrichment analysis
For each module that was significantly associated with infection status and DEGs in each tissue, we performed a functional enrichment analysis to identify which biological process gene ontology (GO:BP) terms and KEGG pathways were overrepresented in those gene clusters, as compared to the genome-wide GO complements of Anolis carolinensis.We performed GO enrichment analysis for all DEGs that were upregulated or downregulated within each tissue separately.We first converted candidate annotated genes in C. viridis genome to A. carolinensis Ensembl IDs using DAVID (Huang da et al., 2009a(Huang da et al., , 2009b) ) and performed functional enrichment analysis using gProfiler (Raudvere et al., 2019).We used g:SCS method (significance threshold <0.05) for computing multiple testing correction for p-values gained from GO and pathway enrichment analysis to account for non-independence among multiple tests among GO terms (Raudvere et al., 2019).We then created phylogenetic trees and functional networks using ShinyGO (Ge et al., 2020) and Pathview (Luo et al., 2017) to visualize enriched pathways.was significantly associated with skin lesions (Allender et al., 2020).

| Identification of potential ophidiomycosis host response loci
All infected snakes at 20°C were euthanized before the end of the study due to the severity of clinical signs, while only one infected snake at 26°C was prematurely euthanized (62 dpi; Table S1).The remaining two infected snakes at 26°C survived with mild to moderate clinical signs until the end of the study (90 dpi).All uninfected control animals survived throughout the study period (Table S1).
RNA was isolated and sequenced from liver, kidney and skin tissues collected from each infected and uninfected snake after euthanasia.Following sample processing and evaluation of library sequence quality, we retained transcriptome sequence data from 10 livers, 11 kidneys and 10 skin tissues from snakes subjected to different temperatures (26°C vs. 20°C) and infection statuses (infected vs. control) (Figure 1a; Table S1).Our sequencing of the RNA-Seq libraries resulted in the generation of a mean of 31.7 million read pairs (SD = 16.8 million) per sample (Table S2).After adapter trimming and removal of low-quality reads, the overall alignment rate for the remaining filtered reads to the C. viridis reference genome for each sample was 59.3% ± 14.8% (mean ± SD; Table S2) with 8.3% ± 2.3% of all aligned read pairs successfully assigned to the annotated regions of the C. viridis assembly.We measured alignment rate as the percentage of transcripts that mapped uniquely to the reference genome.The alignment percentage to the C. viridis reference genome was lowest in skin tissues (liver = 61.6 ± 6.2%; kidney = 65.0 ± 3.7%; skin = 50.8± 23.3%) and one sample had only 1% of total transcripts aligning (Table S2).
We first visualized how samples clustered based on overall gene expression data using multivariate analyses.PCA plots generated from log-normalized read counts of all assembled transcripts indicated that infected snakes at the lower temperature had an expression profile that differed from all other tissue-treatment combinations for all three tissue types (Figure 1b-d).Both hierarchical clustering and discriminant analysis on principle components (DAPC) clustered infected snakes at 20°C as one group across tissue types with all other samples forming either one or two groups consisting of mixtures of samples (Figure 1b-d; Figure S1).The single exception to this pattern was one infected skin sample at the lower temperature, whose expression profile was similar to infected skin samples at the higher temperature (Figure 1d).Strikingly, infected snakes maintained at 26°C showed a similar overall gene expression profile as the control snakes.

| Fungal gene expression in different tissues
To identify if fungal transcripts were present in our RNA-Seq data from different tissues, we mapped the aligned filtered reads to the O. ophidiicola reference genome.We expected that tissues with more actively growing fungus would have greater expression of fungal transcripts and thus, would have a higher alignment rate to the O. ophidiicola reference genome.We observed that most of the fungal transcripts in our samples were identified in skin tissues at the lower temperature (29.8% ± 37.4%; mean ± SD) whereas all other samples (including skin tissues at the higher temperature treatment) had <0.5% fungal transcripts (Figure S2).We also found <0.03% transcripts in control samples at both temperature conditions mapped to the fungal genome (0.025% ± 0.032%) but were not detected by qPCR of cotton swabs (data not shown).The low levels in control samples suggest that these transcripts do not indicate active infections but rather represent statistical artefacts or conserved regions between snake and fungal genomes.Next, we identified which fungal genes were being expressed on infected snake skins at the lower temperature based on highest transcript counts, that is, highest depth of coverage.For each skin sample with fungal transcripts, we observed peaks of high expression at the same region of the fungal genome (see Figure S3).Since the O. ophidiicola reference genome is not yet annotated, we extracted the reference sequence corresponding to the highest transcript peaks and used BLAST (https:// blast.ncbi.nlm.nih.gov/ ) to identify homologous genes in other fungal genomes.We identified expression of genes than encode myosin class V proteins (identity 80.1%, e-value = 0), chaperone protein DnaK (identity 86.6%, e-value = 0) and beta-glucosidase

Ophidiomyces ophidiicola infection in different tissues
We conducted differential gene expression (DGE) analyses from the counts of all assembled transcripts for each tissue separately (Figure 2a-c).Our analysis design included two fixed effects: temperature ('Temperature'; low = 20°C vs. high = 26°C) and infection status ('Infection'; O. ophidiicola infection vs. control).We also included an interaction term ('Infection × Temperature') to identify differentially expressed genes (DEGs) in infected group due to differences in temperature (Figure 1a).We identified a The number within each plot represents the number of differentially expressed genes.There was no interaction between infection and temperature in the analysis of skin tissue.Panels d-f show the number of upregulated (green) and downregulated (pink) genes unique to each fixed effect.There were more upregulated than downregulated genes for each effect in (d) liver, whereas more genes were downregulated due to infection and temperature in (e) kidney.More genes were upregulated in (f) skin due to infection.identified five modules significantly associated with infection status in liver tissue (Table S3), including genes enriched in the Peroxisome KEGG pathway (Figure S6).Most of the genes within the peroxisome pathway were downregulated in liver tissues from infected snakes (Figure S7).Overall, our DGE analysis suggests that O. ophidiicola infection is associated with stress-induced liver inflammation, as well as lower steroid production and protein metabolism, which potentially leads to lipid and protein accumulation within hepatic cells.
In kidney tissues from infected snakes, most of the overexpressed genes due to O. ophidiicola infection reduce protein metabolism by inhibiting proteolysis, protein degradation and catalytic activity in kidney cells (Figure 4a).Among the downregulated genes in kidneys of infected snakes, most belong to biological pathways that regulate steroid metabolism, localization transmembrane transport, integrin synthesis and ribosome biosynthesis (Figure 4b; Figure S7).Just as in liver, steroid synthesis was also downregulated in kidney, which affects many cell-signalling pathways, including electrolyte balance and calcium ion signalling and transport.Seven gene network modules were significantly associated with infection status (Table S3) and were enriched in Ribosome and N-Glycan biosynthesis KEGG pathways (Figure S8).Taken together, our DGE and functional enrichment analysis in the kidneys indicate that O. ophidiicola infections are likely associated with reduced kidney function, nephrotoxicity, loss of protein metabolism and muscular weakness.
In contrast to our findings in liver and kidney, we found little evidence for differential gene expression in skin tissues between infected and control individuals (Figure 2).The overall gene expression was also lower in skin (Table S2), partly because skin lesions tend to be necrotic in cases of ophidiomycosis (Lorch et al., 2016), and most of the transcripts isolated from the skin of infected snakes had fungal origins (Figure S2; Table S2).Many of the upregulated genes were part of the immune response and stress response pathways, including JCHAIN, which links immunoglobulin antibodies, and PSMB8, which is likely involved in the inflammatory response pathway (Table S4).IGLV5 encodes the variable domain of the immunoglobulin light chains that participates in the antigen recognition.Transcription factors were also upregulated in the skin, such as ATF3, which binds to many genomic regions that contain genes involved in cellular stress responses.Among the genes that were downregulated (Table S4), the WNT10B gene is expressed in epidermal keratinocytes and plays crucial roles in regulating skin development and homeostasis whereas downregulation of PLGC2 is associated with immune system disorders.In skin, three modules were significantly correlated with infection status and contained genes involved in DNA damage repair, cell cycle regulation and their associated metabolic pathways like nucleic acid metabolism (Figure S9).The in infected liver and kidney tissues when compared between higher versus lower temperature.Numbers above bars indicate numbers of genes in each infected tissue that were upregulated (green) and downregulated (pink) at the higher temperature as compared to the lower temperature.(b-e) Biological pathway networks for genes that were (b) upregulated in liver, (c) downregulated in liver, (d) upregulated in kidney and (e) downregulated in kidney of infected snakes at the higher temperature as compared to infected snakes at the lower temperature.Biological pathway networks show relationships between enriched pathways (black nodes) connected by edges based on percentage of genes shared between a pair of GO:BP terms.Thicker edges represent more overlap of genes between a pair of GO terms.

| Effect of temperature on differential gene expression within Ophidiomyces ophidiicola-infected snakes
We wanted to identify if infected snakes showed different patterns of gene expression due to difference in temperature (26°C vs. 20°C).
Since we identified no DEGs due to interaction in skin (Figure 2c), we focused on liver and kidney tissues.A total of 91 genes were differentially expressed in the livers of infected snakes due to the higher temperature when compared to the lower temperature, of which 41 were upregulated (45%) and 50 were downregulated (55%; Figure 5a).Functional enrichment analysis showed fibrinolysis pathways were upregulated in infected snakes at the higher temperature (Figure 5b).Fibrinolysis is the enzymatic breakdown of fibrins in blood clots.This indicates that at the lower temperature, blood coagulation processes were ongoing in the livers of infected snakes, which could potentially lead to fibrosis.Most of the downregulated processes in liver due to interaction were organic hydroxy compounds metabolic processes, possibly playing a role in alcohol metabolism (Figure 5c) and thus, adding to already lowered protein metabolism due to ophidiomycosis (Figure 3).In kidney tissues, we identified 67 DEGs due to interaction between treatment and temperature, of which 56 genes (84%) were upregulated and 11 (16%) genes were downregulated (Figure 5a).Upregulated genes in infected snakes at the higher temperature were involved in response to xenobiotics and metabolism of ether lipids (Figure 5d).
Pro-inflammatory responses like production of macrophage inflammatory proteins and chemokines, and neutrophil and myeloid leucocyte migration were in infected snakes at the higher temperature (Figure 5e).Gene regulation related to blood recirculation and reduced inflammatory responses in infected snakes at the higher temperature possibly means that at higher temperatures, infected snakes showed signs of recovery whereas infected snakes continued to fight infections at the lower temperature.Temperature played a crucial role in overall host response to O. ophidiicola infection (Figure 1) as all the infected snakes at the lower temperature were euthanized due to more severe clinical signs, whereas two of three infected snakes at the higher temperature survived the experiment (Table S1).

| Identification of potential ophidiomycosis host response loci
Lastly, we wanted to explore the sequence variation in DEGs as these variant loci could contribute towards disease resistance or susceptibility (e.g.variation in immune or metabolic genes).Within the 906 DEGs that were differentially regulated in at least one tissue due to infection and the interaction between infection and temperature (Data S2), we identified 7700 exons in the C. viridis genome (total exon length ~ 1.93 Mb).We found 2354 non-synonymous SNPs within the exons of DEGs that were polymorphic within C. viridis genomes, and all analysed genomes had high diversity within these functional loci (mean heterozygosity = 0.187; SD = 0.032; Figure S12).These protein-coding functional variants could underpin variability in host response and thus carry utility as diagnostic markers for future ophidiomycosis susceptibility studies in other vulnerable snake species.

| DISCUSS ION
Emerging fungal pathogens are a major threat to biodiversity, as evident by the devastating impacts of diseases like chytridiomycosis in amphibians and white-nose syndrome (WNS) in hibernating bats.
Multiple studies in each disease system have shown the complexity of host-pathogen relationships, including variability in pathogenicity of different fungal strains (McDonald, Ellison, et al., 2020;Piovia-Scott et al., 2014), effects of coinfection (Davy et al., 2018;McDonald, Longo, et al., 2020), difference in response and/or susceptibility of different host species under similar or dissimilar environmental conditions (Davy et al., 2020;Ellison et al., 2020;Eskew et al., 2018;Poorten & Rosenblum, 2016).Compared to chytridiomycosis and WNS, ophidiomycosis is a more recently identified disease (Lorch et al., 2016) and as such there is a lack of information on how snakes respond to infections by this fungus (Burbrink et al., 2017).
Our study is the first molecular genetic study aimed at understanding how snake hosts respond to O. ophidiicola infection by analysing differences in gene expression between infected and control hosts in a species with known susceptibility in the wild (Allender, Raudabaugh, et al., 2015).These results provide insights into what biological processes underlie O. ophidiicola infection, which can help to understand this fungal disease system.Additionally, we have identified a list of candidate genes with potential roles in host response, and functional variation within those genes could potentially be used to monitor vulnerability of snake populations and species in the wild.

| Insights into ophidiomycosis host pathology from gene expression data
Ophidiomyces ophidiicola infects a snake host when the fungus enters the epidermis through an area of injured or scarred skin and fungal invasion is generally limited to the epidermis (Lorch et al., 2016).
From the expression data in this study, we also noted that fungal transcripts were only observed in skin tissues of infected snakes maintained at the lower temperature (Figure S2).Fungal genes that were most expressed included highly conserved fungal myosin class V proteins (Tang et al., 2018), which have been described in chytrid fungi (Prostak et al., 2021) and shown to localize in the actively growing fungal hyphae (Renshaw et al., 2020).The presence of these transcripts is corroborated with previous histopathological studies that show growing fungal hyphae in snake dermis and epidermis (Last et al., 2016;Sun et al., 2022).Other fungal proteins that were highly expressed on host skin (i.e.DnaK and beta-glucosidase 4) are known to be active during fungal pathogenesis and damage host skin by breaking down complex macromolecules like keratin, respectively (Burnie et al., 2006;Huang et al., 2021), indicating how O. ophidiicola damages snake skin and likely produces the skin lesions associated with ophidiomycosis.
In terms of host pathology, many wildlife fungal diseases have skin lesions as a characteristic clinical sign of infection (Baker et al., 2019;Berger et al., 1998;Kilpatrick et al., 2010), which suggests a similar biological host response among different disease systems.For example, genes responsible for activating host pathways like inflammatory responses or skin regeneration (keratin synthesis) were upregulated in O. ophidiicola-infected snakes, which has also been reported to be upregulated in animals infected with Bd (Ellison et al., 2014;Eskew et al., 2018;Poorten & Rosenblum, 2016) and P. destructans (Davy et al., 2020).Snakes with ophidiomycosis have been observed to shed more frequently, and necrotic skin due to penetration of fungal hyphal into the epidermis is replaced with new epidermis via ecdysis (Allender et al., 2018;Lorch et al., 2015).
Skin physiology and ecdysis behaviour are highly variable in different snake species and dependent on the environmental properties of their habitat (Tu et al., 2002).Thus, we suspect that differences in skin physiology among different snake species might play a role in vulnerability to ophidiomycosis (Allender et al., 2018).Future research should compare host skin response in vulnerable snakes from different habitats, as snakes in different environments (e.g.aquatic vs. terrestrial) have been described to vary in their response to ophidiomycosis (McKenzie et al., 2018).
Mortality occurs in cases of severe ophidiomycosis, but the cause of host death is speculated to be complications from the fungal infections, rather than the fungus itself (Lorch et al., 2016).For example, mycotic alterations in lungs and granulomatous pneumonia were implicated as the causes of death in O. ophidiicola-infected garter snakes (Dolinski et al., 2014).In support of this idea, we found that genes associated with inflammatory response, cell proliferation and metabolism were upregulated in the liver tissues of our snakes.
These physiological changes are hallmarks of chronic liver diseases that are characterized by fibrosis and cirrhosis of the liver (Younossi et al., 2011).In kidney tissues of infected snakes, downregulation of steroid synthesis pathways may affect normal cell-signalling and lowered protein synthesis gene expression may cause muscle atrophy, impaired growth of new muscle fibres and loss of kidney function (Wang & Mitch, 2014).Chronic infections may have negative impacts on internal organs over time and disrupt normal organ physiology (Verant et al., 2014), which is likely to alter behaviour of the host.Our DGE results from internal tissues suggest that O. ophidiicola infection likely leads to chronic inflammation in internal organs.Ophidiomycosis has been previously suggested to be a systemic disease (Lorch et al., 2016), and biological processes disrupted by long-term O. ophidiicola infection may provide a link to the behavioural and physiological changes such as lethargy (McKenzie et al., 2020) and higher metabolic rates (Agugliaro et al., 2019) documented in free ranging and captive snakes with ophidiomycosis.
Finally, the observation that most of the differential gene expression was identified in internal organs as opposed to the skin (Figure 2) indicates that most of the physiological changes (i.e.immune response and metabolic activity) are likely to occur within the host's body.Ophidiomycosis can be difficult to diagnose in the wild (Davy et al., 2021)  in Bd infections, an early adaptive immune response is likely beneficial to the host (Grogan, Cashins, et al., 2018) and plays a crucial role in identifying disease susceptibility to chytridiomycosis (Ellison et al., 2014;Eskew et al., 2018;Zamudio et al., 2020).Similarly, differences in the timing of immune response activation can determine susceptibility to WNS in different bat species (Lilley et al., 2019).

| Effects of temperature
Environmental temperature plays a crucial yet complex role in overall ophidiomycosis susceptibility as both host response to fungal infection and fungal activity are temperature dependent (Robert & Casadevall, 2009).In our controlled experiments, all infected snakes at the lower temperature (20°C) displayed more severe clinical signs, leading to early euthanasia (Table S1).Some field studies have shown that PCR detection of O. ophidiicola is more prevalent in early spring than fall season indicating that infection might initiate in winters, but hosts might recover from chronic infection during summers (McKenzie et al., 2018).Higher environmental temperature facilitates host immune activation and infected hosts prefer warmer temperatures (behavioural fever) to increase their immune response and survival during infection, as seen in both amphibians (Boltaña et al., 2013) and snakes (Burns et al., 1996;Godwin et al., 2021).Behavioural thermoregulation via basking has been reported in free-ranging snakes affected by ophidiomycosis (McBride et al., 2015) and is possibly a defence mechanism against infection (Lorch et al., 2016).In cases of chytridiomycosis, amphibian host immune responses influence disease outcome as infected hosts show a greater ability to activate acquired immune responses at higher temperature, whereas the host has to rely on innate immunity and inflammatory pathways to resist infection at lower temperature where the mortality rate is higher (Ellison et al., 2020).Our DEGs between infected snakes at different temperatures suggest that biological processes related to chronic inflammation at the lower temperature might have played a role in the higher mortality in infected snakes in the 20°C group.More research is necessary to better address the role of temperature in snakes' response to ophidiomycosis and to understand how fungal load and host immune response vary in different temperature conditions.

| Limitations of present study
We acknowledge that our conclusions are based on statistical inferences from relatively low numbers of biological replicates for each treatment group, as compared to similar analyses in model systems or lab-raised animals (e.g.Harder et al., 2019;Savage et al., 2020;Sun et al., 2021; Figure 1a, Table S1).However, even though statistical inferences are more accurate with larger replicates, descriptive studies designed to generate hypothesis for future testing with larger samples can be conducted with minimum biological replicates (e.g.n = 3 in Franchini et al., 2017;Videvall et al., 2015) required for statistical tests (Conesa et al., 2016;Schurch et al., 2016;Zhu et al., 2019).
Low replicates per group means that the experiment has less statistical power to identify genes with low effect size on differential expression due to treatment (Conesa et al., 2016).However, genes with high effect sizes due to treatments can be detected; provided appropriate adjustments are made to minimize results due to false positives (Zhu et al., 2019).We minimized false-positive DEGs for each experimental effect (infection and temperature) and corrected for high variance in transcript counts due to small sizes by shrinking fold changes and FDR adjusted p-values to define significance (see Section 2).As such, our results provide a conservative assessment of the impact of O. ophidiicola infection on gene expression in snakes, which is valuable given the lack of information about this important disease that impacts a wide range of species.We acknowledge that the genes we have identified are likely a subset of all DEGs responsible for host response and that a larger sample of individuals would allow us to better account for individual variation in these analyses and allow the identification of genes that have more subtle effects on host response to this disease.

| Potential molecular markers for future ophidiomycosis studies
Our differential gene expression analysis generated a list of candidate protein coding genes that potentially mediate ophidiomycosis response in snakes.From a conservation and disease management perspective, comparing variation in sequence diversity and/or expression profiles of genes likely involved in metabolic and immune pathways that affect host physiology among host species could be a potential avenue to identify vulnerable species (Savage et al., 2020).
However, use of these genetic markers to identify potentially vulnerable species comes with important qualifications.First, inferences about biological pathways from differentially expressed genes were based on functions predicted from orthologs in other reptile species.GO terms are only available for small set of model organisms and the functional impact of many species-specific genes not identified in the genome annotation of the model organism are unknown.Relying on functional information of homologous genes in related species (Anolis carolinensis in this study) does yield a comprehensive inventory of genes whose functional involvement in disease response can be explicitly tested independently.Second, we acknowledge that the physiological response may vary among different host snake species, as has been seen in multiple amphibian and bat species over the course of multiple studies in Bd and WNS, respectively (Hoyt et al., 2021;Zamudio et al., 2020).However, we suspect that the biological mechanisms underlying host response (e.g.immune activation, metabolism, skin regeneration) are similar among different snake species and our curated set of genes (Data S1 and S2) can be used to compare how changes in expression and/ or certain polymorphisms within these genes could be associated with ophidiomycosis response.For example, the link between lower acquired immune gene expression and Bd survival has been observed in many amphibian species (McDonald, Longo, et al., 2020;Savage et al., 2020) and specific genotypes can be associated with disease resistance (McDonald, Longo, et al., 2020;Savage & Zamudio, 2011).

| CON CLUS ION
Ophidiomycosis (snake fungal disease; SFD) is a recently detected and emerging fungal disease and a potential threat to ecologically and phylogenetically diverse snake species.Our transcriptomic analysis of liver, kidney and skin tissues between infected and control groups at different temperatures indicates that even though fun- that included two fixed effects and the interaction between the fixed effects (Figure1a).The two fixed effects were temperature('Temperature'; low [20°C]    vs. high [26°C]) and infection status ('Infection'; O. ophidiicola infected vs. control).The interaction (Infection × Temperature) to cluster highly correlated genes into different modules and to then quantify the relationship of different modules with each other and to the temperature and infection treatments following the workflow described by the authors (https:// horva th.genet ics.ucla.edu/ html/ Coexp ressi onNet work/ Rpack ages/ WGCNA/ Tutor ials/ index.html; last accessed 03/15/2023).We first identified modules of co-expressed genes by calculating pairwise Pearson correlations between each pair of genes with non-zero expression in at least eight of our samples in each tissue separately (N genes_liver = 45,021; N genes_kidney = 47,522;

Finally
Figure S12 for details).

total of 776
Figure2b) and only 17 DEGs in skin (Infection = 15; Temperature = 2; Figure2c).We did not identify any DEGs due to the interaction between infection and temperature in skin.In terms of changes in expression, the majority of genes in liver (457/776; 58.9%) and skin (12/17; 70.6%) were upregulated, whereas majority of the genes in kidney (452/705; 64.1%) were downregulated.To identify DEGs due to the impact of a single fixed effect, we isolated DEGs that did not overlap with the two other fixed effects.These results are shown in Figure 2d-f and Figure S4.In all three tissue types, O. ophidiicola infection alone caused more changes in gene expression than either temperature or the interaction between temperature and infection (liver = 423; kidney = 429; skin = 15; Biological processes that were differentially regulated in liver due to Ophidiomyces ophidiicola infection.Hierarchical clustering tree summarizing the correlation among significant biological process gene ontology terms (GO:BP) for genes that were (a) upregulated or (b) downregulated in liver.Numbers indicate p-values after FDR correction for multiple testing.Each process terms were categorized as subclass of specific parent terms as represented by coloured boxes.Biological processes that were differentially regulated in kidney due to Ophidiomyces ophidiicola infection.Hierarchical clustering tree summarizing the correlation among significant biological process gene ontology terms (GO:BP) for genes that were (a) upregulated or (b) downregulated in kidney.Numbers indicate p-values after FDR correction for multiple testing.Each process terms were categorized as subclass of specific parent terms as represented by coloured boxes.
modules were enriched in DNA replication and DNA mismatch repair KEGG pathways (Figures S10 and S11).Expressed in the nucleus during DNA synthesis and mismatch repair, these pathways are essential F I G U R E 5 Differential gene expression in infected snakes due to difference in temperature.(a) Differentially expressed genes (DEGs) division and proliferation.Overall, our DEG analysis of skin tissues hints towards the activation of the host immune response due to infection-induced stress and rapid skin regeneration through upregulation of DNA replication and DNA repair mechanisms.
gal growth was likely localized to the skin, most differential gene expression occurred in the internal organs.Temperature played an important role as infected snakes maintained at the lower temperature had more severe clinical signs.Our assessment of gene expression suggests that ophidiomycosis may be a systemic disease where chronic fungal infections might lead to the most profound physiological changes.Future research should use diverse -omic approaches to determine how different snake species and populations vary in their response to ophidiomycosis over time.Our reported set of candidate genes and putatively functional loci could be used in future studies to improve our understanding of host-pathogen dynamics to help mitigate ophidiomycosis in wild and captive snakes.AUTH O R CO NTR I B UTI O N S EH, MCA, SM and HLG designed the research; MCA performed controlled exposure experiments; HLG and MCA contributed essential resources; SM performed bioinformatics research and analysed the data; SM and HLG wrote the initial manuscript; MCA and EH revised and edited the final manuscript.
(McKenzie et al., 2018)host biology and likely long-term exposure of the fungus, we suspect that physiological impacts of ophidiomycosis may start internally before clinical signs manifest on the skin(Verant et al., 2014).For example, a field study of O. ophidiicola in 15 wild snake species found that 23.9%-42.3% of snakes tested were qPCR-positive for O. ophidiicola in the spring season without exhibiting any clinical signs(McKenzie et al., 2018).
A potential avenue for future research would be to compare host gene expression between early and late stages of O. ophidiicola infection to understand ophidiomycosis susceptibility.For example,