Genetic diversity among populations of Antarctic springtails (Collembola) within the Mackay Glacier ecotone

Complete List of Authors: Beet, Clare; University of Waikato, School of Science Hogg, Ian; University of Waikato, School of Science Collins, Gemma; University of Waikato, School of Science Cowan, Don; University of Pretoria, Centre for Microbial Ecology and Genomics, and Genomics Research Institute Wall, Diana; Colorado State University, Evolutionary Ecology Laboratories, Department of Biology Adams, Byron; Brigham Young University, Evolutionary Ecology Laboratories, Department of Biology, Monte L. Bean Life Sciences Museum


Introduction
Antarctica's terrestrial biota is exposed to cold temperatures, low moisture levels, and steep chemical gradients (Convey et al. 2014;Hogg et al. 2014;Velasco-Castrillón et al. 2014a).These conditions, coupled with repeated glacial cycles over the past 80 million years, have shaped a terrestrial landscape characterised by low biodiversity (Adams et al. 2006;Convey 2011).Furthermore, this limited diversity is spread across the highly fragmented 0.32% of Antarctica that is currently considered ice-free (Pugh and Convey 2008).Accelerating rates of climate change, anthropogenic activity, and a heightened risk of invasive species all have the potential to severely disrupt present-day Antarctic biodiversity and thereby ecosystem functioning (Wall 2005;Hogg and Wall 2011;Chown et al. 2012).To understand and monitor changing biodiversity patterns, Terrestrial Observation Networks have been proposed to track environmental changes and consequences for the biota (e.g., Levy et al. 2013).Such networks will require a range of suitable locations to adequately characterise any biotic changes.
Ecotones, transitions from one biological community to another, may be ideal sites to study shifts in species distributions as a consequence of global climate changes (Gosz 1993;Risser 1993).The Mackay Glacier in southern Victoria Land is situated immediately to the north of the McMurdo Dry Valleys.Here, the landscape changes from one dominated by the expansive Dry Valleys to one of small rocky outcrops isolated by glaciers.These changes in habitat are also reflected in the communities of invertebrates, particularly springtails and mites.For example, only one species of springtail (Gomphiocephalus hodgsoni) is known from the McMurdo Dry Valleys, whereas an additional two species (Antarcticinella monoculata and Cryptopygus nivicolis) are found north of the Dry Valleys (Wise 1967;Hogg et al. 2014).Accordingly, the Mackay Glacier vicinity can be considered an ecotone or transitional zone (sensu Risser 1993).
Within the Ross Sea region, three biogeographic zones have been designated (northern Victoria Land, southern Victoria Land, and the Queen Maud Mountains), each with three, and in one case four, reported species of springtail for a total of 10 species within the region (Salmon 1965;Wise 1967;Adams et al. 2006;Sinclair and Stevens 2006).However, the basis of this zonation was derived from early work utilizing classical morphological taxonomy and sampling which had been undertaken in proximity to historical bases or camps.Vast areas were still largely unexplored (e.g., see Peat et al. 2007;McGaughran et al. 2011).
Springtails are prominent features of Antarctic terrestrial ecosystems as they are the largest terrestrial invertebrates and are also highly sensitive to environmental disturbances, making them ideal biological indicators of climate change (Hopkin 1997;Hogg et al. 2014;Collins and Hogg 2016).Of the three currently recognized species of springtail within the wider Mackay Glacier region (Gomphiocephalus hodgsoni (Carpenter 1921), Antarcticinella monoculata (Salmon 1965), and Cryptopygus nivicolus (Salmon 1965)), G. hodgsoni is relatively widespread and well-studied, while C. nivicolus and A. monoculata have received comparatively less attention (Adams et al. 2006;Stevens et al. 2006;Hogg et al. 2014).For example, A. monoculata was previously known only from two extremely restricted locations (Salmon 1965;Bennett et al. 2016).
The aim of this study was to assess springtail distribution and genetic (mitochondrial cytochrome c oxidase subunit I, COI) diversity within the vicinity of Mackay Glacier.These data will establish the current distributional ranges of species and baseline levels of genetic diversity, against which any future changes can be determined.Such knowledge can then be used to inform, and contribute to, developing Terrestrial Observation Networks (Levy et al. 2013).

Materials and methods
Individuals of three species of springtail were collected in January 2015 from 12 sites in the vicinity of the Mackay Glacier covering a north-south range of approximately 120 km (Table S1 2 ; Fig. 1).Animals were collected directly from rocks using modified aspirators (Stevens and Hogg 2002) or from soil samples.At each of the 12 sites, 3-7 soil samples were collected within a 1 km 2 area (via random trowelling within a metre squared area at each sub-sampling location).Soil samples weighed between 300-500 g and were returned to the laboratory at McMurdo Station with individuals removed via floatation of samples in a sucrose solution (Jenkins 1964).All individuals were preserved in 95% ethanol for subsequent genetic analyses.
Whole individuals were placed into 95-well plates and sent to the Canadian Centre for DNA barcoding (CCDB), University of Guelph, for sequencing of the barcode region of the COI gene (sensu Hebert et al. 2003).Genomic DNA was extracted via the AcroPrep TM PALL Glass Fibre plate method (Ivanova et al. 2006).A 658-bp fragment of the mitochondrial COI gene was then amplified in accordance with standard CCDB protocols (see Ivanova et al. 2006;Ivanova andGrainger 2007a, 2007b) using the primer cocktails C_LepFolF (5=-ATTCAACCAATCATAAAGATA-TTGG-3=) and C_LepFolR (5=-TAAACTTCTGGATGTCCA-AAAAATCA-3=) (Folmer et al. 1994;Hebert et al. 2004;Ivanova et al. 2006).Successfully amplified products were then cleaned using Sephadex ® before being sequenced in both directions on an ABI 3730xl DNA analyser.Of the 139 springtails extracted, 95 produced complete 658-bp sequences.All photographs, collection information, and sequence data have been added to Barcode of Life Data Systems (www.boldsystems.org;Ratnasingham and Hebert 2007) and are available from the dataset DS-ANTSP (dx.doi.org/10.5883/DS-ANTSP)and cross referenced to GenBank (accession numbers: KX232683-KX232856).
Sequences were confirmed as the target species via comparison with BOLD sequences and translated to amino acids to check for the presence of stop codons.Sequences were aligned using MUSCLE in Geneious 7.1.9and all trimmed to 527 bp to incorporate 69 sequences from a previous study by Bennett et al. 2016 and three Symphypleona outgroups (Sminthurides aquaticus IHCO046-03, Sminthurides malmgreni IHCO047-03, and Sminthurinus el-egans COONT067-08).This alignment was then used in all subsequent analyses except for the haplotype network construction, which did not include the outgroup taxa.Chi square ( 2 ) tests conducted in PAUP* 4.0 (Swofford 2002) were used to determine whether base frequencies were equal among all sites, identify parsimony-informative sites, and designate first, second, or third codon positions.The most appropriate model of evolution was determined using jModelTest 2.1.1 under the lowest AIC criterion (Posada 2008).Bayesian trees were generated using BEAST v2 (Bouckaert et al. 2014).A strict clock model and speciation yule process as the tree prior were employed in BEAUTI v2, with the Markov chain Monte Carlo (MCMC) set at 50 000 000 generations, sampling trees every 5000 generations.The Bayesian analysis was run in BEAST, with the quality of the results evaluated in TRACER v1.5 (Drummond et al. 2010).The burn-in inputted into Tree Annotator v2 was calculated by plotting log-likelihood values against the generation time in TRACER (Rambaut and Drummond 2007) and determining the point at which normalization occurred.The final tree was visualized in FigTree v1.4.0.Neighbour Joining (NJ) and Maximum Likelihood (ML) analyses were conducted in MEGA v5.05 (Tamura et al. 2011).ML and NJ settings both included 1000 bootstrap replicates, with GTR+I+G used as the model of evolution for ML and Tamura-Nei in NJ.Tamura-Nei was used as it allows unequal base frequencies and multiple substitution types (Simon et al. 2006).All other settings were set to the default options in MEGA.MEGA was also utilized to create a pairwise distance matrix to calculate intraspecific divergences (uncorrected pair-wise distances), while Barcode Index Numbers (BINs) were assigned by BOLD (Ratnasingham and Hebert 2013) and used as a measure of Molecular Operational Taxonomic Units (MOTUs).An analysis of molecular variance (AMOVA) was performed using Arlequin v3.5 to further evaluate genetic variation.Uncorrected pairwise distances were calculated for haplotypes and analysed with 1600 permutations.Genetic structure was enforced with haplotypes grouped hierarchically into reported species and then BINs.Haplotype networks were constructed in TCS v.1.21using all 164 sequences.PAUP* was used to determine whether the data were evolving in a clock-like manner.A likelihood ratio test (-2logL = 2(logL0 -logL1), n-2 df where n is the number of terminal taxa) was conducted using the likelihood scores of trees with and without a molecular clock enforced (L0, L1).Subsequently, a molecular clock was calibrated by dividing uncorrected p-distances by mutation rates of 1.5%-2.3%sequence divergence per million years (sensu Knowlton et al. 1993;Brower 1994;Stevens et al. 2006) and then used to estimate the timing of divergence events.A Bayesian tree with a strict clock and rate set as 0.0115 to simulate 2.3% sequence divergence per million years, with all other settings the same as in the previous Bayesian analysis, is provided (Fig. S2 2 ).We caution that if dispersal in springtails is male-biased, the maternally inherited COI marker would underestimate gene flow between populations.However, previous analyses of Antarctic springtails (e.g., Stevens and Hogg 2003) have suggested that COI data are congruent with allozyme analyses and that COI is a reliable indicator of genetic differences within and among populations.
Phylogenetic analyses revealed three previously unreported genetically distinct MOTUs (>2% sequence divergence), one for each of the three recognized species within the region.NJ, ML, and Bayesian analyses all produced concordant results, with all nodes at and above the BIN level having high support (Fig. 2).Maximum intraspecific (within morphologically recognized species) divergences ranged from 5% to 11.2%, resulting in a total of seven Barcode Index Numbers (BINs) for the wider Mackay Glacier region.Based on the AMOVA, 64.96% of the genetic variation was found among the three recognized species, 32.81% was found among the distinct intraspecific populations (BINs), and the remaining 2.24% was within the divergent populations (BINs).Two of the BINs, one each for G. hodgsoni and C. nivicolus, were found only at Towle Glacier, indicating a potentially important site for arthropod diversity within the region.Antarcticinella monoculata was found at two new locations (Benson Glacier, Pegtop Mountain) and C. nivicolus at an additional two locations (Tiger Island, Towle Glacier).The widespread G. hodgsoni, previously reported as far north as Mt.Murray (Wise 1967), was not detected or collected at this location.
We examined haplotype networks to visualize the fine-scale geographic relationships between these "new" populations and other Mackay Glacier locations.From the 121 G. hodgsoni sequences, we found 27 unique haplotypes.Overall, these haplotypes exhibited high levels of genetic connectivity among populations throughout the study area, albeit with some subtle geographic structuring (Fig. 3).Four of the 27 haplotypes occurred in more than one location, with the most common haplotype (Gh1) found at four sites up to an estimated 55 km apart.However, two populations of G. hodgsoni were found to be highly divergent with the newly found population at Towle Glacier (Gh23, 24), 28 mutational steps (7.3% sequence divergence) from Mackay Glacier populations (>45 km).In contrast, the other divergent population (Gh25), previously identified at Mt. Gran (within the Mackay Glacier vicinity) by Bennett et al. 2016, was 32 mutational steps and 7.3% divergent from other nearby populations.
Cryptopygus nivicolus haplotypes exhibited broadly similar patterns to that of G. hodgsoni, although on a smaller scale owing to their more restricted distribution (Fig. 3).A total of seven haplotypes and three BINs were identified from the 21 C. nivicolus sequences obtained.The main haplotype (Cn1) was found in three locations separated by distances of up to 30 km, with the distinct Towle Glacier population (Cn5) being 24 mutational steps (5% sequence divergence) away from Mackay haplotypes.
Antarcticinella monoculata were found in very low numbers at all sites where they were found, with the exceptions of Mt.Murray and Cliff Nunatak (Fig. 3).Of the 23 individuals we sequenced from these two nearby locations (approximately 6 km apart), all had the identical haplotype (Am1).This haplotype was also 50 mutational steps away, and had sequence divergences of 11.2%, from A. monoculata collected from the more southern Mackay Glacier locations.Both G. hodgsoni and C. nivicolus had high levels of genetic variability within the immediate Mackay Glacier area (77°S), including highly divergent haplotypes for each species (Mt.Gran haplotype Gh25 at 37 steps, Springtail Point haplotypes Cn6 and Cn7 at 17 steps).In contrast, A. monoculata populations found throughout the Mackay Glacier area (excluding northern haplotype) were all similar, separated by a maximum of six mutational steps.
The application of standard molecular clock calibrations of 1.5%-2.3%sequence divergence per million years (sensu Knowlton et al. 1993;Brower 1994;Stevens et al. 2006) suggested that individuals from three of the populations from Towle Glacier and Mt.Murray (Gh23, Gh24, Cn5, Am1) diverged within the last 5 million years.The C. nivicolus haplotypes diverged most recently (2.2-3.3 million years ago), followed by G. hodgsoni at 3.2-4.9 million years ago, with the most divergent group of A. monoculata estimated to have been isolated between 4.9 and 7.5 million years ago.

Discussion
The mtDNA COI sequences revealed a considerably greater diversity of collembolans than previously known in the Ross Sea region of Antarctica.In particular, we found several highly divergent populations of springtails including three previously unrecognized BINs for the three putative species, for a total of seven BINs within the wider Mackay Glacier area.The high levels of divergence suggest the presence of potentially cryptic species and that diversity in Antarctica, while comparatively depauperate in the global sense, is much more diverse than once thought (Griffiths 2010;Chown et al. 2015).The use of molecular techniques such as DNA barcoding to study diversity has been particularly helpful in detecting genetically divergent populations (e.g., Stevens andHogg 2003, 2006;Hawes et al. 2010;Stevens and D'Haese 2014).A study by McGaughran et al. (2010) found distinct MOTUs within the springtail Cryptopygus antarcticus, while Mortimer et al. (2011) identified two divergent groups within Western Antarctic Peninsula populations of the mite Alaskozetes antarcticus.Similarly, Torricelli et al. (2010aTorricelli et al. ( , 2010b) ) found high levels of COI divergence (14.4%-17.2%)between Victoria Land and Antarctic Peninsula populations of the "pan-Antarctic" springtail Frisea grisea.Such high levels of variability among springtails and other Antarctic biota could be the result of a number of factors including genetic drift, environmental heterogeneity, and potential diversifying selection in these extreme and variable environmental conditions (Fanciulli et al. 2001).
Habitat fragmentation restricts gene flow, leading to allopatric speciation, and is a likely explanation for the different BINs found only at the Towle Glacier site.The genetic uniformity of all of the northern A. monoculata individuals could indicate that the population has also been derived allopatrically, with genetic drift leading to a reduction in diversity and fixation of the single haplotype (e.g., Costa et al. 2013).Alternatively, the population could be comparatively new, resulting from a founder effect from a more genetically diverse population nearby.Although there was a general trend of increasing divergence with distance, similar to observations of Fanciulli et al. (2001), there were a number of G. hodgsoni haplotypes that were shared among locations up to 55 km apart, demonstrating that populations currently exist in both sympatry and allopatry.
The genetically distinct populations (BINs) were all estimated to have diverged within the last 5 million years (range: 2.2-7.5 million years).While there is controversy surrounding rates of molecular evolution, a number of studies have similarly reported clades of Antarctic invertebrates diverging within the Pliocene (Hawes et al. 2010;McGaughran et al. 2010;Bennett et al. 2016).It was during this time that the Western Antarctic Ice Sheet (WAIS) was thought to have completely collapsed (Pollard and Deconto 2009).This would have resulted in sea levels rising and increased dispersal opportunities for Collembola via rafting in meltwater streams and open seaways (sensu Hawes 2011), before the WAIS eventually reformed and glaciers isolated newly established populations (e.g., McGaughran et al. 2010).This could explain why some G. hodgsoni haplotypes are shared at sites up to 55 km apart despite being currently isolated by glaciers.Although springtails have been collected in air traps over small-scale distances (<8 km; see Hawes et al. 2007;McGaughran et al. 2011), long-distance wind dispersal is considered unlikely.Further, springtails are prone to desiccation due to their permeable cuticle and do not have the anhydrobiotic capacity such as that of nematodes, tardigrades, and rotifers which would facilitate Aeolian transport (Freckman and Virginia 1998;Hogg et al. 2014;Velasco-Castrillon et al. 2014a).
The region immediately to the north of the Mackay Glacier (<77°S) has received limited sampling attention, with this study providing the first records of species for many of the sites visited.The sampling of A. monoculata at Mt. Murray and Cliff Nunatak represented the first time they have been collected there since the original collections in 1958 (Salmon 1965;Wise 1967).Interestingly, G. hodgsoni was not collected from the previously recorded northern location (Mt.Murray).This absence could suggest that species have altered their distributions within the past 50 years (e.g., Stevens and Hogg 2002).It is possible that some springtails have been essentially living at their environmental limit and that recent climatic changes may have led to detectable range shifts.However, their absence from Mt. Murray could also be attributed to stochastic events at local sites (e.g., natural disturbance).The distribution of springtails is largely governed by the paucity of ice-free habitat, with abiotic factors such as salinity and the availability of liquid water further limiting available habitat (Hogg et al. 2006;Sinclair and Stevens 2006;Velasco-Castrillon et al. 2014b).
Historical glacial cycles appear to have influenced the current distribution of springtail populations, with habitat fragmentation and population dynamics refining these distributions.Improving our understanding of genetic variability is fundamental to preserving maximal springtail diversity and thereby the evolutionary potential of natural systems and populations to respond to changing environmental conditions (Adams et al. 2006;Hogg et al. 2006Hogg et al. , 2014)).Given the presence of several geographically restricted BINs within the Mackay Glacier region, any future changes in species' distributions can be easily tracked through monitoring of DNA barcodes within a site.This finer level of resolution will thereby enhance our capacity to detect potentially subtle biological responses resulting from gradual climate changes.

Fig. 1 .
Fig. 1. (A) General map of Antarctica showing the study area (black box) within the Ross Sea region of the Ross Dependency.Also pictured are the Western and Eastern Antarctic Ice Sheets (WAIS, EAIS) in addition to northern and southern Victoria Land (nVL, sVL) and the Queen Maud Mountains.Map adapted from the Antarctic Digital Database v6.0,British Antarctic Survey (http://www.add.scar.org/home/add6).(B) Map showing the location of the study area around the Mackay Glacier in relation to Ross Island.Map data: Google, DigitalGlobe.(C) Map of exact collection sites and distribution of the three springtail species.1, Mount Murray/Cliff Nunatak; 2, Towle Glacier; 3, Benson Glacier; 4, Tiger Island; 5, Mount Gran; 6, 7, Mount Seuss; 8, The Flatiron; 9, Sperm Bluff; 10, Springtail Pt; 11, St Johns Range.Map data: Google, DigitalGlobe.[Colour online.]

Fig. 2 .Fig. 3 .
Fig. 2. Collated phylogenetic tree of sequences representing 40 unique haplotypes for three springtail species including three Symphypleona outgroups (Sminthurides aquaticus IHCO046-03, Sminthurides malmgreni IHCO047-03, and Sminthurinus elegans COONT067-08).Maximum likelihood base tree with support values over 90/0.9 displayed in order of ML bootstrap values/Bayesian posterior probabilities/NJ bootstrap values.* indicates that the support values were either unavailable or under 90/0.9.Tree is coloured according to the seven BINs present, with bars indicating the location where specific haplotypes were collected.Numbers in brackets refer to specific locations pictured in Fig. 1C.For full Bayesian and NJ trees, see the supplementary data, Fig. S2 2 .

Zealand
Antarctic Research Institute (NZARI) for financial support.C.R.B. was the recipient of a New Zealand Post Antarctic Scholarship, a University of Waikato Environmental Research Institute Scholarship, and a Waikato Graduate Women Educational Trust Scholarship.Sequencing at the Canadian Centre for DNA Barcoding, University of Guelph, was supported through funding to the International Barcode of Life Project (iBOL) from the Ontario Genomics Institute (2008-OGI-ICI-03), Genome Canada, the Ontario Ministry of Research and Innovation, and the Natural Sciences and Engineering Research Council of Canada.B.J.A. and D.H.W. were supported by McMurdo LTER NSF OPP grant 1115245.D.A.C. gratefully acknowledges the financial support of the South African NRF SANAP program.