Brain transcriptome analysis reveals gene expression differences associated with dispersal behaviour between range-front and range-core populations of invasive cane toads in Australia.

Understanding the mechanisms allowing invasive species to adapt to novel environments is a challenge in invasion biology. Many invaders demonstrate rapid evolution of behavioural traits involved in range expansion such as locomotor activity, exploration and risk-taking. However, the molecular mechanisms that underpin these changes are poorly understood. In 86 years, invasive cane toads (Rhinella marina) in Australia have drastically expanded their geographic range westward from coastal Queensland to Western Australia. During their range expansion, toads have undergone extensive phenotypic changes, particularly in behaviours that enhance the toads' dispersal ability. Common-garden experiments have shown that some changes in behavioural traits related to dispersal are heritable. At the molecular level, it is currently unknown whether these changes in dispersal-related behaviour are underlain by small or large differences in gene expression, nor is known the biological function of genes showing differential expression. Here, we used RNA-seq to gain a better understanding of the molecular mechanisms underlying dispersal-related behavioural changes. We compared the brain transcriptomes of toads from the Hawai'ian source population, as well as three distinct populations from across the Australian invasive range. We found markedly different gene expression profiles between the source population and Australian toads. By contrast, toads from across the Australian invasive range had very similar transcriptomic profiles. Yet, key genes with functions putatively related to dispersal behaviour showed differential expression between populations located at each end of the invasive range. These genes could play an important role in the behavioural changes characteristic of range expansion in Australian cane toads.

INTRODUCTION

Invasive species are known for their high fecundity, phenotypic plasticity and excellent dispersal ability, which are all traits that enable them to adapt to and expand across novel environments (Prentis et al., 2008). As a result, invaders can pose significant threats to native species, ecosystem services, public health, agriculture and global economies (Clavero & Garcia‐Berthou, 2005; Paini et al., 2016; Pimentel et al., 2005; Walsh et al., 2016). Behavioural traits are key to promoting range expansion of invasive species (Duckworth & Badyaev, 2007; Myles‐Gonzalez et al., 2015; Rehage & Sih, 2004). Across taxa, dispersal success is associated with a suite of correlated behaviours. Individuals that are more active and more prone to taking risks (e.g., exploring unfamiliar habitats) have a greater probability of dispersing further and for longer periods of time (Duckworth & Badyaev, 2007; Myles‐Gonzalez et al., 2015; Réale et al., 2007). Additionally, individuals with a greater propensity or aptitude for dispersal are more likely to be located on the expanding edge of populations than at the range‐core, further increasing the overall dispersal rate (Duckworth & Badyaev, 2007; Phillips et al., 2006; Shine et al., 2011; Sol et al., 2002). The rate and extent of range expansion is influenced by the interplay between genetic factors, environmental factors, and phenotypic plasticity (Holway & Suarez, 1999).

Phenotypic variation between individuals is the substrate for selection on heritable traits, and is thus essential for evolution. Phenotypic variation in a population can result from genetic diversity and/or phenotypic plasticity. Because invasive populations often have low genetic diversity as a direct result of founder effects during introduction (Rollins et al., 2013), phenotypic plasticity may be especially important in these populations. Plasticity can affect behaviour, physiology and morphology (Piersma & Drent, 2003), and can evoke trait variation that could accelerate evolution (Robinson & Dukas, 1999). For example, increased plasticity at the leading edge of an expanding population could increase the expression of advantageous traits (e.g., locomotor ability and exploratory behaviour), which in turn can enhance the rate of dispersal. Phenotypic plasticity may thus be a key driver behind the rapid rates of range expansion common during invasions (Chuang & Peterson, 2016; Sexton et al., 2009).

The cane toad (Rhinella marina), native to South America, is a notorious invasive species now found across many parts of the world (Lever, 2001). In Australia, cane toads were introduced from a source population in Hawai'i to multiple sites along 1200 km of Queensland coastline in 1935, in a failed attempt to control sugar cane beetles (Easteal, 1981; Shine, 2014; Shine et al., 2020). Since its introduction, the cane toad has expanded southward to New South Wales and westward through the Northern Territory, reaching Western Australia by 2010 (Shine, 2010; Urban et al., 2008). Across the western half of that invasion trajectory, climatic conditions are hotter and seasonally much more arid compared to both the toad's native range and its Queensland introduction sites (Kearney et al., 2008; Kosmala et al., 2020).

Genetic variation in invasive Australian cane toads is low (Lillie et al., 2017; Selechnik, Richardson, Shine, DeVore, et al., 2019; Slade & Moritz, 1998). Interestingly, although genome‐wide genetic diversity in Australia is overall greatly reduced as compared to the native range, diversity is increased at some loci that appear to be under selection and that are associated with climatic conditions (Selechnik, Richardson, Shine, DeVore, et al., 2019).

Moreover, substantial variation in morphology, physiology and behaviour linked to dispersal ability have been documented across the toads’ westward invasion in Australia (Shine, 2010). For example, toads located on the invasion front have longer legs (Phillips et al., 2006), wider forelimbs, narrower hindlimbs, more compact skulls (Hudson, McCurry, et al., 2016), smaller relative head widths (Hudson et al., 2018), greater endurance (Llewelyn et al., 2010), move more frequently (Alford et al., 2009) and further in a given period of time (Lindstrom et al., 2013), and are more exploratory and more likely to exhibit risk‐taking behaviour in a novel environment than are individuals from range‐core populations (Gruber et al., 2017a). These phenotypic changes have culminated in increased dispersal ability in toads from range‐front populations (Urban et al., 2008).

Common‐garden experiments have shown that some of the phenotypic traits showing extensive variation across Australia are heritable (e.g., morphological, Hudson, Brown, et al., 2016; Hudson et al., 2018; physiological, Brown, Phillips, et al., 2015; Kosmala et al., 2018; or behavioural, Gruber et al., 2017b; Stuart et al., 2019). This indicates a genetic basis underlying the significant phenotypic variation seen in cane toads across the Australian invasive range.

With respect to plasticity in Australian cane toads, studies have yielded mixed results. Less plasticity in growth and development was found in individuals from the western range‐front vs. those from the range‐core (Ducatez et al., 2016). When juvenile toads were exposed to high vs. low exercise regimes, body size in range‐front toads was less plastic compared to that of range‐core toads, but this trend was reversed in toads whose diets were supplemented with calcium (Stuart et al., 2019). Greater plasticity with respect to temperature tolerance was found in toads from the southern range‐front vs. those from the range‐core (Kolbe et al., 2010; McCann et al., 2014). These studies demonstrate that patterns of plasticity in morphological and physiological traits are complex in toads sampled across the Australian invasive range. To date, little is known about behavioural plasticity in these populations.

Furthermore, it is poorly understood whether the large phenotypic response in morphology, physiology and behaviour is mirrored at the molecular level. Muscle and spleen transcriptome analyses found hundreds of differentially expressed genes (hereafter, DEGs) between range‐core and range‐front populations. These genes are putatively involved in metabolism, cellular repair and immune function, which might indicate a molecular response to greater environmental stressors at the range edge (Rollins et al., 2015; Selechnik, Richardson, Shine, Brown, et al., 2019). The magnitude of gene expression differences underlying changes in dispersal‐related behaviour remains unknown. Behaviour is a particularly complex phenotype at the interplay between genetic and environmental factors. It is thought that specific gene regulatory networks interact with neuronal networks within particular brain regions to orchestrate behavioural responses to both internal (e.g., hormonal) and external (e.g., environmental) stimuli (Sinha et al., 2020). Thus, investigating brain gene regulation is fundamental to our understanding of behaviour (Sinha et al., 2020). Yet, no study has investigated brain gene expression differences between range‐core and range‐front populations.

Here, we used RNA‐seq to analyse brain transcriptomes of toads sampled from Hawai'i, the area from which toads were taken to Australia, as well as three geographically distinct populations from across the Australian invasive range. We hypothesised that there would be large differences in brain gene expression between the Hawai'ian source population and the Australian invasive population, and that there would be a decrease in gene expression variability in Australian toads compared to Hawai'ian toads. Further, considering the extent of phenotypic changes documented within the invasive Australian range, we hypothesised that there would also be at least moderate differences in brain gene expression and/or variability between different Australian populations. We finally hypothesised that changes in brain gene expression patterns across the Australian range would be linked to the timeline of cane toad invasion, as has been found for other dispersal‐related traits (Brown, Kelehear, et al., 2015; Selechnik, Richardson, Shine, Brown, et al., 2019). This sampling design allowed us to identify candidate genes underlying range expansion traits in Australian cane toads and to explore differences in variation in gene expression (i.e., plasticity) across these sampling sites.

MATERIALS AND METHODS

Sample collection

We collected brain tissue from a total of 54 wild, adult female cane toads during April and May of 2014 and 2015 from 13 locations (N = 4/5 individuals per location) representing four populations in Hawai'i and Australia (Figure 1 and Table S1). Populations corresponded to distinct genetic clusters across the invasive range following Selechnik, Richardson, Shine, DeVore, et al. (2019): source (Hawai'i), range‐core (coastal Queensland), intermediate (central Queensland and Northern Territory) and range‐front (Western Australia) (Figure 1). Immediately following capture, we humanely euthanized toads, extracted their brains and stored them in RNAlater (Qiagen) at 4°C to preserve tissue integrity until samples were transported to the lab, where they were stored at −80°C prior to RNA extraction. All experimental procedures involving live toads were approved by the University of Sydney Animal Care and Ethics Committee (2014/562) and the Deakin University Animal Ethics Committee (AEX04‐2014).

FIGURE 1

Location of samples. QLD, Queensland; NT, Northern Territory; WA, Western Australia. The shaded area represents the cane toad’s Australian invasive range

RNA extraction and sequencing

We extracted total RNA from whole brains using Qiagen RNeasy Lipid Tissue Mini Kits (Qiagen), following the manufacturer's protocol. We homogenised tissues using a Fast Prep‐24 Classic homogeniser (MP Biomedicals) and 1 mm Zirconia/Silica beads (Daintree Scientific) for 1 min at 6 m/s. We digested genomic DNA using a Qiagen RNase‐Free DNase set on column during extraction. We quantified extracted RNA using a Qubit RNA HS assay kit on a Qubit 3.0 Fluorometer (Life Technologies). Library preparation and sequencing was conducted commercially at Macrogen (South Korea). Library preparation followed the TruSeq mRNA 2 (Illumina) protocol and libraries were sequenced on an Illumina HiSeq 2500 platform (two lanes of 125 bp paired‐end sequencing), generating 872 million reads. Average RNA Integrity Number values for all samples were 8.0 ± 1.0 (mean ± SD).

Data preprocessing, alignment and gene expression quantification

We used fastqc 0.11.7 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) to check the quality of the raw data. trimmomatic 0.38 (Bolger et al., 2014) was used to remove adaptor sequences and trim low quality reads with the following parameters: ILLUMINACLIP: path/to/TruSeq3‐PE.fa:2:30:10:4 HEADCROP:13 AVGQUAL:30MINLEN:36. We mapped trimmed reads to the multitissue reference cane toad transcriptome (Richardson et al., 2018) using star 2.7.2b (Dobin et al., 2013) in two‐pass mode with default parameters. We used resultant BAM files to quantify gene expression using salmon 1.2.1 (Patro et al., 2017).

The number of raw reads, the number of post‐trimming reads and the number of mapped reads were similar between all populations (source: respectively 11.2–14.9 million reads, 7.2–9.9 million reads and 2.3–3.8 million reads; range‐core: respectively 11.2–13.2 million reads, 7.5–9.3 million reads and 3.2–4.0 million reads; intermediate: respectively 10.6–14.1 million reads, 7.1–9.9 million reads and 2.5–4.1 million reads; range‐front: respectively 10.9–14.1 million reads, 7.4–9.6 million reads and 3.1–4.0 million reads; Table S2).

Differential expression analysis

We used edger 3.32.1 (Robinson et al., 2010) to filter out genes using the filterByExpr function with default parameters in r 4.0.4 (R Core Team, 2021). We further filtered out genes that had <10 counts per million in at least 10 samples. To assess the presence of outliers, we normalised and rlog‐transformed counts before computing pairwise correlation for all the samples. We then used the resultant correlation matrix to plot a heatmap with pheatmap 1.0.12 (Kolde, 2019). This revealed two outliers (Figure S1), B24 and B31, that we excluded from subsequent analyses. We performed differential expression analysis using deseq2 1.30.1 (Love et al., 2014). We performed differential expression analysis between the source population and each population of the Australian range, and second between all three populations of the Australian range. We considered all genes with Benjamini‐Hochberg adjusted p‐values <.05 (Benjamini & Hochberg, 1995) to be significantly differentially expressed between any pairwise comparison. We further performed differential expression analysis using a likelihood ratio test with deseq2 between all nine Australian locations ordered along an east‐to‐west transect. This transect reflects the timeline of the cane toad invasion and was aimed at testing whether differences in brain gene expression followed a continuum across the Australian range rather than population‐specific changes. p‐values were obtained by comparing the full model (i.e., with each toad assigned to a geographic location) to a reduced model (i.e., without the geographic location factor). We then used the degPatterns function in deseq2 on all DEGs across this transect to identify clusters of genes with similar expression profiles across the Australian range. This tool performs a hierarchical clustering based on gene pairwise Kendall correlations.

Differential variability analysis

We used mdseq 1.0.5 (Ran & Daye, 2017) to test whether brain genes differed in their expression variability (hereafter, dispersion) both between the source population and the Australian population as a whole, and across the Australian range. Gene expression variability is being recognised as an important driver of phenotypic differences (Ecker et al., 2018). We combined all three Australian populations (i.e., range‐core, intermediate and range‐front) to compare gene dispersion with the source population. We then conducted differential dispersion analysis separately between the three Australian populations. For both dispersion analyses, we randomly selected two locations (Croydon and Mataranka) of the intermediate Australian population to balance sample sizes with the range‐core and range‐front populations. We normalised gene counts using the trimmed mean of M‐values (TMM) method (Robinson & Oshlack, 2010) in edger. Any gene with a Benjamini‐Hochberg adjusted p‐value <.05 was considered to be significantly differentially dispersed between any pairwise comparison.

Functional analysis

We performed gene ontology (GO) analysis using goseq 1.26.0 (Young et al., 2010) to infer biological function of all differentially expressed or differentially dispersed gene sets. We conducted GO analysis using the probability weighting function to adjust for transcript length bias and the Wallenius approximation to test for over‐representation. p‐Values were adjusted with the Benjamini‐Hochberg method. We visualised GO results using enrichplot (Yu, 2021).

RESULTS

Cane toads from the Australian range showed extensive differences in brain gene expression compared to toads from the Hawai'ian source population. Out of 15,622 filtered genes, there were 6529 (41.8%) DEGs between range‐core and source populations, 7769 (49.7%) DEGs between intermediate and source populations, and 6770 (43.3%) DEGs between range‐front and source populations (Figure 2a–f). The number of upregulated genes was higher in each Australian population compared to the source population (Wilcoxon test: W = 9, N = 6, p = .040). A core set of 4903 genes were commonly differentially expressed between the source population and each Australian population (Figure 2g). GO analysis of these 4903 DEGs revealed significant enrichment for processes such as translation and mitochondrial function (Figure 2h).

FIGURE 2

(a–c) Heatmap of normalised gene expression values for all DEGs between (a) range‐core and source populations, (b) intermediate and source populations, and (c) range‐front and source populations. Columns correspond to samples. Rows correspond to genes. Colour depicts Z‐score normalised gene expression value. (d–f) Volcano plots of significantly DEGs between (d) range‐core and source populations, (e) intermediate and source populations, and (f) range‐front and source populations. Nonsignificant genes are represented in grey. (g) Overlap of DEGs across each pairwise comparison. (h) GO analysis of the core set of DEGs overlapping across each pairwise comparison. The size of each circle is proportional to the number of genes being significantly enriched, while the colour of each circle is proportional to its FDR‐corrected p‐value. Gene ratio corresponds to the proportion of genes being enriched out of the total number of genes in that GO category. (i) Volcano plot of significantly differentially dispersed genes between Australian toads and Hawai'ian toads. Nonsignificant genes are represented in grey

There was a modest reduction in brain gene expression variability in toads from the Australian range compared to toads from the Hawai'ian source population. Out of 63 differentially dispersed genes between Australian toads and Hawai'ian toads, 52 (82.5%) genes were over‐dispersed in Hawai'ian toads, while only 11 (17.5%) genes were over‐dispersed in Australian toads (Figure 2i). GO enrichment analysis failed to find any GO term being significantly over‐represented in both populations.

By contrast, cane toads within the Australian range showed few geographic differences in brain gene expression. Out of 15,023 filtered genes, there were only 59 (0.4%) DEGs between intermediate and range‐core populations (Table 1 and Figure 3a,c), and 21 (0.1%) DEGs between range‐front and range‐core populations (Table 2 and Figure 3b,d). Out of 71 DEGs in both pairwise comparisons, nine (12.7%) genes were commonly differentially expressed in intermediate vs. range‐core populations and range‐front vs. range‐core populations (Figure 3e). Specifically, F12, FANCD2 and POL were upregulated in intermediate and range‐front populations compared to range‐core populations, while POMC, L1RE1, MAN2B1 and POL4 were downregulated in intermediate and range‐front populations compared to range‐core populations. There was no significant enrichment for GO terms in any pairwise comparison. Range‐front and intermediate populations showed no significant differences in brain gene expression.

TABLE 1

DEGs between intermediate and range‐core populations. Genes respectively up‐ and downregulated in intermediate vs. range‐core populations are indicated by log2 fold change values respectively >0 and <0

Gene	Protein	Log2 fold change	FDR‐corrected p‐value
Genes upregulated in intermediate vs. range‐core populations
TTC3	E3 ubiquitin‐protein ligase TTC3 isoform X4	2.00	.0163
F12	Coagulation factor XII	3.07	.0192
TOP1 (isoform 1)	DNA topoisomerase 1	1.67	.0193
TOP1 (isoform 2)	DNA topoisomerase 1	0.91	.0193
FANCD2 (isoform 1)	Fanconi anaemia group D2 protein	0.80	.0193
POL (isoform 1)	Pol polyprotein	0.72	.0193
NUP155	Nuclear pore complex protein Nup155	0.47	.0193
LKAAEAR1	Protein LKAAEAR1	2.49	.0217
NUCB2	Nucleobindin‐2	0.87	.0217
Rm3131d8773450t8	Unknown	2.46	.0269
PLPBP	Proline synthase cotranscribed bacterial homologue protein	0.43	.0277
Rm63737d1184914t1	Unknown	1.74	.0291
FANCD2 (isoform 2)	Fanconi anaemia group D2 protein	1.23	.0355
SLC6A17	Sodium‐dependent neutral amino acid transporter SLC6A17	0.54	.0355
FCRL3	Fc receptor‐like protein 3	0.64	.0361
POL (isoform 2)	Pol polyprotein	0.65	.0367
Rm44170d1321510t1	Unknown	0.79	.0385
OTX2	Orthodenticle homeobox protein OTX2	0.96	.0386
LMX1A	LIM homeobox transcription factor 1‐alpha	0.72	.0386
PRPF6	Pre‐mRNA‐processing factor 6	0.37	.0386
SHB	SH2 domain‐containing adapter protein B	0.48	.0449
PGRP‐SC2	Peptidoglycan‐recognition protein SC2	1.87	.0497
SRXN1	Sulfiredoxin‐1	0.59	.0497
IL13RA1	Interleukin‐13 receptor subunit alpha‐1	0.59	.0497
TSEN15	tRNA‐splicing endonuclease subunit Sen15	0.46	.0497
HSP90AA1	Heat shock protein HSP 90‐alpha	0.36	.0497
Genes downregulated in intermediate vs. range‐core populations
POMC	Pro‐opiomelanocortin	−6.80	<.00001
MAN2B1	Lysosomal alpha‐mannosidase	−2.43	.0002
L1RE1	LINE‐1 retrotransposable element ORF2 protein	−5.42	.0018
CRYM	Ketimine reductase mu‐crystallin	−0.97	.0050
FAM168A	Protein FAM168A	−1.00	.0193
NCOR2	Nuclear receptor corepressor 2	−0.29	.0217
TRIM25	E3 ubiquitin/ISG15 ligase TRIM25	−1.29	.0247
CALB1	Calbindin	−3.00	.0247
Rm74787t4	Unknown	−3.28	.0247
NRP2	Neuropilin‐2	−0.46	.0263
IFIT5 (isoform 1)	Interferon‐induced protein with tetratricopeptide repeats 5	−0.75	.0265
Rm2405t2	Unknown	−0.81	.0268
SPACA6	Sperm acrosome membrane‐associated protein 6	−1.38	.0269
LEMD2	LEM domain‐containing protein 2	−0.55	.0291
HPS5	Hermansky‐Pudlak syndrome 5 protein homologue	−0.62	.0293
IFIT5 (isoform 2)	Interferon‐induced protein with tetratricopeptide repeats 5	−0.73	.0293
IFIT5 (isoform 3)	Interferon‐induced protein with tetratricopeptide repeats 5	−0.83	.0293
MCOLN3	Mucolipin‐3	−1.96	.0299
ASMTL	N‐acetylserotonin O‐methyltransferase‐like protein	−0.52	.0361
GREB1L	GREB1‐like protein	−0.96	.0361
MPO	Myeloperoxidase	−1.92	.0361
DOCK4	Dedicator of cytokinesis protein 4	−0.45	.0367
TBX2	T‐box transcription factor TBX2	−1.47	.0367
AIM1L	Absent in melanoma 1‐like protein	−2.26	.0367
NAPA	Alpha‐soluble NSF attachment protein	−0.27	.0386
IFIT5 (isoform 4)	Interferon‐induced protein with tetratricopeptide repeats 5	−0.87	.0392
NGEF	Ephexin‐1	−0.48	.0401
TRIM8	Probable E3 ubiquitin‐protein ligase TRIM8	−1.22	.0433
POL4	Retrovirus‐related Pol polyprotein from transposon 412	−1.61	.0438
ATP2B1	Plasma membrane calcium‐transporting ATPase 1	−0.84	.0438
LRRK2	Leucine‐rich repeat serine/threonine‐protein kinase 2	−1.83	.0438
NAB2	NGFI‐A‐binding protein 2	−0.57	.0497
CYLD	Ubiquitin carboxyl‐terminal hydrolase CYLD	−0.97	.0497

FIGURE 3

(a,b) Heatmap of normalised gene expression values for all DEGs between (a) intermediate and range‐core populations, and (b) range‐front and range‐core populations. Columns correspond to samples. Rows correspond to genes. Colour depicts Z‐score normalised gene expression value. (c,d) Volcano plots of significantly DEGs between (c) intermediate and range‐core populations, and (d) range‐front and range‐core populations. Nonsignificant genes are represented in grey. (e) Overlap of DEGs across each pairwise comparison

TABLE 2

DEGs between range‐front and range‐core populations. Genes respectively up‐ and downregulated in range‐front vs. range‐core populations are indicated by log2 fold change values respectively >0 and <0

Gene	Protein	Log2 fold change	FDR‐corrected p‐value
Genes upregulated in range‐front vs. range‐core populations
F12	Coagulation factor XII	3.98	.0022
POL (isoform 1)	Pol polyprotein	1.72	.0036
APOA4	Apolipoprotein A‐IV	1.52	.0073
ZNF16	Zinc finger protein 16	2.33	.0185
FANCD2 (isoform 1)	Fanconi anaemia group D2 protein	0.90	.0227
POL (isoform 2)	Pol polyprotein	0.82	.0227
POL (isoform 3)	Pol polyprotein	1.70	.0325
FOXJ1B	Forkhead box protein J1‐B	1.27	.0330
POL (isoform 4)	Pol polyprotein	0.79	.0330
FZD10	Frizzled‐10‐A	1.40	.0353
FANCD2 (isoform 2)	Fanconi anemia group D2 protein	1.42	.0379
Genes downregulated in range‐front vs. range‐core populations
POMC	Pro‐opiomelanocortin	−6.24	.0008
L1RE1 (isoform 1)	LINE‐1 retrotransposable element ORF2 protein	−6.26	.0033
L1RE1 (isoform 2)	LINE‐1 retrotransposable element ORF2 protein	−0.81	.0036
MAN2B1	Lysosomal alpha‐mannosidase	−2.40	.0062
SLC35E2	Solute carrier family 35 member E2	−2.34	.0090
EPB41L4A	Erythrocyte membrane protein band 4.1‐like protein 4A	−0.92	.0227
TUBB4B	Tubulin beta‐4B chain	−0.39	.0227
POL4	Retrovirus‐related Pol polyprotein from transposon 412	−2.17	.0290
PCDHGB1	Protocadherin gamma‐B1	−2.05	.0290
SLC43A2	Large neutral amino acids transporter small subunit 4	−0.44	.0330

There were 228 DEGs across all nine Australian locations ordered along an east‐to‐west transect (top 50 DEGs are shown in Table 3). Of these, 31 genes had been previously identified in this study as showing significant differences in expression between populations, including POMC, MAN2B1, CALB1, NUCB2 and LRRK2. Among these 228 DEGs, we identified four clusters of genes (together comprising 113, i.e., 49.6%, genes) showing similar expression profiles across the Australian range. These four clusters showed curvilinear patterns of gene expression whereby gene expression in intermediate areas differed from that in range‐core and range‐front areas (Figure 4a–d). Clusters 2 and 3 showed significant enrichment for GO terms such as RNA‐mediated transposition, cyclooxygenase pathway and oxidation‐reduction process (Figure 4e,f).

TABLE 3

Top 50 genes showing significant expression changes across the whole Australian range. Genes in bold also show significant expression differences between populations

Gene	Protein	FDR‐corrected p‐value
POMC	Pro‐opiomelanocortin	<.00001
MAN2B1	Lysosomal alpha‐mannosidase	<.00001
Rm16162d1119677t2	Unknown	<.00001
EXT2	Extensin‐2	<.00001
MPO	Myeloperoxidase	<.00001
AKAP10	A‐kinase anchoring protein 10	.00003
Rm73203d1719430t2	Unknown	.00004
JUND	Transcription factor jun‐D	.00006
CBR1 (isoform 1)	Carbonyl reductase [NADPH] 1	.00011
LOC120994811	Uncharacterized LOC120994811	.00012
ABCF2	ATP‐binding cassette subfamily F member 2	.00016
LOC120989780	Cell wall protein IFF6‐like	.00019
TRIM8 (isoform 1)	Probable E3 ubiquitin‐protein ligase TRIM8	.00040
Rm21589d67083t7	Unknown	.00047
MFN1	Mitofusin‐1	.00047
LOC121007588	Uncharacterized LOC121007588	.00049
TTC3	E3 ubiquitin‐protein ligase TTC3 isoform X4	.00049
DNAH7	Dynein heavy chain 7, axonemal	.00051
TOP1	DNA topoisomerase 1	.00062
DAP3	28S ribosomal protein S29, mitochondrial	.00067
TMED8	Protein TMED8	.00068
AN3	Putative ATP‐dependent RNA helicase an3	.00068
TRIM25	E3 ubiquitin/ISG15 ligase TRIM25	.00070
CBR1 (isoform 2)	Carbonyl reductase [NADPH] 1	.00071
TRIM8 (isoform 2)	Probable E3 ubiquitin‐protein ligase TRIM8	.00078
SARM1	Sterile alpha and TIR motif‐containing protein 1	.00083
TLE5	Amino‐terminal enhancer of split	.0013
L1RE1 (isoform 1)	LINE‐1 retrotransposable element ORF2 protein	.0016
NOVH	Protein NOV homologue	.0016
NSUN2	tRNA (cytosine(34)‐C(5))‐methyltransferase	.0019
L1RE1 (isoform 2)	LINE‐1 retrotransposable element ORF2 protein	.0019
DOCK4	Dedicator of cytokinesis protein 4	.0024
ANGEL1	Protein angel homologue 1	.0024
DNAJB5	DnaJ homologue subfamily B member 5	.0024
LEMD2	LEM domain‐containing protein 2	.0024
SNX29	Sorting nexin‐29	.0026
CRYM	Ketimine reductase mu‐crystallin	.0028
DZIP1	Zinc finger protein DZIP1	.0028
SLC35E2	Solute carrier family 35 member E2	.0028
IFT52	Intraflagellar transport protein 52 homologue	.0029
CAD	CAD protein	.0032
CBR1 (isoform 3)	Carbonyl reductase [NADPH] 1	.0033
RND2	Rho‐related GTP‐binding protein RhoN	.0036
PPA2	Inorganic pyrophosphatase 2, mitochondrial	.0036
SLC4A1	Band 3 anion transport protein	.0037
NR1I3	Nuclear receptor subfamily 1 group I member 3	.0040
ZRANB3	DNA annealing helicase and endonuclease ZRANB3	.0040
L1RE1 (isoform 3)	LINE‐1 retrotransposable element ORF2 protein	.0044
APOA1	Apolipoprotein A‐I	.0044
ARHGAP35	Rho GTPase‐activating protein 35	.0044

FIGURE 4

(a–d) Z‐score abundance of gene expression of DEGs across the Australian range. Locations on the x‐axis are ordered from east to west. Genes showing similar patterns of gene expression are grouped together. Box plots represent median, interquartile range and 95% confidence interval. Black lines represent the trend in expression change. Colours correspond to populations (red, range‐core; blue, intermediate; green, range‐front). (e,f) GO analysis of DEGs belonging to clusters 2 (e) and 3 (f). The size of each circle is proportional to the number of genes being significantly enriched, while the colour of each circle is proportional to its FDR‐corrected p‐value. Gene ratio corresponds to the proportion of genes being enriched out of the total number of genes in that GO category

Toads within the Australian range also showed moderate geographic differences in the variability of brain gene expression. There were 17 differentially dispersed genes between intermediate and range‐core populations (Figure 5a and Table S3), 14 differentially dispersed genes between range‐front and range‐core populations (Figure 5b and Table S4), and 29 differentially dispersed genes between range‐front and intermediate populations (Figure 5c and Table S5). The number of over‐dispersed genes was similar in both populations in each pairwise comparison (Figure 5a–c and Tables S3–S5). Thus, changes in gene expression variability were not characteristic of any population within the Australian range. Five genes were commonly differentially dispersed in intermediate vs. range‐core populations and range‐front vs. range‐core populations (Figure 5d), namely WRB, MAN2B1, POL, PAK2 and L1RE1. Four genes were commonly differentially dispersed in range‐front vs. intermediate populations and intermediate vs. range‐core populations (Figure 5d), namely MCOLN3, LRRK2 and two uncharacterised genes. One gene, NPW, was commonly differentially dispersed in range‐front vs. intermediate populations and range‐front vs. range‐core populations (Figure 5d). No gene was commonly differentially dispersed in all three pairwise comparisons. There was no significant enrichment for any GO term in any pairwise comparison. Six genes, MCOLN3, MAN2B1, LRRK2, EPB41L4A, NUCB2 and L1RE1, were both differentially expressed and differentially dispersed between at least two populations (Figure 5e).

FIGURE 5

(a–c) Volcano plots of significantly differentially dispersed genes between (a) intermediate and range‐core populations, (b) range‐front and range‐core populations, and (c) range‐front and intermediate populations. Nonsignificant genes are represented in grey. (d) Overlap of differentially dispersed genes across each pairwise comparison. (e) Overlap of differentially expressed genes and differentially dispersed genes across all pairwise comparisons

DISCUSSION

In this study, we set out to investigate transcriptomic changes in cane toads’ brains associated with expansion across its Australian invasive range. We found approximately five thousand DEGs between the Hawai'ian source population and the Australian invasive population. Therefore, extensive transcriptomic differences exist between the source and Australian population, despite a general similarity in genetic composition between source and range‐core populations (F ST = 0.04; Selechnik, Richardson, Shine, DeVore, et al., 2019). We also found a decrease in brain gene expression variability in Australian toads compared to Hawai'ian toads. These differences in gene expression patterns probably reflect abiotic and biotic differences between the Hawai'ian source and Australian populations. Interestingly, more genes were upregulated than downregulated in the Australian invasive population than in the Hawai'ian source population. A general increase in brain gene expression in the Australian invasive population, particularly in genes involved in metabolism and mitochondrial function, might indicate a response to a novel and potentially more stressful environment (Kosmala et al., 2020; Tingley et al., 2014).

In contrast, the brain transcriptome of toads was remarkably similar across the entire Australian invasive range. Intermediate and range‐front populations showed only a few dozen DEGs compared to the range‐core population, and no significant difference in gene expression with each other. The lack of gene expression differences observed between intermediate and frontal populations is compatible with their genetic structure, because these two populations are part of the same genetic cluster (Selechnik, Richardson, Shine, DeVore, et al., 2019). Looking at the overall transcriptomic response across all locations ordered along an east‐to‐west transect showed essentially the same result. While a few more DEGs were identified, the overall gene expression pattern remained consistent along the entire transect. This pattern confirms that the general similarity in brain transcriptomes identified between range‐core, intermediate and range‐front populations is not an artefact of arbitrarily defined populations, but rather a biological feature of the Australian cane toad invasion.

Focusing on DEGs, we found that a few genes showed parallel regulatory changes in intermediate and range‐front vs. range‐core populations. However, most genes being differentially expressed in intermediate and range‐front populations compared to range‐core populations were nonoverlapping, indicating that intermediate and range‐front populations do show distinct transcriptomic profiles. This was also evident when looking at gene expression across all nine Australian locations ordered along an east‐to‐west transect, where about half of DEGs showed distinct expression patterns in intermediate vs. range‐front and range‐core populations. These genes showed a curvilinear pattern of gene expression along the range‐core to range‐front transect, with gene up‐ or downregulation in intermediate populations compared to both range‐core and range‐front populations. Curvilinear patterns across the invasive range have been found for cane toads’ spleen gene expression (Selechnik, Richardson, Shine, Brown, et al., 2019), spleen mass, fat body mass, lungworm infection (Brown, Kelehear, et al., 2015), and limb length (Hudson, Brown, et al., 2016; Stuart et al., 2019). These curvilinear relationships with invasion history could be partly underlain by a “travelling wave” density pattern, where higher population densities in intermediate populations compared to range‐front and range‐core populations cause a change in selection on dispersal‐related traits (Brown, Kelehear, et al., 2015). If so, this could explain some of the changes in brain gene expression that differentiate intermediate populations from range‐front and range‐core populations.

The variability in brain gene expression was also very limited across the Australian invasive range, with only 14–29 differentially dispersed genes between any two populations. Further, we failed to detect any sign of a change in gene expression variability from range‐core to range‐front populations, and therefore find no genome‐wide support for changes to plasticity across this range. Although variability in gene expression can affect phenotypic differences (Ecker et al., 2018), this factor does not seem to have been important in cane toads. Future studies comparing gene expression variability in toads from across the Australian range and experimentally subjected to various biotic and abiotic stressors are necessary to confirm this finding.

The overall conservatism in brain gene expression seems to be at odds with the magnitude of phenotypic differences observed within the Australian invasive range (Brown, Kelehear, et al., 2015; Brown et al., 2013; Gruber et al., 2017a; Hudson, Brown, et al., 2016; Hudson et al., 2018; Llewelyn et al., 2010; Phillips et al., 2006; Pizzatto et al., 2017; Tingley et al., 2012; Urban et al., 2008). Nonetheless, invasiveness may sometimes be promoted by a small number of genes (Bock et al., 2015), and minor changes in brain transcriptomes can underlie substantial differences in behaviour (Renn & Schumer, 2013; Saul et al., 2017). Moreover, other studies have reported low divergence in brain gene expression between populations with substantial morphological, behavioural and physiological differences, for example, between domesticated and wild mammals (Albert et al., 2012), between domesticated and wild zebrafish (Drew et al., 2012), or between African and European fruit flies (Catalan et al., 2012). In other words, the extensive, and partially heritable, phenotypic differences seen between range‐core, intermediate and range‐front cane toad populations could be underlain by only a few key genes.

In favour of this hypothesis, some of the genes showing differential expression and/or dispersion between range‐core, intermediate and range‐front populations have putative functions related to behaviour. Calbindin (CALB1), lysosomal alpha‐mannosidase (MAN2B1), leucine‐rich repeat serine/threonine‐protein kinase 2 (LRRK2) and mucolipin‐3 (MCOLN3) were all downregulated in intermediate populations (and range‐front populations for MAN2B1) compared to range‐core populations. CALB1 is involved in locomotor behaviour (Barski et al., 2003). CALB1 ko mice have lower anxiety‐like behaviour, increased exploratory behaviour, and are less prone to exhibiting freezing behaviour (Harris et al., 2016). MAN2B1 is involved in neurocognitive functions such as learning and memory, and plays a role in motor function (Damme et al., 2011; D'Hooge et al., 2005). LRRK2 is involved in exploration behaviour and is linked with Parkinson's Disease (Melrose et al., 2010). MCOLN3 is involved in locomotor behaviour, mutant mice showing erratic circling behaviour (Di Palma et al., 2002). MAN2B1, LRRK2 and MCOLN3 were further underdispersed in intermediate populations (and range‐front populations for MAN2B1) compared to range‐core populations. Both LRRK2 and MCOLN3 were also underdispersed in intermediate vs. range‐front populations. Sodium‐dependent neutral amino acid transporter (SLC6A17) and LIM homeobox transcription factor 1‐alpha (LMX1A) were upregulated in intermediate vs. range‐core populations. SLC6A17 mutations cause behavioural problems in humans (Iqbal et al., 2015). LMX1A is involved in memory and locomotor and olfactory behaviour, with mutations linked with Parkinson's Disease (Laguna et al., 2015). Calcipressin‐2 (RCAN2) and relaxin‐3 (RLN3) were underdispersed in range‐front vs. intermediate populations. RCAN2 is involved in locomotor behaviour, stress responses and memory (Miyakawa et al., 2003). RLN3 is involved in exploratory behaviour, stress responses and the regulation of feeding behaviour (Smith et al., 2009). Therefore, changes in expression and/or dispersion in the above genes may contribute to the behavioural shift in dispersal‐related behaviour in toads from range‐front and intermediate populations vs. range‐core populations.

Genes involved in the regulation of feeding behaviour (other than the above‐mentioned RLN3) also showed differential levels and/or variability of expression between range‐core, intermediate and range‐front populations. These genes could also play a role in the phenotypic changes associated with range expansion in Australian cane toads. Indeed, toads from range‐front areas have higher feeding rates, larger fat bodies, better body condition and faster growth (Brown et al., 2013). Pro‐opiomelanocortin (POMC) was downregulated in intermediate and range‐front populations compared to range‐core populations. This gene plays a key role in the regulation of stress responses (Harno et al., 2018) and is also involved in the regulation of feeding behaviour (Millington, 2007). Nucleobindin‐2 (NUCB2) was upregulated in intermediate vs. range‐core populations, and was under‐dispersed in range‐front vs. range‐core populations. Neuropeptide W (NPW) was underdispersed in range‐front populations compared to range‐core and intermediate populations. Agouti‐related protein (AGRP) was over‐dispersed in range‐front vs. intermediate populations. NUCB2, NPW and AGRP have all been shown to play a role in the regulation of feeding behaviour (Dore et al., 2017; Mondal et al., 2003; Ollmann et al., 1997).

We must consider some potential limitations in our study that might explain our findings. Behaviour is a complex phenotype orchestrated by specific brain regions (Sinha et al., 2020). The low magnitude of gene expression changes that we observed in the Australian invasive population could thus be partially obscured by concomitant and opposing changes in gene expression across various brain regions (Nadler et al., 2006). Furthermore, heterogeneity in physiological and/or environmental conditions between individuals can confound any transcriptomic analysis, especially for individuals sampled in the wild (as in the present study). It is possible that by sampling adult toads, we missed the critical developmental window during which key genes underlying behavioural plasticity are differentially expressed (Aubin‐Horth & Renn, 2009). Populations might also differ in their ability to show plastic responses to short vs. long‐term stressors. Finally, additional molecular mechanisms inaccessible with RNA‐seq (such as post‐translational modifications) might play an important role in behavioural changes (Cash et al., 2005).

The observed geographic structuring of transcriptomic differences suggests two phases of divergence in Australian cane toads’ gene expression: (i) A first divergence across the Australian Great Dividing Range that occurred early during the invasion process, and that was potentially driven by differences in climatic conditions; and (ii) a second divergence in intermediate populations, potentially driven by a change in selection on dispersal‐related traits following an increase in population density. This phenomenon might be an example of rapid evolution, or a case of environmentally‐induced variation (phenotypic plasticity). We note, however, that phenotypic plasticity can also lead to adaptive evolution (Ghalambor et al., 2007), for example, through genetic assimilation (Pigliucci & Murren, 2003; West‐Eberhard, 2003). If evolution is at play, does it act through shifts in selective regimes (adaptive evolution), or through nonadaptive processes (e.g., drift, spatial sorting, admixture), or through a combination of factors?

In conclusion, the modest differences that we have documented in brain gene expression along the invasive range of cane toads within Australia might trigger significant changes in the toads’ phenotypic traits, in particular in relation to dispersal behaviour. Common‐garden experiments monitoring dispersal‐related behaviour together with gene expression in the identified genes would provide valuable insight into the heritability, and thus evolvability, of these phenotypic changes.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Mark F. Richardson and Lee A. Rollins conducted fieldwork. Andrea J. West, Mark F. Richardson and Daniel Selechnik conducted laboratory work. Boris Yagound, Andrea J. West and Mark F. Richardson conducted analyses. Boris Yagound, Andrea J. West, Mark F. Richardson, Daniel Selechnik, Richard Shine and Lee A. Rollins contributed to interpretation and writing.

Fig S1

Click here for additional data file.

Table S1‐S5

Click here for additional data file.