Immune and environment‐driven gene expression during invasion: An eco‐immunological application of RNA‐Seq

Abstract Host–pathogen associations change rapidly during a biological invasion and are predicted to impose strong selection on immune function. It has been proposed that the invader may experience an abrupt reduction in pathogen‐mediated selection (“enemy release”), thereby favoring decreased investment into “costly” immune responses. Across plants and animals, there is mixed support for this prediction. Pathogens are not the only form of selection imposed on invaders; differences in abiotic environmental conditions between native and introduced ranges are also expected to drive rapid evolution. Here, we use RNA‐Seq to assess the expression patterns of immune and environmentally associated genes in the cane toad (Rhinella marina) across its invasive Australian range. Transcripts encoding mediators of costly immune responses (inflammation, cytotoxicity) showed a curvilinear relationship with invasion history, with highest expression in toads from oldest and newest colonized areas. This pattern is surprising given theoretical expectations of density dynamics in invasive species and may be because density influences both intraspecific competition and parasite transmission, generating conflicting effects on the strength of immune responses. Alternatively, this expression pattern may be the result of other evolutionary forces, such as spatial sorting and genetic drift, working simultaneously with natural selection. Our findings do not support predictions about immune function based on the enemy release hypothesis and suggest instead that the effects of enemy release are difficult to isolate in wild populations, especially in the absence of information regarding parasite and pathogen infection. Additionally, expression patterns of genes underlying putatively environmentally associated traits are consistent with previous genetic studies, providing further support that Australian cane toads have adapted to novel abiotic challenges.


| INTRODUC TI ON
Invasive species pose a massive threat to biodiversity (Bax, Williamson, Aguero, Gonzalez, & Geeves, 2003;Clavero, Brotonsa, Pons, & Sol, 2009). The potential for pathogens to limit the impact of invaders, or to exacerbate that impact, makes it critical to understand how the process of range expansion alters invader immunity.
In equilibrial systems, pathogen-mediated selection (PMS) favors host individuals with traits that enable them to resist or tolerate infections (Spurgin & Richardson, 2010), but it is unclear whether the same traits remain favorable during invasion. Furthermore, some invasive traits follow patterns that are not explained by selection (Berthouly-Salazar, Rensburg, Roux, Vuuren, & Hui, 2012;Lowe, Muhlfeld, & Allendorf, 2015). This may be because the dispersive tendencies of invaders give rise to additional evolutionary forces: (a) As an invasive population expands, genetic drift may reduce genetic diversity across the range (Rollins, Woolnough, Wilton, Sinclair, & Sherwin, 2009) and modify phenotypic traits. (b) An expanding invasion front is dominated by individuals with the highest rates of dispersal simply because the fastest arrive at new areas first and can only breed with each other (spatial sorting; Shine, Brown, & Phillips, 2011). Thus, a geographic separation of phenotypes occurs; traits that enhance individual fitness are favored in established populations, and traits that enhance dispersal rate are common in expanding populations (Hudson, Brown, & Shine, 2016;Shine et al., 2011).
(c) Admixture between individuals from different introductions or sources, as well as hybridization, may also drive evolutionary change in some invasive populations (Mader, Castro, Bonatto, & Freitas, 2016). Alternatively, some traits differ simply due to heightened environmental variability experienced by invaders and do not reflect evolutionary change. Given the complex possibilities of adaptive and nonadaptive changes imposed on immune systems during range expansion, identifying their underlying basis is an important task.
We examined the effects of range expansion on expression of immune and environmentally associated genes in the invasive Australian cane toad (Rhinella marina) using RNA-Seq data from whole spleen tissue from individuals collected from long-established areas in Queensland (QLD, the "range core"), geographically "intermediate" areas in the Northern Territory (NT), and the leading edge of the range expansion in Western Australia (WA, the "invasion front"; Figure 1). The invasive range of cane toads in Australia includes highly varied environments; climatic conditions in the range core are similar to those in the native range (Central and South America), but intermediate areas and the invasion front receive much less annual rainfall (2,000-3,000 mm in QLD, 400-1,000 mm in NT and WA) and have higher annual mean temperatures (21-24°C in QLD, 24-27°C in NT and WA;Bureau of Meteorology A.G., 2018). Toads cluster genetically based on these environmental patterns: Toads from the range core are genetically distinct from those from intermediate areas and the invasion front . Furthermore, loci putatively under selection are involved in tolerance of temperature extremes and dehydration . Traits such as locomotor performance at high temperatures also follow this pattern (Kosmala, Brown, Christian, Hudson, & Shine, 2018), but others do not. For example, behavioral propensity for exploration increases with distance from the introduction site (Gruber, Brown, Whiting, & Shine, 2017). Traits such as leg length (Hudson et al., 2016), spleen size, and fat body mass (Brown, Kelehear, Shilton, Phillips, & Shine, 2015b) follow a U-shaped (curvilinear) pattern F I G U R E 1 Geographic distribution of the cane toad in Australia (dark gray region). Since arriving in Queensland in 1935, cane toads have further expanded their range through New South Wales, the Northern Territory, and into Western Australia. Black diamonds indicate our toad collection sites (from east to west): QLD (Gordonvale and Daintree, N = 5 each), NT (Cape Crawford and Timber Creek, N = 4 each), and WA (Caroline Pool and Durack River, N = 5 each). Map adapted from Tingley, et al. (2017) across the range, in which they are smallest in toads from intermediate areas of the range and larger at either end.
One process likely to be important in shaping immune systems during range expansion is loss of parasites and pathogens. The enemy release hypothesis (ERH) predicts that the processes of introduction and range expansion decrease rates of infections with coevolved pathogens and parasites in invasive hosts due to the former's inability to spread effectively and persist in novel environmental conditions (Colautti, Ricciardi, Grigorovich, & MacIsaac, 2004). Because of this, Lee and Klasing (2004) predict that invaders may downregulate powerful immune responses such as systemic inflammation due to a decreased need (Cornet, Brouat, Diagne, & Charbonnel, 2016;Lee & Klasing, 2004;Martin, Hopkins, Mydlarz, & Rohr, 2010). Such immune responses are also costly due to energetic expenditure (the reduction of nutrients available for partitioning across tissues due to their use in mounting immune responses; Klasing & Leshchinsky, 1999) and to the potential for collateral damage (tissue injury due to the effects of the immune response; Martin et al., 2010). Reduced energetic investment into these immune responses may enhance invasion success (McKean & Lazzaro, 2011). Nonetheless, loss of immunocompetence (the ability to mount a normal immune response after exposure to an antigen; Janeway, Travers, & Walport, 2001) could render invaders susceptible to infection by novel pathogens and parasites in their introduced range (Cornet et al., 2016;Lee & Klasing, 2004). Thus, invaders are predicted to exhibit lower investment in costly (but not all) immune responses than are seen in their native ranges (Cornet et al., 2016;Lee & Klasing, 2004). Enemy release is likely to be relevant to the ecoimmunology of cane toads in Australia.
Consistent with the ERH, many species of bacteria, protozoan, and metazoan parasites of cane toads have been left behind in the native range (Selechnik, Rollins, Brown, Kelehear, & Shine, 2017a). A major parasite (lungworm Rhabdias pseudosphaerocephala) from the native range that infects toad populations in the Australian range core is absent from toads at the invasion front (Phillips et al., 2010); this lungworm has been reported to increase mortality by 17% in metamorph cane toads (Kelehear, Webb, & Shine, 2009) and 6% in adult cane toads (Finnerty, Shine, & Brown, 2017). Conversely, the Australian soil bacterium Brucella (Ochrobactrum) anthropi causes spinal spondylosis in toads primarily at the invasion front (Brown, Shilton, Phillips, & Shine, 2007), which may represent a novel infection that forces invaders to remain immunocompetent. Furthermore, Rhimavirus A has only been detected in transcriptomes of toads from areas relatively close to the invasion front (Russo et al., 2018). Invasion history has complex effects on toad immunity (Brown, Phillips, Dubey, & Shine, 2015c;Brown & Shine, 2014;Selechnik, West, et al., 2017b).
The loss of pathogens underlying the ERH depends on a decline in pathogen transmission, which likely occurs when host densities are lower. The densities of many invasive populations follow a "traveling wave," in which population density is low at recently colonized areas (e.g., the invasion front), high in areas that have been colonized for several years (e.g., intermediate areas), and low at longcolonized areas (e.g., the range core; Hilker, Lewis, Seno, Langlais, & Malchow, 2005;Simberloff & Gibbons, 2004). Although absolute population densities of cane toads across Australia are unknown, toads appear to follow this trend as well (Brown, Kelehear, & Shine, 2013;Freeland, Delvinquier, & Bonnin, 1986), as does at least one of their major parasites (Rhabdias pseudosphaerocephala). This parasitic lungworm is absent from toads at the invasion front and is most prevalent in toads from intermediate areas (Brown, Kelehear, et al., 2015b;Phillips et al., 2010). Therefore, we predicted that the expression of immune genes may follow a curvilinear pattern, in which those encoding mediators of costly immune responses (e.g., inflammation) may be more highly expressed in toads from intermediate areas than in toads from the core or front. In terms of environmentally associated genes, we predicted that expression patterns would depend on the functional roles of the genes; for example, genes involved in aridity tolerance may be differentially expressed between toads from the (moist) range core and toads throughout the rest of the (more arid) range. . We euthanized toads using 150 mg/ kg sodium pentobarbital, decapitated them as soon as they became unresponsive, and excised their spleens immediately following decapitation. We selected spleen tissue for our investigation because mature immune cells travel to secondary lymphoid tissue (spleen, lymph nodes, and mucosa-associated lymphoid tissue) for activation through pathogenic encounter (Janeway et al., 2001). Thus, the cellular compositions of these tissues should reflect the host's immune functioning. Each spleen was initially preserved in RNAlater (Qiagen), kept at 4°C for <1 week, and then drained and transferred to a −80°C freezer for long-term storage.

| Sample collection and RNA extraction
Prior to RNA extraction, we flash-froze all spleens individually in liquid nitrogen and ground them with a mortar and pestle to lyse preserved tissue. We carried out RNA extractions using the RNeasy Lipid Tissue Mini Kit (QIAGEN) following the manufacturer's instructions, with an additional genomic DNA removal step using on-column RNase-free DNase treatment (Qiagen). We quantified the total RNA extracted using a Qubit RNA HS assay on a Qubit 3.0 fluorometer (Life Technologies). Extracts were then stored at −80°C until sequencing was performed.

| Sequencing
Prior to sequencing, we added 4 µl of either mix 1 or mix 2 of External RNA Controls Consortium (ERCC; Thermo Fisher Science) spike-in solutions diluted 1:100 to 2 µg of RNA to examine the technical performance of sequencing (Table S1). Macrogen (Macrogen Inc., ROK) constructed mRNA libraries using the TruSeq mRNA v2 sample kit (Illumina Inc.), which included a 300 bp selection step.
All samples from the core and the front were individually barcoded and sequenced across two lanes of Illumina HiSeq 2,500 (Illumina Inc.); samples from intermediate areas were sequenced in a separate batch on a single lane of Illumina HiSeq 2500 (also individually bardcoded). Capture of mRNA was performed using the oligo dT method, and size selection parameter choices were made according to the HiSeq2500 manufacturer's protocol. Overall, this generated 678 million paired-end 2 × 125 bp reads. Raw sequence reads are available as FASTQ files in the NCBI short read archive (SRA) under the BioProject Accession PRJNA395127.

| Data preprocessing, alignment, and expression quantification
First, we examined per base raw sequence read quality (Phred scores) and GC content, and checked for the presence of adapter sequences for each sample using FastQC v0.11.5 (Andrews, 2010).
We conducted per sample alignments of our trimmed FASTQ files to this reference using STAR v2.5.0a (Patro, Duggal, Love, Irizarry, & Kingsford, 2017) in basic two-pass mode with default parameters, a runRNGseed of 777, and specifying BAM alignment outputs. We used the BAM outputs to quantify transcript expression using Salmon v0.8.1 (Patro et al., 2017) in alignment mode with libtype = IU, thus producing count files.

| Count filtering and log-ratio transformations
Most methods for analyzing RNA-Seq expression data assume that raw read counts represent absolute abundances (Quinn, Richardson, Lovell, & Crowley, 2017). However, RNA-Seq count data are not absolute and instead represent relative abundances as a type of compositional count data (Quinn, Erb, Richardson, & Crowley, 2018c;Quinn, Richardson, et al., 2017). Using methods that assume absolute values is invalid for compositional data (without first including a transformation) because the total number of reads (library size) generated from each sample varies based on factors such as sequencing performance, making comparisons of the actual count values between samples difficult (Fernandes et al., 2014;Quinn, Erb, et al., 2018c). As such, relationships within RNA-Seq count data make more sense as ratios, either when compared to a reference or to another feature within the dataset. Hence, we analyzed our count data (from Salmon) taking the compositional nature into account using the logratio transformation (Aitchison & Egozcue, 2005;Erb & Notredame, 2016;Lovell, Pawlowsky-Glahn, Egozcue, Marguerat, & Bahler, 2015;Quinn, Erb, et al., 2018b;Quinn, Richardson, et al., 2017). Our total number of expressed transcripts across all toads was 22,930.
To filter out transcripts with low expression, we removed transcripts that did not have at least 10 counts in 10 samples. This reduced our list of expressed transcripts to 18,945. We then used the R (Team, 2016) package ALDEx2 v1.6.0 (Fernandes, Macklaim, Linn, Reid, & Gloor, 2013) to perform an interquantile log-ratio (iqlr) transformation of the transcripts' counts as the denominator for the geometric mean calculation (rather than centered log-ratio transformation) because it removes the bias of transcripts with very high and low expression that may skew the geometric mean (Quinn, Richardson, et al., 2017  mixes (Ambion) to assess whether or not there was a batch effect due to sequencing run. Each set contains four groups of sequences with different ratios between the two mixes, representing "known" differences in abundance (mix 1 vs. mix 2 fold changes: 4:1, 1:1, 1:1.5, 1:2).

| Technical and diagnostic performance
We used the R (Team, 2016) package erccdashboard v1.10.0 (Munro et al., 2014) to analyze the counts of these sequences and generate receiver operator characteristic (ROC) curves and the area under the curve (AUC) statistic, lower limit of DE detection estimates (LODR), and expression ratio variability and bias measures based on these sequence abundance ratios ( Figure S1). In the ROC curves, the AUC for two sets of true-positive ERCC sequences (4:1 and 1:2) was 1.0, indicating perfect diagnostic performance, and the AUC for the third (1:1.5) was 0.89, indicating good diagnostic performance ( Figure   S1b). The MA plot shows that the measured ratios in our ERCC sequences converge around the r m corrected ratios, indicating low variability ( Figure S1c). Finally, the LODR plot indicated that DE p-values were lower for ERCC sequences with wider ratios (i.e., 4:1 has the lowest p-values, then 1:2, then 1:1.5), which is expected because the most pronounced fold change differences should yield the highest DE significance ( Figure S1d). From these results, we inferred that the observed relative abundances between mixes of each set of ERCC sequences were close to the known relative abundances, and thus, batch effects do not appear to have occurred.
We also examined the counts of the invariant (1:1) group across all samples. Seven invariant ERCC transcripts remained after count filtering (same as used for the DE testing); we generated a boxplot of their counts, normalized by library size ( Figure S2). The consistency of the boxplot distributions of the seven invariant ERCC sequences further indicated that there was no batch effect. Because the erccdashboard package indicated that the sets of true-positive ERCC sequences (4:1, 1:2, 1:1.5) existed in observed ratios close to the known ratios, and because the invariant sequences (1:1) exhibited consistency across samples, we proceeded with downstream analyses.

| Differential gene expression in discrete phases of the invasion
After applying a log-ratio transformation to the count data, we were able to implement statistical tests that would otherwise be invalid for relative data. We grouped populations by phase (Daintree and

| Spatial gene expression patterns across the range
The DE testing performed in ALDEx2 generates differences between discrete groups; however, our data are sampled across a continuous variable: space. So, to visualize expression patterns across the toad's Australian range, we performed soft (fuzzy c-means) clustering on our log-transformed count data (with samples grouped by collection site, and sites ordered from east to west) using the R package Mfuzz v2.34.0 (Kumar & Futschik, 2007). The fuzzy c-means algorithm groups transcript together based on similar expression patterns (using a fuzzifier parameter, m) across conditions to identify prominent, recurring patterns (clusters). Each transcript within a cluster is assigned a membership value, indicating how closely its expression pattern aligns with that of the cluster to which it belongs.
To prevent random data from being clustered together, we used the mestimate command in the Mfuzz package to determine the optimal fuzzifier parameter value using a relation proposed for fuzzy cmeans clustering (Schwämmle & Jensen, 2010). We then used the cselection and Dmin commands to determine the optimal number of clusters, c, to generate. The results of both tools suggested using four clusters (c = 4); however, these tools need to be used with caution because automatic determination of the optimal value of c is difficult, and it is advised to review the data before choosing (Kumar & Futschik, 2007).

| Environmentally influenced gene expression
Genes affected by natural selection may have expression levels that are associated with environmental variables. We downloaded climatic data from the BioClim database (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) using the raster package (Hijmans, 2015) in R.
Because different areas of Australia vary in aridity, we downloaded data on rainfall during the driest quarter and maximum temperature in the warmest month; these data are averages of annual statistics over the period of 1970 to 2000. We then used the lfmm v2.0 package (Frichot & Francois, 2015) in R to perform a latent factor mixed model (LFMM) to test the association between the log-transformed count values of every expressed transcript and these two environmental variables. We applied a Benjamini-Hochberg correction to all p-values from the LFMM.

| Coordination in gene expression
To identify genes with coordinated (coassociated) expression, we calculated proportionality (ρ) between all pairs of transcripts in our dataset using the propr package (Quinn, Richardson, et al., 2017).

A full description of this analysis is available in the Supplementary
Information.

| Annotation and gene ontology enrichment
We performed gene ontology (GO) enrichment analysis to identify the most enriched (recurring) biological functions in which our transcripts were involved. This allowed us to determine whether certain functional categories from the GO database were overrepresented in our DE, fuzzy clustering (spatial expression), and proportionality (coordination and differential coordination) datasets more than would be expected by chance. We used the Bioconductor tool (Huber et al., 2015)

| Identification of immune genes
To identify additional transcripts that have known functions in the immune system (outside of those with the largest DE effect size or cluster membership), we cross-matched the output transcript lists from all of our analyses with several lists of known immune genes within the database InnateDB (Breuer et al., 2013)

| Isolation by distance
To assess the effect of geographic distance on divergence in gene expression (thereby testing for isolation by distance), we performed a Mantel test using the ade4 package (Thioulouse & Dray, 2007). A full description of this analysis is available in the Appendix.

| Identification of differentially expressed genes and their expression patterns
Overall, our DE analysis revealed 1,151 transcripts that were differentially expressed between invasion phases. These consisted  Table 1.

| Immune genes
The most common immune functions that we found among our differentially expressed transcripts were activation (list of transcripts and their expression patterns/clusters in Table 2a) and suppression (Table 2b) of inflammatory pathways, and cytotoxicity (Table 2c).
Immune genes with other roles were too few to allow clear inferences to be drawn. Most proinflammatory transcripts were downregulated in toads from intermediate areas relative to toads from the range core and invasion front. This trend is reflected in the spatial expression data; in the fifth cluster (low expression in intermediate areas, high expression on either end of the range), approximately onequarter of transcripts were identified as relating to immune function (the highest proportion of any cluster), and this is the only cluster in which a GO term directly related to immunity was among the most significant ( Figure 3, Table 1) Table 1). Nonetheless, a few pro-and anti-inflammatory transcripts were upregulated in intermediate toads (Table 2).

| Climate-influenced gene expression
Our LFMM revealed eleven transcripts with expression levels associated with maximum temperature during the hottest month, rainfall during the driest quarter, or both (list of transcripts in Appendix S2). Three of these transcripts followed the expression pattern of the first cluster (low expression at the range core, high expression throughout the rest of the range): Two are involved in cell adhesion and platelet activity, and function of the third is unknown. Two other transcripts were also downregulated at the core (but not a member F I G U R E 3 REVIGO plots displaying gene ontology (GO) terms (concepts/classes used to characterize gene function) depicting biological processes associated with transcripts following six major expression patterns in cane toads (Rhinella marina) across their Australian range (Figure 1). RNA-Seq data from spleens were used to identify differentially expressed transcripts between invasion phases, then soft clustering was performed to visualize the expression patterns that these transcripts follow (Figure 2

| Coordination in gene expression
We tested whether some transcripts were coassociated across invasion phases by examining the expected value of ρ metric (Quinn, Richardson, et al., 2017), but only found large groups of coassociated transcripts involved in fundamental cellular processes such as translation ( Figure S3; list of proportional transcripts in Appendix S3).

| Isolation by distance
Our Mantel test revealed a significant relationship between geographic distance and gene expression distance (p = 0.005, R 2 = 0.02; Figure 4).

| D ISCUSS I ON
The expression patterns of our differentially expressed genes suggest that multiple evolutionary forces may be at work. The curvilinearity in expression of many of our differentially expressed genes resembles that of other phenotypic traits affected by spatial sorting, such as leg length (Hudson et al., 2016).
Physical activity can have transgenerational effects on gene expression (Barres & Zierath, 2016;Murashov et al., 2016), so dispersal may directly affect expression patterns. Additionally, the significant result of our Mantel test could be due to either genetic drift or a balance between geographically varying selection and gene flow (Endler, 1977). Although this result may be influenced by environmental gradients, it is unlikely that selection driven by environmental factors would act on a genome-wide level (all expressed transcripts were used in our calculations for gene expression distance). This result suggests that genetic drift may be responsible for isolation by distance in Australian cane toads.
Because spatial sorting and genetic drift drive nonadaptive variation, their effects may obscure adaptive variation, particularly when they act on the same traits as selection (i.e., physical activity is also linked to expression of inflammatory genes; Baynard, Vieira-Potter, Valentine, & Woods, 2012; Gjevestad, Holven, & Ulven, 2015).
We also searched for evidence of adaptation to environment across the Australian invasive range, which is more similar to the native range in the east of Australia, and warmer and drier in the west. The expression pattern depicted in our first cluster (low expression in the core, high throughout the rest of the range) matches the trend across the Australian range in temperature and is opposite to the trend in rainfall. Population genetic structure in toads may be driven by these environmental variables .
Three of the eleven transcripts with expression levels associated with temperature or rainfall belonged to this first cluster; the rest were not members of any of the six clusters. Two out of these three are involved in blood clotting, as are most of the transcripts in the first cluster. Proportionality analysis revealed two small groups of clotting-related transcripts with coordinated expression. Blood clotting is affected by hydration levels (El-Sabban, Fahim, Al Homsi, & Singh, 1996), and excessive blood clotting can impair health (PubmedHealth, 2014). Changes to the rate of blood clotting may reflect hydric-related adaptations of intermediate and frontal toads living in drier conditions. Furthermore, four candidate SNP loci (with outlier F ST values and associations with rainfall, indicating they may TA B L E 1 Most common functions of transcripts in each of the six major expression patterns in cane toads (Rhinella marina) across their Australian range (Figure 1). RNA-Seq data from spleens were used to identify differentially expressed transcripts between invasion phases; then, soft clustering was performed to visualize the expression patterns that these transcripts follow (Figure 2)  be under selection) have previously been found in a gene involved in blood clotting . However, if lowering blood clotting rates is indeed an adaptation to aridity in toads from intermediate areas and the front, then one might expect that transcripts promoting blood clotting would be downregulated in these individuals (yet in our study, they are upregulated). This may be because toads were collected from the wild and the environment was not controlled for. A common-garden experiment with toads from across the range may reveal whether modification of blooding clotting rates is an evolved adaptation to aridity or simply a physiological reaction to different environments. Nonetheless, the differences in the expression of genes underlying these traits supports the hypothesis that cane toads are adapting to their abiotic environment.

Most significantly enriched GO term(s) Less common biological functions
The conceptual scheme of Lee and Klasing (2004) was not assayed, manipulated nor controlled for. Thus, spatial heterogeneity in pathogen pressure or environmental conditions may play a role in gene expression; however, parasite prevalence is generally higher in intermediate-age populations than in younger and older populations (Freeland et al., 1986), inconsistent with the general pattern for downregulation of immune transcripts in toads from intermediate-age populations. Our results may be driven by other environmental variables or may have a genetic basis through polymorphism in promoter regions, which are not sequenced when using RNA-Seq (Wang, Gerstein, & Snyder, 2009).
Our findings are consistent with previous phenotypic data on toads: Spleen sizes and fat body masses also follow a curvilinear pattern in which they are lowest in toads from intermediate parts of the range (Brown, Kelehear, et al., 2015b). If toads indeed follow a "traveling wave" density pattern, then the increased densities that occur several years postcolonization may not just heighten pathogen-mediated selection, but also intraspecific competition; reduced access to food resources may thus constrain the toad's ability to increase investment in immunity (Brown, Kelehear, et al., 2015b). This may explain smaller spleen sizes and lower expression of transcripts involved in inflammatory immune responses in toads from intermediate areas. Availability of energetic resources was not considered in the predictions of Lee and Klasing (2004), and although individuals at an invasion front may have a reduced benefit from energetic investment into immunity, they may invest anyway due to reduced intraspecific competition for food. Thus, the strength of immune functioning in invasion front individuals may depend on factors other than enemy release, such as prey abundance and interspecific competition.

| CON CLUS IONS
Although natural selection causes adaptive change in immune function in response to infection, this can co-occur with selection by abiotic environmental factors, as well as nonadaptive variation driven by forces such as spatial sorting and genetic drift, particularly in invasive species.
Because these forces can exert different effects on the same traits, predicting trait variation based solely on selection may not be an accurate approach. Nonetheless, methods such as RNA-Seq remain a powerful tool for uncovering the diverse and sometimes opposing forces that underpin rapid evolution in invasive species. The expression patterns that we observed in pro-and anti-inflammatory transcripts generally do not support the predictions of Lee and Klasing (2004) based on the ERH. Notably, our study suggests that the predicted consequences of enemy release are difficult to study in the wild, especially in the absence of infection data, because host-parasite associations are affected by many factors. Common-garden studies could be used to test these alternative hypotheses, but it is worth noting that keeping animals in captivity may also have functional impacts, sometimes making it difficult to interpret these results in the context of wild populations. These complicated dynamics are likely to explain why support for the predictions of ecological hypotheses, including ERH, is often inconsistent.
F I G U R E 4 Correlation between geographic distance and gene expression distance of invasive cane toad (Rhinella marina) populations across their Australian range (Figure 1). Euclidean distances in geographic location and gene expression between populations were calculated using the dist function in R. A mantel test (performed with the ade4 package) confirmed that these were significantly correlated (p = 0.005)