Relationships between immune gene expression and circulating cytokine levels in wild house mice

Abstract Quantitative PCR (qPCR) has been commonly used to measure gene expression in a number of research contexts, but the measured RNA concentrations do not always represent the concentrations of active proteins which they encode. This can be due to transcriptional regulation or post‐translational modifications, or localization of immune environments, as can occur during infection. However, in studies using free‐living non‐model species, such as in ecoimmunological research, qPCR may be the only available option to measure a parameter of interest, and so understanding the quantitative link between gene expression and associated effector protein levels is vital. Here, we use qPCR to measure concentrations of RNA from mesenteric lymph node (MLN) and spleen tissue, and multiplex ELISA of blood serum to measure circulating cytokine concentrations in a wild population of a model species, Mus musculus domesticus. Few significant correlations were found between gene expression levels and circulating cytokines of the same immune genes or proteins, or related functional groups. Where significant correlations were observed, these were most frequently within the measured tissue (i.e., the expression levels of genes measured from spleen tissue were more likely to correlate with each other rather than with genes measured from MLN tissue, or with cytokine concentrations measured from blood). Potential reasons for discrepancies between measures including differences in decay rates and transcriptional regulation networks are discussed. We highlight the relative usefulness of different measures under different research questions and consider what might be inferred from immune assays.


| INTRODUC TI ON
Since the field of ecoimmunology emerged over two decades ago (Sheldon & Verhulst, 1996), measuring the immune function of wild animals has remained challenging (Pedersen & Babayan, 2011).
To date, the vast majority of ecoimmunological studies have used non-model species (Jackson, 2015) and so have been limited in their scope or by the availability of appropriate assays or reagents (Fassbinder-Orth, 2014;Zimmerman et al., 2014). Because of this, many early studies relied on less specific immunological assays, such as leukocyte counts, hemagglutination, responses after stimulation with phytohemagglutinin, and bactericidal assays (Demas et al., 2011;Pedersen & Babayan, 2011;Zimmerman et al., 2014).
As modern sequencing technologies have become more accessible and affordable, it has become possible to measure immune responses in non-model organisms using more specific assays (Jackson, 2015). For example, reverse transcription quantitative real-time PCR (RT-qPCR, often abbreviated to qPCR) can be used to measure the expression level of immune genes, such as cytokines, via mRNA concentrations (Adams, 2020). This gene expression-based approach has been widely used to investigate immune function in wild populations of a range of taxa (Fassbinder-Orth, 2014).
Gene expression is regularly used as a proxy for measuring functional protein products (Fassbinder-Orth, 2014), and it is often assumed that the levels of mRNA expression will correlate with the levels of the corresponding proteins (Fassbinder-Orth, 2014;Maier et al., 2009;Vogel & Marcotte, 2012). Indeed, many researchers have focused solely on qPCR and measuring levels of gene expression (Koussounadis et al., 2015). However, the correlation between mRNA concentration and corresponding protein concentration is often dependent on cell type and state (Silva & Vogel, 2016). Indeed, in bacteria and eukaryotes only around 40% of variation in cellular protein levels can be predicted by mRNA concentrations (de Sousa Abreu et al., 2009;Vogel & Marcotte, 2012). Therefore, mRNA concentrations may not always accurately represent the expression levels of the bioactive protein (Jackson, 2015;Pradet-Balade et al., 2001) Further, the immune response may be localized and tissue specific, for example, in response to localized infections or tissue damage (Hu & Pasare, 2013).

Although measures of immune gene expression in verte-
brates may indicate an induced immune response in the host (Jackson, 2015), it should not necessarily be assumed that levels of gene expression completely align with their corresponding effector proteins. In addition to localization of responses, there may be levels of regulation occurring between genes, or delays between transcription and translation and protein modification (Payne, 2015;Pradet-Balade et al., 2001). Indeed, studies have found that the extent to which mRNA expression correlates with the final gene products can vary with transcriptional and post-translational processes, and that the rates of decay between the two can significantly differ (Munsky et al., 2012;Schwanhäusser et al., 2011). Therefore, making robust measures of molecular phenotype requires the measurement of active proteins using antibody reagents against key immunological biomolecules (Jackson, 2015). This is particularly important when dealing with natural populations in which variation is intrinsically higher (Abolins et al., 2011(Abolins et al., , 2017. Protein levels can only reliably be analyzed using sophisticated assays that are typically only available for model species, such as house mice (Mus musculus). For these species, commercial off-theshelf reagents are readily available. Due to structural variation in immune molecules such as cytokines, however, the transferability and cross-reactivity of these assays to non-model species are extremely limited (Jackson, 2015). Developing specific antibodies is typically expensive and time-consuming (Bradley & Jackson, 2008;Friberg et al., 2010;Jackson, 2015;Oko et al., 2006;Pedersen & Babayan, 2011). Therefore, it is informative to understand the precise relationships between measured gene expression and the effector proteins that will actually interact at the immune interface with a pathogen or parasite, particularly in wild populations.
Despite the inherent difficulties associated with applying sophisticated molecular techniques to wild animals, a number of studies have investigated the immune function of wild populations of model species, primarily house mice. These have generally employed specific measures of immune function, such as antibody concentrations and avidity (Abolins et al., 2011;Lochmiller et al., 1991), or circulating cytokine concentrations (Abolins et al., 2017).
Here, using the extensively studied natural population of M. m.
domesticus on the Isle of May (for examples see Berry et al., 1990;Goertz et al., 2019;Taylor et al., 2019;Triggs, 1991), we used qPCR to measure the expression levels of key cytokine genes representing four arms of the adaptive immune response (Th1, Th2, Treg, and Th17), as well as the innate immune response. These were measured from two key immunological tissues; the spleen and mesenteric lymph nodes (MLN). In addition, we measured the circulating levels of the corresponding secreted cytokines from serum using multiplex bead assay. Further, the relative mRNA expression of several important immunity-related transcription factors was measured using qPCR from spleen and MLN tissue.
We aimed to assess the relationship between the expression of these key cytokine genes and the corresponding circulating protein levels in our population of M. m. domesticus. Further, by using both MLN and spleen tissue, we aim to explore the relationship between immune gene expression from two distinct but immunologically important tissue types.
Due to post-translational processes and protein degradation, we hypothesized that there would not be widespread correlation between gene expression and circulating cytokine levels. The expression level of genes from MLN tissue, for example, are unlikely to relate to the expression level of genes from spleen tissue or concentrations of circulating cytokines as these represent a localized, gastrointestinal response. Further, we hypothesized that there may be antagonistic relationships between some response types. For example, the Th1-type response would not be expected to correlate with the Th2-type response due to their antagonistic relationship (Kaiko et al., 2008). Note: "Gene" gives the gene name and official gene symbol (in italics). "Accession number" gives the gene on which primers were based. "Amplicon size" indicates the size of the product in base pairs. "Single exon" shows where the primers cover a single exon ("Yes"), span an intron ("No"), or where these data are not determined ("ND"). "Tm" shows the primer melting temperature. "Function" summarizes the immunological role of the target gene. Gene expression was measured from either spleen or mesenteric lymph node (MLN) tissue.

| Sample collection
Mice were live-trapped from a wild population on the Isle of May Longworth traps (Longworth Scientific Instrument Co.) were primarily used, along with small numbers of Ugglan (Granhab) and homemade "Jordan" traps (Perrow & Jowitt, 1995). Mice were euthanized in a CO 2 chamber with a rising gradient of CO 2 . The eye and foot reflexes were tested, and death was confirmed through exsanguination by cardiac puncture. Immediately following the confirmation of death, a sample of blood was removed to a 1.5 ml microcentrifuge tube. This was stored at room temperature for ~60 min and then stored at 4°C for a further ~60 min. The resulting clot was carefully removed, and the remaining serum was spun at 3,000 g for 10 min.
The serum was removed and snap frozen using liquid nitrogen, before being stored at −80°C until its use in multiplex bead assays to measure cytokine levels.
The spleen and MLN were removed and placed in RNAlater solution (Life Technologies). Following the manufacturer's recommendations, tissues were sectioned into smaller pieces (≤ 5 mm) prior to immersion in the solution. Samples were kept at 4°C for 24 hr, then the supernatant was removed, and samples were stored at −80°C until RNA extraction.

| RNA extraction & cDNA synthesis
RNA was extracted from up to 30 mg spleen and MLN tissue stored in RNAlater using the NucleoSpin RNA kit (Machery-Nagel), following the manufacturers' protocol. A DNase treatment step is included in this protocol, and RNA was eluted in 60 μL of nucleasefree water. The purity and concentration of RNA were assessed on a NanoDrop 1,000 spectrophotometer (Thermo Scientific), with a desired 260/280 absorbance ration of >1.80 (Robertson et al., 2016).
To check for contamination, and to assess the purity of RNA, 6 μL of each sample was mixed with 4 μL of Orange G gel loading dye (Sigma Aldrich). This was incubated at 65°C for 10 min, followed by visualization on a 2% agarose gel stained with ethidium bromide run at 90V for 30 min.
Synthesis of cDNA was performed on up to 2 µg of total RNA using the nanoScript2 Reverse Transcription kit (Primerdesign), using a combination of oligo-dT and random nonamer primers (0.5 µl of each). Nuclease-free water was used to make a total reaction volume of 10 µl. This was incubated at 65°C for five minutes, before being cooled on ice. To this, 10 µl of extension-mix (5.0 µl of nanoScript2 4X Buffer, 1.0 µl of 10 mM dNTP mix, 1.0 µl of nanoScript2 enzyme and 3.0 µl of nuclease-free water) was added. This was incubated at 25°C for five minutes then 42°C for 20 min, before being heat inactivated at 75°C for ten minutes. No-enzyme controls were used to monitor samples for contamination by genomic DNA. All cDNA samples were diluted 1:10 with nuclease-free water and then stored at −20°C for future analysis by qPCR.

| qPCR conditions
qPCR was carried out to measure normalized mRNA expression in spleen and MLN tissue of a suite of genes reflecting different functional arms of the immune system. mRNA levels of seven cytokine genes of interest were measured from both spleen and MLN tissue including Tumor necrosis factor-α (TNF-α), Interferon-γ (IFNγ), Interleukin (IL)-1β, IL-6, IL-10, IL-12β, IL-13, and IL-17. Expression levels of IL-5 were also measured from spleen tissue, but not MLN (Table 1). Further, the relative mRNA expression of several important immunity-related transcription factors was measured using qPCR, from both spleen and MLN tissue. This included expression levels of GATA binding protein 3 (Gata3), Interferon regulatory factor 5 (Irf5) and T-box 21/Tbet (Tbx21) ( Table 1). All primers used were designed and validated by Primerdesign (Southampton). The genes analyzed, along with their corresponding GenBank accession numbers and primer sequence information, are summarized in Table 1.
Each qPCR reaction contained 0.5 µl of the relevant primer mix (

| Quantification of gene expression
The immunological target genes analyzed here ( (also including Actb, Gapdh, Rn18s, and Rpl13a) using the geNorm algorithm (Vandesompele et al., 2002). The stability of expression was analyzed using qbase+ (Biogazelle). In this geNorm assay, 15 randomly selected spleen or MLN cDNA samples were included.
The expression of each target in each sample relative to expression in the reference cDNA was calculated using the 2− ΔΔC T method (Livak & Schmittgen, 2001;Pfaffl, 2001;Winer et al., 1999), standardized against the geometric mean C T of the two reference genes for each sample (Vandesompele et al., 2002).

| Multiplex ELISA bead assay
To determine serum cytokine levels, diluted serum samples were processed in duplicate using a custom Bio-Rad Bio-Plex Pro mouse cytokine multiplex kit (magnetic bead-based multiplex immunoassay) according to manufacturer's protocol (Bio-Rad). Included in the custom multiplex panel were antibodies for the detection of the following mouse cytokines: IFN-γ, IL-1β, IL-5, IL-6, IL-10, IL-12β, IL-13, IL-17 and TNF-α (see Table 1 for functions).
Following incubations, the reaction mixture was analyzed using a Bio-Plex 200 Luminex-based multiplex analysis system (Bio-Rad).
Unknown cytokine concentrations in samples were determined using Bio-Plex Manager Software and standard curves derived from recombinant cytokine standards-data were expressed as fold change from control (Al Gadban et al., 2012). Data that were below the assay's range of detection were assigned values of 0.001 (Abolins et al., 2017).

| Data analysis
All data were corrected to account for interplate variation by correction against the mixed reference sample. All analyses were conducted in R version 3.6.3 (R Core Team, 2020).
To test for correlation between gene expression levels of all genes expressed by spleen and MLN tissue, and the circulating serum cytokine levels, the nonparametric Spearman's correlation coefficient was used. Spearman's rank correlation coefficients and associated significance levels were generated using the rcorr() function in the "Hmisc" package (Harrell Jr., 2020). To control for multiple testing, adjusted p-values were obtained following the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995), using the p.adjust() function. Correlations were deemed to be significant when the adjusted p-value <.05.  Figure S1). Spearman's correlation coefficient was again used to test for correlations between the pro-and anti-inflammatory gene and circulating protein PC scores.
Following Robertson et al. (2016), genes and proteins were then further divided into more specific response types (Table 2). Where appropriate, PCA was again used to group relevant gene expression or circulating cytokine levels. For all response types in all tissues or circulating protein levels, PC1 explained a high proportion of variance (>68.21%; see Figure S2).

| Sample sizes
The total sample size (i.e., number of individual mice) available varied across the genes and circulating cytokines investigated due to missing samples and assay failures (Table 3; Figure S3).

| Correlations between all genes and circulating cytokines
The strongest correlations came within tissue types for gene ex-  Table S1).

| Correlations between inflammation responses
When grouped into "inflammatory" or "anti-inflammatory" genes and proteins, the strongest relationships are again within tissue types No other relationships were significant or strongly correlated (r s ranging from −0.173 to 0.016, adjusted p > .051; Table S2).

| Correlations between response type
Gene expression and circulating cytokine levels were further split into specific response types (e.g., innate, Th1 etc., see Table 2), following Jackson et al. (2011) and Robertson et al. (2016). The strongest relationships between these functional groupings came within

TA B L E 3
The sample size (n) of expression levels of immune genes expressed by spleen or mesenteric lymph node (MLN) tissue, or the circulating levels of cytokines as measured by multiplex bead assay. Also showing corresponding gene expression level (ΔΔC T ) and circulating cytokine concentrations (pg/ml), ±Standard Deviation. Samples were obtained from wild Mus musculus domesticus on the Isle of May the tissue type from which gene expression was measured, or within circulating cytokine levels (26/38 significant correlations were within tissue type; Figure 3). While there were some significant correlations (adjusted p < .049) between different response types as measured by gene expression from different tissue types, or between gene expression and circulating cytokine levels (e.g., between Th17 from spleen and Th2 from MLN), these were often weak relationships (r s ranging from −0.191 to 0.263; Table S3).
All other correlations were weak and not significant (r s ranging from −0.168 to 0.132, adjusted p > .056). The majority of nonsignificant correlations were between gene expression levels from different tissues, or between circulating cytokine levels and gene expression. However, the innate response measured by gene expression from spleen tissue showed no significant relationship with either a Treg (r s = −0.056, adjusted p = .499), Th17 (r s = −0.048, adjusted p = .598), or Th2-type response (r s = 0.107, adjusted p = .146), as measured by gene expression from spleen tissue. Further, the relationship between the innate response and Th2-type response, both measured by gene expression from MLN tissue, was not significant (r s = 0.132, adjusted p = .056) (see Table S3).

| D ISCUSS I ON
Here, we measured the expression levels of key immune genes from two immunologically important tissues (spleen and MLN), and circulating levels of key cytokines in blood, in wild house mice. We found a lack of correlation between the expression levels of immune genes

| Temporal effects
Cross-sectional sampling techniques, as were used here, are not ideal when considering delays between transcription and translation, and differences in degradation rates between proteins and RNA. Cell culturing and stimulatory cellular assays can be used to measure the potential immune response to future challenges F I G U R E 2 Correlation matrix showing the relationship between the inflammatory and anti-inflammatory immune response in wild Mus muculus domesticus from the Isle of May. Labeled in the matrix are the Spearman's correlation coefficients of significant relationships (adjusted p < .05). Immune response was measured using gene expression from spleen (S) or mesenteric lymph node (M) tissue and circulating levels of inflammatory and antiinflammatory cytokines measured by multiplex bead assay (C)

F I G U R E 3 Correlation matrix
showing the relationship between immune response type in wild Mus musculus domesticus from the Isle of May. Labeled in the matrix are the Spearman's correlation coefficients of significant relationships (adjusted p < .05). Immune response was measured by gene expression from spleen (S) or mesenteric lymph node (M) tissue, or by circulating cytokine (C) levels measured by multiplex bead assay (Bradley & Jackson, 2008;Robertson et al., 2016), but were not used in this study due to logistical restraints.

| Downstream effects
Discrepancies between measures may be indicative of the complexities of the molecular pathways involved, as many of the cytokines analyzed have key roles in inducing or suppressing the expression of different cytokines. A lack of correlation between immune measures in humans was largely explained by post-transcriptional and post-translational regulation (Kozak, 2007;Mehra et al., 2003;Nie et al., 2006;Shebl et al., 2010;Shen & Pili, 2008 (Zhu & Kanneganti, 2017).
Measurements of these cytokines may be less informative when trying to assess the immune state in a holistic manner, as the level of immune stimulation affecting the host is unknown (Bartoccioni et al., 1994;Nilsson et al., 1995). Edfors et al. (2016) used an RNA-to-protein conversion factor index to increase predictability of protein copy numbers from RNA levels. This, however, is calculated using a whole transcriptome analysis, and so is not appropriate for targeted qPCR studies such as this one. In addition, it is highly unlikely that the weak correlations between gene expression and protein levels seen here would become strong, even with post-translational or post-transcriptional correction. Further, this study was investigating the correlation of genes expressed by two different tissue types with secreted proteins, so universally strong correlations were not necessarily expected.

| Localization of responses
The gut and the spleen are subject to different immunological challenges. The spleen responds to infections in the blood (Mebius & Kraal, 2005), and can therefore respond to infections throughout the body if immunoactive elements enter the blood and are filtered through the spleen (Bronte & Pittet, 2013 (Finney et al., 2007). Increased expression of MHC genes in response to helminth infection might not be detectable in the spleen, but would be detected in tissues local to the infection, that is, intestinal tissues or MLNs (Schwensow et al., 2011).

| CON CLUS IONS
Most studies of wild immunology have focused on non-model species and so have been restricted to measuring immune gene expression (Jackson et al., , 2014. Protein levels can only be reliably analyzed using sophisticated assays, generally only available for model species such as laboratory mice, M. musculus. Although multiplex ELISA assays can be developed for non-model species, it can be a costly and time-consuming process. Therefore, gene expression, measured by qPCR, can be a useful tool when it is not possible to measure circulating protein levels. Indeed, while RNA concentrations and protein levels may not correlate, gene expression may provide insight into the types and quantities of immune cell present in a tissue. Regardless, the discrepancy between RNA levels and corresponding proteins should be a consideration that is addressed when reporting results.
Other means of comprehensively assessing immune expression do exist and have potential for applications in wild populations.
RNA-seq provides a whole transcriptome approach, and so any relevant changes in gene expression outside of those target genes can be identified. qPCR, however, has the benefit of being comparatively cheap, thus allowing a higher number of replicates and greater statistical power. Further, the genes of interest selected were based on a priori knowledge of immune function, providing a much more targeted approach.
While the results here show some correlations between different measures of immune response, they highlight that care should be taken when interpreting either gene expression or secreted protein-level data. Dependent on the precise questions being asked during an experiment, it might be more appropriate to focus on one measure or tissue type than others. For instance, one might be interested in the immune response to a parasite or pathogen localized to a specific tissue, or want to ascertain how immune responses are regulated by upstream gene expression rather than looking at downstream effector proteins. Selecting the measure of immune function should be handled with great care, particularly in wild populations, as factors such as coinfection, seasonal changes, and resource availability can have far-reaching and overlooked impacts upon immune function. As such, it is advisable to obtain the most comprehensive picture of immunity as is possible in most circumstances.