Evidence for Cryptic Diversity in the “Pan-Antarctic” Springtail Friesea antarctica and the Description of Two New Species

The invertebrate terrestrial fauna of Antarctica is being investigated with increasing interest to discover how life interacts with the extreme polar environment and how millions of years of evolution have shaped their biodiversity. Classical taxonomic approaches, complemented by molecular tools, are improving our understanding of the systematic relationships of some species, changing the nomenclature of taxa and challenging the taxonomic status of others. The springtail Friesea grisea has previously been described as the only species with a “pan-Antarctic” distribution. However, recent genetic comparisons have pointed to another scenario. The latest morphological study has confined F. grisea to the sub-Antarctic island of South Georgia, from which it was originally described, and resurrected F. antarctica as a congeneric species occurring on the continental mainland. Molecular data demonstrate that populations of this taxon, ostensibly occurring across Maritime and Continental Antarctica, as well as on some offshore islands, are evolutionarily isolated and divergent and cannot be included within a single species. The present study, combining morphological with molecular data, attempts to validate this hypothesis and challenges the taxonomic status of F. antarctica, suggesting that two additional new species, described here as Friesea gretae sp. nov. and Friesea propria sp. nov., are present in Continental Antarctica.


Introduction
Human-driven environmental change, of which the long-term consequences are only now being appreciated, currently threatens biodiversity, the stability of ecosystems and the services they provide. Knowledge of biological processes that sustain endangered natural environments is relatively incomplete, especially so for remote and isolated ecosystems, such as those of polar regions. Actions to counteract the impact of global environmental change on threatened ecosystems rely on evaluation of their resilience to disturbance [1]. Simulations to predict the impact of changing abiotic factors have already been designed, although their ability to model future trajectories is impeded by geographical heterogeneity and distinctive biological diversity. In Antarctica, poor knowledge of species diversity and of their dispersal ability/resilience are considered limiting factors that prevent the implementation of sustainable conservation plans [2][3][4][5]. The dispersal abilities of species [6,7] plaques (absent in F. antarctica) and a different number and shape of both clavate tenant hairs and anal spines [3]. Recently, specimens from near the Russian Molodyozhnaya Station, recorded by Wise [32] as F. grisea and Bulavintsev [33] as Friesea sp., were found to be another new species, F. eureka [34].
The anomalous genetic divergence observed between populations of F. antarctica from the Antarctic Peninsula (AP) and Victoria Land (VL) sites (e.g., an overall genetic distance of 21.7% between the mitochondrial genomes from specimens of the two regions [35]) contrasts with an identical morphology [26,35]. The presence of cryptic species hidden within the F. antarctica complex was hypothesized and an in-depth morphological revision of the taxon suggested [35]. No molecular data are so far available from the type locality of F. grisea to allow comparison. Within F. antarctica, population genetic and phylogenetic studies have identified well-defined evolutionary lineages, here referred to as AP (Antarctic Peninsula), VL (Victoria Land) and Cape Hallett, clearly separated from each other in geographical terms. However, as yet, no morphological analyses are available to support a distinction between these three lineages. This may be due to either phylogenetic niche conservatism (e.g., [36]; see Discussion section) or morphological diversity that has been previously overlooked. The current study evaluates their morphology and genetic differentiation through a multidisciplinary comparative analysis. We tested the efficacy of an integrative taxonomic approach to assess whether or not cryptic speciation took place, while continuing morphological uniformity has mainly been selected for by the extreme environmental conditions (i.e., phylogenetic niche conservatism). Morphological revision of F. antarctica specimens from both the Antarctic Peninsula and Victoria Land, with specimens from Cape Hallett further compared against both regions, was performed. Morphological data were re-evaluated in the light of molecular data. Different bioinformatic tools for species delimitation, based both on genetic distances and phylogenetic reconstruction, were used to analyze nuclear and mitochondrial markers. Descriptions of two new species, previously included within F. antarctica and emerging from both morphological and molecular analyses, are provided.

Sample Collection
Specimens initially identified as F. antarctica (known at that time as F. grisea) were collected from both the Antarctic Peninsula and Victoria Land (sampling sites and respective abbreviations given in Figure 1 and Table 1). Specimens from Victoria Land were collected during the Italian National Antarctic Program (PNRA) expeditions (2005/2006, 2017/2018 and 2018/2019). Individuals of F. antarctica from the Antarctic Peninsula were sampled during the 2002/2003 expedition through a collaboration between the British Antarctic Survey (BAS) and the PNRA. After their collection from field habitats using a suction device, specimens were subjected to preliminary morphological identification and stored at −80 • C until further morphological and molecular analyses.

Morphological Analyses
Between five and 20 specimens of F. antarctica, from the localities highlighted by asterisks in Table 1, were incubated in lactic acid (at 35-45 • C) and cleared by short immersion (2-10 min) in 10% KOH. Samples were mounted on slides and observed under a Leitz Laborlux S microscope. Some other specimens were mounted on slides after clearing in Nesbitt to improve the observation of chetae and spines on Ab. V-VI, as needed. When possible, body length, the cephalic diagonal and antennae, as well as the length and width of each thoracic and abdominal segment, were measured. Since data distributions were not normal, the nonparametric Mann-Whitney test with Bonferroni correction for multiple comparisons was applied to assess whether or not there were significant differences in these parameters, as well as in the number of setae in the male and female genital openings. Calculations were performed using R 3.3.2 [37], comparing AP vs. VL, AP vs. Cape Hallett and VL vs. Cape Hallett specimens (see Results).
Some samples were also prepared for scanning electron microscopy. Firstly, they were dehydrated in absolute ethanol followed by critical-point drying in a Balzer CPD 030. Samples were coated Insects 2020, 11, 141 4 of 28 with gold in a Balzer MED 010. Scanning electron microscopy observations were performed using a Philips XL20.
Some samples were also prepared for scanning electron microscopy. Firstly, they were dehydrated in absolute ethanol followed by critical-point drying in a Balzer CPD 030. Samples were coated with gold in a Balzer MED 010. Scanning electron microscopy observations were performed using a Philips XL20.  Table  1 for collecting site acronyms and geographical locations.

Molecular Dataset
The datasets assembled by Torricelli et al. [26] and Carapelli et al. [2] were enlarged, adding specimens of F. antarctica from Danco Island and Detaille Island (Table 1), as well as four outgroup species: the congeneric F. topo (from Alexander I.) and F. salmoni, as well as Bilobella aurantiaca and Folsomotoma octooculata. All specimens were screened not only for two previously obtained mitochondrial markers (cox1 and atp6 [2,26]) but also for two further nuclear markers (28S and EF-1α; large ribosomal subunit and elongation factor subunit 1α, respectively).
Whole genomic DNA was extracted from F. antarctica and outgroup specimens using the Wizard ® SV genomic DNA Purification System (Promega, Madison, WI, USA) and eluted in 50 μl ddH2O. The nuclear markers analyzed in the present study (28S and EF-1α) were amplified and sequenced with universal primer pairs.   Table 1 for collecting site acronyms and geographical locations.

Molecular Dataset
The datasets assembled by Torricelli et al. [26] and Carapelli et al. [2] were enlarged, adding specimens of F. antarctica from Danco Island and Detaille Island (Table 1), as well as four outgroup species: the congeneric F. topo (from Alexander I.) and F. salmoni, as well as Bilobella aurantiaca Insects 2020, 11, 141 5 of 28 and Folsomotoma octooculata. All specimens were screened not only for two previously obtained mitochondrial markers (cox1 and atp6 [2,26]) but also for two further nuclear markers (28S and EF-1α; large ribosomal subunit and elongation factor subunit 1α, respectively).
Whole genomic DNA was extracted from F. antarctica and outgroup specimens using the Wizard ® SV genomic DNA Purification System (Promega, Madison, WI, USA) and eluted in 50 µl ddH 2 O. The nuclear markers analyzed in the present study (28S and EF-1α) were amplified and sequenced with universal primer pairs. The oligonucleotides D1a [5 CCCSCGTAAYTTAAGCATAT 3 ] and D7b [5 GACTTCCCTTACCTACAT 3 ], as well as the nested primers D3a [5 GACCCGTCTTGAAACACGGA 3 ] and D5b1 [5 ACACACTCCTTAGCGGA 3 ] were used to amplify and sequence the large ribosomal subunit (28S rDNA). Specific primers (sequence available on request), matching the flanking portions of the deleted regions of the 28S ribosomal RNA (rRNA) pseudogene, were tested for their ability to amplify this DNA fragment in all samples. The primers M3 [5 CACATYAACATTGTCGTSATYGG 3 ] and rcM4 [5 ACAGCVACKGTYTGYCTCATRTC 3 ], as well as the nested M44.1 [5 GCTGAGCGYGARCGTGGTATCAC 3 ] and rcM51.1 [5 CATRTTGTCKCCGTGCCAKCC 3 ], were applied to amplify and sequence the EF-1α. PCRs were performed in a 25 µL reaction volume consisting of: 2.5 µL of whole genomic DNA from each sample, 1.25 µL of both forward and reverse primers (10 µM), 12.5 µL of GoTaq ® Long PCR Master Mix (Promega, Madison, WI, USA) and 7.5 µL of ddH 2 O. Amplifications were run on a GeneAmp ® PCR System 2700 (Applied Biosystems, Foster City, CA, USA) thermal cycler. The initial denaturation step was set at 95 • C for 5 min, followed by 35 cycles at: 94 • C for 1 min, 50-55 • C (depending on marker) for 1 min and 60 • C for 90 s. The final elongation step was set at 72 • C for 7 min. PCR products were then purified with the kit Wizard ® SV Gel and PCR Clean-up (Promega, Madison, WI, USA) and sequenced on both strands with a DNA Analyzer ABI 3730 at the core facility of the Biofab Research Lab (Rome, Italy). The sequences were then manually corrected and assembled using the software Sequencher 4.4.2 (Gene Codes Corporation, Ann Arbor, MI, USA). For each of the nuclear and mitochondrial markers, the dataset was aligned using the online tool Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) and manually corrected in Mesquite [38]. This latter program was further applied to concatenate the four loci into a multi-locus dataset. The four single-and the multi-locus alignments were then used for phylogenetic and species delimitation analyses.
For both phylogenetic optimization criteria, the best models of evolution were selected before the tree search, considering different dataset partitions: 1st; 2nd and 3rd codon positions for the protein-encoding genes (cox1, atp6 and EF-1α) and one single partition was considered for the 28S. The EF-1α dataset was further partitioned into coding and noncoding regions (i.e., the three codon positions, as well as the introns) (models of evolution applied in the present study are summarized in Table 2).
For ML analyses, the best models of evolution were detected through ModelFinder [41] based on the Bayesian information criterion (BIC) and a greedy strategy. The ML tree searches were performed with various settings of perturbation strength (-pers command) of 0.2; 0.5; 0.7 and stop condition (-nstop command) of 100, 200 and 300. The support values were obtained with 1000 replicates of ultrafast bootstrap approximation [42].
Before the BI analyses, the best models of evolution were selected by means of the software PartitionFinder v2.1 [43] based on the Akaike's information criterion (AIC) and a greedy strategy. The models of evolution that best fitted our datasets were then applied to MrBayes 3.2, run with four chains (three hot and one cold) for 10 6 generations with a sampling frequency of one tree every 1000 iterations and the first 25% of the tree topologies discarded (burn-in step).

Species Delimitation Analyses
Given the high genetic differentiation, previously detected through mitogenomic and population genetic studies [26,35], and the morphological uniformity of the F. antarctica specimens, several methods of species delimitation based on molecular markers were applied. Both the single-(cox1, atp6, 28S and EF-1α genes, individually aligned) and multi-locus (concatenation of cox1, atp6, EF-1α and 28S loci) datasets were tested, starting with species discovery approaches (see Sections 2.5.1-2.5.3 below). Then, putative species clusters detected were further investigated through a validation method of species delimitation (see Section 2.5.4 below). Both distance-and phylogeny-based approaches were used.

Distance-Based Species Delimitation
The cox1 single-locus dataset was processed with the online tool automatic barcode gap discovery (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html [25]). The analysis was performed with default settings (relative gap width of 1.5; intraspecific divergence: between 0.001 and 0.1), as well as applying the Kimura two-parameter (K2P) model when calculating the distribution of pairwise distances.

Phylogeny-based methods: the Poisson Tree Process
The identification of species boundaries was also carried out using phylogenetic methods, based either on the number of substitutions (bPTP [44]) or on divergence time estimations (GMYC [45]; see Section 2.5.3). Since both branch length (that may reflect the number of nucleotide substitutions) and the topology depend on the optimization criteria applied to the inference analysis, the four single-locus datasets were used for both BI and ML tree search, as described above. The two topologies obtained for each of the four molecular markers were then analyzed through the Poisson Tree Process (PTP) on the web server bPTP (https://species.h-its.org/ptp/ [44]). The analyses were run for 500,000 Markov chain Monte Carlo (MCMC) generations and with a burn-in value of 0.25, both including and excluding the outgroup species.

Phylogeny-Based Methods: The Generalized Mixed Yule Coalescent Model
In order to define species boundaries according to the divergence time estimations, the cox1 dataset was run with the software BEAST 2.4.8 [46], thus obtaining an ultrametric tree. The software PartitionFinder v2.1 was applied to detect the models of evolution that better fitted our single-locus dataset, organized into only one charset. A strict molecular clock was applied to the cox1 alignment, defining the clock rate based on the average mutation rate per million years identified in Brower [47] and Papadopoulou [48], with a coalescent model of constant population size as tree prior. Two independent MCMC runs were performed, constituting of 10 6 generations with parameters sampled every 1000 iterations and a burn-in of 0.25. The two runs were combined through the BEAST package LogCombiner Insects 2020, 11, 141 7 of 28 2.4.8 [46]. The convergence of MCMC chains was assessed with Tracer v 1.7 [49] and the tree visualized through the software FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).
The ultrametric topology, obtained as outlined before, was then used as an input for the single-threshold GMYC analysis [45]. This latter was conducted in R 3.3.2 [37] with the use of the splits package [50].

Validation Approach of Species Delimitation
The hypothetical species boundaries detected were then tested and validated under the multispecies coalescent model, by means of the software BP & P (Bayesian Phylogenetic & Phylogeographic [51]). Both the single-(cox1, atp6, 28S and EF-1α genes individually aligned) and the multi-locus datasets were used for the validation analyses of species delimitation. For each analysis, joint species delimitation and species-tree inference were performed (i.e., speciesdelimitation = 1 and speciestree = 1; A11 [51]). The algorithm 0 and the default settings for fine-tuning parameters were used (ε = 5), as well as the species model prior 1 (uniform probability for rooted tree). The species delimitation analyses were run with the following combination of gamma distributions of the parameters θ (ancestral population size) and τ (root age): (1) θ:G (1:10), τ:G (1:10); (2) θ:G (2:100), τ:G (2:500) and (3) θ:G (2:1000), τ:G (2:200). The analyses were run for 100,000 MCMC generations, with a sampling frequency of 50 and a burn-in of 1000 generations; each analysis was run twice in order to confirm the consistency of the results.

Systematics
Characters that have been considered species-specific in the past for the genus are: number of ocelli, number of clavate tenent hairs, degree of reduction and form of the furca, number of teeth on the retinaculum, number and position of s chaetae on antennal segment IV and number of anal spines [3,52].
All specimens of the F. grisea complex are identical, as far as number of ocelli and form of furca are concerned. We have found the following additional characters useful in distinguishing species: the relative lengths and form of the chaetae on the body, particularly on abdomen V and VI; whether spines, spine-like or normal chaetae are present on abdomen VI; the relative lengths of s chaetae and ordinary and macrochaetae on abdomen V, as well as the presence of dorsal cuticular granulation, body size and colour of adults and development of clavate tenent hairs and clavate macrochaetae. In spite of some intra-population variability, the reliability of these characters is reinforced by the number of specimens we have examined and by the observed consistency with the results of molecular analyses. Line drawings are presented for each species (F. antarctica: Figures 2 and 3; F. propria: Figure 4; F. gretae: Figures 5 and 6), scanning electron microscope details are provided below (Figures 7 and 8).

Morphological Description
Friesea antarctica (Willem, 1901). = Achorutoides antarcticus (Willem, 1901 and1902 [53,54] described A. antarcticus using materials from Harry Island, Gerlache Straits, erecting a new genus for it, Achorutoides. It was later synonymized by [55] with F. grisea. It was re-erected by Greenslade [3] for specimens from the Antarctic Peninsula and South Shetland Islands. Willem's description and figures [54] are detailed and accurate except for minute details of the mouthparts. We add descriptions of the chaetotaxy and give sizes and ratios of chaetal lengths here. The original type material cannot be now found. It is probable that Willem's original collection was deposited at Ghent University, Belgium, although the author did not indicate where his material was deposited, as it was not customary at that time. Nevertheless, efforts to trace it failed, and the type material probably no longer exists. Thus, according to Article 75.3.4 of the International Code of Zoological Nomenclature, a neotype is erected here from a location as near as possible to the original type locality, Brabant Island. The original specimen was collected on Gerlache's Belgian expedition in the ship Belgica, who explored the Gerlache Straits and its islands in January 1898.
Neotype  Table 1 for geographical details of collecting sites. Specimens will be conserved in the South Australia Museum and in the Collembola Collection of the Department of Life Sciences in the University of Siena.
Habitus size and colour: Neotype 1.8 mm and Paraneotype 1.9 mm, both types on same slide. Average length 1.55 mm, same in males and females, range from 1.16 mm to 2.00 mm. Original body length in Willem = 1.4 mm [54].
Body color dark blue, grey, sometimes lighter on the ventral side, especially in specimens preserved for a longer period.

Comments
Willem's 1902 figures show well the important characters of his species. He illustrates its squat and broad habitus and dark black color, unlike F. grisea from South Georgia, which is more elongate and a lighter grey, and the tibiotarsus with 8 clavate tenent hairs, although he states 6 are present in the text. The furca, tenaculum, ocelli patch and dorsal abdomen VI are also well-figured. However, there are some slight errors in that he does not show the tubercular insertion of L on the labium, his ventral segment numbers are incorrect, and so he places the furca ventrally on abdomen V and the maxilla lacks one lamella. Additionally, he illustrates the mandible with 11 teeth, while we find the number to be 8.
Friesea propria sp. nov.: Greenslade (Table 3). The shape of the body is similar to this latter species, as is the color, which appears dark blue (sometimes brown dark), with some specimens lighter after longer preservation in alcohol.
Head and Antennae: Antennae shorter than the cephalic diagonal (ratio 1.41); Ant. III and IV dorsally fused. Ant. I and II with 7 and 12-14 chaetae, respectively; Ant. III with 17-20 chaetae; the sensory field includes the AOIII constituted by two curved sensilla placed in a cuticular pit, two subcylindrical guard sensilla and a ventral ms. Chaetotaxy of Ant. IV with six cylindrical sensillae (S1-4 and S7-8), a ms between S7 and S8, a subapical organite and an apical vesicle, usually trilobate. Ocelli 8 + 8, each group as 5 + 3 ocelli with 3 chaetae between them (Figures 4A and 7E). PAO absent. Ventral line of the head with 2 + 2 chaetae; labium with 4 + 4 chaetae in the bmf, 6 + 6 chaetae in the blf and 4 + 4 proximal chaetae and papillated chaetae L ( Figure 8C); mandible with 3-4 distal and 3 basal teeth; maxilla as in F. gretae ( Figure 6J). Labrum with 3, 5 and 4 chaetae and a0 and m0 much longer than a1 and m1, respectively; m2 and a2 very lateral and 2 + 2 preclypeal chaetae.  Habitus Size and Color: The body length of the holotype is 0.95 mm; the average length of specimens studied is 1.28 mm, with values ranging from 0.95 to 1.60 mm. The size of F. gretae sp. nov. is significantly smaller than F. antarctica and larger than of F. propria (Table 3); the shape of the body (especially in the specimens preserved in alcohol before the preparation of the slides) appears cylindrical, slenderer than those in F. antarctica and F. propria.
Cuticle: Integument with uniformly distributed medium-small size granules, without any plaque or tubercles along the body.
Chaetotaxy: F. gretae n. sp. also shows a paurochaetose chaetotaxy both dorsally and ventrally ( Figures 5A,B,D), constituted by sensory, meso-and macro-chaetae, more developed in the last abdominal segments (i.e., V and VI Abd.); as in the previous species, some asymmetries and small individual differences were observed. Thoracic tergum I with 4 + 4 chaetae and thoracic tergum II with 12 + 12 chaetae, Di with a2, De with S chaetae in p4, Dl with 1 ms and 1 S chaeta. Thoracic tergum III as the previous one but without a2. Abdominal terga I-IV Di with 3 chaetae and De with 4 + S Thorax: Chaetotaxy of the legs: subcoxae 1 of legs I, II and III with 1, 2 and 2 chaetae; subcoxae 2 of legs I, II and III, usually with 0, 2 and 2 chaetae. Coxae of legs I, II and III with 3, 8 and 8 chaetae; trochantera I, II and III with 6 chaetae each ( Figure 4E). Femora I, II and III, usually with 13-14, 12 and 10-11 chaetae, respectively; variations sometimes observed. Tibiotarsi of legs I, II and III with 19, 19 and 18 chaetae, of which 5-6 clavate, usually 3-4 external and 1-2 internal ( Figure 4F,G); distal whorl with 11 chaetae; also, in this species, as in F. antarctica, variability in the distal end shape of clavate chaetae was observed; claw untoothed and empodial appendage absent ( Figure 4F,G).
Abdomen VI anal spines two, short and broad at p1, me chaeta at a1; p0 short chaeta ( Figure 4H).    Table 1 for geographical details of collecting sites. Specimens will be conserved in the South Australia Museum and in the Collembola Collection of the Department of Life Sciences in the University of Siena.
Etymology: Named after Greta Tintin Eleonora Ernman Thunberg, young advocate for action on climate change.
Habitus Size and Color: The body length of the holotype is 0.95 mm; the average length of specimens studied is 1.28 mm, with values ranging from 0.95 to 1.60 mm. The size of F. gretae sp. nov. is significantly smaller than F. antarctica and larger than of F. propria (Table 3); the shape of the body (especially in the specimens preserved in alcohol before the preparation of the slides) appears cylindrical, slenderer than those in F. antarctica and F. propria.
Cuticle: Integument with uniformly distributed medium-small size granules, without any plaque or tubercles along the body.

Phylogenetic Analyses
The single locus phylogenetic trees obtained through both BI and ML showed a clear distinction between three main F. antarctica lineages (AP, VL and Cape Hallett), with moderate to high statistical support (PP = 0.92-1.00; BS = 58-79; Figures 9 and 10). In all single locus BI and ML topologies, except that based on the cox1 dataset, F. topo (specimens collected on Alexander I.) is basal to VL and Cape Hallett samples and not to F. antarctica from the AP, as simple geography would suggest. The same phylogenetic relationships are obtained when the four-locus matrix is used for comparison ( Figure 11). In this latter instance, all the nodes were strongly supported, further suggesting that the three F. antarctica clusters may be well-defined different evolutionary entities. In addition, when the amplification of the 28S rRNA-encoding gene was performed using the primer pair D1a/D7b, one single amplification product (of~2100 bp) was obtained from almost all samples, with the exception of those of F. antarctica from VL (seven specimens from as many localities but not from Cape Hallett; Table 1), where two bands of different size were coamplified and sequenced. Of these, the largest band was clearly homologous to the single amplified product that was obtained from the AP and Cape Hallett samples and corresponds to the complete and expected gene-fragment size; the smaller probably corresponds to a 28S rRNA pseudogene representing a copy missing more than 400 bp. Attempts to amplify the smaller band in specimens collected elsewhere from Victoria Land and in those from Cape Hallett were unsuccessful, suggesting that the pseudogene copy may have evolved only in the ancestor of this F. antarctica group, for which it could be considered a synapomorphic character.

Species Delimitation Analyses
Regardless of the species delimitation method applied, all the discovery approaches tested suggested that F. antarctica should be considered a species complex. The lineages corresponding to AP, VL and Cape Hallett were identified as different taxa by most of the discovery methods applied to the single-locus datasets. Overall, the statistical support obtained for each of these hypotheses was moderate to high (Figures 9 and 10). Few exceptions to this pattern were observed. First, when the PTP method was applied to the BI topology of the cox1 dataset and the ultrametric tree was used in the GMYC analysis, the AP lineage was split into two clusters: one including only the haplotype P6 and the other including all the remaining cox1 sequences. Second, the PTP method applied to the ML topology of the EF-1α dataset suggested that the VL and Cape Hallett lineages should be considered as a single species. However, this analysis provided low statistical support to the single species hypothesis for the continental Antarctic (i.e., VL + Cape Hallett, PP = 0.44) and for the AP (PP = 0.09; Figure 11). Third, the application of the PTP method to the 28S dataset proposed VL and Cape Hallett lineages as a single species with high posterior probability values ( Figure 10A). This outcome was observed applying both the BI and ML topology.
The validation method of species delimitation was employed giving as prior information the presence of three different species corresponding to the AP, VL and Cape Hallett lineages, since most of the discovery analyses supported this hypothesis. During the validation analyses, all the possible combinations of species were further tested, and the posterior probability of a given delimitation was calculated.
In all instances (single-and multi-locus datasets, as well as the three combinations of θ and τ parameters), the presence of three distinct species was robustly supported (PP between 0.70 and 1.00; . When the existence of one species (AP + VL + Cape Hallett) or two (AP vs. VL + Cape Hallett) was suggested, the posterior probability values were very low (all PP < 0.01; data not shown). The only ambiguous output was obtained with a combination of θ:G (2:1000), τ:G (2:200) parameters applied to the cox1 and atp6 single-locus datasets; this output suggests the existence of a single species including all three F. antarctica lineages plus the outgroups, a suggestion that is obviously not tenable (data not shown).   In all instances (single-and multi-locus datasets, as well as the three combinations of θ and τ parameters), the presence of three distinct species was robustly supported (PP between 0.70 and 1.00; . When the existence of one species (AP + VL + Cape Hallett) or two (AP vs. VL + Cape Hallett) was suggested, the posterior probability values were very low (all PP < 0.01; data not shown). The only ambiguous output was obtained with a combination of θ:G (2:1000), τ:G (2:200) parameters applied to the cox1 and atp6 single-locus datasets; this output suggests the existence of a single species including all three F. antarctica lineages plus the outgroups, a suggestion that is obviously not tenable (data not shown).

Discussion
The application of an integrated taxonomy approach, here tested for springtail species, appears to be useful in delimiting species boundaries in cases in which morphological diagnostic characters are ambiguous [56,57]. A similar method has recently been applied to Antarctic terrestrial nematodes and was successful in detecting nematode lineages in the Maritime Antarctic [58]. The present study further confirms the efficacy of an integrated taxonomic approach to better define and characterize different ecosystem biodiversity.
Morphological analyses have thus far struggled to find diagnostic characters among different local forms of F. antarctica. On the other hand, the morphological reassessment performed herein in light of the outcome of the molecular analysis led to the identification of some characters that, though marginal, are consistent with the molecular results. Specifically, one character (setae in subcoxa 2) is fixed as 0-3-3 in populations from the Antarctic Peninsula and as 0-2-2 in populations from Victoria Land (inclusive of Cape Hallett). Three other characters (body length, setae on the male and female genital opening) are quantitative in nature but display a statistically significant difference among groups (Tables 3 and 4).
The analysis of mitochondrial, as well as nuclear markers, provided consistent results, as well as supporting the efficacy of the barcode cox1 gene in the identification of species boundaries. The only exception to this outcome was represented by the 28S, which, at variance, supported the

Discussion
The application of an integrated taxonomy approach, here tested for springtail species, appears to be useful in delimiting species boundaries in cases in which morphological diagnostic characters are ambiguous [56,57]. A similar method has recently been applied to Antarctic terrestrial nematodes and was successful in detecting nematode lineages in the Maritime Antarctic [58]. The present study further confirms the efficacy of an integrated taxonomic approach to better define and characterize different ecosystem biodiversity.
Morphological analyses have thus far struggled to find diagnostic characters among different local forms of F. antarctica. On the other hand, the morphological reassessment performed herein in light of the outcome of the molecular analysis led to the identification of some characters that, though marginal, are consistent with the molecular results. Specifically, one character (setae in subcoxa 2) is fixed as 0-3-3 in populations from the Antarctic Peninsula and as 0-2-2 in populations from Victoria Land (inclusive of Cape Hallett). Three other characters (body length, setae on the male and female genital opening) are quantitative in nature but display a statistically significant difference among groups (Tables 3 and 4).
The analysis of mitochondrial, as well as nuclear markers, provided consistent results, as well as supporting the efficacy of the barcode cox1 gene in the identification of species boundaries. The only exception to this outcome was represented by the 28S, which, at variance, supported the presence of a single species from the Continental Antarctic (VL + Cape Hallett; Figure 10A). We can hypothesize that the latter two taxa (i.e., VL and Cape Hallett) started to diverge more recently, and therefore, that the large ribosomal subunit may have not yet accumulated sufficient mutations to discriminate between the two lineages. Nevertheless, during the amplification of the 28S gene (portion between D1 and D7), the presence of two copies of the large ribosomal subunit was observed which differ by about 470 bp in length. This feature was present only in the specimens of F. antarctica sampled in Victoria Land, with the exception of those from Cape Hallett. These latter, together with the specimens from AP, possessed only the longer copy of the 28S. All phylogenetic reconstructions support F. topo as basal to an assemblage of F. antarctica populations from Victoria Land (inclusive of Cape Hallett), with the exclusion of F. antarctica populations from the Antarctic Peninsula. This would make F. antarctica as currently defined polyphyletic, while it is consistent with the proposed subdivision of F. antarctica into different lineages.
The sharp genetic distinction between the three F. antarctica lineages is not striking per se, as recent studies already put forward the hypothesis that F. antarctica may represent a cluster of multiple cryptic species [4,35]. In fact, striking was the previous recognition of the presence over Antarctica of a single "pan-Antarctic" springtail species, especially considering that most of the sampling sites correspond to physically remote geographic regions (Antarctic Peninsula and Victoria Land) and also taking into account the vicariant origin of Antarctic springtails in general [59,60] and their limited dispersal capabilities. This latter characteristic is likely to underlie the differentiation of a third evolutionary lineage in Cape Hallett, the northern-most sampling locality along the Victoria Land coastline, which is separated from all the other sites by the major barrier of the Tucker Glacier ( Figure 1 and Table 1). This geographic barrier seems to be insurmountable for F. antarctica, as has already been shown to be for other springtail species endemic to Victoria Land (K. klovstadi and C. cisantarcticus) that show high levels of genetic divergence between populations from either side of the Tucker Glacier (e.g., [4]). This level of fragmentation raises significant conservation issues. This is particularly the case for F. gretae, which is currently known only from its type locality at Cape Hallett and is now one of three terrestrial micro-arthropods whose type locality is this location (Antarctic Treaty Secretariat, 2015, available at http://documents.ats.aq/recatt%5Catt567_e.pdf). These occur within an Antarctic specially protected area (ASPA No. 106), an area declared both to protect a significant Adélie penguin population and in recognition of its exceptional vegetation development and terrestrial biodiversity by regional standards. This ASPA, however, faces the risk of long-term disturbance to its terrestrial habitats due to human research activities and, in particular, through the very limited options available for erecting and operating field camps. As highlighted by [30], the agreed management plan of this ASPA requires updating, as the primary recommended camping area overlaps with the known distribution of this new species of endemic springtail. Protection of "the type locality or only known habitat of any species" is listed as one of the primary criteria supporting the declaration of areas as ASPAs (Annex V, Article 3(2) of the Environmental Protocol to the Antarctic Treaty) (https://www.environments.aq/informationsummaries/specially-protected-and-managed-areas-in-antarctica-updated/). The highly conserved morphology observed between the three putative species suggests the occurrence of morphological stasis (or phylogenetic niche conservatism), a phenomenon that arises when a species is under particularly strong selective pressures (adversity selection [61]). In a context such as the extreme Antarctic terrestrial environment, the selection of particular characters is greatly influenced by abiotic factors that, over Antarctica, are the major forces shaping terrestrial biota [19,62,63]. Among them, freezing temperatures, drought and exposure to high levels of UV radiation would likely support the selection of darker pigmentation and smaller body size, as well as earlier sexual maturity due to the very short active season [3]. In this context, the significant difference in body size and the reduced number of setae in the genital openings observed between F. antarctica lineages may be considered of substantial biological importance.

Conclusions
The Antarctic Continent and Peninsula environments are predictably severe, a defining characteristic for adversity selection [61,62], which would be here operating on organisms. In this type of selection, morphological stasis is the norm, because the common drivers of morphological adaptation, such as biotic and climatic variability, are low or absent. Competition and predator pressures are low because of extremely low species numbers and the types of habitats present are few. In habitats where adversity selection operates, there is selection against dispersal and for any change in the organism, because it is already highly adapted to severe conditions. Consequently, it is not then surprising that morphological changes in the Friesea species group in Antarctica has been extremely slow over long time periods. The same phenomenon is found in other Antarctic springtail genera, such as the genus Tullbergia [64].
The identification of one single springtail species whose distribution spans areas as distant as Victoria Land and the Antarctic Peninsula, contrasting with all other Antarctic species that are restricted to a single region, together with its fairly high levels of "intra-specific" genetic divergence, has puzzled researchers for over two decades. Nevertheless, hindered by a marked inter-individual population-level variability in the context of an otherwise morphologically uniform taxon, classical morphological taxonomy had thus far failed to find evidence for multiple species.
In this study, a massive sampling effort and multi-locus molecular data and phylogenetic analyses, as well as species delimitation methods, were applied to the study of F. antarctica diversity on a continental scale. These led to the identification of three well-defined and distinct lineages. Based on this initial hypothesis, the morphology of specimens from the entire continent was revised, and a number of subtle morphological characters were identified that provide further support for the three lineages being species-level entities: one corresponds to F. antarctica and two are here formally described as F. gretae nov. sp. and F. propria nov. sp. From a methodological standpoint, we underline the strength of the combined application of molecular, biogeographic and morphological data, an approach that has gained momentum in the last 10-15 years as "integrated taxonomy". From the point-of-view of Antarctic Collembolan diversity, this study resolves the enigma of the uniqueness of a single "pan-Antarctic" species and provides a well-supported scenario leading to a formal taxonomic assessment of the three species.
Funding: This study was funded by the Italian Program of Research in Antarctica (PNRA16_00234). Partial support was also provided by the University of Siena. P.C. is supported by NERC core funding to the British Antarctic Survey's 'Biodiversity, Ecosystems and Adaptation' Team. The paper also contributes to the SCAR 'State of the Antarctic Ecosystem' (AntEco) international program.