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

Background Biological invasions are major anthropogenic changes associated with threats to biodiversity and health. What determines the successful establishment of introduced populations still remains unsolved. Here we explore the appealing assertion that invasion success relies on immune phenotypic traits that would be advantageous in recently invaded sites. Results We compared gene expression proﬁles 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. Transcriptome analyses revealed respectively 364 and 83 differentially expressed genes along the mouse and rat invasion routes. Among them, 20.0% and 10.6% were annotated with functions related to immunity. All immune-related genes detected along the mouse invasion route were over-expressed in recently invaded sites. Genes of the complement activation pathway were over-represented. Results were less straightforward when considering the black rat as no particular immunological process was over-represented. Conclusions We revealed changes in transcriptome profiles along invasion routes. Patterns differed between invaders. They could potentially be driven by increased infection risks in recently invaded sites for the house mouse and stochastic events associated with colonization history for the black rat. These results provide a first step in identifying the immune ecoevolutionary processes potentially involved in invasion success.


Background
Biological invasions are considered as a major component of global change growing threats to biodiversity 1 , ecosystem functioning 2 and human health 3 . They may also lead to huge negative economic impacts 4 . There is hence an urgent need for understanding what determines invasion success (i.e. the successful introduction, establishment and spread of species outside their native range) to design effective prevention, control and management strategies.
From an eco-evolutionary perspective, invasion success may rely on pre-adaptation within the original range 5,6 or on the rapid evolution of phenotypic traits that would be advantageous in newly colonized areas 5,7 . Some supports for this latter process come from the identification of phenotypic variation along invasion gradients, what has been evidenced in various animal species 7 . In this context, a corpus of appealing hypotheses has focused on traits that are related with the ability to cope with novel environments 8 . In particular, important changes of parasite pressures 8,9 can occur on introduced species along invasion gradients (e.g. enemy release, spill-over of native parasites 10-12 ). Thus, invasion success could rely on the efficiency of immune system regulation, either through the reallocation of the energetic resources dedicated to defences towards other life history traits favouring range expansion (including dispersal or reproduction, evolution of increased competitive ability (EICA) hypothesis 13 ), or by the modulation of defence strategies to limit the costs of immunity (EICA-refined hypothesis 14 ). A recent analysis reviewing empirical studies has shown that in most cases, variations in immune responses were observed during invasions although no general pattern could be emphasized to corroborate the EICA hypotheses 15 . One of the problems pointed out by the authors was the lack of studies providing a large array of immune responses to properly assess immunocompetence and estimate immune trade-offs. Unfortunately, immune phenotyping often relies on a small number of immune effectors, because of the low quantity 4 of materials available (blood, tissue), or the difficulty to get immune kits for non-model species.
The recent advent of 'omics' technologies provides opportunities to explore a wide set of genes and pathways involved in the response to novel conditions experienced at invasion fronts. Genome scans have been successfully performed to study some invasion cases (e.g. [16][17][18] ). They have enabled to detect signatures of contemporary adaptation associated with invasion, some of them being related to immunity 19 . Nevertheless, this genomic approach based on the detection of single nucleotide polymorphism outliers does not enable to identify the molecular processes driving the rapid adaptation that might occur within few generations.
It also prevents from catching the importance of phenotypic plasticity in invasion success.
Transcriptomics, the analysis of gene expression at a genome wide scale, is a complementary approach that might allow filling these gaps and deciphering the molecular mechanisms that underlie phenotypic changes 20,21 . Indeed, gene expression is an essential mechanism for rapid acclimatization or adaptation to novel environments [22][23][24][25][26] . Transcriptomics has hence been applied to identify genes and gene regulatory pathways that contribute to phenotypic variation at invasion fronts. Evidence for transcriptomic divergence along invasion gradients have already been reported 18 ,27 , sometimes emphasizing the importance of immune regulation in invasion success 28 . However there is yet no study combining analyses of immune gene expression and potential variations in parasite pressures along invasion routes. These two elements are essential to test the EICA hypotheses and further understand the interactions between immune system regulation and invasion success 15 .
Here, we examined transcriptional profiles along invasion routes of the house mouse Mus musculus domesticus and the black rat Rattus rattus currently invading Senegal (West Africa).
These two rodent species have been colonizing Senegal eastwards from western coastal colonial cities since the beginning of the 20 th century 29,30 . Our previous study based on few 5 immune effectors highlighted higher inflammatory and/or antibody-mediated responses in recently invaded sites compared with anciently invaded ones for both invaders 31 . These immune variations could be considered as responses to novel parasite pressures encountered in recently invaded areas. Surveys of bacterial communities in Senegal seemed to corroborate this prediction 32 as both quantitative and qualitative changes in bacterial composition were observed along the house mouse and black rat routes of invasion. However this pattern could not be generalised to any 'parasite' pressure as a decrease of the overall prevalence in gastrointestinal helminths was also observed, suggesting enemy release for these macroparasites 33 . We therefore developed an RNAseq approach to assess without any a priori differential gene expression patterns between populations from anciently and recently invaded sites. We tested the null prediction that immune gene expression patterns would not differ along invasion routes. We investigated two alternative hypotheses. On one hand, we expected an overall higher immune gene expression of rodent populations in recently invaded sites, as a response to novel parasite pressures encountered 31 . On the other hand, we expected trade-offs between immune pathways in recently invaded sites under the EICA-refined hypothesis. Enemy release of macroparasites and changes in the composition of bacterial communities in recently invaded sites could mediate such trade-offs. More specifically, we investigated the following questions: (i) Are there differences in immune gene expression between anciently invaded and recently invaded populations of mice and rats? (ii) Are there functional categories of immune genes that are over-represented among the differentially expressed genes? iii) Do these results support the EICA or EICA-refined hypotheses?

Results
Qualitative description of expression patterns 6 The 32 transcriptome libraries sequenced (see Table 1 for details) produced an average of 51,6 million read pairs. The numbers of retained reads post-filtering ranged from 37.2 M to 76.7 M, corresponding to 92.82% of the reads. When considering the house mouse data (3 mismatches/read pair allowed), an average of 85.68 % of the quality-filtered reads mapped to the M. musculus domesticus reference genome and 1.58% of these latter mapped to multiple regions. When considering the black rat R. rattus data and the possibility for 9 mismatches/read pair with R. norvegicus reference genome, 81.66% of the quality-filtered reads were mapped to the reference genome of the related R. norvegicus. Among them, 0.98% mapped to multiple regions of the reference genome.
The principal component analysis (PCA) performed on the house mouse standardized read counts separated sampling sites along the second principal component, which accounted for 4.58% of the variance observed (Suppl. Fig. S1a). We observed a slight opposition between sites corresponding to anciently and recently invaded areas, except for Aere Lao that did not group with recently invaded sites. No separation of sites according to invasion-related categories was visible on the PCA performed on the black rat data. On the second principal component that represented 2.98 % of the variance explained, two sites from anciently invaded areas were opposed to the six other ones (Suppl. Fig. S1b).

Differential expression (DE) analyses
For both datasets, a high dispersion of the total number of reads could be observed between biological replicates (Suppl. Fig. S2a,b). Therefore we applied differential expression (DE) analyses using two strategies. The first one named '4vs4 approach' considered a single randomly chosen biological replicate for each site and analysed the 256 comparisons of the relative log expression (RLE)-transformed normalized read counts between anciently and recently invaded sites. This approach is very conservative. The second one named '8vs8 approach' considered the eight possibilities of 'site x replicate' and employed a generalized linear model (GLM) to take into account the fact that each site was represented by two biological replicates in this analysis. This approach is less conservative.
Along the mouse invasion route, over the 17,853 filtered genes (more than 10 occurrences) tested using edgeR and the '4vs4 approach', 18 genes were always differentially expressed between anciently and recently invaded sites. Over the 16,630 filtered genes (more than 10 occurrences) tested using the '8vs8 approach', we detected 593 differentially expressed genes among which 364 were found in at least 85% of the 256 previous '4vs4 approach' comparisons. Also note that only five genes were detected by both '4vs4' and '8vs8' approaches (Hal, Rnf183, Serpina6, Wif1, 9030619P08Rik; see details in Suppl. Table S1a).
Along the black rat invasion route, 54 genes were always differentially expressed between anciently and recently invaded sites among the 13,190 ones (more than 10 occurrences) analysed when using the '4vs4 approach'. We found 268 DE genes among the 12,747 ones (more than 10 occurrences) analysed with the '8vs8 approach' among which 83 were detected in at least 85% of the '4vs4' comparisons. Forty-two differentially expressed genes were common to the '4vs4' and '8vs8' approaches. Results are detailed in Suppl. Table S1b.
We have focused further analyses on the 364 and 83 genes that were 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, which we qualified as "robust" genes. The magnitude of expression differences were significantly different between the mouse and rat, with higher 'log fold changes' observed for mouse genes (Fig. 1, Wilcoxon test, p = 0.01). The log fold change (LogFC) per gene ranged between -9.04 and 2.34 for the mouse and between -5.79 and 4.49 for the rat robust genes. The average coverage varied from -3.16 to 11.94 LogCPM (counts per million) for the mouse and -2.41 to 9.56 for the rat robust genes.

Functional analyses of DE genes
8 Considering the house mouse dataset, STRING tool provided at least one gene ontology (GO) annotation for 345 of the 364 DE genes (Suppl . Table S1). Among the 2,098 GO annotations found, 357 GO terms corresponding to biological processes were significantly enriched (FDR < 0.05) when compared to the house mouse reference genome. In particular, 10 immune-  Table S1). When taking into account redundancies using REVIGO (SimRel semantic similarity measure, similarity cutoff = 0.7), we found 20 main clusters for the category 'Biological Process'. Inflammation and acutephase response were part of the second most represented one (Fig. 2). Using KEGG algorithm, we detected 47 enriched pathways. We found that 29 immune related genes belonged to these over-represented KEGG biological pathways. The most significant one being the immune 'Complement and coagulation cascades' (FDR = 1.56 e-30).
Interestingly, the distribution of average coverage (LogCPM) for these immune-related genes corresponding to enriched GO terms or pathways was shifted towards higher values compared to the 364 robust genes (Wilcoxon test, p = 1.52 x 10 -3 ; Fig. 3, Suppl. Fig. S3). The heatmaps built on the normalized read counts for the immune related genes belonging to the biological processes (as well as pathways) that were found to be over-represented revealed a downregulation of these genes for all anciently invaded sites and an up-regulation for all recently invaded sites except Aere Lao (Fig. 3, Suppl. Fig. S4). In this latter population, genes exhibited similar expression patterns than anciently invaded sites. 9 Finally, using STRING, we found 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), meaning that proteins have more interactions among themselves than what would be expected for a random set of proteins of similar size, drawn from the mouse genome. Two main clusters were detected, one including Cytochrome P450 and UDP glucuronosyltransferase proteins and the other being mainly composed of immune related proteins, including fibrinogen, serin peptidase inhibitor (serpin), apolipoprotein and complement proteins. The PPI interaction network built upon the 73 immune related genes belonging to the biological processes previously found to be over-represented was also statistically significant (p < 1.0e-16, 561 edges found, 13 expected). The main cluster included the fibrinogen and serine peptidase inhibitor (serpin) proteins. Other highly connected proteins included apolipoprotein A and haptoglobin on one hand, and complement proteins on the other hand (Fig. 4).
Considering the R. rattus dataset, GO annotations were available for 75 of the 83 DE genes in black rat (Suppl. Table S2). Eight of these annotated DE genes were related to immune biological processes (i.e. 10.66 % of all DE genes; Suppl. Table S2). We did not find any biological process or pathway significantly enriched, even when considering redundancies using REVIGO. The protein-protein interaction network built from STRING was not significantly enriched (p = 0.09).

Discussion
In this study, we presented a comparative transcriptomics approach based on spleen RNA sequencing for two invasive rodent species studied along their invasion routes in Senegal. We aimed at identifying adaptive divergence between recently and anciently established rodent populations that may at least partly explain invasion success. More specifically, we investigated differential immune related gene expression that could provide support in favour of the EICA hypotheses. We observed differentially expressed genes along the mouse and rat invasion routes. Among them, 20.0% for the mouse and 10.6% for the rat were annotated with functions related to immunity. None of the results found with regard to these genes supported one of the EICA hypotheses. Contrary to what was expected, along the mouse invasion route, all immune-related genes detected were over-expressed in recently invaded sites, and among them, inflammation and complement pathways were over-represented. Considering the black rat, results were less straightforward -changes of expression levels were of lower amplitude than for the mouse genes -and no marked pattern could be drawn. Genes found to be differentially expressed were either up-or down-regulated in recently-invaded sites, but no particular biological process was significantly over-represented. Although bioinformatics limitations could underlie these differences between the house mouse and black rat (transcriptome annotation was less efficient for R. rattus as we used R. norvegicus transcriptome as a reference), our results suggested that variations of immune phenotypes were a less important strategy for R. rattus invasion success. Overall, it is likely that different processes, including colonization history and/or alternative mechanisms by which species adapt to novel environments, may have mediated the invasion success of these two rodent species. In the future, it could be interesting to take into account the potential variability of these patterns with regard to the targeted tissue or organ. We here focused on rodent spleen. Results would perhaps have been different if we had targeted other lymphoid organs, lymphatic tissues or non-immune targets.
Invasion success is supposed to rely on the evolution of phenotypic traits that would be advantageous at invasion fronts, where species are likely to be exposed to novel 11 environmental conditions 7,9,34 . Here, using comparative transcriptomics in wild populations, we have evidenced phenotypic variation along the mouse and rat invasion routes in Senegal.
Indeed, we detected variations of gene expression levels, a phenotypic characterisation that results from a combination of genotype, environment and genotype-environment interactions 35 . This phenotypic variation may result from stochastic changes due to population history, including founder events and genetic drift 7 . Previous population genetic studies suggested serial founder events along invasion routes for both the house mouse 36  Variations of gene expression levels may also result from plastic responses to novel selection pressures (e.g. 37 ) or from deterministic shifts involving natural selection and adaptation to local conditions 38 . Although mouse and rat populations experienced reduced genetic diversity due to founder events, they may have developed adaptive responses to novel selection pressures through high levels of plasticity 39,40 . In the house mouse, the functional categories of genes that were over-represented among the differentially expressed ones gave us insights into the phenotypic traits that may be concerned by deterministic shifts during invasion.
Transcriptome sequencing revealed strong evidence for significantly biased variations in immune gene expression towards an up-regulation of complement and pro-inflammatory cascades in recently invaded sites. These results corroborated our previous functional studies 31 : immune challenges revealed that antibody-mediated (natural antibodies and complement) and inflammatory (haptoglobin) responses were increased in mice sampled in recently 12 invaded sites compared to those sampled in anciently invaded ones. This study was therefore in line with the hypothesis that parasite-mediated selection could contribute to phenotypic differentiation along invasion routes, potentially in interaction with other environmental factors. However, it still did not support the EICA or EICA refined hypotheses 14 . Despite the enemy release patterns detected for mouse helminths in Senegal (decrease of specific richness and of overall prevalence in recently invaded sites, see 33 ), we did not find any evidence of immunity or particular immune pathways being dampened at the expense of other life history traits or of less energetically costly immune defences. The up-regulation of all immune genes found to be differentially expressed in sites recently invaded by the house mouse strongly supported the assumption of an increased overall infection risk in recently invaded sites. The changes in bacterial communities observed between anciently and recently invaded sites introduction of exotic bacterial genera and/or changes in the prevalence of endemic ones, see 32 could be in line with this potential increased risk. Similar rapid changes in immunity have been observed in other invading organisms, but these studies revealed a wide array of patterns that did not enable to conclude on the potential genericity / specificity of the EICA hypotheses 15 . We advocate for further studies using wide ranges of immune effectors, including 'omics methods' 41 , to improve our understanding of the relationships between biological invasion and immunity.
Interestingly, the site of Aere Lao did not provide the expected signature along the mouse invasion route, as it provided similar results than anciently invaded sites. The difference of Aere Lao among recently-invaded sites was already noticeable in some of the effectors (complement; haptoglobin) previously measured by immune challenges 31 . This pattern may not be related to parasitological data, because endemic parasites were found at high levels of prevalence in this site 32,33 . On the contrary, it may reflect a more ancient invasion history of the house mouse in this locality than in the neighbouring ones, for which we have historical 13 data 29 . Indeed, Aere Lao is a bigger city than Dodel, Croisement Boubé or Lougue (P. Handschumacher, pers. com.) and a market place that may have been colonised earlier by the house mouse.
Concerning R. rattus, no particular pathway (including immune related ones) was significantly over-represented. In addition, only eight immune related genes were found to be differentially expressed. They encode for proteins that interact with bacteria (Rnase6, Ptprz1), viruses (Ltc4s, Hspa1b) or parasites (Gfi1b). Some of them were down-regulated in recently invaded sites (e.g. Rnase6) while others were strongly up-regulated (e.g. Ptprz1). So, no specific pattern emerged with regard to potential variations in infection risks. Interestingly, we also previously found from a functional immune study that differences in immune responses between the rat anciently and recently invaded sites were less marked than for the mouse 31 . In addition, we could hardly identify other biological processes that could reflect life history traits adaptation with regard to invasion. Those processes that were found to be associated with differentially expressed genes covered a wide range of functions including clustering of voltage sodium channels, negative regulation of cell-matrix adhesion and protein dephosphorylation and response to cyclic Adenosine Monophosphate (cAMP). This latter function could be interesting with regard to the rat invasion success as cAMP is central to the regulation of insulin and glucagon secretion 42 . Variations in cAMP gene expression could therefore mediate differences in response to stressful situations, including starvation or fightor-flight response 43 . It would be interesting to analyse whether such differences could result in different behavioral phenotypes between rats trapped in anciently invaded sites (expected to exhibit little performance and high stress response in novel environments) and recently invaded ones (expected to exxhibit high response capacity in novel environments) 44 .
Therefore it seems even more important to further perform transcriptomics analyses on other 14 organs and tissues to identify the phenotypic changes and ecoevolutionary processes linked with the black rat invasion success. Brain, where key genes underlying behavioral invasion syndrome are expected to be expressed, could be a relevant organ to target in the future 45 .
In conclusion, our work revealed changes in transcriptomic profiles along invasion routes for both the house mouse and the black rat. The patterns observed could potentially be driven by increased infection risks in recently invaded sites for the house mouse and by stochastic changes for the black rat. Whether differences in expression were driven by phenotypic plasticity or directional selection during or after invasion, or reflected the colonization history of rodents, has yet to be confirmed, and should deserve genomic studies and experimental work 18,21 . Moreover, it would also be important to validate whether these changes in phenotypic traits influenced ecological dynamics (e.g. 46 ) and in turn, invasion success.

Ethical statements
Sampling campaigns within private properties were realized with the prior explicit agreement

Sampling sites, sample collection and RNA extraction
We used some of the spleen samples of M. m. domesticus and R. rattus that were previously analysed for the bacterial metabarcoding study previously described 32 . Briefly, for both invasive species, rodent trapping was performed in four anciently invaded (more than 100 years) and four recently invaded (less than 30 years) sites defined according to historical and longitudinal surveys 29 . Autopsies were realized at approximately the same period of the day for all individuals (11 am -14 pm). Spleens were collected and immediately placed in RNAlater, stored at 4 °C overnight then at −20 °C until further analyses.
Twenty adult rodents per site were considered excepted for five sites where only 14-18 samples were available (see Table 1). Total RNA was extracted from approximately 5mg of spleen for each sample using AllPrep 96 DNA/RNA Kit (Qiagen). The quality and quantity of the purified RNA was assessed by gel electrophoresis and NanoDrop spectrophotometer (Thermo Scientific) before pooling, and Bioanalyzer 2100 (Agilent) after pooling (see below).
For each rodent species, sixteen equimolar pools (2 per sites) that combined 6 to 10 individual samples with equilibrated sex ratio were made.

Transcript assembly
Image analyses and base calling were performed using the Illumina HiSeq Control Software and Real-Time Analysis component. Demultiplexing and production of Fastq files were performed using Illumina's 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 using FASTQC from the Babraham Institute and the Illumina software SAV (Sequencing Analysis Viewer).
A splice junction mapper, TOPHAT 48 (v2.0.9 for the mouse data, v2.0.13 for the black rat data), using BOWTIE 49 (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 gene model 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 having more than 3 mismatches (house mouse) or 9 mismatches (black rat) were discarded. Again we used a different threshold for mouse and rat data as the reference genome available for this latter is from Rattus norvegicus. Read alignment rates were above 83,3 % for all M. m. domesticus libraries and 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. Then, gene counting was performed with HTSEQ count (v0.5.4p5 for the mouse data, v0.6.1p1 for the black rat data) using the union mode 50 . The data is from a strand-specific assay, the read has to be mapped to the opposite strand of the gene.

Analysis of differential gene expression between populations
Differentially expressed (DE) genes were identified using the BIOCONDUCTOR 51 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, 52 .
It enabled to take into account the potentially high variability observed between biological replicates. Data were normalized using the relative log expression Rle normalization factors 53 . Two analysis approaches have been used to identify the DE genes. First, one out of the two locality replicates has been randomly selected, in order to have 4 replicates in each condition ('recently invaded' versus 'anciently invaded' sites, '4vs4 approach'). In total, 256 combinations of samples were possible and the comparisons performed. Second, the eight samples from each condition were kept, and a 'locality' factor was added to the design in order to consider that replicates were paired ('8vs8 approach'). Before statistical analysis, genes with less than 20 reads (first approach) or 40 reads (second approach), cumulating all the 8 or 16 analysed samples per species, were filtered and thus removed. It enabled to limit the number of statistical tests and therefore to lower the impact of multiple tests corrections.
In these two approaches, the p-value threshold was set to 5%, after application of the Benjamini-Hochberg method for multiple testing correction.
The differentially expressed genes identified with the '8vs8 approach' and found in more than 85% of the 256 combinations of the '4vs4 approach' were declared as DE. The threshold of 85% was set according to the visualisation of a barplot representing the number of '4vs4' comparisons in which a gene (highlighted in the '8vs8' approach) is found to be differentially expressed.
We performed a principal component analysis (PCA) on normalized read counts of expressed genes (excluding weakly expressed ones, see above) using R v3.0.1 54

Additional Information (including a Competing Financial Interests Statement)
We have no conflict of interest to declare.

Data availability
All sequence data, including raw reads and assemblies have been deposited on GEO. Gene    Heatmap was built in R using heatmap.2 for a) all DE genes and b) 73 immune-related genes belonging to over-represented biological pathways. The genes (rows) and samples (columns) were clustered using dendrograms built with Ward distance and hierarchical clustering. Table S1. a) Details of the 18 genes found to be differentially expressed between the anciently and recently invaded sites of M. musculus domesticus invasion route using the 4x4 approach. Genes indicated in bold were found to be differentially expressed with both 8x8 and 4x4 approaches. b) Details of the 364 genes found to be differentially expressed between the anciently and recently invaded sites of M. musculus domesticus invasion route using the 8x8 approach and found in 85% of the comparisons made using the 4x4 approach. The 73 genes indicated in red are related with immunity and have over-represented GO annotations.
The 29 genes underlined and in italics correspond to the immune-related genes with over represented Kegg biological pathways. GO ID corresponds to Gene ontology annotations. Table S2. a) Details of the 54 genes found to be differentially expressed between the anciently and recently invaded sites of R. rattus invasion route using the 4x4 approach. Genes indicated in bold were found to be differentially expressed with both 8x8 and 4x4 approaches. b) Details of the 83 genes found to be differentially expressed between the anciently and recently invaded sites of R. rattus invasion route using the 8x8 approach and found in 85% of the comparisons made using the 4x4 approach. Genes indicated in bold were found to be 22 differentially expressed with both 8x8 and 4x4 approaches. Genes indicated in red are related with immunity. GO ID corresponds to Gene ontology annotations.