Genetic analysis reveals spatial structure in an expanding introduced rusa deer population

Study area

The Illawarra region (Fig. 1) is a coastal stretch of land south of Royal National Park adjacent to Sydney. The region is 7000 km² (AdaptNSW NSW Governmnent 2022) and ~320 000 people live there, including ~220 000 in the city of Wollongong (The Australian Bureau of Statistics 2020). The region is divided east–west by the Illawarra escarpment (or Illawarra Range; Fig. 1b), a steep outcrop ranging from 300 m above sea level (asl) in the north to >700 m asl in the south. The cliff-like slopes of the escarpment become more moderate in steepness above the city of Wollongong (Young 1979), where main roads cross them. Above the escarpment (to the west), there is a large plateau covered by continuous native bushlands dominated by tall eucalypt forests surrounding large artificial water reservoirs that supply Sydney and the Illawarra region (NPWS 2002). Public access is limited to maintain high water quality. The steep slopes of the escarpment consist of a mosaic of moist forest and rainforest communities, some unique to the Illawarra region (Ashcroft 2009). The narrow coastal plain between the escarpment and the Tasman Sea is almost completely urbanised north of the City of Wollongong, with high-density roading and housing covering most of the foothills of the escarpment (Fig. 1). Further south, the coastal plain widens into a grassy flood plain between the escarpment and the Lake Illawarra (Young 1979), and is less urbanised than north of the City of Wollongong.

The climate is mild marine, with 1500–1600 mm annual rainfall (NSW Government, Office of Environment and Heritage 2014). In summer, average temperatures range from 20 to 22°C along the coast and 18–20°C on the escarpment. In winter, temperatures range from 12 to 14°C along the coast and 4–6°C near the escarpment (NSW Government, Office of Environment and Heritage 2014).

Rusa deer tissue samples

During 2019 and 2020, ear tissue samples were collected from rusa deer shot during vehicle-based ground shooting operations organised by the IWDMP. The shooting methodology is detailed in Hampton et al. (2022). Operations were conducted at night and throughout the year except on weekends and public holidays. The shooters were asked to collect samples from deer spread throughout the study area, rather than from a few high deer-density sites. Tissue samples were collected using commercially available Allflex Tissue Sampling Units (TSU, Allflex Australia Pty Ltd, Capalaba, Queensland). For each tissue sample, we recorded the collection date, deer age (adult [>18 months] or juvenile [<18 months]), deer sex (based on external genitalia and the presence of antlers) and the GPS coordinates of the sampling site. To maximise coverage of the study area in our genetic data set, we selected samples evenly across all sites, while ensuring an equal sex ratio. Our final sample size of tissues collected for sequencing was 178.

DNA extraction and sequencing

DNA was extracted from ear tissue by Diversity Array Technology Pty Ltd (Canberra, Australia), using a NucleoMag Tissue Extraction Kit (Macherey-Nagel, Duren, Germany). DNA was double-digested using the restriction enzymes Pstl and Sphl, barcodes were added for disaggregation of samples and the fragments sequenced by DArT-seq (Diversity Array Technology Pty Ltd, Canberra) (Kilian et al. 2012; Georges and Gruber 2019). DNA could not be successfully extracted and sequenced for four deer, resulting in genetic data for a total of 174 deer.

Sequence data were analysed on the computational cluster at the University of New South Wales (PVC (Research Infrastructure), UNSW Sydney 2010). We used a reference-based variant calling approach to analyse fastq files using stacks v2.2 (Catchen et al. 2013). We used the program process_radtags in stacks to process the raw data to generate cleaned reads. We used process_radtags to check barcode intactness and for error correction, and then used fastqc v0.11.9 (Andrews 2010) to assess read quality using a sliding window (15% of the read) approach.

There are no reference genomes available for rusa deer, so we mapped reads to the hog deer genome (Axis porcinus, GCA_003798545.1) using the burrows-wheeler aligner (BWA) v0.7.17 software package (Li and Durbin 2009). The hog deer genome was selected based on the phylogeny and topology of the intrafamilial relationships among deer (Pitra et al. 2004). We used bwa align (‘aln’ command) (Li and Durbin 2009) to align the reads to the genome and trimmed the first five nucleotides (‘-B 5’ command) based on the fastqc results. We used the populations program (Catchen et al. 2013) in stacks to call single nucleotide polymorphisms (SNPs). We conducted several filtering steps prior to the initial analysis.

We performed initial quality filtering using the population program in stacks, using a minimum loci log likelihood value of −10 (‘--lnl_lim -10’ command), and we retained one random SNP per locus (‘--write_random_snp’ command). We applied additional filtering parameters with VCFTools v0.1.16 (Danecek et al. 2011) to the variant file produced by the population program in stacks, with minimum genotype quality score of 15 (‘--minGQ 15’ command), a maximum locus depth of 50 (‘--maxDP 50’ command) and minimum locus depth of 2 (‘--minDP 2’ command). Following these filtering steps, there were 62 461 SNPs remaining in the dataset. Using VCFTools, we identified and removed two individuals that had a high level of missingness (over 70%; ‘--missing-indv’ command). Following this step, 172 individuals were retained for all downstream analyses (Supplementary File 1), including 88 females (74 adults and 14 juveniles) and 84 males (64 adults and 20 juveniles). We chose a minimum minor allele frequency (‘--maf’ command) of 1% and maximum missing (‘--max-missingness’ command) of 0.3. This resulted in a data set including 12 949 SNPs.

Analysis

To evaluate population structure across the study site, we conducted an initial admixture analysis with all samples across the study site. admixture v. 1.3.0 (Alexander et al. 2009) produced maximum likelihood estimations of the individual ancestries and representative K-values (a key user-defined parameter of hypothesised number of populations). For all admixture analyses throughout this study, a random seed generated from the current time (‘-s time’), and K-values 1–10 were generated for each admixture analysis. Using admixture’s cross-validation procedure (‘--cv flag’ command), we selected appropriate admixture plots based on examining the K-values’ CV error and selecting the K that corresponded with the lowest CV error. Plots were visualised using the python script distruct v2.3 (Chhatre 2018). In addition to this, we also performed Principal Component Analysis (PCA) using snprelate v1.24.0 (‘snpgdsPCA’) (Zheng et al. 2012) in r v3.6.3 (R Core Team 2020) and visualised the results using ggplot2 v3.4 (Wickham 2016).

We removed outlier individuals identified in the initial admixture and PCA analyses before we conducted downstream analyses, which left a total of 165 individuals. Taking into consideration the current rusa deer distribution in New South Wales (Crittle and Millynn 2020) and number of samples per sampling site and location, we assigned individuals to geographic groups using hierarchical clustering with a maximum distance of 5 km (R-package stats v 4.2.2) (Murtagh and Legendre 2014). No individual was assigned to more than one of these 14 groups (Fig. 1b, Supplementary File 1). We estimated relatedness values for all pairs of individuals using the dartR v1.9.6 (Gruber et al. 2018) function ‘gl.grm’.

To answer Question 1 (Table 1, Do deer in the Illawarra region represent a single, panmictic population?), two predictions were tested. To test Prediction 1A: Large Distance, we investigated whether large geographic distances across the Illawarra region result in major genetic subdivision, and to test Prediction 1B: Isolation By Distance (IBD), we investigated whether there is a relationship between geographical and genetic distances (Table 1). All individuals passing quality filters were included in these analyses, excluding outliers.

To test Prediction 1A: Large Distance, we conducted admixture analysis (Fig. 2A) and produced a PCA (as above; Supplementary File 2, Fig. S1). A pairwise matrix of fixation indices (F_ST) was created using the ‘gl.fst.pop’ function in dartR to estimate levels of genetic differentiation across our study site. For this analysis, individuals were assigned to the same 14 groups described above (Supplementary File 1). For Prediction 1B: IBD, we used a Mantel test performed in dartR function ‘gl.ibd’ with 999 permutations. We visualised patterns of IBD by plotting transformed pairwise differentiation (F_ST/1–F_ST; (Rousset 1997) against log-transformed geographic distance), again using dartR.

Fig. 2.

Location of rusa deer samples used to answer Question 1A and Questions 2A–C (Table 1). The admixture plots below the maps indicate population genetic structure and contain genetic data from individuals included on the corresponding map, with individuals ordered in the plots from north to south (left to right). admixture analysis determines the optimal number of ancestral genetic groups that explain patterns in a genetic data set, with each ancestral genetic group being represented by a separate colour. To visualise this information, each vertical bar within an admixture plot represents an individual, with the proportion of colours within the bar representing the likelihood that that individual was descended from each genetic group; for example, a vertical bar in the admixture plot that is all blue indicates an individual with no mixed ancestry, while a vertical bar that contains blue, red, and purple indicates an individual of mixed ancestry. Ordering individuals based on their sampling location allows us to see broad trends within and across sampling groups. For all panels, the colours and shapes on maps correspond to those indicated below the associated admixture plot, as well as the PCA plots presented in Supplementary File 2, Figs S2 and S3. Map 1A: deer sampled in the north (Group 1A1), central (Group 1A2), and south of the Illawarra region (Group 2A3) for Prediction 1A: Large Distance. Map 2A: Deer sampled above (Group 2A1) and below (Group 2A2) the escarpment (shown in green) used to test Prediction 2A: Escarpment. Map 2B: Deer sampled north of the Unanderra-Moss Vale railway (Group 2B1), between Lake Illawarra and highways/railways (Group 2B2), and south of the Unanderra-Moss Vale railway (Group 2B3) used to test Prediction 2B: Railways and Highways. Map 2C: deer sampled north of the narrow corridor (Group 2C1), within the narrow corridor (Group 2C2) and south of the narrow corridor (Group 2C3) used to test Prediction 2C: Narrow Corridor.

To answer Question 2 (Table 1, Do natural and human-made landscape features influence dispersal?), we tested three predictions including different subgroups of individuals (Table 1, Fig. 2.2A–C). To test Prediction 2A: Escarpment, we investigated whether the steep slopes of the escarpment above the city of Wollongong act as a barrier to dispersal. The individuals included in this analysis were those found above (Group 2A1, Fig. 2.2A) or below (Group 1A2, Fig. 2.1A) the Illawarra escarpment (group assignments described in Supplementary File 1).

To test Prediction 2B: Railways and Highways, we investigated whether the railways and highways act as barriers to dispersal. We assigned individuals to three groups (group assignments described in Supplementary File 1) based on their proximity to the railways and highways near the city of Wollongong and Lake Illawarra (Fig. 2.2B). Groups 2B1 and 2B3 were sampled on each side of the Unanderra–Moss Vale railway line (north and south, respectively), and west of the main highways and railways. The Group 2B2 was sampled east of the main highways and railways, in a densely urbanised area from Lake Illawarra to the Wollongong CBD (Fig. 2.2B). These railways and highways were partially fenced, but fence material and height varied.

To test Prediction 2C: Narrow Corridor, we investigated whether the narrow corridor between the steepest part of the escarpment (north of the city of Wollongong) and the Tasman Sea reduced dispersal of rusa deer (Table 1, Fig. 2.2C). The individuals included in this prediction were those sampled north of the corridor close to Royal National Park (Group 2C1, Fig. 2.2C), those sampled within the corridor (Group 2C2, Fig. 2.2C) and those sampled south of the corridor close to the city of Wollongong (Group 2C3, Fig. 2.2C; group assignments described in Supplementary File 1). Additionally, we calculated genetic diversity of the north and south groups, using a random subsample of 34 individuals from each group and estimating expected heterozygosity, number of private variants (only found in one site) and percentage of polymorphic loci within each site using the populations program in stacks.

To answer Question 3 (Table 1, Is there evidence of sex-biased dispersal?), we selected three groups of individuals that were determined through admixture analysis to be genetically distinct from one another (these groups were 2B3, 2C1, and 2C3, described in Supplementary File 1), which was comprised of 60 adult females and 41 adult males. We used the R package hierfstat (‘sexbias.test’) to calculate the assignment index (AI) statistic (mean (mAIC), variance (vAIC), F_ST, and F_IS were assessed, with the settings alternative = ‘two.sided’, nperm = 100) (Goudet et al. 2002).

Ethical approval

This study did not require approval from an Animal Ethics Committee because it used tissue samples collected from animals killed in pest management programs.

Initial analyses

When all 172 individuals were included, we identified seven individuals in the PCA and admixture analyses that were genetically dissimilar to all remaining samples (n = 165), and to one another (Fig. 3; Supplementary File 2, Fig. S1). Outliers were geographically dispersed (Fig. 3) and included four adult males, one juvenile male, one adult female and one juvenile female.

Fig. 3.

The genetic outliers identified during our initial analysis of rusa deer sampled in the Illawarra region. (a) Filled circles indicate sampling locations of all deer (n = 172) included in this study (outliers are indicated by red, all other samples by black). (b) PCA analysis of genetic differences based on single nucleotide polymorphism data shows outliers in red (n = 7) and all other samples in black (n = 165).

Question 1: Do deer in the Illawarra region represent a single, panmictic population?

Both the admixture and PCA analysis showed evidence of genetic differentiation between individuals in the northern group (Group 1A1) and the middle group (Group 1A2; Fig. 2.1A admixture plot; Supplementary File 2, Figs S2 and S3). The two individuals sampled farthest to the south of our study area (both adult males; Group 1A3) were genetically similar to Group 1A2 (Fig. 2.1A admixture plot; Supplementary File 2, Figs S2 and S3), despite being separated from this group by >25 km. Overall, pairwise relatedness values across our study were low (highest value was 0.36; only 12 of 13 530 pairs had a value that exceeded 0.25). The Mantel test revealed a weak but statistically significant positive relationship between genetic and geographic distance (Mantel statistic r = 0.404, P = 0.01; Fig. 4). Pairwise F_ST values ranged between 0.00 and 0.15, and most were significant (pairwise F_ST and associated P-values are reported in Supplementary File 3).

Fig. 4.

Isolation by distance (IBD) analysis of 165 rusa deer included in this study (excluding seven genetic outliers). Individuals were grouped based on geographical proximity (n = 14 groups, Supplementary File 1). The result of the Mantel test showed a significant positive relationship between geographic and genetic distance (r = 0.404, P = 0.01). Each point in this graph represents a pairwise comparison between two geographic groups. Within a population displaying IBD, we would expect two groups that have a small geographic distance between them to also have a low genetic distance between them (and vice versa). The genetic distance measurement used in this analysis is based on fixation index (F_ST), which is a pairwise comparison that captures the total amount of genetic variance between groups as a proportion of the total genetic variance across both groups.

Question 2: Do natural or human-made landscape features influence dispersal?

The prediction that the Illawarra escarpment acts as a barrier to gene flow was not supported; samples from above (Group 2A1) and below (Group 2A2) the escarpment were genetically indistinct (Fig. 2.2A admixture plot; Supplementary File 2, Figs S2 and S3). However, the Unanderra–Moss Vale railway line was associated with reduced gene flow between the group of samples north of the railway and the group south of the railway (Fig. 2.2B admixture plot; Supplementary File 2, Figs S2 and S3). In contrast, the 10 samples collected close to Lake Illawarra (Group 2B2) clustered with the group south of the Unanderra–Moss Vale railway (Group 2B3), despite being separated by two highways and a busy railway (Fig. 2.2B admixture plot; Supplementary File 2, Figs S2 and S3).

The analysis of samples along the narrow corridor between the escarpment and the Tasman Sea showed evidence of genetic differentiation. We found two distinct genetic clusters in both the PCA and admixture analyses (Fig. 2.2C admixture plot; Supplementary File 2, Figs S2 and S3) for the north (Group 2C1) and south groups (Group 2C3). The 11 samples collected at the narrowest part of the corridor (Group 2C2) contained mixed ancestry from both the northern and southern groups (Fig. 2.2C admixture plot; Supplementary File 2, Figs S2 and S3). There was slightly higher genetic diversity in the northern group (expected heterozygosity = 0.21 and 65.7% polymorphic loci) compared with the southern group (expected heterozygosity = 0.19 and 59.2% polymorphic loci). In addition, the northern group had >2.5 times more private variants than the southern group (1281 and 449, respectively).

Question 3: Is there evidence of sex-biased dispersal?

We found no evidence of sex-biased dispersal: mAIc (t = 0.00, P = 1.00), vAIc (t = 0.63, P = 0.99), F_IS (t = 0.00, P = 0.57) and F_ST (t = 0.00, P = 1.00).

Question 1: Do deer in the Illawarra region represent a single, panmictic population?

Our admixture analyses indicated that there was population genetic structure across the study area, but F_ST values indicated that differentiation between most sampling sites was low. After removing the seven individual outliers, genetic distance between sampling sites was weakly correlated to the geographical distance between them. This accords with what is known of the home range size of rusa deer in Australia (i.e. <8 km²; Moriarty 2004b; Amos et al. 2023), so our study area far exceeds home range size (approximately 85 km from two furthest sample points). Although the pattern of IBD we observed was significant, the large variance in the strength of this relationship indicates that other factors are also important for population structure in this species. A more comprehensive study of landscape features may better explain these patterns. The analyses discussed below indicate that there is finer-scale population structure.

Question 2: Do natural or human-made landscape features influence dispersal?

Surprisingly, we found no evidence that the escarpment was a barrier to gene flow. Steep environmental gradients commonly influence gene flow, and animal movement patterns often follow prominent physical features of the landscape (Pérez-Espona et al. 2008). This includes escarpments in Australia that have been found to impede gene flow in other invasive mammal species such as feral pigs (Loving 2022). However, in New Zealand, rusa deer live on steep forested slopes (Allen et al. 2015). Our results suggest that deer are likely to continue to move across this landscape feature in the southern Illawarra region. We note that most of our samples were collected below the escarpment, and we were unable to collect any samples above the escarpment at its northern extent, where the slope is steepest (Young 1979). Deliberate (illegal) translocation could also explain this result. Additional sampling of deer above the escarpment, especially in the north, may clarify whether this landscape feature affects rusa deer population structure.

Despite close geographic proximity, there was genetic clustering of deer on either side of the Unanderra–Moss Vale railway. This is interesting, because this railway is an infrequently used freight line, and is only partially fenced. Railways and roads may act as significant barriers or facilitators to gene flow in deer populations (Coulon et al. 2006; Long et al. 2010; Frantz et al. 2012; Robinson et al. 2012), but some studies have noted that these anthropological features need to be interpreted carefully because they are often associated with natural landscape features that could affect movement (Pérez-Espona et al. 2008). The genetic differentiation that we observed could be due to deer not wanting to cross the railway line due to the risk of collision and/or because they are scared by the loud noises made by trains (Drolet et al. 2016; Bhardwaj et al. 2022), but another railway in our study did not appear to impede gene flow (see below). The genetic clustering on either side of the Unanderra–Moss Vale railway might also represent larger-scale differences in landscape features affecting deer movement. This railway is associated with a major change in the landscape. North of the railway there is mostly high-density residential housing that encroaches the lower forested slopes of the escarpment. South of the railway, the landscape becomes flatter and more open with extensive industrial areas and pasture paddocks. Rusa deer generally prefer dense vegetation cover as refuge during the day (Moriarty 2004b; Spaggiari and de Garine-Wichatitsky 2006); therefore, the region south of the railway might be a lower quality habitat that deer avoid.

We found no evidence of population structure in samples separated by major highways and a busy railway near Lake Illawarra. Major highways have a varying impact on gene flow in deer, acting as barriers in some studies but not others (Long et al. 2010; Meisingset et al. 2013; Peterson et al. 2017). In our study, there is irregular fencing along these anthropological features, and this may be an additional factor contributing to the permeability of these landscape features (Olsson 2007; Hepenstrick et al. 2012; Huijser and Begley 2022). Reasons for deer moving across human-made structures associated with high mortality risks, such as roads and railways, include searching for food and shelter, and seeking mating opportunities (Ando 2003). It is likely that the rusa deer moving across these structures have learned behaviours (e.g. flight responses) to reduce the risks of collisions with trains and road vehicles, as seen in other mammals (Borda-de-Água et al. 2017; Bhardwaj et al. 2022). Our study investigated the impacts of only three human-made landscape features. Substantially increasing the number of rusa deer samples, and their geographic spread, would help better understand the influence of human-made landscape features on gene flow.

We examined a subset of samples near the original introduction site in Royal National Park and contiguous sites to the south, including a geographically intermediate site along the narrow corridor of land between the escarpment and the coast (Table 1, Groups 2C1–2C3). The deer sampled across these areas exhibited population genetic structure, and it appears that individuals in the narrow corridor are genetically intermediate to those to the north and those to the south. The genetic diversity is higher in the northern sites, near the original introduction site in Royal National Park (Bentley 1998; Moriarty 2004b; Webley et al. 2004). There is less genetic diversity further south (Group 2C3), which is expected if deer in this area are more recently established following dispersal from sites to the north. Our genetic diversity data suggest that more southerly sites resemble patterns expected of range-edge populations (Brussard 1984; Phillips et al. 2008; Langin et al. 2017), supporting the prediction that the rusa deer population is expanding south and west from Royal National Park. Despite lower genetic diversity in the south, we still identified 449 private genetic variants in these samples. Given the results of Hill et al. (2023), these variants are unlikely to be due to hybridisation with other species but could result from admixture following immigration from different populations of rusa deer, from introgression following hybridisation events, or could represent rare alleles that exist in the north but were not sampled.

Question 3: Is there evidence of sex-biased dispersal?

Male-biased dispersal is commonly observed in deer (Spaggiari and de Garine-Wichatitsky 2006; Pérez-Espona et al. 2010), but we did not find any genetic evidence that rusa deer dispersal in our study area was sex-biased. The method we used compares the variance of how well individuals of each sex genetically assign to the group from which they were sampled. In studies like ours, where population differentiation is low, the accuracy of assignment may be lower than in populations with high levels of genetic differentiation. The two adult male deer sampled in the far south of our study area (seen in Fig. 2, Group 1A3) were genetically similar to individuals sampled ~25 km to the north. These individuals could have been dispersers and/or illegal translocations. More samples from this area are needed to determine if they were dispersers.

Management implications

Although we found low overall levels of genetic differentiation in rusa deer sampled across the Illawarra region, we identified three spatially segregated genetic clusters that could inform management (Supplementary File 2, Fig. S4). The clustering indicates reduced mixing of deer in each unit with deer from outside that unit, as identified by genetic structure associated with populations to the north versus the south of the narrow corridor between the escarpment and the sea (Prediction 2C) and genetic structure across one of the corridors of major highway and rail (Prediction 2B). Vehicle-based shooting in peri-urban landscapes is limited by the need to minimise actual and perceived risks to people (e.g. shooting does not occur on weekends and public holidays), and by the cost. Deer killed during management operations can be replaced by in situ births and/or immigration. In the management units identified in our study, immigration is limited (i.e. low genetic mixing with outside populations), and therefore only local recruitment can compensate for population reduction by shooting. Thus, concentrating the culling within a particular management unit is expected to have a greater long-term impact on local deer abundance, compared with spreading the culling effort across the whole region.

Interestingly, we found that approximately 4% of the samples included in our study were genetic outliers and are likely to be hybrids with another deer species (Hill et al. 2023). It would be useful to identify their origin so that future incursions can be limited (e.g. by culling in the population of origin). Management of an invasive population is more cost-effective when implemented at the start of the invasion, during the well-characterised ‘lag period’ that precedes exponential population growth (Sakai et al. 2001), rather than when a population is well-established.

Additionally, it appears that there is ongoing gene flow from outside our study area, because 449 genetic variants we found in samples taken near the city of Wollongong (Fig. 2.2C) were not found in any of the 36 individuals sampled close to Royal National Park (Group 2C1). If the rusa deer sources outside the Illawarra region could be identified, managers might consider culling in the population of origin or take legal actions to prevent further introductions.