Springtail phylogeography highlights biosecurity risks of repeated invasions and intraregional transfers among remote islands.

Human-mediated transport of species outside their natural range is a rapidly growing threat to biodiversity, particularly for island ecosystems that have evolved in isolation. The genetic structure underpinning island populations will largely determine their response to increased transport and thus help to inform biosecurity management. However, this information is severely lacking for some groups, such as the soil fauna. We therefore analysed the phylogeographic structure of an indigenous and an invasive springtail species (Collembola: Poduromorpha), each distributed across multiple remote sub-Antarctic islands, where human activity is currently intensifying. For both species, we generated a genome-wide SNP data set and additionally analysed all available COI barcodes. Genetic differentiation in the indigenous springtail Tullbergia bisetosa is substantial among (and, to a lesser degree, within) islands, reflecting low dispersal and historic population fragmentation, while COI patterns reveal ancestral signatures of postglacial recolonization. This pronounced geographic structure demonstrates the key role of allopatric divergence in shaping the region's diversity and highlights the vulnerability of indigenous populations to genetic homogenization via human transport. For the invasive species Hypogastrura viatica, nuclear genetic structure is much less apparent, particularly for islands linked by regular shipping, while diverged COI haplotypes indicate multiple independent introductions to each island. Thus, human transport has likely facilitated this species' persistence since its initial colonization, through the ongoing introduction and inter-island spread of genetic variation. These findings highlight the different evolutionary consequences of human transport for indigenous and invasive soil species. Crucially, both outcomes demonstrate the need for improved intraregional biosecurity among remote island systems, where the policy focus to date has been on external introductions.

INTRODUCTION

Human transport is eroding the isolation of island ecosystems, with remote islands especially vulnerable to invasion (Moser et al., 2018). Many island biotas now comprise high proportions of introduced species (e.g. Cicconardi et al., 2017; Gaston, Jones, Hänel, & Chown, 2003), which have a range of often substantial impacts (e.g. Clavero, Brotons, Pons, & Sol, 2009; Reaser et al., 2007). Human movements also threaten to merge previously isolated populations of indigenous species, leading to genetic homogenization and subsequent reductions in fitness (Crispo, Moore, Lee‐Yaw, Gray, & Haller, 2011; Olden, 2006), or even reverse speciation (Seehausen, 2006). Ultimately, the response of island communities to increased human transport will be shaped by dispersal, colonization and evolutionary processes unique to both the indigenous and introduced biota (Gillespie, Claridge, & Roderick, 2008). While invasion patterns and dynamics have been well studied for a growing range of island taxa (e.g. Lloret et al., 2005; Moser et al., 2018), those for the soil biota remain poorly known (Coyle et al., 2017; Novo et al., 2015; Ricciardi et al., 2017), despite the prevalence of introduced soil‐dwelling species on islands worldwide (e.g. Cicconardi et al., 2017; Frenot et al., 2005). Changes in soil community composition have far‐reaching effects on ecosystem functioning (e.g. Coyle et al., 2017; Geisen, Wall, & Putten, 2019); therefore, understanding how human transport can alter the structure and evolutionary processes of soil populations is a pressing priority for predicting broader impacts on island systems.

The sub‐Antarctic islands (located in the Southern Ocean) provide an exemplary setting to do so. They are among the world's most remote islands and have substantial conservation value (Chown, Rodrigues, Gremmen, & Gaston, 2001; UNESCO, 2019), with highly specialized soil ecosystems considered vulnerable to invasion (Convey & Lebouvier, 2009; Frenot et al., 2005). Transport across the Southern Ocean is burgeoning as human activities intensify in the Antarctic, increasing the risk of nonindigenous species introductions into the region, as well as intraregional movement of populations between islands or biogeographic zones (Hughes & Convey, 2010; Lee & Chown, 2011). Two key potential consequences have been identified for intraregional transport: the genetic homogenization and biogeographic breakdown of distinct indigenous populations, and the further spread of invasive species established in the region. The genetic structure underlying indigenous and invasive populations will largely determine how this plays out (Crispo et al., 2011; Olden, 2006; Petit, 2004). Thus, an improved understanding of the phylogeography of sub‐Antarctic populations can begin to inform generalities for isolated soil ecosystems under growing transport regimes.

Present understanding of indigenous phylogeographic structure across the terrestrial sub‐Antarctic is based on a limited number of molecular studies (see Moon, Chown, & Fraser, 2017). Instances where few genotypes predominate across the region are typically thought to represent postglacial recolonization of islands following population bottlenecks during glacial maxima (e.g. Fraser, Nikula, Spencer, & Waters, 2009; Wouw, Dijk, & Huiskes, 2008), presumably aided by rafting on marine currents (Fraser, Nikula, Ruzzante, & Waters, 2012; McGaughran, Stevens, Hogg, & Carapelli, 2011). However, many sub‐Antarctic soil populations are now considered genetically isolated with pronounced differentiation between islands (e.g. van Vuuren, Lee, Convey, & Chown, 2018), often indicating allopatric speciation processes (e.g. Allegrucci, Carchini, Convey, & Sbordoni, 2012; Stevens, Greenslade, Hogg, & Sunnucks, 2006). At a finer scale, intra‐island population differentiation has been detected for several species, aligning with historical barriers to gene flow such as glaciers and topography (e.g. Chau et al., 2019; Mortimer, Vuuren, Meiklejohn, & Chown, 2012). Although these studies are limited in geographic and molecular scope (Halanych & Mahon, 2018; Moon et al., 2017), they provide preliminary evidence that the indigenous terrestrial faunas comprise genetically distinct populations vulnerable to homogenization by humans, as is the case for islands elsewhere (e.g. Gillespie et al., 2008; Juan, Crespo, Cowan, Lexer, & Fay, 2004; Vandergast, Gillespie, & Roderick, 2004).

The risk of human activity not only introducing new alien species, but also spreading or reintroducing existing invasives, is also of concern for island systems such as the sub‐Antarctic (Houghton, Terauds, Merritt, Driessen, & Shaw, 2019; Lee & Chown, 2011). Multiple introductions from distinct sources can enhance genetic variation within invasive populations, potentially promoting their establishment and persistence (Dlugosch, Anderson, Braasch, Cang, & Gillette, 2015; Lavergne & Molofsky, 2007; Wagner, Ochocki, Crawford, Compagnoni, & Miller, 2017). The few studies that have employed high‐resolution molecular markers for sub‐Antarctic invasive species have indeed revealed multiple colonizations of a single island (Ouisse, 2016; Piertney et al., 2015), though not always (e.g. Hardouin et al., 2010). The absence of genetic studies on more widespread sub‐Antarctic invasives makes it difficult to determine whether initial introductions have been followed by natural dispersal to other islands, or whether each island represents an additional introduction. Data on the genetic structure of invasive populations are thus essential to better understand invasion dynamics across the region (Estoup & Guillemaud, 2010; Le Roux & Wieczorek, 2009) and to design adequate biosecurity measures for remote island ecosystems (e.g. Abdelkrim, Pascal, Calmet, & Samadi, 2005; Dudaniec, Gardner, Donnellan, & Kleindorfer, 2008).

We therefore conducted the first genome‐wide SNP‐based analysis of phylogeographic structure for an indigenous and an invasive soil species distributed across multiple sub‐Antarctic islands. Given that previous studies in the region primarily examine COI barcodes, we also incorporated COI data to compare mitochondrial outcomes with nuclear genomic patterns. We focussed on Collembola (springtails), which represent a substantial component of sub‐Antarctic soil communities and include a large and concerning number of alien invasive species on these and other islands (Baird, Janion‐Scheepers, Stevens, Leihy, & Chown, 2019; Cicconardi et al., 2017; Greenslade & Convey, 2012).

MATERIALS AND METHODS

Sampling

We sampled two springtail species from the Order Poduromorpha: the indigenous Tullbergia bisetosa Börner, 1902 (Tullbergiidae), and the highly invasive Hypogastrura viatica (Tullberg, 1872) (Hypogastruridae). Both of these species occur on at least four sub‐Antarctic islands and have been confirmed not to comprise cryptic species (T. bisetosa: Louis Deharveng, pers. comm.; H. viatica: Greenslade & Convey, 2012). To explore small‐ and large‐scale genetic structure, multiple sites were sampled from the sub‐Antarctic islands of Macquarie, Heard, Kerguelen, Possession, Marion and South Georgia (see Figure 1) between the years 2004 and 2017 (see Tables S1.1 and S1.2). Given the logistical challenges accessing these islands, the spatial spread of within‐island sampling varied and sample sizes were occasionally limited. Nonetheless, at least two sites were sampled from each island, spaced a minimum of 100 m and a maximum of 21 km apart, while the islands themselves are separated by 1,000–8,000 km of open ocean (see Tables S2.1 and S2.2 for distance matrices). Animals were sampled either by direct aspiration (from an area <30 cm2) or by Tullgren funnel extractions of soil or vegetation cores (<30 cm2 width), and were immediately preserved in ethanol.

Figure 1

Principal coordinate analysis (PCoA) and maximum likelihood (ML) tree outcomes for the indigenous springtail Tullbergia bisetosa (a,b) and the invasive springtail Hypogastrura viatica (c,d), based on 5,680 and 7,275 SNPs, respectively. A map of Antarctica showing the sub‐Antarctic islands from which the two species were sampled is depicted in the centre. Individuals in both PCoA plots and ML trees are colour‐coded by their geographic origin; each island has a main colour, and sites within islands are shades of that colour (see legends for each species). Sites within an island discerned by a letter suffix (1a, b, c..) are separated by <5 km; sites within an island discerned by a number suffix (1, 2…) are separated by >10 km. PCoA plots show the per cent variation explained by each axis. Branch lengths on the ML trees represent the average number of substitutions per site (see scale bars), and nodes show bootstrap support (nodes with <90% support are not shown). Numbers on the ML tree for T. bisetosa show the groupings inferred by fixed allelic differences to be strongly isolated from gene flow. APF, Antarctic Polar Front; ACC, Antarctic Circumpolar Current

DNA extraction and SNP genotyping

Specimens were sent to Diversity Arrays Technology (DArT) Pty. Ltd. (Canberra, Australia) for genomic DNA extraction and sequencing using DArT‐SeqTM technology. DArT‐SeqTM is a reduced representation next‐generation sequencing method, which employs a combination of restriction enzymes to digest genomic DNA, followed by the ligation of adaptors to both restriction enzyme overhangs, to facilitate sample identification and Illumina sequencing (see methodology in Kilian et al., 2012). Raw Illumina data are processed using a proprietary DArT pipeline to provide a quality‐filtered SNP data set, which we filtered even further using the “dartR” package v.1.0.5 (Gruber, Unmack, Berry, & Georges, 2018) in R v.3.4.3 (R Core Team, 2017). Appendix S1 provides more detail on DArT‐Seq SNP library preparation, quality control and filtering measures.

SNP analyses

For the analyses where statistical significance was required for decision‐making (tests for selection, fixed allelic differences and deviations from the Hardy–Weinberg equilibrium), our results are adjusted for multiple tests using the most common approach in the genomic literature (e.g. Bonferroni, false discovery rate). For all other analyses, we only calculated uncorrected p‐values where necessary and generally refrain from heavy emphases on statistical significance, acknowledging growing concern over the appropriate adjustment for multiple testing in genomics (e.g. Sethuraman et al., 2019; White, Ende, & Nichols, 2019) and the over‐interpretation of p‐values in general (Amrhein, Greenland, & McShane, 2019).

The final set of SNPs for each species was tested for evidence of selection in BayeScan 2.1 (Foll & Gaggiotti, 2008), using the F ST outlier approach and a false discovery rate of α = 0.01. Each SNP locus was also tested for deviations from Hardy–Weinberg equilibrium (HWE) in each population, using exact tests with a Bonferroni correction, conducted in dartR.

To examine the genetic diversity within each population, observed and expected heterozygosity and inbreeding coefficients were calculated for each species at each site using GenoDive v.2.0 (Meirmans & Van Tienderen, 2004). The 95% confidence intervals for these indices were generated with 5,000 bootstraps to compare differences in diversity between sites and between islands, within each species.

Population differentiation was initially explored within each species by generating global and pairwise F ST estimates (as per Weir & Cockerham, 1984) in GenoDive, using 5,000 permutations to assess significance. Given the hierarchical spacing of our populations, we employed hierarchical AMOVA (Excoffier, Smouse, & Quattro, 1992) to assess the subdivision of genetic variation at each level of sampling (i.e. among islands and among sites within islands). AMOVA was carried out in GenoDive.

Principal coordinate analysis (PCoA) was used to visualize the clustering of genetically similar groups within each species and was carried out using the “gl.pcoa” command in dartR. Additionally, 3D PCAs (which enabled visualization of a third axis) were created using the “pca3d” package (Weiner, 2017) in R. We also used a maximum likelihood approach to infer unrooted phylogenetic trees, as implemented in IQ‐TREE (Nguyen, Schmidt, Haeseler, & Minh, 2014). Node support was inferred based on 10,000 bootstraps of the data, and output files were visualized with FigTree v.1.4.3 (Rambaut, 2009).

Population structure was further explored using the Bayesian clustering algorithm implemented in fastSTRUCTURE (Raj, Stephens, & Pritchard, 2014). Ten replicates were run for each K value from K = 1 to K = 15, using the simple prior model. The number of populations that best fit the data was determined using the “chooseK.py” function, and these population assignments were visualized using Distruct v.2.3 (Chhatre, 2019).

Finally, each data set was tested for fixed allelic differences (when two populations share no alleles at a given SNP), which can indicate a complete lack of contemporary gene flow (Slatkin, 1987). This approach is less affected by departures from HWE and small sample sizes than the more conventional “Structure” family of clustering software (Georges et al., 2018). Fixed difference analysis was conducted in dartR using the “gl.collapse.recursive” function. For each pairwise population comparison, we used 5,000 simulations to estimate the expected rate of false positives (which arise due to the finite number of individuals in each population), and any two populations deemed to have an insignificant number of fixed differences (based on this false‐positive rate, set at α = 0.01) were amalgamated. This was performed iteratively until a final set of “pseudo‐populations” remained, each discernable by having a significant number of fixed differences from one another.

Mitochondrial DNA analyses

All available COI sequences for T. bisetosa and H. viatica were collated from GenBank (http://www.ncbi.nlm.nih.gov/genbank) and the Barcode of Life Data System (“BOLD”: http://www.barcodinglife.org), by searching for public records within their respective Barcode Index Number (BIN) on BOLD. BINs employ a Refined Single Linkage algorithm to cluster sequences that are highly likely to represent the same species (Ratnasingham & Hebert, 2013), and we conservatively excluded any BIN members that had been morphologically identified as a species other than T. bisetosa or H. viatica. All remaining sequences ≥500 bp within each BIN were downloaded (see Tables S3.1 and S3.2 for accession numbers), manually checked for stop codons or indels that could indicate pseudogenes, and aligned using MUSCLE in MEGA v.7 (Kumar, Stecher, & Tamura, 2016). Haplotype networks were created using the median‐joining option (Bandelt, Forster, & Röhl, 1999) in PopART (popart.otago.ac.nz). Further exploration of the data, such as testing for population expansion, was not attempted given the opportunistic (and therefore geographically nonuniform) sampling of the sequences.

RESULTS

SNP data sets

The DArTsoft14 pipeline delivered an initial 39,847 and 37,803 SNPs for Tullbergia bisetosa and Hypogastrura viatica, respectively. After quality filtering (see Table S4 for the data set size at each filtering step), the final data set for T. bisetosa comprised 5,680 SNPs (see Baird, Moon, Janion‐Scheepers, & Chown, 2019a), and for H. viatica comprised 7,275 SNPs (see Baird, Moon, Janion‐Scheepers, & Chown, 2019b). Less than 5% of SNPs were found to significantly violate HWE in any single population for either species, and the vast majority of these HWE deviations represented heterozygote deficits (>95% for both species). The number of outlier SNPs putatively under selection was also minimal, representing 2% of the data set for both species (123 SNPs for T. bisetosa and 158 SNPs for H. viatica). Thus, outlier SNPs putatively under selection or deviating from HWE were not removed from the main data sets. However, we prepared a conservative data set for each species with all of these SNPs removed (see Table S4) on which we performed a subset of the analyses (genetic diversity statistics, AMOVA and PCoA), to verify that this did not affect our results. These outcomes are almost identical to those from the main data sets (see Figures S1, S2 and Tables S5–S8); we henceforth refer only to results from the full data sets.

Genomic diversity and structure

Tullbergia bisetosa

Populations of T. bisetosa had low levels of observed and expected heterozygosity (global H O = 0.069, global H E  = 0.099; see Table 1 for local population values), in part reflecting the large proportion of SNPs that were fixed for a single allele within each population. Heterozygote deficiencies were evident for all populations by high and positive coefficients of inbreeding (global F IS = 0.304; Table 1), indicating that nearly one third of heterozygosity was lost. Comparisons between populations of T. bisetosa revealed that South Georgia harboured significantly higher levels of genetic diversity (based on observed and expected heterozygosity) than Macquarie, Heard or Marion islands (p < .05 for both H O and H E; Table 1).

Table 1

Sample size and genetic diversity statistics (with 95% CIs in parentheses) for sub‐Antarctic populations of the indigenous springtail Tullbergia bisetosa and the invasive springtail Hypogastrura viatica, based on 5,680 and 7,275 genome‐wide SNPs, respectively

	Site	N	H O	H E	F IS
T. bisetosa
Macquarie	MQ‐1a	17	0.046 (0.043, 0.049)	0.067 (0.063, 0.071)	0.316 (0.295, 0.336)
	MQ‐1b	10	0.051 (0.047, 0.054)	0.069 (0.065, 0.073)	0.263 (0.238, 0.288)
	MQ‐2	19	0.059 (0.055, 0.062)	0.069 (0.066, 0.073)	0.156 (0.137, 0.175)
Heard	HD‐1a	4	0.058 (0.052, 0.059)	0.075 (0.067, 0.076)	0.225 (0.185, 0.264)
	HD‐1b	13	0.052 (0.048, 0.054)	0.069 (0.064, 0.072)	0.250 (0.223, 0.278)
Marion	MR‐1a	16	0.046 (0.043, 0.049)	0.071 (0.067, 0.075)	0.353 (0.332, 0.375)
	MR‐1b	23	0.069 (0.065, 0.072)	0.103 (0.098, 0.107)	0.333 (0.316, 0.350)
	MR‐1c	17	0.056 (0.053, 0.060)	0.093 (0.088, 0.097)	0.392 (0.373, 0.410)
Sth Georgia	SG‐1a	20	0.106 (0.102, 0.110)	0.159 (0.154, 0.164)	0.332 (0.318, 0.346)
	SG‐1b	19	0.111 (0.107, 0.115)	0.153 (0.148, 0.158)	0.274 (0.259, 0.289)
	SG‐2	20	0.110 (0.105, 0.113)	0.159 (0.153, 0.164)	0.312 (0.297, 0.327)
	Global	178	0.069 (0.068, 0.071)	0.099 (0.098, 0.101)	0.304 (0.296, 0.312)
H. viatica
Macquarie	MQ‐H‐1a	16	0.138 (0.135, 0.142)	0.231 (0.226, 0.235)	0.399 (0.389, 0.410)
	MQ‐H‐1b	22	0.151 (0.148, 0.154)	0.229 (0.225, 0.233)	0.341 (0.332, 0.350)
	MQ‐H‐1c	14	0.133 (0.129, 0.136)	0.229 (0.225, 0.234)	0.421 (0.411, 0.431)
	MQ‐H‐2	21	0.144 (0.141, 0.148)	0.217 (0.213, 0.221)	0.334 (0.325, 0.344)
Kerguelen	KR‐1a	19	0.150 (0.147, 0.154)	0.231 (0.227, 0.236)	0.351 (0.342, 0.360)
	KR‐1b	19	0.143 (0.139, 0.146)	0.229 (0.225, 0.234)	0.378 (0.368, 0.387)
Possession	PS‐1a	17	0.147 (0.143, 0.150)	0.223 (0.218, 0.227)	0.341 (0.330, 0.351)
	PS‐1b	19	0.148 (0.144, 0.151)	0.227 (0.222, 0.231)	0.348 (0.339, 0.358)
Sth Georgia	SG‐H‐1	20	0.133 (0.129, 0.136)	0.218 (0.213, 0.223)	0.390 (0.380, 0.400)
	SG‐H‐2	11	0.111 0.107, 0.114)	0.221 (0.216, 0.225)	0.498 (0.486, 0.509)
	Global	178	0.140 (0.137, 0.142)	0.225 (0.222, 0.228)	0.380 (0.374, 0.385)

Sites within an island discerned by a letter suffix (1a, b, c..) are separated by <5 km; sites within an island discerned by a number suffix (1, 2…) are separated by >10 km. Diversity indices for T. bisetosa from South Georgia were found to significantly differ from all other islands (p < .05), as shown in bold.

Abbreviations: H O, observed heterozygosity; H E, expected heterozygosity; F IS, inbreeding coefficient; N, sample size.

Substantial genetic differentiation was found across all populations for T. bisetosa (global F ST = 0.542), with F ST estimates as high as 0.677 between sites on different islands and estimates ranging from 0.027 to 0.346 among sites within an island (Table S9). The largest degree of variation occurred between islands (40%), with a small amount partitioned between sites within islands (10%) and the remaining variance explained by among‐ and within‐individual variation (hierarchical AMOVA; Table 2). The PCoA for T. bisetosa also demonstrated a high degree of clustering by island, which captured 43.6% of variation within the first two axes, with subtle clustering by site also evident within South Georgia and Marion Island (Figure 1a). Genetic variation captured by a third PC axis dropped to 5.5%, with this axis further exaggerating the differences between island and, to some degree, the clustering of sites within South Georgia and Marion Island (Figure S3).

Table 2

Analysis of molecular variance (AMOVA) outcomes for Tullbergia bisetosa and Hypogastrura viatica based on 5,680 and 7,275 genome‐wide SNPs, respectively

	T. bisetosa			H. viatica
	SS	% var	F‐statistic	SS	% var	F‐statistic
Among island	82,165.13	40.0	0.40	29,815.00	10.1	0.10
Among site within island	19,914.76	10.3	0.17	10,565.68	1.9	0.02
Among samples within site	77,553.09	14.4	0.29	186,693.77	29.0	0.33
Within samples	45,497.50	35.3	0.65	99,756.50	59.0	0.41

Abbreviations: % var, per cent of total variation; SS, sum of squares.

The maximum likelihood (ML) tree for T. bisetosa was well supported with most nodes having >90% bootstrap support (Figure 1b). Each island formed a distinct clade, and further clear groupings were associated with intra‐island sites for South Georgia and Marion Island (these islands also had the highest F ST values between their sites; Table S9). Despite having some of the furthest spread intra‐island sites, Macquarie Island only showed weak grouping by site and intra‐island clades were not discernible within Heard Island.

Seven clusters (K = 7) were identified by fastSTRUCTURE to best represent the data for T. bisetosa, although the resultant structure plot showed five clusters to be most explicit (Figure S4). These corresponded to each island, with an extra division separating the furthest‐spaced sites on South Georgia, as was clear on the ML tree also. An additional cluster was assigned to individuals at low frequency and sporadically across all source populations, which may represent admixture with an unsampled island (T. bisetosa is also found on Prince Edward Island and the South Sandwich Islands, which we were unable to sample), or missing data occurring at the same SNPs.

Numerous fixed allelic differences were evident between populations of T. bisetosa (Table S10). The amalgamation of populations without significant fixed differences resulted in seven discrete “pseudo‐populations” separated by an inferred lack of gene flow. These are indicated by numbers 1–7 on the ML Tree (Figure 1b) and correspond closely to the clustering of branches within the tree, representing the four islands plus additional subdivisions within South Georgia and Marion Island. These pseudo‐populations were also largely in keeping with the fastSTRUCTURE clusters.

Hypogastrura viatica

Diversity estimates for each population of H. viatica are provided in Table 1; none of the populations were found to differ significantly from one another in observed or expected heterozygosity (global H O = 0.140, global H E  = 0.225). Populations consistently exhibited high levels of inbreeding, ranging from F IS = 0.334 to F IS = 0.498 (global F IS = 0.380; Table 1) indicating that up to half of the expected heterozygosity was lost.

Genetic differentiation across all populations of H. viatica was modest (global F ST = 0.115; Table S11). The vast majority of genetic variation occurred within sites or individuals (88% total), with a small amount partitioned between islands (10%) and a negligible amount between sites within islands (2%) (hierarchical AMOVA; Table 2). Weak geographic structure was also demonstrated by the PCoA, which, although revealing a clustering by island (with the exception of the Possession and Kerguelen that clustered together), only captured 12.8% of variation within the first two axes (Figure 1c). Visualization of a third axis revealed further differentiation between sites on South Georgia and Macquarie Island, though only explaining an additional 1.3% of variation (Figure S5).

The ML tree for H. viatica was weakly supported at several nodes, and, in line with AMOVA results, branch lengths between individuals within populations were longer than branches between islands (Figure 1d). Three clades were evident, however, corresponding with the three clusters in the PCoA (South Georgia, Macquarie Island and Possession + Kerguelen Islands) and weak intra‐island clades for South Georgia and Macquarie Island. Three clusters (K = 3) were also identified by fastSTRUCTURE (Figure S6), which were in agreement with the main island groupings in the PCoA and ML tree, and which again did not discriminate between Possession and Kerguelen Islands.

There were very few fixed allelic differences between populations of H. viatica (Table S12), and none were above the false‐positive rate. Thus, none of the populations could be distinguished by a complete lack of gene flow; rather, they were all amalgamated into a single pseudo‐population.

Mitochondrial DNA diversity

A total of 70 COI sequences were retrieved for T. bisetosa, originating from South Georgia, Marion, Macquarie and Heard Islands, as well as a single sequence from Chile (Table S3.1). The haplotype network is provided in Figure 2a; most haplotypes from the sub‐Antarctic islands were separated by one or two mutational steps. The nine sequences from Heard Island represented a single haplotype, which was also the most frequent haplotype at Macquarie Island and additionally present on Marion Island. Limited observations can be made for South Georgia, which was only represented by a single sequence, although this haplotype was the closest related to the only continental haplotype (from Chile).

Figure 2

Median‐joining haplotype networks for the indigenous springtail Tullbergia bisetosa (a) and the invasive springtail Hypogastrura viatica (b), based on COI sequences retrieved from GenBank and BOLD. Circle size is proportional to the haplotype frequency, black nodes represent inferred unsampled haplotypes, and hatches show the number of inferred mutations between haplotypes. Sequences are colour‐coded by their country or island of origin

For H. viatica, 94 COI sequences were retrieved originating from South Georgia, Macquarie, Possession and Kerguelen Islands, plus an additional six countries globally (Table S3.2); the haplotype network is shown in Figure 2b. The most common haplotype had a global distribution, occurring on all four sub‐Antarctic islands and Canada, New Zealand and South Africa, and often at relatively high frequency. Additional haplotypes from the sub‐Antarctic islands were separated from this common global haplotype (and from each other) by several mutational steps. Possession and Kerguelen Islands shared a haplotype, which was at least five mutations from the common haplotype. Even further diverged (separated from the common haplotype by at least six mutations) were two distinct haplotypes from South Georgia—one shared with Greenland, the other with Norway—and a haplotype occurring at relatively high frequency on Macquarie Island, which was shared with New Zealand and Denmark.

DISCUSSION

We present some of the first nuclear genomic data for sub‐Antarctic soil fauna, revealing that gene flow between sub‐Antarctic islands is highly restricted in the indigenous springtail Tullbergia bisetosa, with structure also evident within some islands. Numerous COI haplotypes are shared among islands, however, potentially reflecting postglacial colonization. For the invasive springtail Hypogastrura viatica, genomic structure is much less apparent, particularly for islands linked by human transport, though broad island groupings remain discernible and diverged COI haplotypes reflect a diverse invasion history. Bearing in mind that each species was represented by a unique SNP data set, our findings suggest that increased human movement throughout the sub‐Antarctic is likely to have very different evolutionary consequences for indigenous and invasive soil communities, with wider implications for inter‐island biosecurity policy and for island systems globally.

Pronounced phylogeographic structure in Tullbergia bisetosa

Within T. bisetosa populations, genetic diversity was low and substantial heterozygote deficits indicated inbreeding. While null alleles potentially contributed, they are unlikely to generate such high inbreeding coefficients alone, particularly given the SNPs were filtered by read depth and call rate (Andrews, Good, Miller, Luikart, & Hohenlohe, 2016). Inbreeding may have been partly due to assortative mating or facultative parthenogenesis, which have been observed in some springtail species (Chernova, Potapov, Savenkova, & Bokova, 2010; Van der Wurff, Isaaks, Ernsting, & Straalen, 2003). Indeed, nonrandom mating would also contribute to the high levels of genetic isolation between populations and, together with the large degree of fixation, may help explain the low overall diversity (global H E = 0.099) of T. bisetosa populations (Charlesworth, 2003).

Populations of T. bisetosa from different sub‐Antarctic islands share little to no gene flow, as revealed by all methods of genomic analysis. The large number of SNPs fixed for different alleles suggests that island populations are likely on independent evolutionary trajectories and allopatric speciation between islands may already be underway. Indeed, a large proportion of springtails from the Southern Ocean region are endemic to a single island (Baird, Janion‐Scheepers, et al., 2019), a situation likely promoted by allopatric genetic divergence (Losos & Ricklefs, 2009). Lacking flight, springtails disperse passively; therefore, genetic isolation over thousands of kilometres of open ocean is not surprising (see Costa et al., 2013).

Populations of T. bisetosa separated by several kilometres within both Marion Island and South Georgia were also strongly differentiated, as indicated by high F ST values, fixed allelic differences and distinct clades on the ML tree. In contrast, populations within Heard Island and Macquarie Island were not discernible by ML nor by fixed difference analysis, with much lower F ST values. While this may reflect the smaller geographic distance separating Heard Island populations (<1 km), populations on Macquarie Island were separated by up to 21 km, vastly exceeding the distances between differentiated Marion Island populations (1–3 km). Differences in local population structure of T. bisetosa may therefore reflect the unique geological history of each island. South Georgia and Marion Island have undergone numerous pulses of glaciation (and volcanism on the latter), which would have promoted population differentiation over time through physical fragmentation (e.g. Mortimer et al., 2012; Myburgh, Chown, Daniels, & Vuuren, 2007), potentially also explaining the higher genetic diversity on South Georgia. In contrast, Macquarie Island has not experienced glacial or volcanic activity, while Heard Island has been so extensively glaciated that it is unlikely to have been colonized by T. bisetosa until relatively recently, following deglaciation (Hodgson et al., 2014). Signatures of this postglacial colonization may also help to explain the seemingly incongruent mitochondrial results observed for T. bisetosa.

The COI data set for T. bisetosa revealed very different patterns to the nuclear genomic data, demonstrating the tenuousness of COI‐based inferences for sub‐Antarctic fauna to date (also see Halanych & Mahon, 2018). The sharing of a common COI haplotype between Marion, Heard and Macquarie Islands could be inferred to show a high degree of connectedness between populations; considering SNP data however, it is clear that there is presently negligible nuclear gene flow between islands. Rather, this haplotype may reflect the retention of an ancient polymorphism involved in the initial colonization of these islands (Fraser et al., 2012), particularly given its central relationship to most other variants. That it was the only haplotype found on Heard Island (from nine samples covering multiple sites) is indicative of founder effect, in accordance with a more recent, postglacial colonization of this island (see Wouw et al., 2008, for similar observations in plants). Alternatively, it remains possible that the lack of any clear mitochondrial subdivision between islands reflects selection acting on COI or another region of the mitochondrial genome (Ballard & Whitlock, 2004; Dowling, Friberg, & Lindell, 2008).

Our findings for T. bisetosa substantiate concerns that human transport of biota between sub‐Antarctic islands may lead to genetic homogenization of otherwise highly structured species (Hughes & Convey, 2010; Lee & Chown, 2011). The low genetic diversity, high inbreeding coefficients and marked genetic differentiation observed for T. bisetosa may indicate localized adaptation of populations, which could be swamped by the transport of individuals between islands, resulting in decreased fitness (“migration load”: see Bolnick & Nosil, 2007; Lenormand, 2002). Breaching the isolation between islands with human‐mediated gene flow could even disrupt the progress of allopatric speciation and associated biodiversity outcomes (Crispo et al., 2011; Seehausen, 2006). If T. bisetosa populations are, however, reproductively isolated and can no longer introgress, inter‐island exchange of these newly speciated populations remains a threat to the region's biogeography, which boasts a high degree of endemicity (Baird, Janion‐Scheepers, et al., 2019). The ultimate course of these microevolutionary processes depends largely on the degree of selection acting on populations (Bourne et al., 2014; Rieseberg & Burke, 2001) and may play out differently for nuclear versus mitochondrial genomes (see Toews & Brelsford, 2012). Regardless, allopatric population divergence plays a significant role in generating island biodiversity (Losos & Ricklefs, 2009); thus, any disruptions to this process via transport‐driven homogenization could pose a threat to island diversity patterns globally. Whether this is especially the case for soil fauna requires testing across other taxa, though high population structure is indeed an emerging characteristic of indigenous soil biota elsewhere (e.g. Costa et al., 2013; von Saltzwedel, Scheu, & Schaefer, 2016).

Human‐mediated gene flow and multiple invasions in Hypogastrura viatica

Most of the genomic variation within H. viatica was partitioned within rather than among populations, with strong levels of inbreeding indicated by heterozygote deficits. High inbreeding coefficients in H. viatica may reflect similar mechanisms to those discussed for T. bisetosa; in particular, parthenogenesis has previously been proposed to underlie the successful invasion history of H. viatica (Greenslade, 2002). Genetic diversity was not particularly low within H. viatica populations (global H E = 0.225), as would be expected if populations had been seeded by a single colonization of few individuals. Indeed, COI data indicated that multiple introductions have occurred on each island, a process known to maintain genetic diversity in other invasive populations (Kolbe, Larson, Losos, & Queiroz, 2008; Lavergne & Molofsky, 2007; Garnas et al., 2016). All of the islands contained representatives of a common global COI haplotype, as well as additional unique haplotypes shared with individuals from different countries. Most of the haplotypes within each island were diverged by numerous mutations and thus are much more likely to reflect colonization by differentiated H. viatica lineages than in situ diversification (see Piertney et al., 2015). The mixing of lineages from different colonization events helps to explain the large proportion of genomic variation attributed within populations of H. viatica rather than at higher geographic levels, and may also be contributing to the heterozygote deficits observed (via Wahlund effect) if different lineages within islands are yet to completely interbreed (see Kolbe et al., 2008).

Although nuclear genomic structure was not particularly pronounced for H. viatica, PCoA and fastSTRUCTURE plots did discern clearly between three main island groups, with only Kerguelen and Possession island populations relatively indistinguishable. Population structure was also weakly evident between the two sites on South Georgia and the most distal sites on Macquarie Island, revealing that some degree of genetic drift has occurred within and among islands, since the introduction of H. viatica. That this geographic differentiation remained only moderate, however, explaining such a small amount of total genetic variation, likely reflects the limited evolutionary time over which the species has existed in the region, meaning it is unlikely to be at migration–drift equilibrium (Slatkin, 1987). Despite being separated by >1,000 km of open ocean, populations from the islands of Kerguelen and Possession were profoundly more homogenous than populations merely 15 km apart within other islands. Kerguelen and Possession are both governed by France and therefore receive cargo from the same ship annually (Frenot, Gloaguen, Massé, & Lebouvier, 2001), a link that no other islands in our study share. This indicates that shipping has played a role in facilitating genetic exchange of the invasive populations or at least that the two islands have been colonized from the same source population, demonstrating the greater influence of human transport over natural dispersal in structuring H. viatica across the sub‐Antarctic. Very similar findings have also been made for invasive aphids from Kerguelen and Possession Islands (Lebouvier et al., 2011).

One COI haplotype for H. viatica occurred on all four sub‐Antarctic islands and most other sampled countries, indicating potentially rapid global spread of a successful genotype (Le Roux, Wieczorek, Wright, & Tran, 2007; Garnas et al., 2016). Yet, the haplotypes unique to each island provided most insight on the colonization history for H. viatica. The occurrence of two haplotypes on South Georgia shared with Norway and Greenland directly accords with Norway's whaling history on the two islands, which saw substantial shipping transport between them (Tønnessen & Johnsen, 1982). Similarly, the sharing of a Macquarie Island haplotype with New Zealand, the most proximal land mass with historic, science and tourism links to the island, is unsurprising. A connection with Denmark is also in keeping with the long history of Danish vessels chartered to access Macquarie Island (Kriwoken & Williamson, 1993). Finally, the sharing of a unique haplotype solely between Kerguelen and Possession Islands likely reflects their mutual visitation exclusively by France, as for the SNP data. Thus, the unique nature of the sub‐Antarctic islands, each governed and visited by different nations, may influence the population structure of introduced species in the region, with discrete colonizations from different sources more evident than range expansion from a single introduction (also see Baird, Janion‐Scheepers, et al., 2019). It is, however, important to note that H. viatica occurs on several more islands across the Southern Ocean than were sampled in this study, for which the colonization history may well differ.

In the context of increased human activity throughout the sub‐Antarctic, the relatively weak geographic structure and negligible fixed allelic differences for H. viatica imply that movement of these individuals between islands is unlikely to disrupt evolutionary trajectories or biogeographic patterns to the same degree as for T. bisetosa. Rather, intraregional transport may help to spread beneficial alleles introduced by the unique colonizations of each island, an important prerequisite to realizing any advantage from the diversity provided by multiple introductions (Barker et al., 2019; Dlugosch & Parker, 2008). The connectivity among islands linked by shipping indicates that ongoing inter‐island transport has the potential to spread H. viatica beyond what might be achieved through natural range expansion. This is likely to apply to other invasive soil species, as natural dispersal is often limited belowground (Cicconardi et al., 2017; Costa et al., 2013). For islands already invaded by H. viatica, additional introductions may continue to contribute new genetic variation, as evidenced by differentiated COI lineages. While genetic diversity is not always a prerequisite for invasion success (see Dlugosch et al., 2015), it is likely to enhance an invasive's evolutionary persistence and spread (e.g. Ahlroth, Alatalo, Holopainen, Kumpulainen, & Suhonen, 2003; Lavergne & Molofsky, 2007; Wagner et al., 2017). Therefore, the management concern of multiple introductions (see Le Roux & Wieczorek, 2009) now also appears pertinent to remote island soil systems (see also Novo et al., 2015).

Biosecurity implications

Soil is a key pathway by which humans unknowingly transport species beyond their natural range (Hulme, 2015); therefore, insights from soil biota are particularly relevant to biosecurity. Our study provides evidence that indigenous sub‐Antarctic soil fauna may be vulnerable to genetic homogenization by inter‐island transport, with potential consequences for population fitness, biodiversity processes and/or biogeography. While it is evident, at least for our focal species, that such transport is yet to occur, the degree and extent of biotic exchange will likely increase in the future due to expanding tourism and science activities involving multiple islands (e.g. Stewart, Espiner, Liggett, & Taylor, 2017; Walton & Thomas, 2018). The threat of inter‐island homogenization likely applies to remote islands globally—while generally understudied, other island‐based soil faunas also show evidence of naturally high genetic and biogeographic structure (e.g. Caruso, Diega, & Bernini, 2005; Cicconardi, Nardi, Emerson, Frati, & Fanciulli, 2010). Substantial allopatric genetic drift between populations, as demonstrated here, is known to play a key role in driving island diversity worldwide (Losos & Ricklefs, 2009). Therefore, our findings emphasize calls to consider the genetic consequences of human transport, in order to conserve biodiversity processes (Hudson, Viard, Roby, & Rius, 2016; Petit, 2004; Seehausen, 2006).

Our work also indicates that repeat introductions of existing invasives are an important management consideration for remote island systems, with invasive spread likely promoted by intraregional transport. Again, this is in line with recent findings for soil ecosystems elsewhere (e.g. Novo et al., 2015; Porco et al., 2013). While most studies concerning island biosecurity risks have focused on novel, external introductions, inter‐island transfers are now being recognized as a major concern globally (e.g. IUCN, 2018; Lee & Chown, 2011; Petit, Hoddle, Grandgirard, Roderick, & Davies, 2009). For the sub‐Antarctic islands, which are largely World Heritage‐listed, stringent among‐island biosecurity measures should be considered and would serve as an exemplar for conserved island systems elsewhere.

CONFLICT OF INTEREST

None declared.

Click here for additional data file.