Comparison of spleen transcriptomes of two wild rodent species reveals differences in the immune response against Borrelia afzelii

Abstract Different host species often differ considerably in susceptibility to a given pathogen, but the causes of such differences are rarely known. The natural hosts of the tick‐transmitted bacterium Borrelia afzelii, which is one of causative agents of Lyme borreliosis in humans, include a variety of small mammals like voles and mice. Previous studies have shown that B. afzelii‐infected bank voles (Myodes glareolus) have about ten times higher bacterial load than infected yellow‐necked mice (Apodemus flavicollis), indicating that these two species differ in resistance. In this study, we compared the immune response to B. afzelii infection in these host species by using RNA sequencing to quantify gene expression in spleen. Gene set enrichment analysis (GSEA) showed that several immune pathways were down‐regulated in infected animals in both bank voles and yellow‐necked mice. Moreover, IFNα response was up‐regulated in B. afzelii‐infected yellow‐necked mice, while IL6 signaling and the complement pathway were down‐regulated in infected bank voles; differences in regulation of these three pathways between bank voles and yellow‐necked mice could thus contribute to the difference in resistance to B. afzelii between the species. This study provides knowledge of gene expression induced by a zoonotic pathogen in its natural host, and possible species‐specific regulation of immune responses associated with resistance.

malaria Haemoproteus majoris varies among host species (Huang, Ellis, Jönsson, & Bensch, 2018). Interspecific variation in susceptibility could be a result of variation in resistance and/or tolerance (where resistance refers to variation in the ability to control pathogen replication, while tolerance means variation in ability to limit the damage of a given pathogen load; Råberg, Sim, & Read, 2007). Variation in resistance and/or tolerance can, in turn, be a result of numerous different factors, including differences in anatomy, behavior, and metabolism between host species, but divergence in immune defense likely plays the most important role (Mandl et al., 2014). Yet, there is relatively limited information on differences in immune responses to a given pathogen between host species, in particular closely related ones (but see Poorten &Rosenblum, 2016 andPalesch et al., 2018 for recent exceptions).
Since there are no differences in prevalence of different B. afzelii strains between these host species (Råberg et al., 2017), the difference in resistance is most likely caused by interspecific variation in the immune response. This hypothesis is supported by the higher level of Borrelia-specific antibodies found in yellow-necked mice than in bank voles (Kurtenbach et al., 1994).
Studies of the immune response to Borrelia have primarily been performed in laboratory mice and humans. Briefly, signaling by pattern recognition receptors, including but not limited to toll-like receptor 1 and 2 (TLR1 and TLR2; Oosting, Buffen, van der Meer, Netea, & Joosten, 2014), leads to production and secretion of pro-inflammatory cytokines, like interleukin (IL) 1β, IL2, IL6, tumor necrosis factor α (TNFα), and type I interferons (IFNs) (Jones et al., 2008;Müllegger et al., 2000;Petzke et al., 2016). In addition, the complement system help clear bacteria by promoting phagocytosis through opsonization (Kurtenbach et al., 2002). A strong and effective antibody response against B. burgdorferi infection in mice also occurs (Connolly & Benach, 2005;McKisic & Barthold, 2000). Despite the immune response, bacteria typically disseminate from the tick bite site and colonize internal tissues (Zhong et al., 2019). Studies of both B. burgdorferi and B. afzelii have shown that infections may be chronic (persist for life), although the length of the infection may depend on bacterial strain (Humair, Rais, & Gern, 1999;Kurtenbach et al., 2006).
To investigate the immune response to B. afzelii in natural hosts, in particular how it differs among host species, we here compare the immune response during infections in the bank vole and yellow-necked mouse, in the wild. To this end, we generated spleen transcriptomes of bank voles and yellow-necked mice and compared gene expression between B. afzelii-infected and B. afzelii-uninfected animals.

| Study species and field work
The bank vole and the yellow-necked mouse are small rodents (adult body mass 18-28 g and 24-48 g, respectively) native to Europe and Western Asia (Wilson, Mittermeier, & Lacher, 2017). In Southern Sweden, the two species occur in similar habitats, primarily deciduous woodland, and parasites and pathogens are to a large extent shared between species. For example, at our study site, both species are infested with ticks (Ixodes ricinus) and fleas and infected with the vector-borne bacteria B. afzelii, Candidatus Neoehrlichia mikurensis, and Bartonella spp Hellgren, Andersson, & Råberg, 2011), as well as various helminths (Clough & Råberg, 2014;X. Zhong and L. Råberg, unpublished data).
Animals for the present study were trapped at one locality (Stensoffa) in the Revinge area, 20 km east of Lund, southern Sweden, during five days in August and September 2016 using live traps (Ugglan special, Grahn AB). The habitat at the trapping locality is mixed deciduous forest dominated by beech (Fagus sylvatica).
To minimize variation in gene expression related to for example sex, development or reproductive status, we focused on adult males (adult age indicated by the presence of a scrotum). Based on previous mark-recapture studies of the same population (Tschirren et al., 2013;L. Råberg unpublished data), the vast majority of adult bank voles present in August-September was born in early summer the same year (rather than being over-wintering individuals) and were thus about three months old at the time of sampling for the present study. To avoid bias in gene expression due to short-term variation in environmental conditions, we also made sure to collect equal numbers of the two species on each trapping day.
The spleen is an important immune organ, containing large numbers of immune cells like phagocytes and lymphocytes. The spleen is often enlarged during infection due to proliferation of lymphocytes.
Gene expression in the spleen should therefore reflect the overall activity of the immune system in an individual.
Animals were transported to Stensoffa field station (300 m from the trapping site), euthanized by cervical dislocation, and spleens for RNA extraction were immediately collected and stored in RNAlater RNA Stabilization Reagent (Qiagen GmbH). Ear biopsies (one from the inner part of each ear) for determining infection status were collected and stored in 70% ethanol. Animals were weighed, and sex was determined phenotypically.

| DNA extraction and infection status identification
Total DNA was extracted from ear biopsies following the protocol of Laird et al. (1991). To improve the accuracy of infection status, biopsies from left and right ear from each individual were prepared for DNA extraction separately. The quantity and quality of DNA were assessed by measuring the concentration of DNA and the ratio of UV absorption at 260-280 nm by Nanodrop (Thermo Fisher Scientific).
Infection status was determined by quantitative PCR (qPCR) based on the flaB gene of B. afzelii as described in Råberg (2012). Animals were considered infected if at least one flaB gene copy per nanogram of host DNA was detected in the sample from each ear. Besides B. afzelii, we assessed infection status with respect to helminths (by dissecting the intestines and counting all helminths) and tick loads (Ixodes ricinus; by counting the number of ticks on the ears).
De novo transcriptome assemblies for bank vole and yellow-necked mouse were constructed using Trinity version 2.3.2 (Grabherr et al., 2011) with default settings. To assess the completeness of the transcriptome assemblies, the number of reads that could be mapped back to the contigs was quantified. In addition, the representation of full-length reconstructed protein-coding genes was estimated by BLASTx version 2.6.0 (Camacho et al., 2009) of the assembled contigs against a set of manually reviewed rodent proteins (N = 26,560) downloaded from Swiss-Uniprot (www.unipr ot.org, 2017-05-13).
To reduce the redundancy in the Trinity assembly, contigs were clustered with CD-HIT version 4.6.8 (Fu, Niu, Zhu, Wu, & Li, 2012) using a minimum sequence identity of 95% (-c 0.95). To further reduce the redundancy and to enrich for protein-coding sequence contigs, TransDecoder version 5.0.1 (Haas et al., 2013) was used to identify open reading frames (ORF). For each CD-hit-clustered group, the contig with the longest ORF was saved. To reduce potential assembly errors, reads were mapped back to the filtered contigs using RSEM version 1.3.0 (Li & Dewey, 2011), and contigs with TPM (transcripts per million mapped reads) values <1 were removed.
The assembled contigs were annotated by using BLASTx version 2.6.0 to search against house mouse (Mus musculus) proteins (91,244 transcripts from 22,237 unique protein-coding genes annotated in Ensembl 87, www.ensem bl.org), with an e-value cutoff at 1e-10.
When several contigs had a best hit to the same mouse protein, the one with highest bit score was selected. Finally, the completeness of our annotated assemblies was determined by searching for 4,104 single-copy mammalian orthologues using BUSCO v3 (Simão, Waterhouse, Ioannidis, Kriventseva, & Zdobnov, 2015).

| Mapping approaches
To obtain estimates of gene expression that were comparable between species, we wanted to use similar approaches for both species. An annotated reference genome is available for the bank vole (Lundberg, Zhong, Konrad, Olsen, & Råberg, 2020) but not for the yellow-necked mouse, so mapping reads from each species to its reference genome was not possible. Instead, we considered two other mapping approaches: (a) mapping reads from both voles and mice to the house mouse genome, and (b) mapping reads to de novo transcriptomes for each species.
To test the first approach, trimmed bank vole and yellow-necked mouse reads were mapped to the house mouse genome (GRCm38) with STAR version 2.5.3 (Dobin et al., 2013). Four maximum mismatch numbers (10, 12, 25, and 50) were tested for each species.
Regardless of mismatch number, mapping rates were considerably higher in yellow-necked mouse than in bank vole (not shown). To avoid biased results because of differences in evolutionary distance between our two species and the reference genome, we therefore decided to map reads from each species to its de novo transcriptome instead.
We evaluated the second approach by comparing gene expression in bank voles when mapping to the reference genome and the de novo transcriptome. Reference genome mapping was performed with STAR version 2.5.3 followed by read count quantification using Cufflinks version 2.2.1 (Trapnell et al., 2010). Mapping to de novo transcriptomes and quantification was performed with RSEM version 1.3.0. To compare the results from mapping to the reference genome and the de novo transcriptome, we calculated the correlation of gene expression values (FPKM, Fragments Per Kilobase Million) for each individual vole.

| Differential expression analysis within species
Differential expression between B. afzelii-infected and B. afzelii-uninfected animals was calculated using edgeR (Robinson, McCarthy, & Smyth, 2009) for bank voles and yellow-necked mice separately.
Counts generated from RSEM for each gene were used as input for edgeR, which normalizes counts data using trimmed mean of M-values (TMM). Genes that are lowly expressed (<10 counts detected in more than half of the individuals) were discarded to reduce the number of tests carried out in the differential expression analysis. Genes were considered as significantly differentially expressed between groups when the |log 2 fold change| >1 and the false discovery rate (FDR) value <0.05.

| Gene set enrichment analysis
As a complement to the differential expression analyses, we performed gene set enrichment analyses (GSEA; Subramanian et al., 2005). This method is useful when the data are noisy, as in our case, where we compared individuals sampled in the wild, rather than from a controlled experiment. All genes were ranked based on the association between their expression and the class distinction (infected or uninfected), as measured by the sign of the fold change multiplied by the inverse of the p value (both obtained from differential expression analysis in edgeR). GSEA then uses a priori defined gene sets (such as the genes in a particular pathway), and tests whether the members of a gene set are randomly distributed throughout the ranked list of all expressed genes, or enriched at the top or bottom. If genes in a gene set are enriched at either the top (up-regulated) or bottom (down-regulated), the gene set is considered to be related to the phenotypic difference between samples.
In the software, the GSEAPreranked module was used to calculate a running enrichment score for each gene in our ranked gene list.
We used the Hallmark gene sets from The Molecular Signatures Database v7.0 (Liberzon et al., 2015), which summarize and represent specific well-defined biological states or processes, for running GSEA. Nominal p value < .05 and adjusted p value < .25 were considered as statistically significant (Subramanian et al., 2005). We were particularly interested in gene sets, such as IL2_STAT5 signaling, IL6_JAK_STAT3 signaling, TNFα_via_NFκB signaling, interferon alpha (IFNα) response, and inflammatory response, which could be expected to be highly expressed in the spleen and cover immunological pathways known to be involved in the immune response to Borrelia infection based on studies of laboratory mice and humans (Jones et al., 2008;Kurtenbach et al., 2002;Müllegger et al., 2000;Petzke et al., 2016). IL2 is a cytokine involved in T-cell proliferation and differentiation (Murphy & Weaver, 2016). IL6 and TNFα are pro-inflammatory cytokines, and have a wide spectrum of effects involved in the response to infection by pathogens, including extracellular bacteria like Borrelia. IFNα is primarily known as an anti-viral cytokine, but recent studies found it is also involved in responses against bacteria (Kovarik, Castiglia, Ivin, & Ebner, 2016).

| RNA sequencing and quality control
Thirty-five samples (18 bank voles and 17 yellow-necked mice) were sequenced, resulting in total ~1.44 billion reads with between 29.42 and 54.91 million reads per individual. After removing low-quality bases and adaptors, the data consisted of ~1.35 billion trimmed and paired reads.

| De novo assembled transcriptomes
The assembly of the bank vole transcriptome resulted in 842,299 contigs (presumptive transcripts) grouped in 641,425 clusters (presumptive genes), with a mean length of 899 bp and N50 of 1,799 bp.
The yellow-necked mouse de novo assembled transcriptome consisted of 761,841 contigs and 565,794 clusters with a mean length of 976 bp and N50 of 2,022 bp. When mapping back the reads used in the assembly, both of the assemblies showed a high proportion (>95%) of properly aligned read pairs (Table 1). The bank vole and yellow-necked mouse assemblies also showed similar numbers of matched rodent proteins (19,065 and 18,819) (Table 1).

| Removal of redundant contigs
Using a similarity-based clustering with CD-HIT and selecting the longest contig in each cluster, the number of contigs were reduced to 660,968 and 588,508 in the bank vole and yellow-necked mouse, respectively. Selecting the longest ORF for each of the CD-HIT clusters resulted in 117,706 contigs for bank vole and 114,379 contigs for yellow-necked mouse. Following the removal of contigs with low expression levels (TPM < 1), the transcriptome assemblies consisted of 31,944 contigs for bank vole and 32,238 contigs for yellownecked mouse (Table 1).

| Gene annotation
In the filtered transcriptome assemblies, 19,147 contigs in the bank vole and 19,181 contigs in the yellow-necked mouse could be matched to house mouse proteins. Selecting the contig with the highest bit score to a specific mouse gene resulted in the retention of 13,631 contigs and 13,744 contigs in bank vole and yellow-necked mouse, respectively (Table 1). The two nonredundant transcriptome assemblies showed similar levels of complete and partial mammal single-copy orthologues (Table 1).

| Mapping to de novo assembly versus reference genome
When mapping bank vole reads to the reference genome, 78.05% of reads were mapped and assigned to a unique gene, while the alignment rate to the de novo transcriptome was 51.44%. Gene expression values (log 2 FPKM) estimated from the two approaches were highly correlated within each individual (r: 0.91 -0.93). Therefore, we decided to use de novo transcriptome assemblies as reference for mapping reads in both species.

| Infection status
Of the eighteen bank voles, seven were infected and nine were uninfected, while two had ambiguous infection status (very low infection intensity in one ear). Of the seventeen yellow-necked mice, eight were infected and six were uninfected, while three had ambiguous infection status. Samples with ambiguous infection status were excluded from further analyses. The bacterial load of B. afzelii in bank voles was 10fold higher than in yellow-necked mice (F 1,13 = 20.5, p = .004, Figure 1

| Differentially expressed genes within species
To compare gene expression in infected and uninfected voles and mice, we first performed differential expression analyses in each species separately. After filtering out genes with low expression (<10 counts detected in more than half of the individuals), 12,185 genes of bank vole and 12,282 genes of yellow-necked mouse remained. In bank voles, eight genes were differentially expressed between infected and uninfected individuals. Six of these were up-regulated in B. afzelii-infected bank voles, while two genes were down-regulated ( Figure 2a, Table 2). However, clustering of samples based on the similarity of their gene expression pattern did not group infected and uninfected animals in separate clusters (Figure 2a).
In yellow-necked mice, five genes with significantly differential expression were found; three genes had higher expression in B. afzelii-infected animals than in controls, while two genes had lower expression ( Figure 2b, Table 2). The infected and uninfected yellow-necked mice clustered separately according to gene expression patterns.

| Gene set enrichment analysis
To To better understand the immune response against B. afzelii, we focused on the seven immunological Hallmark gene sets. Of these, six were differentially regulated between infected and uninfected individuals in at least one of the species (all but IFNγ response). Each of these six gene set consists of 87-200 genes involved in a specific immune pathway. In bank vole, 73%-86% of the genes in each gene set were represented in the de novo assembly, while 72%-86% were represented in the yellow-necked mouse assembly. In each gene set, 70%-81% of the genes were shared between bank vole and yellow-necked mouse (Table 3).
Three gene sets (IL2_STAT5 signaling, TNFα_via_NFκB signaling and inflammatory response), showed significant negative enrichment scores in both species, indicating that these three gene sets were down-regulated in B. afzelii-infected individuals in both host species ( Figure 5, Table 4). Two gene sets (complement system and IL6_JAK_STAT3 signaling) were down-regulated in B. afzelii-infected bank voles, while there was no difference in yellow-necked mice.
Finally, the gene set IFNα response was up-regulated in B. afzelii-infected yellow-necked mice, but not in bank voles ( Figure 5, Table 4).
Three of the immunological Hallmark gene sets were down-regulated in both species (IL2_STAT5 signaling, TNFα_via_NFκB signaling and inflammatory response). To investigate to what extent genes within each of these gene sets were regulated in the same way in both species, we performed two analyses. First, we calculated the correlation between the rankings of genes in the two species. For all three gene sets, rankings were moderately positively correlated ( Figure 6a). Secondly, we calculated to what extent "core genes" were shared between species. Core genes are the genes at the leading edge of a gene set, that is the subset of genes that contribute most to the enrichment score of gene set. As shown in Figure 6b, for each gene set, the core genes were shared to 27%-30%.

TA B L E 4 Enrichment of six immunological gene sets between phenotypes (Borrelia afzelii-infected and uninfected)
F I G U R E 6 Similarity of regulation of individual genes in gene sets that were down-regulated in Borrelia afzelii-infected individuals in both bank voles and yellow-necked mice. (a) Correlation of gene ranks between bank voles and yellow-necked mice; (b) Venn diagrams showing the number of overlapping core genes from each gene set between two species. Rs, Spearman's correlation coefficient. *p < .05, **p < .01, ***p < .001 of these pro-inflammatory cytokines in skin and blood is up-regulated during the early phase of Borrelia infection (Bouquet et al., 2016;Jones et al., 2008;Müllegger et al., 2000;Zlotnikov et al., 2017). In contrast, expression of for example TNFα in tissues from late-stage Borrelia infection in mice was reduced (Hodzic, Feng, & Barthold, 2013). Moreover, expression of inflammatory cytokines in response to in vitro stimulation with LPS of blood from late-stage borreliosis patients was suppressed (Diterich, Härter, Hassler, Wendel, & Hartung, 2001). Thus, the down-regulation of pro-inflammatory cytokine expression we observed in naturally infected hosts also occurs during the chronic stage of Borrelia infection in non-natural hosts.
Previous studies have found that the bacterial load is about ten times higher in B. afzelii-infected bank voles than infected yellow-necked mice, indicating that mice are more resistant (Råberg, 2012;Zhong et al., 2019). The same pattern was found in the present study, despite a relatively small data set. The GSEA showed that the main differences between voles and mice in the immune response to B. afzelii is that the complement system and IL6 signaling is down-regulated in infected voles, while IFNα response is up-regulated in infected mice. Thus, differences in regulation of these pathways between bank voles and yellow-necked mice potentially contribute to the difference in resistance to B. afzelii between the species.
The complement system plays an important role in resistance against bacteria, including Borrelia (Kurtenbach et al., 2002). The key proteins involved in the complement system (C2, C3, C4 etc) are mainly synthesized in the liver, but some of these genes are highly expressed also in the spleen, in particular the genes encoding the different components of C1 (at least in humans; GTEx Consortium, 2015). Moreover, the Hallmark complement gene set comprises 200 genes, and thus includes a large number of genes besides the ca. 30 traditional complement proteins. Hence, it is not surprising to find differences in expression of complement genes in the spleen.
IL6 is one of the inflammatory cytokines required for the development and maintenance of the T helper 17 subset, which is important in defense against infections with extracellular bacteria, like Borrelia (Murphy & Weaver, 2016). IL-6 was the predominant cytokine in skin lesions ("erythema migrans") of Lyme borreliosis patients and is strongly induced in Borrelia-stimulated human monocytes in vitro (Cervantes et al., 2011;Salazar et al., 2003). Besides the canonical immune gene sets, there were some other gene sets potentially involved in the immune response that were differently regulated between infected and uninfected individuals.
For example, notch signaling (down-regulated in infection animals in both species) regulates T-cell differentiation together with other signaling pathways, such as IL2 and type I IFN (Amsen, Helbig, & Backer, 2015). Moreover, epithelial-mesenchymal transition (down-regulated in infected animals in both species), which leads to formation of multipotent stromal cells, is a common consequence of inflammation and is involved in for example wound healing (Kalluri & Weinberg, 2009).
To conclude, the present study provides the de novo spleen transcriptome assemblies for two non-model species, bank vole, and yellow-necked mouse and reveals some general patterns as well as differences between host species in the immune response against B. afzelii infection in naturally infected hosts in the wild.
Specifically, several immune pathways were down-regulated during B. afzelii infection in both bank voles and yellow-necked mice, something which might contribute to their role as reservoir hosts. Moreover, differences in regulation of the complement system, IL6 signaling, and IFNα response may contribute to differences in resistance to B. afzelii between bank voles and yellow-necked mice. We would like to stress, though, that controlled infection experiments are required to confirm the present analyses and provide a more detailed picture of the response to B. afzelii in these host species. Vinkler for comments on the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.