Temporal dynamics of the bat wing transcriptome: Insight into gene-expression changes that enable protection against pathogen

ABSTRACT Skin acts as a mechanical barrier between the body and its surrounding environment and plays an important role in resistance to pathogens. However, we still know little regarding skin responses to physiological changes, particularly with regard to responses against potential pathogens. We herein executed RNA-seq on the wing of the Rhinolophus ferrumequinum to assess gene-expression variations at four physiological stages: pre-hibernation, hibernation (early-hibernation and late-hibernation), and post-hibernation, as well as the gene-expression patterns of infected and uninfected bats with the Pseudogymnoascus destructans (Pd). Our results showed that a greater number of differentially expressed genes between the more disparate physiological stages. Functional enrichment analysis showed that the down-regulated response pathways in hibernating bats included phosphorus metabolism and immune response, indicating metabolic suppression and decreased whole immune function. We also found up-regulated genes in post-hibernating bats that included C-type lectin receptor signalling, Toll-like receptor signalling pathway, and cell adhesion, suggesting that the immune response and skin integrity of the wing were improved after bats emerged from their hibernation and that this facilitated clearing Pd from the integument. Additionally, we found that the genes involved in cytokine or chemokine activity were up-regulated in late-hibernation compared to early-hibernation and that FOSB regulation of immune cell activation was differentially expressed in bats infected with Pd during late-hibernation, implying that the host’s innate immune function was enhanced during late-hibernation so as to resist pathogenic infection. Our findings highlight the concept that maintenance of intrinsic immunity provides protection against pathogenic infections in highly resistant bats.


Introduction
Hibernation is a seasonal physiological adaptation that allows some mammals to survive in harsh winter circumstances with limited food availability [1]. During hibernation, much of this time is spent in torpor, a physiological stage characterized by inactivity, attenuated body temperature, reduction in heart and metabolic rates, and immune suppression [2][3][4]. While the molecular and genetic bases of hibernation physiology in mammals have been partially evaluated by analysing the differential gene expression in some species such as bears, ground squirrels, primates, and bats [5][6][7][8], the majority of extant studies focused on determining the differentially expressed genes in the brain, heart, and liver -organs involved in carbohydrate and lipid metabolism, detoxification, and molecular transport [9][10][11].
Animal skin, which is the tissue that serves as a portal between the body and external environment, acts as a mechanical barrier in protecting against pathogens. An abundance of evidence now suggests that the skin is actually an active immune organ [12][13][14][15]. Additionally, due to its direct contact with the external surroundings, skin is continuously exposed to large numbers of pathogenic. Thus, animal skin may undergo marked physiological changes that cause alterations in gene expression to combat potential pathogens, such as local inflammation response and the production of cytokines, especially during hibernation when the host's systemic immune function is inhibited. Therefore, exploring gene expression in skin during hibernation is essential to understand how the host copes with potential pathogens and also contributes to understanding the skin's adaptation to environmental changes. However, little is known regarding changes in skin gene expression during hibernation in mammals, particularly with respect to bat species that are at risk of population extinction due to Pseudogymnoascus destructans (Pd) infection.
The emerging infectious disease known as white-nose syndrome (WNS) (also called white-nose disease [WND]) is caused by Pd (previously referred to as Geomyces destructans), and WNS has precipitated a catastrophic decline in the bat fauna of North America [16][17][18]. Pd invades the wing tissues, forming characteristic cup-shaped erosions and ulcerations [19]. Bats awaking from hibernation may also undergo an inflammatory immune-reconstitution syndrome and an acute inflammatory response to cope with Pd infection, potentially leading to their mortality [20]. Intriguingly, bat species are seemingly unequally affected by Pd [21][22][23]. For example, many North American bats have been severely affected, while Chinese bats were also infected with Pd without any mortality. Variations in skin responses to Pd may play a role in the inherently different susceptibilities to Pd between North American and Chinese bat faunas. This variety of responses is similar to that reflected in some previous studies on amphibians where those resistant to Batrachochytrium dendrobatidis exhibited up-regulated genes involved in pathways associated with maintaining structural integrity of skin (e.g. cell adhesion and epidermis development), whereas susceptible amphibians did not [24,25]. Although limited data on skin responses to Pd infection [26][27][28][29][30] as well as on differences in the chemical composition of bat epidermis [31] in susceptible and tolerant bat species have been reported, there are no extant studies on bats that manifest high resistance to Pd in endemic regions. We posit that studying the mechanisms by which highly resistant bats cope with Pd infection would facilitate revealing the adaptive evolutionary processes associated with the occurrence or reduction of fungal infectious-disease severity. This is critical to our understanding of pathogen virulence and bat survival in natural immune states.
In the present study, we used next-generation RNA sequencing on wing tissues from wild greater horseshoe bats (Rhinolophus ferrumequinum) that exhibit high resistance to Pd infection so as to provide a more complete picture of the changes in wing gene expression in hibernating and active bats. We assume that R. ferrumequinum is resistant rather than tolerant due to the fact that resistance refers to host defences that reduce pathogen growth, whereas tolerance refers to host defences that reduce damage experienced by the host without reducing pathogen growth [32]. A previous study showed that although Chinese R. ferrumequinum was similar in colony size and temperatures in hibernation to North American bats, it manifested very low infection intensity and Pd load in China, suggesting host resistance to pathogen [21]. Herein, our aims were to 1) to explore the temporal changes in the bat wing transcriptome and to seek critical responses for coping with potential pathogens, and 2) to compare skin responses to Pd infection during late-hibernation.

Sample collection
We collected 30 R. ferrumequinum adults over a longitudinal time-course that spanned four physiological stages between October of 2020 and May of 2021 from a hibernaculum in Jilin Province of Northeast China where we had detected Pd in previous studies [22,33]. The four physiological stages were prehibernation (five individuals in October), earlyhibernation (five individuals in December), latehibernation (12 individuals in April), and posthibernation (eight individuals in May). We captured active bats with a net (8 m × 3 m) when they flew out of the cave at night, and removed them from the net using sterile latex gloves. After 2 h, we untangled the nets and brought them back to the laboratory to prevent interfering with bat predation. Since only R. ferrumequinum lives in this cave, no other bat species were caught during sampling. We captured hibernating bats from their roosting locations. We swabbed each individual bat five times along its forearm and muzzle using sterile polyester swabs for Pd detection and punched 5-mm biopsy samples from areas of the right and left plagiopatagium (the membranous area between the last digit and the hindlimbs) of the wing. We then stored the samples in RNase-free tubes containing 500 μL of RNAlater (TIANGEN, Beijing, China) and flash-froze the tubes in liquid nitrogen for RNA extraction. Weights were measured and recorded to estimate the relative physical condition of each individual [34], and bats were then released immediately after sampling.

RNA extraction and transcriptome sequencing
We extracted RNA from each punch biopsy sample of bat tissue based on the TRIzol method (TIANGEN, Beijing, China). We quantified and assessed RNA integrity using the Bioanalyzer 2100 System (Agilent Technologies, CA, USA). RIN values greater than 8.0 were confirmed for all samples in this study. RNA sequencing was subsequently executed on the Illumina HiSeq 4000 platform with paired-end 150-bp nucleotides by the Beijing Allwegene Technology Co. Ltd, China.

Mapping analysis and quantification of gene-expression levels
After obtaining raw data, quality control was performed using Trimmomatic v 0.33 according to the following criteria: 1) removing reads with adapter sequences; 2) filtering out reads with uncertain bases (N) greater than 10%; 3) filtering out reads with low quality bases (Q < 20) greater than 50% [37]. We then mapped the clean reads to the reference genome of R. ferrumequinum from Ensembl 100 using STAR v 2.5.2b [38]. Only reads with a perfect match or with one mismatch were further analysed and annotated based on the bat reference genome. We exploited HTSeq v 0.5.4 p3 to count the read numbers mapped to each gene [39] and calculated the values for fragments per kilobase of transcript, per million mapped reads (FPKM) using RSEM v1.3.1 to estimate gene-expression levels [40,41].

Differential expression analysis
DESeq2 v1.32.0 was employed to detect differentially expressed genes (DEGs) for each sample pair [42], and resulted in six sets of DEGs, i.e. early-vs. prehibernation, late-vs. early-hibernation, post-vs. latehibernation, late-vs. pre-hibernation, post-vs. prehibernation, and post-vs. early-hibernation. To explore whether gene expression differed between wing tissues of infected and uninfected Pd bats, we used DESeq 2 to screen for DEGs between the two groups during latehibernation, and P values were corrected based on the Benjamini and Hochberg method [43]. In this study, genes with an adjusted P value <0.05 and |log2 fold-change| > 1 were considered to be differentially expressed. In each pairwise comparison, we defined specific up-and down-regulated genes, with upregulated genes designating expression at higher levels in the former physiological stage than in the latter, while down-regulated genes reflected attenuated expression in the former physiological stage relative to the latter.
To visualize gene-expression patterns across all 30 samples of uninfected bats and bats infected with Pd, PCA and hierarchical clustering with Pearson's correlation were performed within R v. 4.1.0 using dudi.pca() and hclust() functions [44]. ANOSIM (Analysis of Similarity) based on Euclidean distance was implemented using the anosim() function in the vegan package to assess differences between groups. To determine the biological function of DEGs identified in each comparison, functional annotation was carried out by GO and KEGG pathway analyses using the OmicShare tool (https://www.omicshare.com) [45,46]. We achieved multiple-test correction using the Benjamini-Hochberg approach and determined significance with an adjusted P value of <0.05.

Pd-infection status of R. ferrumequinum
While we detected no Pd infection in the bats collected from pre-, early-, and post-hibernation groups, for bats collected during late-hibernation, we detected six bats infected with Pd, with an average fungal load was −5.38, ranging from −6.24 to −4.58 per individual (Ct values ranged from 37.38 to 42.94). We also ascertained that infected and uninfected Pd bats possessed a similar weight, suggesting their similar physical condition (t = 0.81, P = 0.24).

Transcriptomic sequencing and mapping
We obtained 1,378.52 million raw reads from the 30 samples, with an average of 45.95 million (range, 39.11-54.39 million/sample) raw reads per sample (Table S1)

Differential expression at four physiological stages
We conducted six pairwise comparisons for the differential gene-expression analyses across the four physiological stages and identified a total of 8778 DEGs. Our principal component analysis (PCA) based on all DEGs showed four distinct groups, with each corresponding to a distinct physiological stage, and PC1 and PC2 explained 64.15% and 10.24% of the variation, respectively, indicating sizable differences in gene expression among groups (ANOSIM: R = 0.53, P = 0.001; Figure 1a). We discerned similar results from the hierarchical clustering (Figure 1b).
When we statistically analysed the number of DEGs, we uncovered 913, 876, 4790, 2361, 5379, and 5682 DEGs in the six groups of pairwise comparisons, respectively, and the total number of DEGs or specific DEGs was the highest in the four groups of comparisons (two post-hibernation vs. hibernation, post-vs. pre-hibernation, and late-vs. prehibernation groups). This reflected the largest differences in physiological stages, while the number of up-regulated genes was also greater than the number of down-regulated genes (Figure 1c and 1d). Furthermore, the total number of DEGs was also smaller in the comparison between the two groups, with smaller differences in other physiological stages.

GO and KEGG pathway enrichment analyses
GO enrichment analysis showed that a total of 261 GO terms were significantly enriched for the six pairwise comparisons of DEGs, of which up-and downregulated genes were significantly enriched for 180 and 81 GO entries, respectively. We realized that numerous GO terms were enriched by DEGs from the post-hibernation group relative to the other stages and that the up-regulated genes were enriched with more GO terms than were the down-regulated genes. However, fewer terms were significantly enriched for the two hibernation stages vs. pre-hibernation and late-  early-hibernation comparisons ( Table 1).
The DEGs in the hibernating bats relative to the prehibernating bats that were down-regulated in earlyhibernation were involved in metabolism (Figure 2a), with many of the most enriched categories involving phosphorus metabolism (e.g. GO: 0006796, phosphatecontaining compound metabolic process, Padj = 1.87 × 10 − 5). In addition, down-regulated DEGs in latehibernation were not only enriched in some similar GO terms but were also included in immune response (e.g. antigen processing and presentation), adhesion response (e.g. biological adhesion and wound healing), and transport process (Figure 2b).
The functional categories of genes in posthibernation relative to the hibernation groups were enriched in GO terms that included three main categories: immune response, adhesion process, and development process ( Figure 3). Specifically, compared to the two stages of hibernation, the up-regulated genes in the post-hibernation group were not only significantly enriched for similar terms (e.g. the three top GO terms with the smallest Padj values were immune response, immune system process, and regulation of immune response) but were also significantly enriched for other disparate terms (Table S2). Compared with latehibernation, the up-regulated genes in the posthibernation group were significantly enriched for developmentally related terms such as epidermis development and system development, which were not found in the comparison with the early-hibernation group. However, the up-regulated genes in posthibernation relative to early-hibernation groups were significantly enriched in T-cell related terms such as T-helper 17 type immune response and T cell differentiation, which were not found in the comparison with the late-hibernation group. More importantly, for the late-hibernation vs. early-hibernation comparison, we also discovered that many of the DEGs were involved in cytokine and chemokine activities (Table S3).
The results of KEGG enrichment analysis revealed that a total of 99 KEGG pathways were significantly enriched in the six pairwise comparisons of DEGs, of which up-and down-regulated genes were significantly enriched for 59 and 40 KEGG pathways, respectively. When we focused on the pathways that were significantly enriched for DEGs during hibernation and posthibernation (Figure 4), we demonstrated that compared with pre-hibernation, down-regulated genes in the hibernation group were significantly enriched in endocrine system, signal transduction, and immune system pathways. In addition, down-regulated genes were significantly enriched in metabolic pathways between latevs. pre-hibernation groups, including fatty acid metabolism, biosynthesis of unsaturated fatty acids, and cell adhesion-related pathways, while these pathways were not significantly enriched in the early-hibernation group (Table S4).
When we compared them between the two stages of hibernation, the up-regulated genes in the posthibernation group were significantly enriched in pathways related to the immune system, signal transduction, and infection diseases. A further analysis of immune system pathways revealed that up-regulated genes were also significantly enriched in the Toll-like receptor signalling pathway, the C-type lectin receptor signalling pathway, and the cytokine signalling pathway ( Table 2). In addition, two digestion-related pathways (protein digestion and absorption, and vitamin digestion and absorption) were significantly enriched for upregulated genes during post-hibernation (Table S5).

Comparison between uninfected bats and bats infected with Pd in late hibernation
To determine the host response mounted by R. ferrumequinum to Pd during hibernation, we first compared host gene-expression patterns between infected and uninfected bats. Using PCA, we observed that the first two principal components accounted for approximately 66% of the variation among samples (Figure 5a) and that the gene-expression patterns were similar between the two groups (ANOSIM: R = 0.07, P = 0.18). Hierarchical clustering also consistently showed that infected and uninfected samples were mixed together (Figure 5b). Furthermore, compared with uninfected bats, bats infected with Pd showed lower fold-changes in gene expression (Figure 5c). Using DESeq2, we only noted three up-regulated DEGs in Pd-infected bats (Figure 5d); of these, the FOSB gene (an AP-1 transcription factor) was able to regulate immune-cell activation.

Discussion
In the present study, we used the RNA-seq approach to examine the bat wing transcriptome based on changes in gene expression in R. ferrumequinum that exhibited high resistance to Pd. Our data highlighted the dynamic nature of the wing transcriptome, which varied temporally throughout the physiological stages from active to hibernation to active stages. Hibernation is a survival strategy in which animals exhibit reduced body temperature and metabolic rate and experience immune suppression to conserve energy. The reduced metabolic rate of animals during hibernation may lead to changes at the molecular level through regulation of gene expression, just as our results revealed that the genes down-regulated during hibernation compared to pre-hibernation were principally involved in metabolic pathways ( Figure 2). Analogous to our findings, the expression of genes involved in oxidoreductase and glycolytic processes (which are also critical to energy metabolism) was shown to be down-regulated in the liver and brain tissues of R. ferrumequinum during hibernation [7,47], suggesting a generalized inhibition of metabolism in different tissues of bats during this physiological process. While DEGs were not enriched for immune-related GO entries during early-hibernation compared to prehibernation groups in our study, down-regulated genes in late-hibernation were involved in the immune response, implying a gradual and global downregulation of host immune system function with increasing hibernation time. It is worth noting, however, that some up-regulated genes in late-hibernation relative to early-hibernation were enriched in GO terms such as chemokine or cytokine activity (Table S3), including CSF3, CCL19, and CCL20. Actually, although hibernation affects host immune system function, the effects of hibernation on intrinsic and specific immune functions during hibernation are divergent due to the exposure of the host to potentially pathogenic infections. Our results therefore provide a new perspective for the study of changes in immune function during hibernation in mammals. On the basis of our results, it is likely that part of the intrinsic immune function may be increased during late-hibernation of the R. ferrumequinum, and we speculate that this may be due to cope with potential pathogenic infections. Among these pathogenic infections, an example that is particularly profound for bats is Pd, which exhibits seasonal transmission dynamics. Just as we likewise did not detect Pd infection in bats during pre-hibernation, but rather detected it during late-hibernation. Fritze et al [48]. similarly found that the whole immune function of highly tolerant Myotis Myotis was suppressed during hibernation, but that host innate immunity was enhanced. Thus, the appearance of enhanced intrinsic immune function in the wing tissues of the R. ferrumequinum during late-hibernation may be the result of host-pathogen interaction in which bats try to resist Pd.
R.ferrumequinum also appears to activate an immune reaction to remove a pathogenic infection after bats emerge from hibernation and resume normal activity, with subsequent restoration of all components of the immune system. We demonstrated that upregulated genes in the post-hibernation group that were involved in C-type lectin receptor signalling pathway and Toll-like receptor signalling pathway included CLEC7A and TLR4 (Table 2) and that these genes were also up-regulated at both the local and systemic levels in response to Pd in M. lucifugus [26,27]. We thus suggest that after bats emerge from hibernation, responses to pathogenic infection (such as Pd) within their bat wing tissue may be mediated by the innate immune system. We also uncovered genes encoding for pro-inflammatory mediators that characterize the innate immune response (including CSF, IL-23A, and IL-6), with increased transcript levels in posthibernation that potentially mediate the recruitment of monocytes and neutrophils to initiate an adaptive Th17 response and thus provide protection [49]. As a low-level infection is usually contained by an innate immune response, we posit that pathogen clearance may not trigger a robust acute inflammatory response, thus avoiding pathological inflammation.
Although our sampling protocol under native conditions could not eliminate the host's immune response to other pathogens (including bacteria and viruses), our results clearly demonstrated that these host actions were beneficial in their response to Pd infection. We first found that some genes up-regulated in the posthibernation group relative to the hibernation group were enriched in GO terms associated with skin structural integrity, including cell adhesion, response to stimulus, biological adhesion and animal organ development (Table S2). More importantly, compared with late-hibernating animals with Pd infection, the upregulated genes in the post-hibernation group were also enriched in some GO terms, that play important functions in skin structural re-modelling, such as epidermis development and system development. And this was not observed in the post-hibernation group relative to the early-hibernation group without Pd infection. Additionally, KEGG pathway enrichment analysis showed that many up-regulated genes in the posthibernation group were significantly enriched in the cell adhesion molecules (CAMs) pathway. In fact, the disruption of skin function is considered to be a significant cause of host mortality due to Pd infection [50]. Therefore, the increasing expression of genes related to skin structure is significant in host physiological homeostasis and in coping with Pd infection after host emergence from hibernation. This is similar to what has been observed in amphibians, where the genes involved in skin structural integrity and remodelling (e.g. cell-matrix adhesion gene set) were upregulated in resistant species and down-regulated in susceptible species [24,25,51]. Additionally, Pdproduced vitamin B2 and some proteases are the potential virulence factors that erode bat wings, resulting in death [52][53][54]. Our results also showed that differentially expressed up-regulated genes were associated with vitamin digestion and absorption pathway and protein digestion and absorption pathway, which could potentially facilitate the absorption and digestion of vitamin B2 and of proteases produced by Pd. Thus, it may be important for resistant-bat survival to avoid Pdinduced damage to wings during hibernation.
The comparison of gene expression between Pdinfected and uninfected R. ferrumequinum in latehibernation may have revealed potential physiological tolerance as a host-defence mechanism against pathogens. Previous studies have found a range of DEGs associated with immune responses or other responses in susceptible bats following exposure to Pd infection [26,27]. Our results showed that R. ferrumequinum appeared to be relatively unresponsive to Pd infection, which is different from those susceptible bats hibernating in similar ecological conditions [21]. This result suggests that R. ferrumeuqinum exhibits specific patterns that inhibit Pd growth. For example, our recent study showed that the skin microbiota of R. ferrumeuqinum was enriched in particular taxa with antifungal abilities [33]. Additionally, the bats could eradicate the fungus based on different hibernating behaviours, including different arousing times and hibernating length, which should be tested by behavioural experiments and pathological analysis in future study. Despite this, we reported three genes as exhibiting significantly differential expression in infection with Pd ( Figure 5). However (and intriguingly), the FOSB transcript levels among the three DEGs were also augmented in M. lucifugus infected with Pd during both torpor and arousal [27]. Investigators reported in a previous study that FOSB knockdown significantly reduced the mRNA and protein levels of immunerelated chemokines [55], suggesting that bats with high resistance to Pd underwent a transformation so as to produce a trade-off between immune response and energy conservation during hibernation, and this then was able to adjust interleukin-17 signalling to prevent tissue damage [56]. Pd strains in Asia have coexisted longer with bats than those in North America and Europe, and Pd-infected bats during latehibernation showing an approximately 1000-fold lower fungal load than North American species [57]; this therefore lowers risk following the establishment of equilibrium in host-pathogen interaction [23]. Thus, we hypothesize that this balancing mechanism underlying pathogenic resistance in hibernation may have long been selected by Pd pressure, and that it has ultimately been retained throughout co-evolutionary adaptation [58].

Conclusion
Hibernation allows bats to survive in resource-scarce environments by decreasing their metabolism and immune function, but this also makes them vulnerable to psychrophilic pathogens that include Pd. However, we found that a highly resistant species, R. ferrumequinum, actually exhibited increased expression of genes related to intrinsic immune responses and avoids acute inflammatory responses during latehibernation, thus potentially tempering the effects of Pd infections. Our data also suggest that improvements in genes related to immune response and skin integrity will facilitate the clearance of Pd after the bats emerge from hibernation. The genes identified in the present study may provide inspiration for designing effective interventions in susceptible bats. It also provides a basis for elucidating the mechanisms involved in disease susceptibility, tolerance and resistance to other emerging infectious diseases.

Funding
This work was supported by the National Natural Science Foundation of China (grant numbers 32171525, 31961123001 and 31770403).

Data availability statement
All raw sequencing reads generated were deposited in the NCBI Sequence Read Archive (SRA) under submission accession number SRP337443. The datasets generated and analysed during the current study are available in the Figshare r e p o s i t o r y ( a v a i l a b l e a t : h t t p s : / / f i g s h a r e . c o m / s / eb8363e704079baa68fe)