# Differential immune gene expression associated with contemporary range expansion of two invasive rodents in Senegal

Nathalie Charbonnel, Maxime Galan, Caroline Tatard, Anne Loiseau, Christophe Amidi Diagne, Ambroise Dalecky, Hugues Parrinello, Stéphanie Rialle, Dany Severac, Carine Brouat

► To cite this version:

Nathalie Charbonnel, Maxime Galan, Caroline Tatard, Anne Loiseau, Christophe Amidi Diagne, et al.. Differential immune gene expression associated with contemporary range expansion of two invasive rodents in Senegal. Scientific Reports, 2020, Peer Community in Ecology, 10, pp.18257. 10.1101/442160. hal-02790736v2

HAL Id: hal-02790736

https://hal.inrae.fr/hal-02790736v2

Submitted on 30 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Differential immune gene expression associated with contemporary range expansion in two invasive rodents in Senegal

Nathalie Charbonnel<fnref n="1" />✉, Maxime Galan<fnref n="1" />, Caroline Tatard<fnref n="1" />, Anne Loiseau<fnref n="1" />, Christophe Diagne>D<fnref n="1" />,², Ambroise Dalecky³, Hugues Parrinello⁴, Stephanie Rialle⁴, Dany Severac⁴ & Carine Brouat<fnref n="1" />

Biological invasions are major anthropogenic changes associated with threats to biodiversity and health. However, what determines the successful establishment and spread of introduced populations remains unclear. Here, we explore several hypotheses linking invasion success and immune phenotype traits, including those based on the evolution of increased competitive ability concept. We compared gene expression profiles between anciently and recently established populations of two major invading species, the house mouse *Mus musculus domesticus* and the black rat *Rattus rattus*, in Senegal (West Africa). Transcriptome analyses identified differential expression between anciently and recently established populations for 364 mouse genes and 83 rat genes. All immune-related genes displaying differential expression along the mouse invasion route were overexpressed at three of the four recently invaded sites studied. Complement activation pathway genes were overrepresented among these genes. By contrast, no particular immunological process was found to be overrepresented among the differentially expressed genes of black rat. Changes in transcriptome profiles were thus observed along invasion routes, but with different specific patterns between the two invasive species. These changes may be driven by increases in infection risks at sites recently invaded by the house mouse, and by stochastic events associated with colonization history for the black rat. These results constitute a first step toward the identification of immune eco-evolutionary processes potentially involved in the invasion success of these two rodent species.

Biological invasions are considered to be one of the five most important anthropogenic changes, with major impacts on biodiversity<fnref n="1" />,² and human health³,⁴. Considerable progress has been made in determining population histories and the invasive routes of exotic species, but the mechanisms driving the successful establishment and spread of introduced populations remain unclear. In particular, the assertion that invasion success depends on the contemporary evolution of phenotypic traits of advantage at invasion fronts has seldom been tested, particularly in vertebrates (but see⁵,⁶).

From an eco-evolutionary perspective, invasion success may depend on pre-adaptation within the original range⁷,⁸ or on rapid changes in phenotypic traits (through selection on new/standing genetic variation or adaptive plasticity) advantageous in newly colonized areas⁷,⁹. In this context, many hypotheses have been formulated concerning traits related to the ability to cope with new environments¹⁰. In particular, large changes in parasite pressure can occur in introduced species, because of enemy release or the spillover of native parasites, along invasion gradients between the native and invaded areas¹¹⁻¹⁵ or within invaded regions, from introduction areas to invasion fronts¹⁶,¹⁷. Invasion success may therefore depend on immune strategies. Two general hypotheses have been developed, in accordance with the life history theory. The “evolution of increased competitive ability”

<footnote n="1">¹CBGP, INRAE, CIRAD, IRD, Montpellier SupAgro, Univ Montpellier, Montpellier, France. ²Département de Biologie Animale, Faculté des Sciences et Techniques, Université Cheikh Anta Diop (UCAD), Fann, Dakar, Senegal. ³LPED, IRD, Aix Marseille Univ, Marseille, France. ⁴MGX-Montpellier GenomiX, c/o Institut de Génomique Fonctionnelle, Montpellier, France. ✉email: nathalie.charbonnel@inrae.fr</footnote>
(EICA)¹⁸ is based on trade-offs between immune traits and other life history traits. According to this theory, the loss of parasites that often occurs during the course of invasion should lead to a release of energy resources previously dedicated to anti-parasite defenses. The reallocation of this energy to other life history traits, increasing competitive abilities and favoring range expansion (e.g. dispersal, reproduction), would therefore be expected to be favored. According to the “EICA-refined” hypothesis¹⁹, invaders may encounter new infection risks during the course of invasion, due to a spillover of native parasites, for example. In such situations, the allocation of energy to anti-parasite defenses would be expected to increase to a lesser extent, instead being modulated in favor of defense strategies with lower costs. Indeed, different immune pathways seem to have different developmental, maintenance and use costs²⁰. In vertebrates, costly immune pathways, such as inflammation, may be attenuated, whereas less costly pathways, such as antibody-mediated responses, may be favored⁵,¹⁹,²¹,²². As an alternative to both the EICA and EICA-refined hypotheses, it is also possible that invaders display stronger immune responses as a result of exposure to new parasites and their pathologic effects at invasion fronts.

A recent analysis reviewing empirical studies showed that, in most cases, immune responses display varia-
tions during invasions, although no general pattern supporting the EICA hypotheses was detected²¹. The authors
highlighted a problem due to the lack of studies reporting a large array of immune responses for the correct
assessment of immunocompetence and estimation of immune trade-offs. Unfortunately, immune phenotyping
is often based on a small number of immune effectors, due to the availability of only small amounts of material
(blood, tissue), or the difficulty obtaining immune kits for non-model species.

The recent advent of ‘omics’ technologies has made it possible to explore a large set of genes and pathways involved in the response to new conditions experienced at invasion fronts. Genome scans analyzing genome-wide variation with a view to detect loci evolving under positive directional selection have been successfully performed in some case studies of invasion (e.g.,²³⁻²⁵). These scans have identified signatures of contemporary adaptation associated with invasion, some of which are related to immunity²⁶. Nevertheless, this genomic approach, based on the detection of single-nucleotide polymorphism outliers, cannot detect rapid adaptation occurring over a few generations, through phenotypic plasticity, for example. Transcriptomics, the analysis of gene expression at a genome-wide scale, is a complementary approach that could potentially fill in these gaps, making it possible to decipher the molecular mechanisms underlying phenotypic changes²⁷,²⁸. Indeed, gene expression is an essential mechanism for rapid acclimatization or adaptation to new environments²⁹⁻³³. Transcriptomics has thus been used to identify genes and gene regulatory pathways contributing to phenotypic variation at invasion fronts. Evidence for transcriptomic divergence along invasion gradients has already been reported²⁵,³⁴, in some cases highlighting the important role of immune regulation in invasion success³⁵. However, only one study to date has analyzed immune gene expression in the context of established variations of parasite pressures along invasion routes³⁶. Such studies are essential, to test the EICA hypotheses and to increase our understanding of the interactions between immune system regulation and invasion success²¹.

We investigated transcriptional profiles along the invasion routes of the house mouse *Mus musculus domesticus* and the black rat *Rattus rattus*, which are currently invading Senegal (West Africa). Both these rodent species originate from Asia and made use of human migration to expand their distribution range worldwide⁵⁷,³⁸. In Senegal, historical records (see references in³⁹ and molecular analyses⁴⁰,⁴¹) show that these rodents were first brought to Senegal, by sea, by European explorers and settlers. They were reported to be common in coastal villages and towns from the middle of the nineteenth century, and their distribution areas remained restricted to the coastal area until the beginning of the twentieth century. Both taxa then spread further inland, due to increases in human activities and the further development of transport infrastructures. Our previous study, based on a small number of immune effectors, highlighted the occurrence of stronger inflammatory and/or antibody-mediated responses at recently invaded sites than at anciently invaded ones, for both invaders⁴². These immune variations may reflect responses to new parasite pressures encountered in recently invaded areas. Surveys of the bacterial communities present in the spleens of rodents from Senegal were found to be consistent with this prediction⁴³, as changes in bacterial composition were observed along the routes of invasion of both the house mouse and the black rat. Exotic rodents were exposed to a greater diversity of *Mycoplasma* sp. and to higher levels of *Bartonella* sp. at recently invaded sites, due to the high prevalence of these bacteria in the native rodent *Mastomys erythroleucus*⁴³. By contrast, decreases in the pressure due to certain parasites were also observed, as the number of helminths present in the gastrointestinal tract of rodents was smaller at recently invaded sites than at anciently invaded ones, suggesting enemy release for these macroparasites⁴⁴. Decreases were observed in both the overall prevalence and species richness of gastrointestinal helminths, for both the house mouse and the black rat⁴⁴.

We therefore developed a whole-RNA sequencing (i.e., RNAseq) approach to assess differential patterns of gene expression between populations from recently invaded and long-established invasion sites. We focused on the spleen in both species, as this immune-related organ of vertebrates is commonly used for immunological studies. We tested the null hypothesis of an absence of difference in immune gene expression patterns along invasion routes (i.e. recently invaded sites vs. long-established invasion sites). We tested two alternative hypotheses. According to the EICA-refined hypothesis, expression levels would be expected to be lower for genes encoding proteins involved in energetically costly immune pathways (e.g., inflammation) and higher for genes encoding proteins involved in cost-effective immune pathways (e.g., antibody-mediated responses), at recently invaded sites¹⁹. The changes in infection risk at recently invaded sites may mediate such trade-offs between energetically costly and cost-effective immune pathways⁴³,⁴⁴. We also tested the alternative hypothesis that rodents from recently invaded sites display higher overall levels of immune gene expression in response to the new parasite pressures encountered⁵,⁴⁵,⁴⁶.

More specifically, we investigated the following questions: (1) Does immune gene expression differ between anciently and recently invaded populations of mice and rats? (2) Are certain functional categories of immune genes overrepresented among the differentially expressed genes? (3) Do these results support the EICA hypothesis
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sampling sites for the house mouse, *Mus musculus domesticus* (symbols in white) and the black rat, *Rattus rattus* (symbols in black) invasion routes. Triangles and squares correspond, respectively, to anciently and recently invaded sites. DAG Dagathie, MBA Mbakhana, THL Thilene, NDB Ndombo, CRB Croisement Boube, DOD Dodel, AEL Aere Lao, LOU Lougue, DIK Diakene Wolof, DIT Diattacounda, MAR Marsassoum, TOB Tobor, BAN Badi Nieriko, BOU Boutougoufara, KED Kedougou, SOU Soutouta.</figcaption></figure>

or the EICA-refined hypotheses? Or are these results consistent with an overall increase in immune gene expression reflecting the new parasite pressures experienced at the invasion front?

**Materials and methods**

**Ethics statements.** Sampling campaigns were performed on private property with explicit prior agreement from family, or the local authority or institution concerned. None of the rodent species investigated here has protected status (see list of the International Union for Conservation of Nature). All animal-related procedures were performed according to the official ethics guidelines of the American Society of Mammalogists⁴⁷. All protocols presented here were performed with the explicit prior agreement (CBGP: D 34-169-1) of the relevant institutional committee (Regional Head of the Veterinary Service, Hérault, France). They were also performed in accordance with the requirements of Senegalese and French law.

**Sampling sites.** We used the spleens of *M. m. domesticus* and *R. rattus* animals used in a previously described bacterial metabarcoding study⁴³. Briefly, live trapping was performed for both these invasive species at four anciently sites and four recently invaded ones (Fig. 1).

No accurate estimate of the time of introduction was available for these rodents, but we used historical records to differentiate between anciently invaded sites and recently invaded sites. Exotic mice and rats were common in the villages and cities of coastal areas of Senegal in the middle of the nineteenth century, but the available records did not provide precise geographic locations (see refs in³⁹). We therefore considered coastal areas to have been colonized between 300 and 100 years ago. Regular rodent community surveys have been performed across Senegal since the 1970s. These surveys indicate that, before the 1990s, the house mouse was absent from the villages and cities of the middle of the Senegal River valley, and the black rat was absent from western Senegal (data summarized in³⁹). We therefore considered inland sites to have been colonized less than 30 years ago. Previous population genetic studies confirmed that most of the house mouse populations from the north of Senegal belonged to the same genetic group, which colonized northern Senegal from Saint Louis⁴¹. Black rat populations from south-eastern Senegal form a separate genetic group, which originated from the coastal area south of Dakar⁴⁰.

**Sample collection and RNA extraction.** Fieldwork was conducted during the dry season, which runs roughly from October to June in northern Senegal and from November to May in southern Senegal. Mice and rats breed year-round in Senegal, as generally reported in commensal conditions (see refs in⁴⁸), but we chose to define a standardized sampling period, to avoid the possibility of confounding factors associated with sea-
<table><caption>**Table 1.** Sampling sites and their invasion-related categorization as well as rodent sample sizes included in each biological replicate along the invasion routes of the house mouse and of the black rat. N total number of rodents included, M number of males, F number of females.</caption><thead><tr><th rowspan="2">Invasion category</th><th colspan="2">House mouse</th><th colspan="2">Black rat</th></tr><tr><th>Sites</th><th>N (M/F)</th><th>Sites</th><th>N (M/F)</th></tr></thead><tbody><tr><td rowspan="10">Anciently invaded sites (&gt; 100 years)</td><td>Dagathie</td><td></td><td>Diakene Wolof</td><td></td></tr><tr><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>6 (3/3)</td></tr><tr><td>Pool of RNA #2</td><td>8 (4/4)</td><td>Pool of RNA #2</td><td>8 (4/4)</td></tr><tr><td>Mbakhana</td><td></td><td>Diattacounda</td><td></td></tr><tr><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>10 (5/5)</td></tr><tr><td>Pool of RNA #2</td><td>10 (5/5)</td><td>Pool of RNA #2</td><td>10 (5/5)</td></tr><tr><td>Thilene</td><td></td><td>Marsassoum</td><td></td></tr><tr><td>Pool of RNA #1</td><td>8 (4/4)</td><td>Pool of RNA #1</td><td>10 (5/5)</td></tr><tr><td>Pool of RNA #2</td><td>8 (4/4)</td><td>Pool of RNA #2</td><td>10 (5/5)</td></tr><tr><td>Ndombo</td><td></td><td>Tobor</td><td></td></tr><tr><td rowspan="10">Recently invaded sites (&lt; 30 years)</td><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>6 (3/3)</td></tr><tr><td>Pool of RNA #2</td><td>10 (5/5)</td><td>Pool of RNA #2</td><td>8 (4/4)</td></tr><tr><td>Croisement Boube</td><td></td><td>BadiNieriko</td><td></td></tr><tr><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>8 (4/4)</td></tr><tr><td>Pool of RNA #2</td><td>10 (5/5)</td><td>Pool of RNA #2</td><td>8 (4/4)</td></tr><tr><td>Dodel</td><td></td><td>Boutougoufara</td><td></td></tr><tr><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>8 (4/4)</td></tr><tr><td>Pool of RNA #2</td><td>8 (4/4)</td><td>Pool of RNA #2</td><td>8 (4/4)</td></tr><tr><td>Aere Lao</td><td></td><td>Kedougou</td><td></td></tr><tr><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>8 (4/4)</td></tr><tr><td>Pool of RNA #2</td><td>10 (5/5)</td><td>Pool of RNA #2</td><td>8 (4/4)</td></tr><tr><td>Lougue</td><td></td><td>Soutouta</td><td></td></tr><tr><td rowspan="2"></td><td>Pool of RNA #1</td><td>10 (5/5)</td><td>Pool of RNA #1</td><td>10 (5/5)</td></tr><tr><td>Pool of RNA #2</td><td>8 (4/4)</td><td>Pool of RNA #2</td><td>10 (5/5)</td></tr></tbody></table>

sonal variations. Rodents were killed by cervical dislocation at approximately the same period of the day for all individuals (11 a.m.–2 p.m.). Autopsies were performed and the spleens were collected and immediately placed in RNAlater, stored at 4 °C overnight and then at –20 °C until further analyses. Only adult rodents (mass > 9 g and/or length > 72 mm for house mice, mass > 54 g and/or length > 145 mm for black rats, see$^{48}$) were included in the analyses, as individual age (at least juveniles vs adults) is known to have a strong influence on immune responses$^{49}$. None of the trapped rodents displayed any sign of disease.

We considered 20 adult rodents per site for all but five sites at which only 14–18 individuals were available (see Table 1). Total RNA was extracted from approximately 5 mg of spleen from each sample with the AllPrep 96 DNA/RNA Kit (Qiagen). The quality and quantity of the purified RNA was assessed by gel electrophoresis and on a NanoDrop spectrophotometer (Thermo Scientific) before pooling, and with a Bioanalyzer 2100 (Agilent) or Fragment Analyzer (Advanced Analytical) after pooling (see below). The RNA integrity (RIN) or quality (RQN) numbers of each pool are provided in Fig. S1. For each rodent species and each site, we constituted two pools of individuals. These two pools contained similar numbers of individuals (from 6 to 10) and had a balanced sex ratio (Table 1). For a given site, the two pools were considered to be ‘biological replicates’ of the site concerned. We therefore constituted 16 equivalent pools of individuals (Table 1). This design made it possible to attain the number of replicates recommended in gene expression analyses for which high levels of variability between replicates are expected (between 6 and 12 replicates per category, see$^{50}$).

**cDNA library preparation and RNA sequencing.** cDNA library construction, sequencing and the alignment of filtered read sequences were performed at the MGX platform. RNA-Seq libraries were constructed with the Truseq stranded mRNA sample preparation (low-throughput protocol) kit from Illumina. We used 1 µg of total RNA for library construction. We first purified poly-A-containing mRNA molecules with a poly-T oligomer attached to magnetic beads (poly-A-based mRNA enrichment). The mRNA was then fragmented into small pieces with divalent cations at high temperature. The cleaved RNA fragments were used to generate first-strand cDNA with actinomycin D, random hexamer primers and SuperScript II reverse transcriptase for *M. m. domesticus* or SuperScript IV reverse transcriptase for *R. rattus*. The second strand of the cDNA was synthesized by replacing dTTP with dUTP. These cDNA fragments had an additional single ‘A’ base and were subsequently ligated to the adapter. The products were then purified and enriched over 15 cycles of PCR. The final cDNA libraries were validated and quantified with a KAPA qPCR kit.

For *M. m. domesticus*, four cDNA libraries per sequencing line were pooled, in equal proportions, denatured with NaOH and diluted to 7 pM before clustering. Cluster formation was performed on a Flowcell V3, primer
hybridization and 50 cycles of single-end read sequencing were performed on a cBot and a HiSeq2000 (Illumina, San Diego, CA), respectively. For *R. rattus*, four cDNA libraries per sequencing line were pooled, in equal proportions, denatured with NaOH and diluted to 8 pM before clustering. Cluster formation was performed on a Flowcell V4, primer hybridization and 100 cycles of single-end read sequencing were performed on a cBot and a HiSeq2500 (Illumina, San Diego, CA), respectively. We performed 100 cycles of single-end sequencing for *R. rattus*, as the only published genome sequence available for rats is from another species, *Rattus norvegicus*.

**Transcriptome analysis.** *Sequence alignment and gene quantification.* Image analyses and base calling were performed with Illumina HiSeq Control Software and the Real-Time Analysis component. Demultiplexing was performed and Fastq files were generated with Illumina conversion software (CASAVA 1.8.2 for the house mouse data, BCL2FASTQ 2.17 for the black rat data). The quality of the raw data was assessed with FASTQC from the Babraham Institute and the Illumina software SAV (Sequencing Analysis Viewer).

A splice junction mapper, *TOPHAT*<sup>51</sup> (v2.0.9 for the mouse data, v2.0.13 for the black rat data), using *BOWTIE*<sup>52</sup> (v2.1.0 for the mouse data, v2.2.3 for the black rat data), was used to align the reads to the *M. musculus* genome (UCSC mm10) or the *Rattus norvegicus* genome (UCSC rn4) with a set of model gene annotations (genes.gtf downloaded from UCSC on March 6, 2013 for the mouse genome, on July 17, 2015 for the brown rat genome). Final read alignments containing more than three mismatches (house mouse) or nine mismatches (black rat) were discarded. Again we used a different threshold for mouse and rat data, as the reference genome sequence for rat was that of *Rattus norvegicus*, and genomic studies have shown that these two species diverged about 2.9 million years ago<sup>53</sup>. Read alignment rates were above 83.3% for all *M. m. domesticus* libraries and above 74.0% for all *R. rattus* libraries. *SAMTOOLS* (v0.1.18 for the mouse data, v1.2 for the black rat data) was used to sort the alignment files. Genes were then counted with *HTSEQ* count (v0.5.4p5 for the mouse data, v0.6.1p1 for the black rat data), in union mode<sup>54</sup>. The data were obtained in a strand-specific assay, so the read had to be mapped to the opposite strand.

*Analysis of differential gene expression between populations.* Differentially expressed (DE) genes were identified with the *BIOCONDUCTOR*<sup>55</sup> package edgeR v3.4.0, under R 3.0.2 for the mouse data, v3.8.6 under R 3.2.3 for the black rat data<sup>56</sup>, and limma<sup>57</sup>. Data were normalized with relative log expression (RLE)-normalization factors<sup>58</sup>. Genes with fewer than 10 occurrences (reads) were filtered out and removed before the statistical analysis. This made it possible to limit the number of statistical tests and, therefore, to reduce the impact of correction for multiple testing.

Two statistical factors were included in the experimental design: the main factor included was the timing of the invasion (two levels: recent and ancient invasion), and the second factor was the site at which the rodent was caught (eight sites in total: four recently invaded sites and four anciently invaded ones). For each site, we considered two biological replicates (i.e. the two pools of individuals with a balanced sex ratio and containing a similar number of individuals).

We used three analytical strategies to identify the differentially expressed genes. For all strategies, the *p* value threshold was set to 0.05 after Benjamini-Hochberg correction for multiple testing. We first retained only one randomly chosen biological replicate per site and compared the four recently invaded sites with the four anciently invaded ones (the ‘4vs4 approach’, see Suppl. Fig. S2). This strategy made it possible to remove the second factor (site) from the analysis. As there was no particular reason to retain either the first or the second replicate for a given site, we analyzed all 256 possible combinations of replicate selection. We then considered all sites and replicates for each invasion period (recent vs. ancient invasion) together (the ‘8vs8 approach’, see Suppl. Fig. S2). The factor “site” was included in a generalized linear model (GLM) to take into account the use of two biological replicates for each site.

Considerable variability was observed between the replicates of a given site (Suppl. Fig. S3a, b), rendering the ‘4vs4 approach’ highly conservative, but potentially resulting in false positives when the ‘8vs8 approach’ was applied. Combining the results of the ‘4vs4 approach’ and the ‘8vs8 approach’ therefore seemed to be the most prudent choice, providing a maximum robustness. The genes identified with the ‘8vs8 approach’ and found to be differentially expressed in more than 85% of the 256 combinations of the ‘4vs4 approach’ were considered to be DE. The threshold of 85% was fixed based on a bar plot of the number of ‘4vs4’ comparisons in which a gene (highlighted in the ‘8vs8 approach’) was found to be differentially expressed.

Finally, we analyzed gene expression with limma and the *duplicateCorrelation* function<sup>57</sup>. We fitted a generalized linear model to the data, including site as a random effect within invasion categories. This approach corresponds to a random effects model, except that genes are constrained to share the same within-category correlation<sup>57</sup>.

**Functional and enrichment analyses.** *M. m. domesticus* and *R. rattus* genes were functionally annotated with gene ontology (GO) terms with Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) v10.5<sup>59,60</sup>. We assessed gene enrichment by performing Fisher’s exact tests on the sets of DE genes from house mouse and black rat, house mouse *M. musculus* and brown rat *R. norvegicus* GO as the background and STRING software. We corrected *p* values for multiple testing by the Benjamini and Hochberg false discovery rate (FDR) method. GO terms with a FDR *p* value <0.05 were considered to be significant. The redundancies of significantly enriched GO terms were reduced with REVIGO<sup>61</sup>, using a similarity cutoff of 0.7. We also performed enrichment analyses on biological processes and biological pathways supplied by KEGG (https://www.genom.e.jp/kegg/pathway.html), implemented in STRING.

Finally we performed a global protein network analysis based on DE genes, using a STRING database including interaction databases (including both direct (physical) and indirect (functional) associations), genetic
interactions and shared pathway interactions. This analysis revealed the interactions between key components of different pathways. Two PPI networks were constructed by mapping all DE genes and immune-related DE genes, respectively, to the STRING database, with confidence scores > 0.7 (high level). PPI networks were visualized and analyzed in CYTOSCAPE software, via the web interface of STRING. Enrichment in these protein-protein interactions was assessed by performing Fisher's exact tests with correction for multiple testing (FDR).

## Results

**Qualitative description of expression patterns.** Thirty-two transcriptome libraries were sequenced (see Table 1 for details). These libraries produced a mean of 51.6 million paired reads. The numbers of reads retained after filtering ranged from 37.2 to 76.7 M, corresponding to 92.82% of all reads. For the house mouse (3 mismatches/paired read allowed), 85.68% of the reads passing the quality filter mapped to the *M. musculus domesticus* reference genome, and 1.58% of these reads mapped to multiple regions. For the black rat (*R. rattus*, for which we allowed 9 mismatches/paired read with the *R. norvegicus* reference genome), 81.66% of the reads passing the quality filter were mapped to the reference genome of the related *R. norvegicus*, and 0.98% of these reads mapped to multiple regions of the reference genome.

**Differential expression (DE) analyses.** Along the mouse invasion route, 18 of the 17,853 filtered genes tested with edgeR and the '4vs4 approach' were systematically found to be differentially expressed between recently and anciently invaded sites. For the 16,630 filtered genes tested with the '8vs8 approach', we detected 593 differentially expressed genes, 364 of which were found in at least 85% of the 256 '4vs4 approach' comparisons. Five genes were detected by both the '4vs4' and '8vs8' approaches (*Hal*, *Rnf183*, *Serpina6*, *Wif1*, *9030619P08Rik*; see details in Suppl. Table S1a).

Along the black rat invasion route, 54 of the 13,190 filtered genes were systematically differentially expressed between recently and anciently invaded sites with the '4vs4 approach'. We found 268 DE genes among the 12,747 filtered genes analyzed with the '8vs8 approach', 83 of which were detected in at least 85% of the '4vs4' comparisons. Forty-two differentially expressed genes were detected by both the '4vs4' and '8vs8' approaches. Detailed results are shown in Suppl. Table S1b.

We focused additional analyses on these 364 and 83 genes identified as DE in the '8vs8' approach and detected in at least 85% of the '4vs4' comparisons, for the house mouse and black rat, respectively. The magnitude of expression differences differed significantly between mouse and rat, with higher 'log fold changes' observed for mouse genes (Fig. 2, Wilcoxon test, *p* = 0.01). The log fold change (LogFC) per gene ranged from -9.04 to 2.34 for the mouse and from -5.79 to 4.49 for the rat DE genes. Mean coverage ranged from -3.16 to 11.94 LogCPM (counts per million) for the mouse and -2.41 to 9.56 LogCPM for the rat DE genes.

When site was included as a random effect in the generalized linear model fitted with limma, we found no gene differentially expressed between the anciently and recently invaded sites along the house mouse invasion route. Along the rat invasion route, 12 genes were differentially expressed between recently and anciently invaded sites. All of these sites were also identified with edgeR ('8vs8 approach'), and 11 were also found in at least 85% of the 256 '4vs4 approach' comparisons.

**Functional analyses of DE genes.** For the house mouse dataset, STRING provided at least one gene ontology (GO) annotation for 345 of the 364 DE genes (Suppl. Table S1a). Of the 2098 GO annotations identified, 357 GO terms corresponding to biological processes displayed significant enrichment (FDR < 0.05) relative to the house mouse reference genome. In particular, 10 immune-related functions were significantly enriched, including inflammatory (GO:0006954), humoral (GO:0002455, GO:0006959), and innate (GO:0045087, GO:0034097) responses or activation/regulation of immunity (GO:0006952, GO:0002253, GO:0050776). We also detected enrichment for two biological processes relating to the response to wounding (GO:0009611, GO:1903035). These enriched GO terms relating to immunity corresponded to 73 DE genes, 20% of all the DE genes (see Suppl. Table S1a). The REVIGO analysis yielded 20 main clusters for the category 'Biological Process', taking redundancies into account. Inflammation and acute-phase response were part of the second most represented cluster (Suppl. Fig. S4). Using the KEGG algorithm, we identified 47 pathways as enriched. We found that 29 immune-related genes belonged to these over-represented KEGG biological pathways, the most significant of which was immune 'Complement and coagulation cascades' (FDR = 1.56e-30). Note that considering differentially expressed genes as the background for these analyses did not affect these results, as we also observed several significantly enriched GO terms corresponding to immune responses, acute-phase and humoral responses being the second main clusters for the category 'Biological Process' while taking into account redundancies with REVIGO and 'Complement and coagulation cascades' being the most significantly over represented biological pathways using KEGG algorithm.

Interestingly, the distribution of log fold-change (LogFC) for these immune-related genes corresponding to enriched GO terms or pathways was shifted towards higher values than for the 364 DE genes (Wilcoxon test, *p* = 1.52 × 10⁻³; Suppl. Fig. S5). The heatmaps built on the normalized read counts for the immune-related genes belonging to the biological processes (and pathways) found to be overrepresented revealed a downregulation of these genes at all anciently invaded sites, and an upregulation for all recently invaded sites except Aere Lao (Fig. 3, Suppl. Fig. S6). In this last population, the expression pattern observed was similar to that at anciently invaded sites.

Finally, using STRING, we showed that the protein-protein interaction (PPI) network based on the 364 DE genes was significantly enriched (*p* < 1.0e-16, 731 edges found, 54 expected), indicating a larger number of interactions between these proteins than would be expected for a set of proteins of similar size drawn at random from the sequences of the mouse genome. Two main clusters were detected, one including cytochrome P₄₅₀ and
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Significantly differentially expressed genes (orange dots) and immune related genes (red dots) between anciently and recently invaded sites along the invasion routes of (a) the house mouse and (b) the black rat. Only the 73 immune related genes belonging to biological pathways found to be significantly enriched are represented for the house mouse. Vertical lines indicate 1 log fold change (Log FC). The x-axis indicates genes that are down- (negative values) and up-regulated (positive values) in anciently compared to recently invaded areas.</figcaption></figure>

UDP glucuronosyltransferase proteins, and the other consisting mostly of immune-related proteins, including fibrinogen, serin peptidase inhibitor (serpin), apolipoprotein and complement proteins. The PPI interaction network built upon the restricted dataset including the 73 immune-related genes belonging to the biological
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Heatmap of the differentially expressed (DE) genes between the anciently and recently invaded sites of the house mouse (*M. musculus domesticus*). The normalized read counts for the expressed genes are shown. For clarity, the heatmap was built in R using heatmap.2 for 29 immune related genes belonging to over-represented KEGG biological pathways. The genes (rows) and samples (columns) were clustered using dendrograms built with Ward distance and hierarchical clustering. Anciently and recently invaded sites are indicated using (AI) and (RI), respectively.</figcaption></figure>

processes shown to be overrepresented also displayed significant enrichment ($p$ < 1.0e-16, 561 edges found, 13 expected). The network had a hierarchical structure, with the main cluster including the fibrinogen and serine peptidase inhibitor (serpin) proteins, and two other clusters of highly connected proteins, one including apolipoprotein A and haptoglobin, and the other, complement proteins (Fig. 4).

For the *R. rattus* dataset, GO annotations were available for 75 of the 83 DE genes identified (Suppl. Table S1b). Eight of these annotated DE genes related to immune biological processes (i.e., 10.66% of all DE genes; Suppl. Table S1b). We identified no significantly enriched biological process or pathway, even when redundancies were considered with REVIGO. The protein-protein interaction network built from STRING displayed no significant enrichment ($p$ = 0.09). Similar conclusions were obtained when the dataset was restricted to the 12 genes found to be differentially expressed, with site considered as a random effect within invasion category.

Discussion

We present here a comparative transcriptomics approach based on the sequencing of spleen RNA, for two invasive rodent species studied along their invasion routes in Senegal. We aimed to detect adaptive divergence between recently and long-established rodent populations that could at least partly account for invasion success. More specifically, we investigated the differential expression of immune-related genes, to determine whether the patterns of expression observed supported either the EICA or the EICA-refined hypotheses.

Invasion success is thought to be dependent on a shift in phenotypic traits of advantage at invasion fronts, where species are likely to be exposed to new environmental conditions⁹,¹⁰,⁶². Using comparative transcriptomics in wild populations, we show here that a phenotypic shift⁶³ has occurred along the mouse and rat invasion routes in Senegal. Indeed, we detected variations in gene expression levels between recently and long-established rodent populations. Climatic variations and changes in habitat types covary with invasion categories along the invasion routes of both mouse and rat⁴⁴, and all the processes that may explain the phenotypic shift may be influenced by environmental factors. This phenotypic shift may reflect stochastic changes due to population history, including founder events and genetic drift⁹, or non-adaptive responses to environmental cues. Alternatively, this phenotypic shift may result from adaptive processes, including phenotypic plasticity⁶⁴ or natural selection⁶⁵. Invasive populations have a lower genetic diversity than their source populations, due to founder events, but may still evolve rapidly through natural selection from the remaining standing genetic variation or de novo mutations⁶⁶. Several experimental immune challenges have been performed, in controlled conditions, on cane toads from parents collected from long-established populations and invasion front populations from Australia. These experiments showed a rapid evolution of immunity in response to the selection pressures occurring during biological invasions⁶⁷,⁶⁸. Variations of gene expression at invasion front would also be mediated by phenotypic plasticity⁶⁹,⁷⁰.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Protein-protein interaction network drawn from the 73 immune genes found to belong to enriched biological processes related to immunity. Nodes correspond to proteins, the thickness of edge network indicates the strength of data support, and the minimum required interaction score was set to 0.7 (high confidence level). The plain circle indicates the alipoprotein A-aloa—and haptoglobin-Hp-protein. The dashed circle indicates the complement proteins (C8). Fibrinogen (Fg) and serine peptidase inhibitor (serpin) proteins are distributed within the highly connected central part of the network.</figcaption></figure>

Statistical analyses of transcriptomic data for natural populations are not trivial. We therefore used a combination of three analytical strategies to identify differentially expressed genes. It was tempting to consider site as a random effect within invasion category, to account for the inherently nested structure of our experimental design. However, the strong variability of gene expression between sites made it difficult to detect differentially expressed genes with this strategy. Several lines of evidence suggested that sites within invasion categories could not really be considered as true replicates and that timing of invasion might be variable, according to the history of connections between the site concerned and colonial cities or areas previously invaded by exotic rodents. For instance, colonial settlements were established at some anciently invaded sites (NDB, MBA), but not at others (DAG, THL). These sites without colonial settlements may have been invaded a little later by the house mouse. Historical data also suggest that the Aere Lao site may have been colonized by the house mouse earlier than other sites of the recently invaded category. Indeed, Aere Lao was not connected to the road network until the 1970–1980s, when a weekly rural market was established at this site, creating commercial links with cities already invaded by the house mouse (e.g. Saint Louis, Dakar). This may have favored the introduction of the house mouse at this site 10 or 20 years before this rodent reached the neighboring villages⁷¹, but several dozens of years after long-standing invasion sites. This history may explain why immune phenotypes in Aere Lao were more similar to those observed in long-established populations from western Senegal for both analysis strategies, and in terms of immune effectors (complement; haptoglobin), as previously shown by immune challenges⁴².

In the house mouse, 20.0% of the genes found to be differentially expressed with a combination of the 4vs4 and 8vs8 strategies were annotated with functions relating to immunity. We also found that some functional categories of genes were overrepresented among these differentially expressed genes. In particular, we obtained strong evidence for a significantly biased variation of immune gene expression towards an upregulation of
serpins, alipoproteins, complement and pro-inflammatory cascades at recently invaded sites. These results were consistent with our previous functional studies<sup>42</sup>, in which immune challenges revealed that antibody-mediated (natural antibodies and complement) and inflammatory (haptoglobin) responses were stronger in mice sampled from recently invaded sites than in those sampled in anciently invaded ones. Our findings therefore support the hypothesis that parasite-mediated selection could contribute to phenotypic differentiation along invasion routes, potentially in interaction with other environmental factors. However, they are not consistent with the EICA hypothesis or the EICA-refined hypothesis<sup>19</sup>. Despite the enemy release patterns detected for mouse helminths in Senegal from sites of long-term invasion to recently invaded sites (decrease in specific richness and overall prevalence, see<sup>44</sup>), we found no evidence of a global decrease in immune-related gene expression at recently invaded sites (the pattern expected according to the EICA hypothesis). We also observed no change in immune-related gene expression that could be attributed to the relative energy costs of immune pathways (the pattern expected according to the EICA-refined hypothesis). Indeed, both inflammation (considered costly in terms of energy) and complement (considered cost-effective in terms of energy) pathways were upregulated<sup>20</sup>. The upregulation of all the immune genes found to be differentially expressed at sites recently invaded by the house mouse therefore strongly supported the alternative hypothesis of an increase in overall infection risk at recently invaded sites. More specifically, the biological processes and pathways displaying significant enrichment, and the protein-protein interaction network, highlighted innate immune responses, which are essential for front-line defense against infections, including bacterial infections. Complement is a protein complex activated immediately after the detection of infection that can induce bacterial lysis<sup>72</sup>. Serine protease inhibitors (serpin) are involved in modulating many physiological functions, including infection and inflammation<sup>73,74</sup>. Apolipoprotein A is involved in lipid transport, and in diverse antimicrobial activities, including the inhibition of bacterial growth and the suppression of viral replication<sup>75</sup>. As such, an increase in expression of the genes encoding these proteins should make it possible for the invading organism to fight effectively against the new bacteria encountered at recently invaded sites<sup>43</sup>. Interestingly, some of these proteins may also have anti-inflammatory activity, through the secretion of particular cytokines (e.g. haptoglobin, serpinG1, C1q, C3a<sup>72,73,76,77</sup>). This function may be particularly important at recently invaded sites, for decreasing the immunopathological effects of inflammation in response to exposure to new parasites<sup>78</sup>. Moreover, the overexpression at invasion fronts of genes encoding anti-inflammatory proteins may counteract the helminth release observed at recently invaded sites along the invasion route of the house mouse<sup>44</sup>. Indeed, as several helminth infections are known to activate regulatory responses that reduce inflammation<sup>79</sup>, helminth release could exacerbate the risk of diseases associated with inflammatory responses.

Similar rapid changes in immunity have been observed in other invading organisms, but these studies revealed a wide array of patterns from which it was not possible to draw any firm conclusions about the potential generic/specific nature of EICA hypotheses<sup>21</sup>. Further studies are required, on a wide range of immune effectors, based on 'omics methods'<sup>80</sup>, to improve our understanding of the relationships between biological invasion and immunity.

Concerning *R. rattus*, the observed changes in gene expression levels were smaller than those observed in the mouse, and no particular pathway (including immune-related pathways) was significantly overrepresented. The absence of a clear pattern suggested that stochastic events might have strongly shaped the phenotypic variations observed along the route of *R. rattus* invasion. However, caution is required in interpreting our results, which may be due to the limitations of bioinformatics methods (transcriptome annotation was less efficient for *R. rattus* than for *M. musculus domesticus*, as we used the *R. norvegicus* transcriptome as a reference). The eight immune-related genes found to be differentially expressed encode proteins that interact with bacteria (*Rnase6*, *Ptprz1*), viruses (*Ltc4s*, *Hspa1b*) or parasites (*Gfi1b*). Some were downregulated at recently invaded sites (e.g., *Rnase6*), whereas others were strongly upregulated (e.g., *Ptprz1*). No specific pattern has yet emerged concerning potential variations in infection risks. Variations of immune phenotypes, as reflected in gene expression data for the spleen, were, therefore, unlikely to be the main driver of invasion success in the black rat (assuming that invasion success depends on a rapid shift in phenotypic traits in newly colonized areas rather than pre-adaptation in the range of long-standing invasion). Interestingly, we also found, in a previous functional immunity study, that differences in immune responses between anciently and recently invaded sites were less marked for the rat than for the mouse<sup>42</sup>. It was also very difficult to identify any other biological processes that might reflect life history trait adaptation with regard to invasion. The processes found to be associated with differentially expressed genes covered a wide range of functions, including the clustering of sodium voltage channels, the downregulation of cell-matrix adhesion and protein dephosphorylation and the response to cyclic adenosine monophosphate (cAMP). This last function is of potential interest with respect to rat invasion success, because cAMP plays a key role in regulating insulin and glucagon secretion<sup>81</sup>. Variations in cAMP gene expression may, therefore, mediate differences in response to stressful situations, including starvation or fight-or-flight response<sup>82</sup>. It would be interesting to analyze whether such differences could result in different behavioral phenotypes between rats trapped at anciently invaded sites (expected to have a poor performance and high stress response in new environments) and those caught at sites of recent invasion (expected to have a high response capacity in new environments)<sup>83</sup>. It therefore appears particularly important to perform additional transcriptomics analyses on other organs and tissues, to identify the phenotypic changes and ecoevolutionary processes linked to the invasion success of the black rat. The brain, in which key genes underlying behavioral invasion syndrome would be expected to be expressed, may be a relevant organ for future studies<sup>84</sup>.

In conclusion, our work revealed changes in transcriptomic profiles along invasion routes for both the house mouse and the black rat. It is likely that different processes, including colonization history and/or alternative mechanisms by which species adapt to new environments, mediated the invasion success of these two rodent species. The patterns observed could potentially be driven by an increase in infection risks at recently invaded sites for the house mouse and stochastic changes for the black rat. It would be interesting to take into account the potential variability of these patterns with regard to the tissue or organ targeted. Different results might have
been obtained if we had targeted other lymphoid organs, lymphatic tissues or non-immune targets. Additional genomic studies and experimental work<sup>25,28,67,68,85</sup> are required to determine whether the observed differences in gene expression were driven by phenotypic plasticity or directional selection during or after invasion, or whether they reflect the colonization history of the rodents concerned. Moreover, it is also important to determine whether these changes in phenotypic traits influence ecological dynamics (e.g.<sup>86</sup>) and, in turn, invasion success.

## Data availability

All sequence data, including raw read sequences and assemblies have been deposited on GEO. Gene ontologies for DE genes are provided in Supplementary Tables 1 and 2. Differential expression data are available from the corresponding author on request.

Received: 21 May 2019; Accepted: 28 June 2020

Published online: 26 October 2020

## References

1. Kamenova, S. et al. in *Networks of Invasion: A Synthesis of Concepts* Vol. 56 *Adv. Ecol. Res.* (eds Bohan, D. A., Dumbrell, A. J., & Massol, F.) 85–182 (2017).

2. Kumschick, S. et al. Ecological impacts of alien species: quantification, scope, caveats, and recommendations. *Bioscience* **65**, 55–63. https://doi.org/10.1093/biosci/biu193 (2014).

3. Daszak, P., Cunningham, A. A. & Hyatt, A. D. Emerging infectious diseases of wildlife—threats to biodiversity and human health. *Science* **287**, 443–449. https://doi.org/10.1126/science.287.5452.443 (2000).

4. Bradshaw, C. J. A. et al. Massive yet grossly underestimated global costs of invasive insects. *Nat. Commun.* https://doi.org/10.1038/s41487-019-00211-x (2016).

5. Martin, L. B., Coon, C. A. C., Liebl, A. L. & Schrey, A. W. Surveillance for microbes and range expansion in house sparrows. *Proc. R. Soc. Lond. B* **281**, 20132690. https://doi.org/10.1098/rspb.2013.2690 (2013).

6. Lindstrom, T., Brown, G. P., Sisson, S. A., Phillips, B. L. & Shine, R. Rapid shifts in dispersal behavior on an expanding range edge. *Proc. Natl. Acad. Sci. USA* **110**, 13452–13456. https://doi.org/10.1073/pnas.1303157110 (2013).

7. Façon, B. et al. A general eco-evolutionary framework for understanding bioinvasions. *Trends Ecol. Evol.* **21**, 130–135. https://doi.org/10.1016/j.tree.2005.10.012 (2006).

8. Hufbauer, R. A. et al. Anthropogenically induced adaptation to invade (AIAI): contemporary adaptation to human-altered habitats within the native range can promote invasions. *Evol. Appl.* **5**, 89–101. https://doi.org/10.1111/j.1752-4571.2011.00211.x (2012).

9. Colautti, R. I. & Lau, J. A. Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation. *Mol. Ecol.* **24**, 1999–2017. https://doi.org/10.1111/mec.13162 (2015).

10. Chuang, A. & Peterson, C. R. Expanding population edges: theories, traits, and trade-offs. *Glob. Change Biol.* **22**, 494–512. https://doi.org/10.1111/gcb.13107 (2016).

11. Selechnik, D., Rollins, L. A., Brown, G. P., Kelehear, C. & Shine, R. The things they carried: the pathogenic effects of old and new parasites following the intercontinental invasion of the Australian cane toad (*Rhinella marina*). *Int. J. Parasitol.* **6**, 375–385. https://doi.org/10.1016/j.ijppaw.2016.12.001 (2017).

12. Torchin, M. E., Lafferty, K. D., Dobson, A. P., McKenzie, V. J. & Kuris, A. M. Introduced species and their missing parasites. *Nature* **421**, 628–630. https://doi.org/10.1038/nature01346 (2003).

13. Torchin, M. E. & Mitchell, C. E. Parasites, pathogens, and invasions by plants and animals. *Front. Ecol. Environ.* **2**, 183–190. https://doi.org/10.2307/3868313 (2004).

14. Dunn, A. M. in *Advances in Parasitology, Vol 68: Natural History of Host-Parasite Interactions* Vol. 68 (ed Webster, J. P.) 161–184 (2009).

15. O'Brien, V. A. et al. An enzootic vector-borne virus is amplified at epizootic levels by an invasive avian host. *Proc. R. Soc. Lond. B* **278**, 239–246 (2011).

16. Perkins, S. E., White, T. A., Pascoe, E. L. & Gillingham, E. L. Parasite community dynamics in an invasive vole—from focal introduction to wave front. *Int. J. Parasitol.* **6**, 412–419. https://doi.org/10.1016/j.ijppaw.2017.07.005 (2017).

17. Yang, C. C. et al. Loss of microbial (pathogen) infections associated with recent invasions of the red imported fire ant *Solenopsis invicta*. *Biol. Invasions* **12**, 3307–3318. https://doi.org/10.1007/s10530-010-9724-9 (2010).

18. Blossey, B. & Notzold, R. Evolution of increased competitive ability in invasive nonindigenous plants—a hypothesis. *J. Ecol.* **83**, 887–889. https://doi.org/10.2307/2261425 (1995).

19. Lee, K. A. & Klasing, K. C. A role for immunology in invasion biology. *Trends Ecol. Evol.* **19**, 523–529. https://doi.org/10.1016/j.tree.2004.07.012 (2004).

20. Lee, K. A. Linking immune defenses and life history at the levels of the individual and the species. *Integr. Comp. Biol.* **46**, 1000–1015. https://doi.org/10.1093/icb/icl049 (2006).

21. Cornet, S., Brouat, C., Diagne, C. A. & Charbonnel, N. EcoImmunology and bioinvasion: revisiting the EICA hypotheses. *Evol. Appl.* **9**, 952–962. https://doi.org/10.1111/eva.12406 (2016).

22. Martin, L. B., Alam, J. L., Imboma, T. & Liebl, A. L. Variation in inflammation as a correlate of range expansion in Kenyan house sparrows. *Oecologia* **164**, 339–347. https://doi.org/10.1007/s00442-010-1654-9 (2010).

23. Bernardi, G., Azzurro, E., Golani, D. & Miller, M. R. Genomic signatures of rapid adaptive evolution in the bluespotted cornetfish, a Mediterranean Lessepsian invader. *Mol. Ecol* **25**, 3384–3396. https://doi.org/10.1111/mec.13682 (2016).

24. Vera, M., Diez-del-Molino, D. & Garcia-Marin, J. L. Genomic survey provides insights into the evolutionary changes that occurred during European expansion of the invasive mosquito fish (*Gambusia holbrooki*). *Mol. Ecol* **25**, 1089–1105. https://doi.org/10.1111/mec.13545 (2016).

25. Hodgins, K. A., Lai, Z., Nurkowski, K., Huang, J. & Rieseberg, L. H. The molecular basis of invasiveness: differences in gene expression of native and introduced common ragweed (*Ambrosia artemisiifolia*) in stressful and benign environments. *Mol. Ecol* **22**, 2496–2510. https://doi.org/10.1111/mec.12179 (2013).

26. White, T. A., Perkins, S. E., Heckel, G. & Searle, J. B. Adaptive evolution during an ongoing range expansion: the invasive bank vole (*Myodes glareolus*) in Ireland. *Mol. Ecol* **22**, 2971–2985. https://doi.org/10.1111/mec.12343 (2013).

27. Todd, E. V., Black, M. A. & Gemmell, N. J. The power and promise of RNA-seq in ecology and evolution. *Mol. Ecol* **25**, 1224–1241. https://doi.org/10.1111/mec.13526 (2016).

28. Alvarez, M., Schrey, A. W. & Richards, C. L. Ten years of transcriptomics in wild populations: what have we learned about their ecology and evolution? *Mol. Ecol* **24**, 710–725. https://doi.org/10.1111/mec.13055 (2015).

29. Rius, M. & Darling, J. A. How important is intraspecific genetic admixture to the success of colonising populations? *Trends Ecol. Evol.* **29**, 233–242. https://doi.org/10.1016/j.tree.2014.02.003 (2014).

30. Fraser, B. A., Ramnarine, I. W. & Neff, B. D. Temporal variation at the Mhc class IIB in wild populations of the guppy (*Poecilia reticulata*). *Evolution* **64**, 2086–2096. https://doi.org/10.1111/j.1558-5646.2010.00965.x (2010).
31. Oleksiak, M. F., Churchill, G. A. & Crawford, D. L. Variation in gene expression within and among natural populations. *Nat. Genet.* **32**, 261–266. https://doi.org/10.1038/ng983 (2002).

32. Whitehead, A. & Crawford, D. L. Variation within and among species in gene expression: raw material for evolution. *Mol. Ecol.* **15**, 1197–1211. https://doi.org/10.1111/j.1365-294X.2006.02868.x (2006).

33. Whitehead, A. & Crawford, D. L. Neutral and adaptive variation in gene expression. *Proc. Natl. Acad. Sci. USA* **103**, 5425–5430. https://doi.org/10.1073/pnas.0507648103 (2006).

34. Rollins, L. A., Richardson, M. F. & Shine, R. A genetic perspective on rapid evolution in cane toads (*Rhinella marina*). *Mol. Ecol.* **24**, 2264–2276. https://doi.org/10.1111/mec.13184 (2015).

35. Vogel, H. J., Schmidtberg, H. & Vilcinskas, A. Comparative transcriptomics in three ladybird species supports a role for immunity in invasion biology. *Dev. Comp. Immunol.* **67**, 452–456. https://doi.org/10.1016/j.dci.2016.09.015 (2017).

36. Selechnik, D., Richardson, M. F., Shine, R., Brown, G. P. & Rollins, L. A. Immune and environment-driven gene expression during invasion: an eco-immunological application of RNA-Seq. *Ecol. Evol.* **9**, 6708–6721. https://doi.org/10.1002/ece3.5249 (2019).

37. Aplin, K. P. et al. Multiple geographic origins of commensalism and complex dispersal history of black rats. *PLoS ONE* **6**, e26357. https://doi.org/10.1371/journal.pone.0026357 (2011).

38. Bonhomme, F. et al. Genetic differentiation of the house mouse around the Mediterranean basin: matrilinear footprints of early and late colonization. *Proc. R. Soc. Lond. B* **278**, 1034–1043. https://doi.org/10.1098/rspb.2010.1228 (2011).

39. Dalecky, A. et al. Range expansion of the invasive house mouse *Mus musculus domesticus* in Senegal, Western Africa: a three decades synthesis of trapping data, 1983–2014. *Mammal Rev.* **45**, 176–190. https://doi.org/10.1111/mam.12043 (2015).

40. Konecny, A. et al. Invasion genetics of the introduced black rat (*Rattus rattus*) in Senegal, West Africa. *Mol. Ecol.* **22**, 286–300. https://doi.org/10.1111/mec.12112 (2013).

41. Lippens, C. et al. Genetic structure and invasion history of the house mouse (*Mus musculus domesticus*) in Senegal, West Africa: a legacy of colonial and contemporary times. *Heredity* **119**, 64–75. https://doi.org/10.1038/hdy.2017.18 (2017).

42. Diagne, C. A. et al. Contemporary variations of immune responsiveness during range expansion of two invasive rodents in Senegal. *Oikos* **126**, 435–446. https://doi.org/10.1111/oik.03470 (2017).

43. Diagne, C. A. et al. Ecological and sanitary impacts of bacterial communities associated to biological invasions in African commensal rodent communities. *Nat. Sci. Rep.* **7**, 14995. https://doi.org/10.1038/s41598-017-14880-1 (2017).

44. Diagne, C. A. et al. Parasites and invasions: changes in gastrointestinal helminth assemblages in invasive and native rodents in Senegal. *Int. J. Parasitol.* **46**, 857–869. https://doi.org/10.1016/j.ijpara.2016.07.007 (2016).

45. Frank, S. A. Immune response to parasitic attack: evolution of a pulsed character. *J. Theor. Biol.* **219**, 281–290. https://doi.org/10.1006/jtbi.2002.3122 (2002).

46. Brown, G. P., Shilton, C., Phillips, B. L. & Shine, R. Invasion, stress, and spinal arthritis in cane toads. *Proc. Natl. Acad. Sci. USA* **104**, 17698–17700. https://doi.org/10.1644/10-MAMM-F-355.1 (2011).

47. Sikes, R. S. & Gannon, W. L. Guidelines of the American Society of Mammalogists for the use of wild mammals in research. *J. Mammal.* **92**, 235–253. https://doi.org/10.1644/10-MAMM-F-355.1 (2011).

48. Granjon, L. & Duplantier, J. M. *Les rongeurs de l'Afrique sahéro-soudanienne* (Publications scientifiques du Muséum, 2009).

49. Nussey, D. H., Watt, K., Pilkington, J. G., Zamoyska, R. & McNeilly, T. N. Age-related variation in immunity in a wild mammal population. *Aging Cell* **11**, 178–180. https://doi.org/10.1111/j.1474-9726.2011.00771.x (2012).

50. Schurch, N. J. et al. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? *RNA* **22**, 839–851. https://doi.org/10.1261/rna.053959.115 (2016).

51. Kim, D. et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. *Genome Biol.* **14**, R36. https://doi.org/10.1186/gb-2013-14-4-r36 (2013).

52. Langmead, B. & Salzberg, S. Fast gapped-read alignment with Bowtie 2. *Nat. Methods* **9**, 357–359. https://doi.org/10.1038/nmeth.1923 (2012).

53. Robins, J. H. et al. Dating of divergences within the Rattus genus phylogeny using whole mitochondrial genomes. *Mol. Phylogenet. Evol.* **49**, 460–466. https://doi.org/10.1016/j.jmpev.2008.08.001 (2008).

54. Anders, S., Pyl, T. P. & Huber, W. HTSeq—a Python framework to work with high-throughput sequencing data. *Bioinformatics* **31**, 166–169. https://doi.org/10.1101/002824 (2015).

55. Gentleman, R. C. et al. Bioconductor: open software development for computational biology and bioinformatics. *Genome Biol.* **5**, R80. https://doi.org/10.1186/gb-2004-5-10-r80 (2004).

56. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. A bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* **26**, 139–140. https://doi.org/10.1093/bioinformatics/btp616 (2010).

57. Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. *Nucl. Acids Res.* **43**, e47. https://doi.org/10.1093/nar/gkv007 (2015).

58. Simon, A. & Huber, W. Differential expression analysis for sequence count data. *Genome Biol.* **11**, R106. https://doi.org/10.1186/gb-2010-11-10-r106 (2010).

59. Szklarczyk, D. et al. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. *Nucl. Acids Res.* **45**, D362–D368. https://doi.org/10.1093/nar/gkw937 (2017).

60. Szklarczyk, D. et al. STRING v10: protein–protein interaction networks, integrated over the tree of life. *Nucl. Acids Res.* **43**, D447–D452. https://doi.org/10.1093/nar/gku1003 (2015).

61. Supek, F., Bošnjak, M., Škunca, N. & Smuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. *PLoS ONE* **6**, e21800. https://doi.org/10.1371/journal.pone.0021800 (2011).

62. Merrick, M. J. & Koprowski, J. L. Altered natal dispersal at the range periphery: the role of behavior, resources, and maternal condition. *Ecol. Evol.* **7**, 58–72. https://doi.org/10.1002/ece3.2612 (2017).

63. DeBiasse, M. B. & Kelly, M. W. Plastic and evolved responses to global change: what can we learn from comparative transcriptomics? *J. Heredity* **107**, 71–81. https://doi.org/10.1093/jhered/esv073 (2016).

64. Whitehead, A., Triant, D. A., Champlin, D. & Nacci, D. Comparative transcriptomics implicates mechanisms of evolved pollution tolerance in a killifish population. *Mol. Ecol.* **19**, 5186–5203. https://doi.org/10.1111/j.1365-294X.2010.04829.x (2010).

65. Barrett, S. C. H. Foundations of invasion genetics: the Baker and Stebbins legacy. *Mol. Ecol.* **24**, 1927–1941. https://doi.org/10.1111/mec.13014 (2015).

66. Dlugosch, K. M., Anderson, S. R., Braasch, J., Cang, F. A. & Gillette, H. D. The devil is in the details: genetic variation in introduced populations and its contributions to invasion. *Mol. Ecol.* **24**, 2095–2111. https://doi.org/10.1111/mec.13183 (2015).

67. Llewellyn, D., Thompson, M. B., Brown, G. P., Phillips, B. L. & Shine, R. Reduced investment in immune function in invasion-front populations of the cane toad (*Rhinella marina*) in Australia. *Biol. Invasions* **14**, 999–1008. https://doi.org/10.1007/s10530-011-0135-3 (2012).

68. Brown, G. P., Phillips, B. L., Dubey, S. & Shine, R. Invader immunology: invasion history alters immune system function in cane toads (*Rhinella marina*) in tropical Australia. *Ecol. Lett.* **18**, 57–65. https://doi.org/10.1111/ele.12390 (2015).

69. Dlugosch, K. M. & Parker, I. M. Invading populations of an ornamental shrub show rapid life history evolution despite genetic bottlenecks. *Ecol. Lett.* **11**, 701–709. https://doi.org/10.1111/j.1461-0248.2008.01181.x (2008).

70. Chevin, L. M. & Lande, R. Adaptation to marginal habitats by evolution of increased phenotypic plasticity. *J. Evol. Biol.* **24**, 1462–1476. https://doi.org/10.1111/j.1420-9101.2011.02279.x (2011).
71. Ninot, O. *Vie de relations, organisation de l'espace et développement en Afrique de l'Ouest : la région de Tambacounda au Sénégal*. Ph.D. thesis, Rouen University (2003).

72. Merle, N. S., Noe, R., Halbwachs-Mecarelli, L., Fremeaux-Bacchi, V. & Roumenina, L. T. Complement system part II: role in immunity. *Front. Immunol.* **6**, 257. https://doi.org/10.3389/fimmu.2015.00257 (2015).

73. Bao, J. et al. Serpin functions in host-pathogen interactions. *PeerJ* **6**, e4557. https://doi.org/10.7717/peerj.4557 (2018).

74. Mangan, M. S. J., Kaiserman, D. & Bird, P. I. The role of serpins in vertebrate immunity. *Tissue Antigens* **72**, 1–10. https://doi.org/10.1111/j.1399-0039.2008.01059.x (2008).

75. Wang, W., Qu, Q. & Chen, J. Identification, expression analysis, and antibacterial activity of Apolipoprotein A-I from amphioxus (*Branchiostoma belcheri*). *Comp. Biochem. Physiol.* **238**, 110329. https://doi.org/10.1016/j.cbpb.2019.110329 (2019).

76. Huntoon, K. M. et al. The acute phase protein haptoglobin regulates host immunity. *J. Leukoc. Biol.* **84**, 170–181. https://doi.org/10.1189/jlb.0208100 (2008).

77. Burger, D. & Dayer, J. M. High-density lipoprotein-associated apolipoprotein AI: the missing link between infection and chronic inflammation?. *Autoimmun. Rev.* **1**, 111–117. https://doi.org/10.1016/S1568-9972(01)00018-0 (2002).

78. Sorci, G. & Faivre, B. Inflammation and oxidative stress in vertebrate host-parasite systems. *Philos. Trans. R. Soc. B* **364**, 71–83. https://doi.org/10.1098/rstb.2008.0151 (2009).

79. McKay, D. M. The beneficial helminth parasite?. *Parasitology* **132**, 1–12. https://doi.org/10.1017/S003118200500884X (2006).

80. Robertson, S., Bradley, J. E. & MacColl, A. D. Measuring the immune system of the three-spined stickleback—investigating natural variation by quantifying immune expression in the laboratory and the wild. *Mol. Ecol. Res.* **16**, 701–713. https://doi.org/10.1111/1755-0998.12497 (2016).

81. Tengholm, A. & Gylfe, E. cAMP signalling in insulin and glucagon secretion. *Diabetes Obes. Metabol.* **19**, 42–53. https://doi.org/10.1111/dom.12993 (2017).

82. Jones, B. J., Tan, T. & Bloom, S. R. Minireview: glucagon in stress and energy homeostasis. *Endocrinology* **153**, 1049–1054. https://doi.org/10.1210/en.2011-1979 (2012).

83. Sih, A., Cote, J., Evans, M. R., Fogarty, S. & Pruitt, J. Ecological implications of behavioural syndromes. *Ecol. Lett.* **15**, 278–289. https://doi.org/10.1111/j.1461-0248.2011.01731.x (2012).

84. Bengston, S. E. et al. Genomic tools for behavioural ecologists to understand repeatable individual differences in behaviour. *Nat. Ecol. Evol.* **2**, 944–955. https://doi.org/10.1038/s41559-017-0411-4 (2018).

85. Martin, L. B. et al. Costs of immunity and their role in the range expansion of the house sparrow in Kenya. *J. Exp. Biol.* **220**, 2228–2235. https://doi.org/10.1242/jeb.154716 (2017).

86. Schoener, T. W. The newest synthesis: understanding the interplay of evolutionary and ecological dynamics. *Science* **331**, 426–429. https://doi.org/10.1126/science.1193954 (2011).

## Acknowledgements

We thank Khalilou Bâ, Mamoudou Diallo, Mamadou Kane, Laurent Granjon, Youssoupha Niang, Aliou Sow and Philippe Gauthier for their help in collecting rodents, and J.-M. Duplantier for invaluable advice regarding field work and rodent invasion history. We are indebted to all the Senegalese people who allowed us to trap rodents in their homes. This work was supported by the ANR ENEMI project (ANR-11-JSV7-0006), and by the French Embassy in Senegal, which provided the funding for PhD scholarships. This preprint has been reviewed and recommended by Peer Community in Ecology (https://dx.doi.org/10.24072/pci.ecology.100011).

## Author contributions

C.B., M.G. and N.C. conceived the study. A.D., C.B., C.A.D., N.C., participated in the population sampling (other people involved have already co-authored other published papers). M.G., C.T. and A.L. performed the RNA extraction. H.P. and D.S. constructed the RNA libraries and performed the HiSeq sequencing. H.P. and S.R. performed bioinformatics and statistical analyses for D.E. gene detection. N.C. conducted functional and enrichment analyses. N.C. wrote the manuscript, with contribution gathered from all co-authors.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41598-020-75060-2.

Correspondence and requests for materials should be addressed to N.C.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020