Introduction

How does the white-footed deermouse Peromyscus leucopus continue to thrive while sustaining infections with disease agents it serves as reservoir for ⁽⁴⁾? The diverse tickborne pathogens (and diseases) for humans include the extracellular bacterium Borreliella burgdorferi (Lyme disease), the intracellular bacterium Anaplasma phagocytophilum (anaplasmosis), the protozoan Babesia microti (babesiosis), and the Powassan flavivirus (viral encephalitis). Most deermice remain persistently infected but display scant inflammation in affected tissues^{(23, 52, 64)}, and without apparent consequence for fitness ^{(83, 92)}.

A related question—conceivably with the same answer--is what accounts for the two-to-three fold longer life span for P. leucopus than for the house mouse, Mus musculus ^{(48, 79)}? The abundance of P. leucopus across much of North America ^{(34, 66)} and its adaptation to a variety of environments, including urban areas and toxic waste sites ^{(11, 50, 67)}, indicates successful adjustment to changing landscapes and climate. Peromyscus species, including the hantavirus reservoir P. maniculatus ⁽⁶⁵⁾, are more closely related to hamsters and voles in family Cricetidae than to mice and rats of family Muridae ⁽¹⁵⁾.

As a species native to North America, P. leucopus is an advantageous alternative to the Eurasian-origin house mouse for study of natural variation in populations that are readily accessible ^{(9, 53)}. A disadvantage for the study of any Peromyscus species is the limited reagents and genetic tools of the sorts that are applied for mouse studies. As an alternative, we study P. leucopus with a non-reductionist approach that is comparative in design and agnostic in assumptions ⁽³⁾. The genome-wide expression comparison for P. leucopus is with M. musculus and, added here, the brown rat Rattus norvegicus. Given the wide range of pathogenic microbes that deermice tolerate, we use the bacterial endotoxin lipopolysaccharide (LPS) as the primary experimental treatment because the inflammation it elicits within a few hours has features common to different kinds of serious infections, not to mention severe burns and critical injuries ⁽⁹⁶⁾.

We previously reported that a few hours after injection of LPS, P. leucopus and M. musculus had distinguishing profiles of differentially expressed genes (DEG) in the blood, spleen, and liver ⁽³⁾. In brief, the inflammation phenotype of deermice was consistent with an “alternatively activated” or M2-type macrophage polarization phenotype instead of the expected “classically activated” or M1-type polarization phenotype that was observed for M. musculus ⁽⁶⁹⁾. The deermice also differed from mice in displaying evidence of greater neutrophil activation and degranulation after LPS exposure. The potentially damaging action from neutrophil proteases and reactive oxygen species appeared to be mitigated in part in P. leucopus by proteins like secretory leukocyte protease inhibitor (Slpi) and superoxide dismutase 2 (Sod2).

Here, we first address whether the heightened transcription of neutrophil-associated genes in P. leucopus is attributable to differences in numbers of white cells in the blood. To better match for genetic diversity, we substituted outbred M. musculus for the inbred BALB/c mouse of the previous study. We retained the experimental protocol of short-term responses to LPS. This main experiment was supplemented by a study of rats under the similar conditions, by an investigation of a different dose of LPS and duration of exposure in another group of deermice, and by analysis of deermice infected with a bacterium lacking LPS. The focus was on the blood of these animals, not only because the distinctions between species in their transcriptional profiles were nearly as numerous for this specimen as for spleen and liver ⁽³⁾, but also because for ecological and immunological studies of natural populations of Peromyscus species blood is obtainable from captured-released animals without their sacrifice.

The results inform future studies of Peromyscus species, not only with respect to microbial infections and innate immunity, but conceivably also determinants of longevity and resilience in the face of other stressors, such as toxic substances in the environment. The findings pertain as well to the phenomenon of infection tolerance broadly documented in other reservoirs for human disease agents, such as betacoronaviruses and bats ⁽⁵⁷⁾. Less directly, the results provide for insights about maladaptive responses among humans to microbes, from systemic inflammatory response syndrome (SIRS) to post-infection fatigue syndromes.

Results

LPS experiment and hematology studies

Twenty adult animals each for P. leucopus and M. musculus and equally divided between sexes received by intraperitoneal injection either purified E. coli LPS at a dose of 10 µg per g body mass or saline alone (Table 1). Within 2 h LPS-treated animals of both species displayed piloerection and sickness behavior, i.e. reduced activity, hunched posture, and huddling. By the experiment’s termination at 4 h, 8 of 10 M. musculus treated with LPS had tachypnea, while only one of ten LPS-treated P. leucopus displayed this sign of the sepsis state (p = 0.005).

Table 1.

Within a given species there was little difference between LPS-treated and control animals in values for erythrocytes. But overall the deermice had lower mean (95% confidence interval) hematocrit at 42 ^(36-48)%, hemoglobin concentration at 13.8 g/dL (12.1-15.5), and mean corpuscular volume for erythrocytes at 49 fL ^(47-51) than M. musculus with respective values of 56 ^(51-62)%, 16.1 g/dL (14.6-17.7), and 60 fL ^(58-62) (p < 0.01). These hematology values for adult CD-1 M. musculus and LL stock P. leucopus in this study were close to what had been reported for these colony populations ^{(18, 94)}.

In contrast to red blood cells, the mean numbers of white blood cells in the LPS groups in both species were lower than those of control groups (Figure 1). Controls had a mean 4.9 (3.5-6.4) x 10³ white cells per µl among M. musculus and 5.8 (4.2-7.4) x 10³ white cells per µl among P. leucopus (p = 0.41). For the LPS-treated animals the values were 2.1 (1.5-2.7) x 10³ for mice and 3.1 (0.9-5.4) x 10³ for deermice (p = 0.39). However, there was difference between species among LPS-treated animals in the proportions of neutrophils and lymphocytes in the white cell population. The ratios of neutrophils to lymphocytes were 0.25 (0.14-0.45) and 0.20 (0.13-0.31) for control M. musculus and P. leucopus, respectively (p = 0.53). But under the LPS condition. the neutrophil-to-lymphocyte ratio was 0.18 (0.11-0.28) for mice and 0.64 (0.42-0.97) for deermice (p = 0.0006). The regression curves for plots of neutrophils and lymphocytes for LPS-treated and control P. leucopus and LPS-treated M. musculus had similar slopes, but the y-intercept was shifted upwards towards a higher ratio of neutrophils to lymphocytes for blood from the LPS group of deermice. Control group mice and deermice and LPS-treated mice had similar percentages (∼5%) of monocytes in their blood; the mean monocyte percentage rose to 10% in LPS treated deermice (p = 0.12). Eosinophil percentages tended to be higher in deermice at a mean 3.4 (2.1-4.7)% than mice at 1.2 (0.5-1.9)% under either condition (p = 0.004).

Figure 1.

In the P. leucopus experiment with a 10-fold lower dose of LPS and a 12 hour duration, the mean (95% confidence interval) white blood cell count (x 10³) at termination 3.5 (2.5-4.5) in controls and 7.9 (6.0-9.7) in the LPS-treated (p = 0.01). Even with the higher overall white blood cell count the increase white cells was proportionately higher for neutrophils than for lymphocytes, as was seen in the deermice in the higher dose LPS experiment. The ratio of neutrophils-to-lymphocytes was 0.20 (0.07-0.32) in the controls and 0.38 (0.26-0.50) in the LPS-treated (p = 0.10).

The higher neutrophil to lymphocyte ratio in the deermice exposed to LPS was consistent with the greater neutrophil activation noted by transcriptional analysis ⁽³⁾. But many individual genes that constitute this and related gene ontology (GO) terms had transcription levels in the deermice that far exceeded a three-fold difference in neutrophil counts. For some genes the differences were a hundred or more fold, which suggested that the distinctive LPS transcriptional response profile between species was not attributable solely to neutrophil counts.

Genome-wide expression in blood of deermice and mice

We used the respective transcript sets from the reference genomes for P. leucopus and M. musculus for deep coverage RNA-seq with paired-end ∼150 nt reads (Dryad Tables D1 and D2). Principle component analyses (PCA) of the P. leucopus data and M. musculus data revealed that untreated controls had coherent profiles within each species (Figure 2). With the exception of one mouse, the LPS-treated M. musculus were also in a tight PCA cluster. In contrast, the LPS-treated deermice displayed a diversity of genome-wide transcription profiles and limited clustering.

Figure 2.

For both species the number of genes with higher expression with LPS exposure exceeded those with lower or unchanged expression (Dryad Tables D3 and D4). For P. leucopus and M. musculus the mean fold-changes were 1.32 (1.29-1.35) and 1.30 (1.24-1.36), respectively (p = 0.31). For GO term analysis the absolute fold-change criterion was ≥ 2. Because of the ∼3-fold greater number of transcripts for the M. musculus reference set than the P. leucopus reference set, application of the same false-discovery rate (FDR) threshold for both datasets would favor the labeling of transcripts as DEGs in P. leucopus. Accordingly, the FDR p values were arbitrarily set at <5 x 10^-5 for P. leucopus and <3 x 10^-3 for M. musculus to provide approximately the same number of DEGs for P. leucopus (1154 DEGs) and M. musculus (1266 DEGs) for the GO term comparison.

Figure 3 shows the GO terms for the top 20 clusters by ascending p-value for up-regulated and down-regulated in P. leucopus and the corresponding categories for M. musculus. The up-regulated gene profile for P. leucopus featured terms associated with “neutrophil degranulation”, “myeloid leukocyte activation”, “leukocyte migration”, and “response to molecule of bacterial origin”. Other sets of up-regulated genes for the deermice were “negative regulation of cytokine production” and “regulation of reactive oxygen species metabolic process”. None of these were among the top 20 up-regulated clusters for M. musculus. Indeed, “leukocyte activation” and “leukocyte migration” were GO terms for down-regulated DEGs in M. musculus. Distinctive GO terms for up-regulated genes distinguishing mice from deermice were “response to virus”, “response to interferon-beta”, “response to interferon-gamma”, “response to protozoan”, and “type II interferon signaling”.

By arbitrary criterion of 100 for the top DEGs by ascending p value for each species, 24 genes were shared between species (Supplementary Materials Table S1). These included up-regulated Bcl3, Ccl3, Cxcl1, Cxcl2, Cxcl3, Cxcl10, Il1rn, and Sod2. Among the 100 mouse DEGs, 20 were constituents of GO terms “response to virus” or “response to interferon-beta” and only 3 were members of GO term sets “response to molecule of bacterial origin” or “response to lipopolysaccharide”. In contrast, among the top 100 deermouse DEGs, there were only 2 associated with the virus or type 1 interferon GO terms, but 12 were associated with either or both of the bacterial molecule GO terms.

We confirmed the sex identification for each sample with sex-specific transcripts of Xist for females and Ddx3y for males ⁽³⁾. For female and male P. leucopus there were 5012 transcripts out of 54,466 in the reference set for which there were TPM values of ≥10 in at least one animal in each of the sexes under either condition (Dryad Table D5). The comparison of females to males by fold-changes between LPS-treated and control animals revealed transcripts that were differentially expressed between sexes under LPS treatment (Table S2 and Figure 4). Some were down-regulated in one sex while unchanged in expression in the opposite sex. Of note in this category were different isoforms or variants of Lilra6 (leukocyte immunoglobulin-like receptor, subfamily A, member 6), one of a family of orphan receptors of myeloid cells ⁽⁸⁾. The opposite case was exemplified by the Dnajc15 and Hspa8 genes for two chaperones: DnaJ heat shock protein family (Hsp40) member C15 and heat shock protein 8, respectively. These were substantially lower in transcription in the LPS-treated females than in untreated animals, but little changed in LPS-treated males. Coordinates for some other genes, e.g. Saa5 and Cxcl2, fell outside the prediction limits at the extreme end of up-regulation, but their vectors were within 20-25° of each other. While these and other sex-associated differences merit attention for future studies, overall they were not of sufficient number or magnitude in our view to warrant division by sex for the subsequent analyses, which had the aim of identifying differences applicable for both females and males.

Figure 3.

Figure 4.

Targeted RNA seq analysis

The emerging picture was of P. leucopus generally responding to LPS exposure as if infected with an extracellular bacterial pathogen, including with activated neutrophils. While M. musculus animals of both sexes shared with P. leucopus some features of an antibacterial response, they also displayed type 1 and type 2 interferon type response profiles associated with infections with viruses and intracellular bacteria and parasites.

Going forward, the challenge for a cross-species RNA-seq was commensurability between annotated transcripts of reference sets. Orthologous genes can be identified, but mRNA isoforms and their 5’ and 3’ untranslated regions may not fully correspond. Accordingly, we limited targeted RNA-seq to protein coding sequences of mRNAs for the corresponding sets of P. leucopus and M. musculus sequences.

The 113 mRNA coding sequences, which are listed in Methods, were drawn from the identified DEGs for P. leucopus and M. musculus from the genome-wide RNA-seq. For cross-species normalization we first evaluated three methods: (1) normalization using the ratio of mean total reads for all samples to total reads for a given sample, (2) the ratio of reads mapping to a given target gene (i.e. the numerator) to the reads to the mitochondrial 12S rRNA gene (i.e. the denominator), or (3) when the denominator instead was the myeloid cell marker CD45 (protein tyrosine phosphatase, receptor type C) encoded by the Ptprc gene. In humans, mice, and hamsters, Ptprc is expressed by nucleated hematopoietic cells, and CD45 is commonly used as a white cell marker for flow cytometry ⁽⁸²⁾. The coefficients of determination (R²) between comparison pairs (e.g. normalized total reads vs. Ptprc ratio) within a species were ≥ 0.95 (Supplementary Materials Figure S1 and Table S3). There was also little difference between the choice of Ptprc or 12S rRNA transcripts as denominator with respect to cross-species comparisons of LPS-treated to control fold changes. Given the widespread adoption of CD45 for flow cytometry, we chose Ptprc as denominator and as an adjustment for white cell numbers in the samples. Pearson correlation between log-transformed total white blood cell counts and normalized reads for Ptprc across 40 animals representing both species, sexes, and treatments was 0.40 (p = 0.01).

Figure 5 comprises plots of the log-transformed mean ratios for the 10 P. leucopus controls and 10 M. musculus controls and for the 10 P. leucopus and 10 M. musculus treated with LPS (Table S4). For untreated animals (left panel) there was high correlation and a regression coefficient of ∼1 between the paired data for deermice and mice. The mitochondrial cytochrome oxidase 1 gene (MT-Co1) and S100a9, a subunit of calprotectin, were comparably transcribed. But, there were other coding sequences that stood out for either their greater or lesser transcription in untreated deermice than mice. Two examples of greater expression were Arg1 and MX dynaminin-like GTPase 2 (Mx2), an ISG, while two examples of lesser expression were matrix metalloprotease 8 (Mmp8) and Slpi. There was low to undetectable transcription of Nos2 and interferon-gamma (Ifng) in the blood of controls of both species.

Figure 5.

For the LPS-treated animals (right panel Figure 5) there was, as expected for this selected set, higher expression of the majority genes and greater heterogeneity among P. leucopus and M. musculus animals in their responses for represented genes. In contrast to the findings with controls, Ifng and Nos2 had higher transcription in treated mice. In deermice the magnitude of difference in the transcription between controls and LPS-treated was less. A comparatively restrained transcriptional response in deermice was also noted for Mx2. On the other hand, there were greater fold-change from baseline in P. leucopus than in M. musculus for interleukin-1 beta (Il1b), Mmp8, Slpi, and S100a9.

Table 2 lists all the selected targets with the means and confidence intervals for the normalized values for controls and LPS-treated M. musculus and controls and LPS-treated P. leucopus (Table S4). Box plots for a selected 54 of these targets organized by functional characteristics are provided in Figure S2. The fold-changes within each species and between treatments across species are given. The final column is the ratio of the fold-change between LPS to control in P. leucopus to the corresponding value for M. musculus. This along with the derived heat-map of these ratios, presented in the second column, indicates the genes for which there was little difference between species in their responses to LPS--either up or down--as well as those that were comparatively greater or lesser in one species or the other. Several of these genes are considered in other specific contexts below. Of note here are the places of Nos2 and Ifng at the bottom of the table, and Il1b near the top at position 20.

Table 2.

“Alternatively-activated” macrophages and “nonclassical” monocytes in P. leucopus

While we could not type single cells using protein markers, we could assess relative transcription of established indicators of different white cell subpopulations in whole blood. The present study, which incorporated outbred M. musculus instead of an inbred strain, confirmed the previous finding of differences in Nos2, the gene for inducible nitric oxide synthase and Arg1, the gene for arginase 1, expression between M. musculus and P. leucopus (Figure 5; Table 2). Results similar to the RNA-seq findings were obtained with specific RT-qPCR assays for Nos2 and Arg1 transcripts for P. musculus and M. musculus (Table 3).

Table 3.

Low transcription of Nos2 in both in controls and LPS-treated P. leucopus and an increase in Arg1 with LPS was also observed in another experiment for present study where the dose of LPS was 1 µg/g body mass instead of 10 µg/g and the interval between injection and assessment was 12 h instead of 4 h (Table 4).

Table 4.

In addition to the differences in Nos2 and Arg1 expression for typing macrophage and monocyte subpopulations, there are also the relative expressions of three other pairs of genes: (1) IL12 and IL10, where a lower IL12/IL10 ratio is more characteristic of alternatively activated or M2 type ^{(68, 90)}; (2) the proto-oncogene kinases Akt1 and Akt2, where the associations are Akt1 with M2-type and Akt2 with M1-type macrophages ^{(1, 91)}; and (3) CD14 and CD16 (low affinity immunoglobulin gamma Fc region receptor III; Fcgr3), where low expression of CD14 and high expression of CD16/Fcgr3 is associated with “non-classical” monocytes ⁽⁷⁰⁾. There is evidence that nonclassical monocytes can change to M2-type macrophages ⁽⁴²⁾.

These four relationships, which are presented as log-transformed ratios of Nos2/Arg1, IL12/IL10, Akt1/Akt2, and CD14/CD16, are shown in Figure 6. We confirmed the difference between P. leucopus and M. musculus in the ratios of Nos2/Arg1 and IL12/IL10 (3) with outbred mice and normalization for white cells. In both species the Akt1/Akt2 ratio declined in LPS-treated animals, but for P. leucopus the ratio remained > 1.0 even among LPS-treated animals, while in the blood of M. musculus the ratio was < 1.0 at baseline and declined further in the LPS-treated animals.

Figure 6.

An ortholog of Ly6C ⁽¹³⁾, a protein used for typing mouse monocytes and other white cells, has not been identified in Peromyscus or other Cricetidae family members. Therefore, for this study the comparison with Cd14 is with Cd16 or Fcgr3, which deermice and other cricetines do have. In mice the Cd14/Cd16 ratio increased from baseline in the LPS group. In the deermice the ratio in control animals was midway between the two groups of mice but there was a marked decrease in the LPS-treated deermice (Figure 6). This was not associated with a fall in the absolute numbers or percentages of monocytes in the blood of these animals.

Taken together, the Nos2-Arg1, Il12-Il10, Akt1-Akt2, and CD14-CD16 relationships document a disposition toward alternatively-activated macrophages and nonclassical monocytes in P. leucopus both before and after exposure to LPS. This contrasts with profiles consistent with a predominance of classically-activated macrophages and classical monocytes in mice.

Interferon-gamma and interleukin-1 beta dichotomy between deermice and murids

For mice the Ifng transcript was one of the top ranked DEGs by both fold-change and adjusted p value by genome-wide RNA-seq (Table 2). In contrast, for P. leucopus Ifng was far down the list, and the comparably ranked DEG instead was Il1b. This inversion of relationships between two pro-inflammatory cytokines was confirmed by analysis of the individual animals of both species (Figure 7). There was little or no detectable transcription of Ifng in the blood of deermice in which Il1b expression was high. There was also scant to no transcription of Ifng in the blood of P. leucopus 12 h after injection of LPS (Table 4).

Figure 7.

The up-regulation of Ifng within 4 hours of exposure to LPS was not limited to the species M. musculus. In an experiment with the rat R. norvegicus, we used two different LPS doses (5 µg/g and 20 µg/g), but the same 4 h endpoint and whole blood as the sample. Both groups of LPS-treated rats had lowered total white blood cells and, like the mice, lower neutrophil-to-lymphocyte ratios compared to controls (Table 5). There were also elevations of interferon-gamma, interleukin-6 and interleukin-10 proteins from undetectable levels in the blood of the treated rats (Dryad Table D5). The values for rats receiving 5 µg/g or 20 µg/g doses were similar (Figure S3), so these groups were combined. By targeted RNA-seq there were 24x fold-changes between the LPS-treated rats and control rats for Ifng and Nos2 but only ∼3x fold-change for Il1b (Table 5).

Table 5.

We asked why the Infg response observed in CD-1 mice and rats here was not as pronounced in BALB/c mice (3). Accordingly, we used the RNA-seq reads from the prior study in combination with the reads of the present study and carried out targeted RNA-seq (Figure S4). The BALB/c inbred mice had, like the CD-1 mice, modest elevations of Il1b transcription. Ifng expression was also elevated in the BALB/c animals but not to the degree noted in CD-1 mice or rats. One explanation is an inherent difference of BALB/c mice from other strains in their lower interferon-gamma response _{to LPS} ^{(47, 84)}.

Interferon-gamma and inducible nitric oxide synthase

Interferon-gamma is a determinant of Nos2 expression ^{(55, 81)}. So, the scant transcription of Ifng in P. leucopus conceivably accounted for the low expression of Nos2 in that species. The analysis shown in upper left panel of Figure 8 shows a tight correlation between the levels of transcription of Ifng and Nos2 for both species and both experimental conditions (Figure S5). A significant correlation was also observed for the combined set of animals between the ratios of Nos2 to Arg1 and Ifng to Il1b (upper right panel), an indication of co-variation between Ifng expression and macrophage polarization.

Figure 8.

The plausible sources of Ifng mRNA in whole blood are T-cells, Natural Killer cells, and Type 1 Innate Lymphoid Cells ⁽⁷⁵⁾. A DEG for M. musculus by both genome-wide and targeted RNA-seq (Table 2; Table S1) was Cd69, a C-type lectin protein and an early activation antigen for these cells(38). In P. leucopus transcription Cd69 occurred in the blood of control P. leucopus, but it was the same or only marginally different for the LPS-treated animals (lower left panel). In contrast, in M. musculus the baseline transcription of Cd69 was below that of P. leucopus, but in the LPS-treated mice it was many fold higher. In mice transcripts for Cd69 correlated tightly with Ifng transcription, while in the deermice there was little correlation between Cd69 and Ifng expression at those low levels (lower right panel).

The findings are consistent with CD69-positive cells being a source of Ifng in mice. Cd69 transcription was comparatively higher in control deermice than in control mice, so we presume that deermice have CD69-positive cells at baseline. One explanation then for the comparatively few Ifng transcripts in the deermice after LPS is a diminished responsiveness of these cells. Tlr4 expression increased ∼3 -fold more in P. leucopus than in M. musculus after LPS (Table 2), but the magnitude of the decline in expression of Cd14 in deermice than mice was even greater. CD14 is required for LPS-stimulated signaling through surface TLR4 ⁽⁵⁸⁾, and, as such, its decreased availability for this signaling pathway is a possible explanation for the moderated response to LPS in P. leucopus.

Interferon-stimulated genes and RIG-I-like receptors

As noted, GO terms differentiating mice from deermice included “response to interferon-beta” and “response to virus” (Figure 3). There was also the example of Mx2, an ISG with anti-viral activity on its own, that showed a greater fold-change from baseline in mice than in deermice (Figure 5). Five other ISGs--guanylate binding protein 4 (Gbp4), interferon-induced protein with tetratricopeptide repeat (Ifit1), interferon regulatory factor 7 (Irf7), ubiquitin-type modifier ISG15 (Isg15), and 2’-5’ oligoadenylate synthase 1A (Oas1a)—had higher transcription in all the LPS-treated animals. But the magnitude of fold change was less in the deermice, ranging from 6-25% of what it was in the LPS group of mice (Table 2; Figure S6).

The up-regulation of these ISGs was evidence of an interferon effect, but transcripts for interferon-1 beta (Ifnb) or -alpha (Ifna) themselves were scarcely detectable in deermice or mice in the blood under either condition. We then considered pattern recognition receptors (PRR) that might be part of a signaling pathway leading to ISG expression. Among the DEGs from the genome-wide analyses were four cytoplasmic PRRs: (1) Rigi (formerly Ddx58), which encodes the RNA helicase retinoic acid-inducible I (RIG-I); (2) Ifih1, which encodes interferon induced with helicase C domain 1, also known as MDA5 and a RIG-I-like receptor; (3) Dhx58, also known LGP2 and another RIG-I-like receptor; and (4) cGAS (cyclic GMP-AMP synthase), which is part of the cGAS-STING sensing pathway (Table S4).

All four cytoplasmic PRRs were upregulated in the blood of LPS-treated mice and deermice (Table 2). But, again, for each of them the magnitude of fold change was less by 50-90% in treated P. leucopus than in M. musculus. The coefficients of determination for the 6 ISGs and the 4 PRRs are provided in Figure 9. For most of the pairs there was evidence of covariation across all 40 animals. When the correlation was low across all the data, e.g. between the ISG Mx2 and the PRR Rigi or the ISGs Mx2 and Gbp4, it was high within a species.

Figure 9.

These findings were evidence that pathways in P. leucopus for PRR signaling and ISG expression functioned similarly to those in M. musculus but differed under these experimental conditions in magnitude of the changes, being more moderate in the deermice.

Endogenous retroviruses in deermice, mice, and rats after LPS exposure

The six ISGs are nonexclusive consequences of activity of type 1 interferons. What we could document was the association of transcription of the gene for the cytoplasmic PPRs, including RIG-I, and the ISGs in both species, as well as the distinction between deermice in the magnitude of the responses of both PRRs and ISGs. These findings led us to ask could be a pathogen-associated molecular pattern (PAMP) for signaling pathways leading to expression of type 1 interferons.

One of these is endogenous retroviruses (ERV). The activity of these diverse, abundant, and pervasive elements have been recognized as one of the drivers of innate immune responses to a microbe ^{(41, 51, 76)}. Our attention was drawn to ERVs by finding in the genome-wide RNA-seq of LPS-treated and control rats. Two of the three highest scoring DEGs by FDR p value and fold-change were a gag-pol polyprotein of a leukemia virus with 131x fold-change from controls and a mouse leukmia virus (MLV) envelope (Env) protein with 62x fold-change (Dryad Table D5).

We returned to the mouse and deermouse data. There were four MLV or other ERV Env proteins among the 1266 genome-wide RNA-seq DEGs for M. musculus. But, there was no ERV Env protein identified as such among the 1154 DEGs identified for P. leucopus (Table S1). One possible explanation for the difference was an incomplete annotation of the P. leucopus genome. We took three approaches to rectify this. The first was to examine the DEGs for P. leucopus that encoded a polypeptide ≥ 200 amino acids and was annotated for the genome as “uncharacterized”. A search with these candidates of both the virus and rodent proteins databases identified two that were homologous with gag-pol polyproteins of ERVs, mainly leukemia viruses, of mammals.

For a second approach we carried out a de novo transcript assembly of mRNA reads from blood of LPS-treated and control P. leucopus and used the resultant contigs as the reference set for RNA-seq analysis. This identified two contigs that were measurably transcribed in the blood, differentially expressed between conditions, and were homologous to ERV sequences. One was revealed to be an Env protein that was identical to a P. leucopus protein annotated as “MLV-related proviral Env protein” (XP_037065362). The second was a gag-pol protein. The latter was near-identical to the gag-pol protein identified by the first approach.

The third approach was to scan the P. leucopus genome for nonredundant sequences, defined as <95% identity, that were homologous with ERV gag-pol sequences, which are not typically annotated because of masking for repetitive sequences. This analysis yielded 615 unique sequences. These were used in turn as a reference set for RNA-seq. There were 4 sequences that met the criterion of FDR p value <0.01. Three were transcribed at 5-to 40-fold higher levels in LPS-treated deermice than in controls. But all three, as well as the fourth, a down-regulated DEG, were ERV relics with truncations, frame shifts, and in-frame stop codons. These were assessed as non-coding RNAs and not further pursued in this study.

To represent P. leucopus in a targeted RNA-seq comparison with mice and rats we settled on the Env protein and gag-pol coding sequences identified present in the blood mRNA and as DEGs. Representing M. musculus were highest ranked MLV Env and gag-pol protein sequences among the DEGs. For rats we chose Env and gag-pol proteins that were second and third ranked DEGs identified in the genome-wide RNA-seq. Because of length differences for the coding sequences, the unit used for cross-species analysis was reads per kilobase before normalization for Ptprc transcription.

The left panel of Figure 10 shows the striking transcriptional fold-change in LPS-treated rats of both Env and gag-pol sequences over controls. Of lesser magnitude but no less significant was the fold-change observed M. musculus for both Env and gag-pol sequences. In both mice and rats Env and gag-pol read values were highly correlated across conditions. In contrast, in P. leucopus the magnitudes of fold-change upwards for gag-pol was less than in mice or rats, and transcription of the Env protein sequence was actually lower in LPS-treated animals than in controls. While there was a tight association between Env protein and the PRR Rigi transcription in the M. musculus, this was not observed in P. leucopus. Rigi transcription was moderately higher at the time that Env protein’s transcription was lower in the LPS group.

Figure 10.

Borrelia hermsii infection of P. leucopus

The phenomena reported so far were consequences of exposures to a particular PAMP--bacterial lipopolysaccharide with its hallmark lipid A moiety--recognized by a particular PRR, TLR4. While the focus was primarily on events downstream from that initial signaling, we asked in a concluding study whether the profile observed in P. leucopus applied in circumstances when the PAMP or PAMPs did not include LPS. This question is germane, given P. leucopus’ role as a natural host for B. burgdorferi. This organism and other members of the spirochete family Borreliaceae do not have LPS ^{(5, 88)}, but they have abundant lipoproteins, which are agonists for TLR2 in a heterodimer with TLR1 ⁽⁸⁰⁾. B. burgdorferi is transiently blood-borne at low densities in P. leucopus, but in its life cycle B. burgdorferi is mainly tissue-associated in vertebrate hosts ⁽⁶⁾. We previously observed that the blood of B. burgdorferi-infected P. leucopus manifested few DEGs in comparison to skin ⁽⁵²⁾. More comparable to the LPS experimental model is infection of P. leucopus with a relapsing fever Borrelia species, which commonly achieve high densities in the blood. P. leucopus is a reservoir for Borrelia miyamotoi, which causes hard tick-borne relapsing fever ⁽⁶⁾, and the related P. maniculatus is a natural host for the soft tick-borne relapsing fever agent B. hermsii ⁽⁴³⁾.

Accordingly, we used blood RNA-seq reads, which were taken from a prior study of B. hermsii infection of P. leucopus ⁽³⁾, for targeted analysis with the same reference set employed for the LPS analyses (Table 6 and Table S5). The blood samples were taken from infected and uninfected animals on day 5, when bacteremia was at its peak, as documented by microscopy of the blood, qPCR of the spleen, and transcripts of a B. hermsii plasmid in the RNA extracts of the blood. As expected for B. hermsii infection ⁽²⁴⁾, the spleen was enlarged in infected animals.

Table 6.

Similarities in the profiles for the LPS-treated and B. hermsii-infected deermice were as follows: (1) low levels of transcription of Nos2 and Ifng that contrasted with the high levels for Arg1 and Il1b in the same animals, (2) maintenance of Akt1/Akt2 > 1.0 under both conditions, (3) reduction of the Cd14/Cd16 ratio, (4) decreased transcription of Cd69, and (5) stable, low transcription of ERV Env and gag-pol loci with only marginal increases in transcription of ISGs and RIG-I-like receptors. Other equivalences under the two experimental conditions included increases in expression of genes for superoxide dismutase 2, low-affinity Fc gamma receptors, and secretory leukocyte protease inhibitor. Thus, the responses that distinguish deermice are not confined to the singular case of LPS as the elicitor.

Discussion

Study limitations

The approach was forward and unbiased, looking for differences between species broadly across their transcriptomes. The findings lead to hypotheses, but reverse genetics in service of that testing was not applied here. In selective cases we could point to supporting evidence in the literature on M. musculus and the phenotypes of relevant gene knockouts, but there are no such resources for Peromyscus as yet. The resource constraint also applies to the availability of antibodies for use with Peromyscus for immunoassays for specific proteins, e.g. interferon-gamma, in serum, or for cell markers, e.g. CD69, for flow cytometry of white blood cells.

While a strength of the study was use of an outbred population of M. musculus to approximate the genetic diversity of the P. leucopus in the study, this meant that some genes of potential relevance might have gone undetected, i.e. from type II error. The variances for a sample of genetically diverse outbred animals, like the LL stock of P. leucopus ^{(52, 53)}, would be expected to be greater than for the same sized sample of inbred animals. For some traits, especially ones that are complex or under balancing selection, even sample sizes of 10 in each group may not have provided sufficient power for discrimination between deermice and mice. For the same reason differences between sexes of a species in their responses might have been undetected. The interpretations applied to mixed-sex groups of deermice and mice. Expression strongly associated with female or male sex could have yielded an average fold change for the whole group that fell below the screen’s threshold.

The parameters for the experiment of LPS dose, the route, and duration of experiment each might have had different values under another design. Those particular choices were based on past studies of deermice and mice ^{(3, 49)}. In another experiment we found that with doses twice or half those given the deermice the responses by rats to the different doses were indistinguishable by hematology, cytokine assays, and RNA-seq. Thus, there seems to be some latitude in the dose and still achieving replication. We obtained similar results for P. leucopus when we looked at a replicate of the experiment with the same conditions ⁽³⁾, when the dose was lower and duration lengthened to 12 h (this study). The analysis here of the B. hermsii infection experiment also indicated that the phenomenon observed in P. leucopus was not limited to a TLR4 agonist.

While the rodents in these experiments were housed in the same facility and ate the same diet, we cannot exclude inherent differences in gastrointestinal microbiota between species and individual outbred animals as co-variables for the experimental outcomes. We reported differences between the LL stock P. leucopus and BALB/c M. musculus of the same age and diet in their microbiomes by metagenomic analysis and microbiologic means ⁽⁶¹⁾. This included a commensal Tritrichomonas sp. in P. leucopus but not in the M. musculus in the study. The presence of these protozoa affects innate and adaptive immune responses in the gastrointestinal tract ^{(19, 29)}, but it is not clear whether there are systemic consequences of colonization by this flagellate.

LPS, ERVs, and interferons

The results confirm previous reports of heightened transcription of ERV sequences in mice or mouse cells after exposure to LPS ^{(35, 44, 87)}. Here we add the example of the rat. The LPS was administered in solution and not by means of membrane vesicles. The sensing PRR presumably was surface-displayed, membrane-anchored TLR4 ⁽⁵⁸⁾. It follows that a second, indirect of LPS on the mouse is through its provocation of increased ERV transcription intracellularly. ERV-origin RNA, cDNA and/or protein would then be recognized by a cytoplasmic PRR. RIG-I was one associated with ERV transcription in this study. Kong et al. reported that LPS stimulated expression of Rigi in a mouse macrophages but did not investigate ERVs for an intermediary function in this phenomenon ⁽⁴⁶⁾. As was demonstrated for LINE type retrotransposons in human fibroblasts, intracellular PRR signaling can trigger a type 1 interferon response ⁽²⁵⁾. The combination of these two signaling events, i.e. one through surface TLR4 by LPS itself and another through intracellular PPR(s) by to-be-defined ERV products, manifested in mice and rats as a response profile that had features of both a response to a virus with type 1 interferon and ISGs and a response to a bacterial PAMP like LPS with acute phase reactants such as calprotectin and serum amyloid.

This or a similar phenomenon has been observed under other circumstances. In humans there was heightened transcription of retrotransposons in patients with septic shock ⁽⁶³⁾, as well as in peripheral blood mononuclear cells from human subjects experimentally injected with LPS ⁽⁷⁴⁾. Bacteria like Staphylococcus epidermidis that express TLR2 agonists, such as lipoteichoic acid, promoted expression of ERVs, which in turn modulated host immune responses ⁽⁵¹⁾. A synthetic analog of a B. burgdorferi lipoprotein activated human monocytic cells and promoted replication of the latent HIV virus in cells that were persistently infected ⁽⁷²⁾.

P. leucopus does not fit well with this model. Instead of the prominent interferon-gamma response observed in mice and rats, there were prominent responses of interleukin-1 beta and genes associated with neutrophil activation. Instead of the much heightened expression of ISGs, like Mx2 and Isg15, in mice treated with LPS, the deermice under the same condition had a more subdued ISG transcription profile. Instead of increased expression of ERV Env protein sequences in blood of mice and rats treated with LPS, there was decreased transcription of the homologous MLV Env gene in like-treated P. leucopus.

This suppression in the deermice may be attributable to defensive adaptations of Peromyscus to repeated invasions of endogenous retroviruses, as Gozashti et al. has proposed for P. maniculatus ⁽³²⁾. This includes expanding the repertoire of silencing mechanisms, such as Kruppel-associated box (KRAB) domain-containing zinc finger proteins ⁽⁹⁷⁾. Like P. maniculatus, P. leucopus has an abundance of Long Terminal Repeat retrotransposons, several named for their endogenous retrovirus heritages ⁽⁵²⁾. Our initial analysis of the P. leucopus genome reported a depletion of KRAB domains compared to Muridae ⁽⁵²⁾. But a subsequent annotation round identified several genes for KRAB domain zinc finger proteins in P. leucopus, including Zfp809 (XP_006982432), which initiates ERV silencing ⁽⁹⁵⁾, and Zfp997 (XP_037067826), which suppresses ERV expression ⁽⁸⁹⁾. Another possible adaptation in P. leucopus is the higher baseline expression of some ISGs as noted here (Figure 9 and Figure S6).

Reducing differences between P. leucopus and murids M. musculus and R. norvegicus to a single attribute, such as the documented inactivation of the Fcgr1 gene in P. leucopus ⁽⁷⁾, may be fruitless. But the feature that may best distinguish the deermouse from the mouse and rat is its predominantly anti-inflammatory quality. This characteristic likely has a complex, polygenic basis, with environmental (including microbiota) and epigenetic influences. An individual’s placement is on a spectrum or, more likely, a landscape rather than in one or another binary or Mendelian category.

One argument against a purely anti-inflammatory characterization is the greater neutrophil numbers and activity in P. leucopus compared to M. musculus in the LPS experiment. The neutrophil activation, migration, and phagocytosis would be appropriate early defenses against a pyogenic pathogen. But if not contained, they bring local and systemic risks for the host. This damage would not likely be from nitric oxide and reactive nitrogen species, given the minimal Nos2 transcription. But deermice showed heightened expression of proteases, such as Mmp8, enzymes for reactive oxygen species, such as NADPH oxidase 1 (Nox1), and facilitators of neutrophil extracellular traps, such as PAD4 (Padi4) (Table 2). We had previously identified possible mitigators, such as the protease inhibitor Slpi and superoxide dismutase 2 ⁽³⁾. These findings were replicated here. The topic of neutrophil activation and these and other possible counters is considered in more detail elsewhere.

An anti-inflammatory disposition but at what cost?

An assignment of infection tolerance to a host and pathogen pairing assumes sufficient immunity against the microbe to keep it in check if elimination fails. P. leucopus and P. maniculatus, are in this sense “immunocompetent” with respect to the microbes they host and with which they may share a long history ⁽⁴⁰⁾. Yet, has this balance of resistance and tolerance for certain host-associated microbes been achieved in a trade-off that entails vulnerabilities to other types of agents?

The selection of LPS as the experimental model was meant to cover this contingency, at least for the common denominator of acute inflammation many types of infections elicit. But LPS studies revealed potential weaknesses that some pathogens might exploit. One of these is the low expression of inducible nitric oxide. Although Nos2 gene knockouts in M. musculus had lower LPS-induced mortality than their wild-type counterparts, the mutants were more susceptible to the protozoan Leishmania major and the facultative intracellular bacterium Listeria monocytogenes ^{(56, 93)}. While there are no known studies of either of these pathogens in P. leucopus, the related species P. yucatanicus is the main reservoir for Leishmania mexicana in Mexico ⁽¹⁷⁾. Compared with M. musculus, which suffer a high fatality rate from experimental infections with L. mexicana, P. yucatanicus infections are commonly asymptomatic ⁽⁵⁴⁾.

Given the restrained interferon and ISG response shown by P. leucopus, another plausible vulnerability would be viral infections. But other studies suggest that neither RNA nor DNA viruses pose an inordinately high risk for Peromyscus. Both tolerance of and resistance to the tickborne encephalitis flavivirus Powassan virus by P. leucopus were demonstrated in an experimental model in which mice, by contrast, gwere severely affected ⁽⁶²⁾. P. maniculatus has been successfully infected with the SARS-CoV-2 virus by the respiratory route, but the infected animals displayed only mild pathology, manifested little if any disability, and recovered within a few days ^{(30, 33)}. Among natural populations and in the laboratory P. maniculatus is noted for its tolerance of hantavirus, which commonly is fatal for infected humans ^{(14, 20)}. P. maniculatus was permissive of infection with monkeypox virus, but the infection was mild and transient ⁽²⁶⁾.

A distinguishing P. leucopus characteristic, which was not expressly examined here, is its aforementioned 2-3 fold greater life span than that of M. musculus. While deermice may not be in the same longevity league as the naked mole-rat (Heterocephalus glaber), which can live for over 30 years ⁽⁷³⁾, some features of naked mole-rat immunology are intriguingly similar to what we have observed for P. leucopus. These include macrophages and blood myeloid cells with low to absent transcription of Nos2 or production of nitric oxide in response to LPS, even in the presence of added inteferon-gamma ⁽³¹⁾. Like P. leucopus and in distinction to M. musculus, naked mole-rats showed an increase in the proportion of neutrophils in the blood 4 hours after intraperitoneal injection of LPS ⁽³⁹⁾. In another comparative study, the hematopoietic stem and progenitor cells of these rodents had a lower type 1 interferon response than mice to a TLR3 agonist ⁽²⁸⁾.

In summary, if there is a vulnerability that Peromyscus accepts in return for relief from inflammation (and perhaps a longer life), it has not been identified yet. However, potential threats and stressors are many, and the number assessed either in the field or laboratory has been limited to date.

Implications for Lyme disease and other zoonoses

Our studies of P. leucopus began with a natural population and documented a >80% prevalence of infection and high incidence of re-infections by B. burgdorferi in the area’s white-footed deermouse, the most abundant mammal there ⁽¹⁶⁾. This was a Lyme disease endemic area ⁽²⁷⁾, where residents frequently presented for medical care for a variety of clinical manifestations, from mild to serious, of B. burgdorferi infection ⁽⁸⁵⁾. Subclinical infections in humans occur, but most of those who become infected have a definable illness ⁽⁸⁶⁾. The localized or systemic presence of the microbe is a necessary condition for Lyme disease, but the majority of the symptoms and signs are attributable to inflammation elicited by the organism’s presence and not from virulence properties per se or the hijacking of host cells ⁽²¹⁾. Since humans are transmission dead-ends for B. burgdorferi and many other zoonotic agents in their life cycles, it is not surprising that human infections are generally more debilitating if not fatal than what adapted natural hosts experience.

It is in the space between the asymptomatic natural host and symptomatic inadvertent host where there may be insights with basic and translational application. With this goal, we consider the ways the results inform studies of the pathogenesis of Lyme disease, where “disease” includes lingering disorders akin to “long Covid” ⁽⁷¹⁾, and where “pathogenesis” includes both microbial and host contributions. Plausibly-germane deermouse-mouse differences identified in our studies are summarized in Figure 11. Two are highlighted here.

Figure 11.

The first is macrophage polarization ⁽⁶⁸⁾. By the criteria summarized above, the response to LPS by P. leucopus is consistent with the alternatively-activated or M2 type, rather than the expected classical or M1 type. But it was not only LPS-treated deermice that had this attribute, the blood of untreated animals also displayed M2 type polarization features. This included a comparatively high Arg1 expression level and a Akt1/Akt2 ratio of more than 1 at baseline. This suggests that studies of other mammals, including humans, need not administer LPS or other TLR agonist to assess disposition toward M1 or M2-type polarization. This reading could serve as a prognostic indicator of the inflammatory response to infection with B. burgdorferi or other pathogen and the long-term outcome.

The second difference we highlight is the activation of ERV transcription that was prominent in the LPS-treated mice and rats but not in similarly-treated deermice. A paradoxical enlistment of antiviral defenses, including type 1 and type 2 interferons, for an infection with an extracellular bacterium, like B. burgdorferi, may bring about more harm than benefit, especially if the resultant inflammation persists after antibiotic therapy. There are various ways to assess ERV activation in the blood, including assays for RNA, protein, and reverse transcriptase activity. A xenotropic MLV-related retrovirus has been discounted as a cause of chronic fatigue syndrome ⁽⁶⁰⁾. However, production of whole virions need not occur for there to be PRR signaling in response to cytoplasmic Env protein, single stranded RNA, or cDNA ⁽⁷⁸⁾.

Methods

Animals

The study was carried out in accordance with the Guide for the Care and Use of Laboratory Animals: Eighth Edition of the National Academy of Sciences. The protocols AUP-18-020 and AUP-21-007 were approved by the Institutional Animal Care and Use Committee of the University of California Irvine.

Peromyscus leucopus, here also referred to as “deermice”, were of the outbred LL stock, which originated with 38 animals captured near Linville, NC and thereafter comprised a closed colony without sib-sib matings at the Peromyscus Genetic Stock Center at the University of South Carolina ⁽⁴⁵⁾. The LL stock animals for this study were bred and raised at the vivarium of University of California Irvine, an AAALAC approved facility. Outbred Mus musculus breed CD-1, specifically Crl:CD1(ICR) IGS, and here also referred to as “mice”, were obtained from Charles River. Fischer F344 strain inbred Rattus norvegicus, specifically F344/NHsd, and here also referred to as “rats”, were obtained from Charles River.

For the combined P. leucopus-M. musculus experiment the 20 P. leucopus were of a mean (95% confidence interval) 158 (156-159) days of age and had a mean 21 (19-22) g body mass. The 20 M. musculus were all 149 days of age and had a mean body mass of 47 (43-50) g. The ratio of average male to average female body mass was 1.04 for P. leucopus and 1.03 for M. musculus. The 6 female P. leucopus for the 12 h duration experiment were of a mean 401 (266-535) days of age and mean body mass of 20 (17-23) g. The 16 adult 10-12 week old female R. norvegicus had a mean 139 (137-141) g body mass. The 7 male P. leucopus for the infection study were of a mean 107 (80-134) days and mean body mass of 21 (18-24) g.

Animals were housed in Techniplast-ventilated cages in vivarium rooms with a 16 h-8 h light-dark cycle, an ambient temperature of 22 °C, and on ad libitum water and a diet of 2020X Teklad global soy protein-free extruded rodent chow with 6% fat content (Envigo, Placentia, CA).

For injections the rodents were anesthetized with inhaled isoflurane. The rodents were euthanized by carbon dioxide overdose and intracardiac exsanguination at the termination of the experiment. No animals died or became moribund before the 4 hour or 12 h termination time points in the LPS experiments or before the 5 d termination point of infection study.

LPS experiments

For the P. leucopus and M. musculus combined experiment, sample sizes replicated the specifications of the previous study, in which there were 20 P. leucopus and 20 M. musculus, equally divided between females and males and equally allotted between conditions ⁽³⁾. The treatments were administered in the morning of a single day. At 15 min intervals and alternating between species, sex, and treatments, animals were intraperitoneally (ip) injected 50 µl volumes of either ion-exchange chromatography-purified Escherichia coli O111:B4 LPS (Sigma-Aldrich L3024) in a dose of 10 µg per g body mass or the diluent alone: sterile-filtered, endotoxin-tested, 0.9% sodium chloride (Sigma-Aldrich). The animals were visually monitored in separate cages continuously for the duration of the experiment. We recorded whether there was reduced activity by criterion of huddling with little or movement for > 5 min, ruffled fur or piloerection), or rapid respiration rate or tachypnea. At 4.0 h time after injection animals were euthanized as described above, and sterile dissection was carried out immediately.

Lower dose and longer duration experiment. In an experiment with 6 P. leucopus, the animals were administered the same single dose of LPS but at 1.0 µg/g and the same control solution. The animals were euthanized 12 h after the injection the following day.

Rat LPS experiment. The same experimental design was used for the rats as for the combined deermice-mice experiment, with the exception that the formulation of the E. coli O111:B4 LPS was “cell culture grade” (Sigma-Aldrich L4391), and the groups were sterile saline alone (n = 5), 5 µg LPS per g body mass (n = 6), or 20 µg LPS per g (n = 5).

Experimental infection

The infection of a group of P. leucopus LL stock with the relapsing fever agent B. hermsii and the processing of blood and tissues for RNA extraction 5 days into the infection were described previously ⁽³⁾. In brief, animals were infected intraperitoneally on day 0 with either phosphate-buffered saline alone or 10³ cells of B. hermsii MTW, a strain that is infectious for Peromyscus species ⁽⁴³⁾. Bacteremia was confirmed by microscopy on day 4, and the animals were euthanized on day 5. For that prior study the RNA-seq analysis was limited to the genome-wide transcript reference set. For the present study we used the original fastq format files for targeted RNA-seq as described below.

Hematology and plasma analyte assays

For the combined P. leucopus-M. musculus experiment automated complete blood counts with differentials were performed at Antech Diagnostics, Fountain Valley, CA on a Siemens ADVIA 2120i with Autoslide hematology instrument with manual review of blood smears by a veterinary pathologist. For the 12 h duration P. leucopus experiment hematologic parameters were analyzed on an ABCVet Hemalyzer automated cell counter instrument at U.C. Irvine. For the rat experiment complete blood counts with differentials were performed at the Comparative Pathology Laboratory of the University of California Davis. Multiplex bead-based cytokine protein assay of the plasma of the rats was performed at Charles River Laboratories using selected options of the Millipore MILLIPLEX MAP rat cytokine/chemokine panel.

RNA extraction of blood

After the chest cavity was exposed, cardiac puncture was performed through a 25 gauge needle into a sterile 1 ml polypropylene syringe. After the needle was removed, the blood was expelled into Becton-Dickinson K2E Microtainer Tubes, which contained potassium EDTA. Anticoagulated blood was split into a sample that was placed on ice for same-day delivery to the veterinary hematology laboratory and a sample intended for RNA extraction which was transferred to an Invitrogen RiboPure tube with DNA/RNA Later and this suspension was stored at −20 °C. RNA was isolated using the Invitrogen Mouse RiboPure-Blood RNA Isolation Kit]. RNA concentration was determined on a NanoDrop microvolume spectrophotometer (ThermoFisher) and quality was assessed on an Agilent Bioanalyzer 2100.

RNA-seq of blood

The chosen sample sizes and coverage for the bulk RNA-seq were based on empirical data from the prior study ⁽³⁾, which indicated that with 10 animals per group and a two-sided two sample t-test we could detect with a power of ≥ 0.80 and at a significance level of 0.05 a ≥1.5 fold difference in transcription between groups for a given gene. We also were guided by the simulations calculations of Hart et al. ⁽³⁶⁾, which indicated for a biological coefficient of variation of 0.4 within a group, a minimum depth of coverage of ≥ 10, and a target of ≥ 2x fold change that a sample size of 7-8 was sufficient for 80% power and type I error of 5%. For the P. leucopus and M. musculus samples production of cDNA libraries was with the Illumina TruSeq Stranded mRNA kit. After normalization and multiplexing, the libraries were sequenced at the University of California Irvine’s Genomic High Throughput Facility on a Illumina NovaSeq 6000 instrument with paired-end chemistry and 150 cycles to achieve ∼100 million reads per sample for the combined P. leucopus-M. musculus experiment. The same method for producing cDNA libraries was used for the R. norvegicus RNA and the P. leucopus in the infection study, but these were sequenced on a Illumina HiSeq 4000 instrument with paired-end chemistry and 100 cycles. The quality of sequencing reads was analyzed using FastQC (Babraham Bioinformatics). The reads were trimmed of low-quality reads (Phred score of <15) and adapter sequences, and corrected for poor-quality bases using Trimmomatic ⁽¹²⁾.

For the combined species experiment the mean (95% CI) number of PE150 reads per animal after trimming for quality was 1.1 (1.0-1.2) x 10⁸ for P. leucopus and 1.1 (1.0-1.2) x 10⁸ for M. musculus (p = 0.91). For P. leucopus of this experiment a mean of 83% of the reads mapped to the genome transcript reference set of 54,466; mean coverages for all transcripts and for the mean 62% of reference transcripts with ≥ 1x coverage were 97x and 157x, respectively. For M. musculus of this experiment a mean 91% of the reads mapped to the genome transcript reference set of 130,329; mean coverages for all transcripts and for the mean of 21% of reference transcripts with ≥ 1x coverage were 103x and 568x, respectively. For the lower dose-longer duration experiment with P. leucopus the mean number of PE150 reads was 2.5 (2.3-2.6) x 10⁷. For the rat experiment the mean number of PE100 reads was 2.4 (2.2-2.5) x 10⁷. For the B. hermsii infection experiment the mean number of PE100 reads was 4.9 (4.5-5.3) x 10⁷.

Batched fastq files were subjected to analysis with CLC Genomics Workbench version 23 (Qiagen). Library size normalization was done by the TMM (trimmed mean of M values) method of Robinson and Oshlack ⁽⁷⁷⁾. The reference genome transcript sets on GenBank were the following: GCF_004664715.2 for P. leucopus LL stock, GCF_000001635.27_GRCm39 for M. musculus C57Bl/6, and GCF_015227675.2_mRatBN7.2 for R. norvegicus. The settings for PE150 reads were as follows for both strands: length fraction of 0.35, similarity fraction of 0.9, and costs for mismatch, insertion, or deletion of 3. For PE100 reads the settings were the same except for length fraction of 0.4. Principal Component Analysis was carried with the “PCA for RNA-Seq” module of the suite of programs.

For the P. leucopus RNA-seq analysis there were 54,466 reference transcripts, of which 48,164 (88%) were mRNAs with protein coding sequences, and 6,302 were identified as non-coding RNAs (ncRNA). Of the 48,164 coding sequences, 40,247 (84%) had matching reads for at least one of the samples. The five most highly represented P. leucopus coding sequences among the matched transcripts of whole blood among treated and control animals were hemoglobin subunits alpha and beta, the calprotectin subunits S100A8 and S100A9, and ferritin heavy chain. For the M. musculus analysis there were available 130,329 reference transcripts: 92,486 (71%) mRNAs with protein coding sequences and 37,843 ncRNAs. Of the coding sequences, 59,239 (64%) were detectably transcribed in one or both groups by the same criterion. The five most highly represented coding sequences of mRNAs of identified genes for M. musculus were hemoglobin subunits alpha and beta, aminolevulinic synthase 2, ferritin light polypeptide 1, and thymosin beta. For R. norvegicus there were 99,126 reference transcripts, of which 74,742 (75%) were mRNAs. The five most highly represented coding sequences of mRNAs of identified genes for R. norvegicus were hemoglobin subunits alpha and beta, beta-2 microglobulin, ferritin heavy chain, and S100a9.

Genome-wide differential gene expression

Differential expression between experimental conditions was assessed with an assumption of a negative binomial distribution for expression level and a separate Generalized Linear Model for each ⁽⁵⁹⁾. Fold changes in TPM (transcripts per million) were log2-transformed. The False Discovery Rate (FDR) with corrected p value was estimated by the method of Benjamini and Hochberg ⁽¹⁰⁾. To assess the limit of detection for differentially expressed genes between 10 animals treated with LPS and 10 with saline alone, we took the data for 4650 reference transcripts for which the mean TPM across 20 P. leucopus was >10 and randomly permuted the data to achieve another 9 sets and calculated the fold-change of sub-groups of 10 and 10 with one random group serving as the proxy of the experimental treatment and other second as the control for each of the sets. The expectation was that mean fold-change of the 9 permuted sets and the 4650 reference sequences would be ∼1. The result was a mean and median of 1.08 with a 99.9% asymmetric confidence interval for the mean of 0.81-1.49. This was an indication that the choices for sample sizes were realistic for achieving detection of ≥ 1.5x fold changes.

Gene Ontology term analysis

M. musculus was selected as the closest reference for the P. leucopus data. The analysis was implemented for data for differentially expressed genes meeting the criteria of a FDR p-value ≤ 0.01 and fold-change of ≥ 1.5. The analysis was implemented with the tools of Metascape (https://metascape.org) ⁽⁹⁹⁾. Functional enrichment analysis was carried out first with the hypergeometric test and FDR p-value correction ⁽¹⁰⁾. Then pairwise similarities between any two enriched terms were computed based on a Kappa-test score ⁽²²⁾. Similarity matrices were then hierarchically clustered and a 0.3 similarity threshold was applied to trim resultant trees into separate clusters. The lower the p-value, the less the likelihood the observed enrichment is due to randomness ⁽⁹⁸⁾. The lowest p-value term represented each cluster shown in the horizontal bar graph. Besides the terms beginning with “GO” and referring to the Gene Ontology resource (http://geneontology.org) ⁽²⁾, others refer to Kegg Pathway database (https://www.kegg.jp) for “mmu..” designations, WikiPathways database (https://www.wikipathways.org) for “WP…” designations, and Reactome database (https://reactome.org) for “R-MMU…” designations.

Targeted RNA-seq

RNA-seq of selected set of protein coding sequences (CDS), which are listed below, was carried out using CLC Genomics Workbench v. 23 (Qiagen). Paired-end reads were mapped with a length fraction of 0.35 for ∼150 nt reads and 0.40 for ∼100 nt reads, a similarity fraction of 0.9, and costs of 3 for mismatch, insertion, or deletion to the CDS of sets of corresponding orthologous mRNAs of P. leucopus, M. musculus, and R. norvegicus. Preliminary expression values were unique reads normalized for total reads across all the samples without adjustment for reference sequence length, as described ⁽³⁾. Exceptions were the endogenous retrovirus coding sequences which differed in lengths between species. For within-and cross-species comparisons we initially normalized three different ways after quality filtering and removing vector and linker sequence: for total reads for the given sample, for unique reads for 12S ribosomal RNA for the mitochondria of nucleated cells in the blood, and for unique reads for the gene Ptprc, which encodes CD45, a marker for both granulocytes and mononuclear cells in the blood. This is described in more detail in Results. Following the recommendation of Hedges et al. we used the natural logarithm (ln) of ratios ⁽³⁷⁾.

The target CDS were as follows: Acod1, Akt1, Akt2, Arg1, Bcl3, Camp, Ccl2, Ccl3, Ccl4, Cd14, Cd177, Cd3d, Cd4, Cd69, Cd8, Cfb, Cgas, Csf1, Csf1r, Csf2, Csf3, Csf3r, Cx3cr1, Cxcl1, Cxcl10, Cxcl2, Cxcl3, Dhx58, Fcer2, Fcgr2a, Fcgr2b, Fcgr3 (CD16), Fgr, Fos, Fpr2, Gapdh, Gbp4, Glrx, Gzmb, Hif1a, Hk3, Hmox1, Ibsp, Icam, Ifih1, Ifit1, Ifng, Il10, Il12, Il18, Il1b, Il1rn, Il2ra, Il4ra, Il6, Il7r, Irf7, Isg15, Itgam, Jak1, Jak2, Jun, Lcn2, Lpo, Lrg, Lrrk2, Ltf, Mapk1, Mmp8, Mmp9, Mpo, MT-Co1, Mt2, Mtor, Mx2, Myc, Myd88, Ncf4, Nfkb1, Ngp, Nos2, Nox1, Nr3c1, Oas1, Olfm4, Padi4, Pbib, Pkm, Ptx, Ptprc (CD45), Retn, Rigi (Ddx58), S100a9, Saa3, Serpine1, Slc11a1 (Nramp), Slpi, Socs1, Socs3, Sod2, Stat1, Stat2, Stat4, Steap1, Sting, Tgfb, Thy1, Timp1, Tlr1, Tlr2, Tlr4, Tnf, Tnfrsf1a, and Tnfrsf9. The sources for these coding sequences were the reference genome transcript sets for P. leucopus, M. musculus, and R. norvegicus listed above. If there were two or more isoforms of the mRNAs and the amino acid sequences differed, the default selection for the coding sequence was the first listed isoform. The lengths of the orthologous pairs of P. leucopus and M. musculus coding sequences were either identical or within 2% of the other. Fcgr1, the gene for high affinity Fc gamma receptor I or CD64, was not included in the comparison, because in P. leucopus it is an untranscribed pseudogene ⁽⁷⁾. For the targeted RNA-seq of blood of deermice infected with B. hermsii (Table 6), the reference sequence was the cp6.5 plasmid of B. hermsii (NZ_CP015335).

Quantitative PCR assays

Reverse transcriptase (RT)-qPCR assays and the corresponding primers for measurement of transcripts of genes for nitric oxide synthase 2 (Nos2) and glyceraldehyde 3-phosphate dehydrogenase (Gapdh) were those described previously ⁽³⁾. These primers worked for M. musculus as well as P. leucopus using modified cycling conditions. For the arginase 1 transcript assays different primer sets were used for each species. The forward and reverse primer sets for the 352 bp Arg1 product for P. leucopus were 5’-TCCGCTGACAACCAACTCTG and 5’-GACAGGTGTGCCAGTAGATG, respectively. The corresponding primer pairs for a 348 bp Arg1 of M. musculus were 5’-TGTGAAGAACCCACGGTCTG and 5’-ACGTCTCGCAAGCCAATGTA. cDNA synthesis and qPCR were achieved with a Power Sybr Green RNA-to-Ct 1-Step Kit (Applied Biosystems) in 96 MicroAmp Fast Reaction Tubes using an Applied Biosystems StepOne Plus real-time PCR instrument. The initial steps for all assays were 48 °C for 30 min and 95 °C for 10 min. For Arg1 and Nos2 assays, this was followed by 40 cycles of a 2-step PCR of, first, 95°C for 15 s and then, second, annealing and extension at 60 °C for 1 min. The cycling conditions for Gapdh were 40 cycles of 95 °C for 15 s followed by 60 °C for 30 s. Quantitation of genome copies of B. hermsii in extracted DNA was carried out by probe-based qPCR as described ⁽⁶⁾.

Additional statistics

Means are presented with asymmetrical 95% confidence intervals (CI) to accommodate data that was not normally distributed. Parametric (t test) and non-parametric (Mann-Whitney) tests of significance were 2-tailed. Unless otherwise stated, the t test p value is given. Categorical variables were assessed by 2-tailed Fisher’s exact test. FDR correction of p values for multiple testing was by the Benjamini-Hochberg method ⁽¹⁰⁾, as implemented in CLC Genomics Workbench (see above), or False Discovery Rate Online Calculator (https://tools.carbocation.com/FDR). Discriminant Analysis, linear regression, correlation, coefficient of determination, and General Linear Model analyses were performed with SYSTAT v. 13.1 software (Systat Software, Inc.). Box plots with whiskers display the minimum, first quartile, median, third quartile, and maximum.

Data availability

The experiments reported here are associated with National Center for Biotechnology Information (https://ncbi.nlm.nih.gov) BioProjects PRJNA975149 for the combined P. leucopus and M. musculus experiment, PRJNA874306 for the 12 h-lower dose P. leucopus experiment, PRJNA973677 for the R. norvegicus LPS experiment, and PRJNA508222 for the B. hermsii infection of P. leucopus experiment. PRJNA975149 includes 40 BioSamples (SAMN35347136-SAMN35347175) and 40 sets of Sequence Read Archive (SRA) fastq files of Illumina paired-end chemistry reads 1 and 2 (SRR24733648-SRR24733687). PRJNA874306 includes 6 BioSamples (SAMN30561752-SAMN30561757) and the corresponding SRA files SRR24451178-SRR24451180, SRR24451183, SRR24451194, and SRR24451195. PRJNA973677 includes 16 BioSamples (SAMN35351370-SAMN35351385) and the corresponding SRA files SRR24731678-SRR24731693. PRJNA508222 includes 7 BioSamples (SAMN10522571-SAMN10522573; SAMN10522575-SAMN10522578) and SRA files (SRR8283809; SRR8283811-SRR8283816). The datasets of the genome-wide RNA-seq and the differentially-expressed gene analysis of RNA-seq data for the 4 h LPS experiments for P. leucopus, M. musculus, and R. norvegicus are available under the title “Differentially-expressed genes in blood in response to lipopolysaccharide in three rodent species” at the public data repository Dryad (https://datadryad.org). The link for Dryad Tables D1-D6 and their descriptions is https://doi.org/10.5061/dryad.zw3r228dh. A Transcriptome Shotgun Assembly (TSA) for stranded mRNA of blood of the 20 female and male P. leucopus of the 4 h LPS experiment comprised 14,982 contigs of ≥ 400 nt and ≥ 20x coverage and have TSA accession number GKOE00000000.

Acknowledgements

We thank Hanjuan Shao for technical assistance, Vanessa Cook for her participation in the R. norvegicus experiment, Anthony Long for bioinformatic advice and contributions, and Brianna Craver-Hoover and Gajalakshmi Ramanathan for assistance with blood cell counts performed at U.C. Irvine. The experimental studies reported here were supported by National Institutes of Health (NIH) grants AI-157513 and AI-136523. The services of the Genomics Research and Technology Hub were administratively supported in part by NIH Cancer Center Support Grant P30 CA-062203 and NIH shared instrumentation grants RR-025496, OD-010794, and OD-021718.