Specimens of P. steinenii were collected from freshwater bodies within three recognized biogeographic regions ( Figure 1 , for more detail see Supplementary Table S1 ): 1) Magellanic sub-Antarctic Ecoregion (MSE), in particular in Navarino Island, within the Cape Horn Biosphere Reserve; 2) the sub–Antarctic island of South Georgia (SG), including individuals from Bird Island (BI) and the South Georgian mainland (TP, STR and HB); and 3) the Maritime Antarctic (MA) specifically King George Island (KGI), Deception Island (DCI) and Livingston Island (LVI) in the South Shetland Islands (SSI). Adult flies were sampled during the austral summer from 2017 to 2023 from their natural habitats around lakes with the use of an entomological aspirator. Additionally, larvae and adults were also manually collected on Navarino Island from aquatic mosses in a lake near the summit of Bandera Hill, located close to the north coast of the island. All individuals were preserved in 95% ethanol and their geolocation was later included in GBIF (see Gañan et al., 2020). Additional searches were conducted in the Falkland/Malvinas Islands and South Orkney Islands, but no P. steinenii were observed. Taxonomic identity of the specimens was confirmed following Edwards and Usher (1985).

DNA was extracted using two main protocols, QIAGEN DNeasy Blood & Tissue, following (Maturana et al., 2022), and the QiAMP Extraction Kit, following the manufacturer’s protocols. All individuals were fully submerged in proteinase K+ATL buffer solution for 4 h at 56°C or overnight at 40°C; adults were not crushed to retain them as complete as possible for morphological examination, while larvae were partially or fully crushed before extraction.

Two DNA loci were amplified, two segments of the mitochondrial cytochrome c oxidase subunit I (cox1), and one segment of the nuclear DNA 28S rRNA gene. For the two segments of cox1, we used the universal primers LCO1490 and HCO2198 (Folmer et al., 1994) and UEA5 (AGTTTTAGCAGGAGCAATTACTAT) and UEA10 (TCCAATGCACTAATCTGCCATATTA) (Lunt et al., 1996), while for the 28S rRNA we used rD1.2a (CCCSSGTAATTTAAGCATATTA) and rD5a (GGYGTTGGTTGCTTAAGACAG) (Whiting, 2002). Detailed descriptions of the protocol and techniques used in the Polymerase Chain Reaction are given in Supplementary Material (S2) . The products obtained were commercially sequenced by LGC Genomics (Germany) and Macrogen (South Korea).

For the phylogenetic reconstructions we incorporated the sequences of P. steinenii into a comprehensive dataset depicting the main lineages of the genus Parochlus Enderlein1912 ( Supplementary Table S3 ). Homologous DNA sequences were independently aligned for each genetic locus using MAFFT v7.505 (Katoh and Standley, 2013) using the web server (Katoh et al., 2019). We conducted traditional phylogenetic estimations using both Maximum Likelihood (ML) and Bayesian Inference (BI) algorithms. The appropriate sequence evolution model for each genetic locus was determined using bModelTest v1.2.1 (Bouckaert and Drummond, 2017) and ModelFinder (Kalyaanamoorthy et al., 2017). The GTR+Γ model was used for cox1 and the HKY+I+Γ model for 28S rRNA. We prepared a phylogenetic tree for each locus to check congruence between both datasets ( Supplementary Figure S4 ). Subsequently, the sequences were concatenated, and we performed posterior analyses on a partitioned dataset, considering each genetic locus. ML reconstruction was carried out using RAxML v8.2.12 (Stamatakis, 2014) and the node support was assessed through a bootstrap analysis with 1000 pseudo-replicates. The tree was rooted with Zavrelimyia sp. following Cranston et al. (2010). Bayesian reconstruction was performed in MrBayes3 v2.6.7 (Ronquist et al., 2012). The BI analysis was run three times for 100 million generations each time, with tree sampling every 10,000 generations. The consensus tree considered a burn-in period of 25%. For the intraspecific phylogenetic tree using the complete cox1 dataset, we make a midpoint rooted Bayesian tree reconstruction based on GTR+CAT model for ML and GTR+Γ for BI, respectively. Convergence was assessed by checking that split frequencies had an average standard deviation below 0.01 and all parameters had effective sample sizes (ESS) > 200 using Tracer v.1.6 (Rambaut et al., 2014).

Levels of genetic polymorphism were determined for each locus using standard diversity indices including number of haplotypes/alleles (K), number of segregation sites (S), haplotype/allelic diversity (H), average number of pairwise differences (Π) and nucleotide diversity (π) for each of the sampling areas, including the entire MA (KGI, DCI and LVI), sub-Antarctic (TP, BI), and the MSE using DnaSP v6.12.3 (Rozas et al., 2017). We estimated the number of private alleles for each sampling area and region as a proxy to assess the potential existence of refugia (Maggs et al., 2008). Neutrality tests (Tajima’s D and Fu’s F _S) were performed to assess deviation from the neutral model. Significant negative values for these tests provide evidence of excess of rare polymorphism in a population, indicating either recent demographic expansion or positive selection. Genealogical relationships were estimated for mitochondrial sequences by constructing a median-joining haplotype network (Bandelt et al., 1999) using PopArt v1.7 (Leigh et al., 2015).

We estimated levels of genetic differentiation among the recognized groups (MSE, SG and SSI), among populations within those groups, and within populations, through mean pairwise differences (Φ_ST, using Kimura-2P genetic distances) and haplotype frequencies (F_ST) in ARLEQUIN v3.5.2.2 (Excoffier and Heckel, 2006). To test statistical significance of differentiation, we performed a permutation test (20,000 iterations). The p-value for pairwise Φ_ST and F_ST between populations was corrected using the false discovery rate correction (FDR; Benjamini et al., 2005). We used two different clustering methods to determine the spatial genetic structure and the spatial boundaries among them in P. steinenii. First, we used the GENELAND v4.0.7 package (Guillot et al., 2005) in the software R v3.6.1 (R Development Core Team, 2019). This Bayesian approach uses the Markov Chain Monte Carlo (MCMC) procedure to estimate the optimal clustering of samples based on geographic information. The analysis considers that spatially structured clusters are more probable than randomly distributed clusters in space. The most probable number of populations (K) was identified through ten independent MCMC analyses, each consisting of 10 million iterations with a thinning interval of 1000 iterations using the correlated frequency model. The range was limited between K = 1 and K = 5, with a burn-in of 25%. Second, we conducted an analysis of molecular variance (AMOVA) in Arlequin v3.5.2.2 to estimate the proportion of genetic variation explained by partitioning sampling sites into the different demographic units detected. This analysis characterizes spatial genetic structure by partitioning the variance within populations, among populations within groups and among groups.

We examined historical demography of the different previously identified genetic clusters (GENELAND and AMOVA), based on the cox1 marker. First, we compared observed mismatch distributions (Rogers and Harpending, 1992) implemented in DnaSP. We calculated Harpending’s Raggedness index (rH) in Arlequin v3.5.2.2 to test the unimodality of observed data. Second, we estimated population growth trends over time, through Bayesian Skyline Plots (BSP, Drummond and Bouckaert, 2014) implemented in BEAST v2.5 (Bouckaert et al., 2019). The running conditions considered 100 million iterations with parameter sampling every 10,000 steps, discarding the initial 25% of steps in the analysis. We used a population substitution rate of tenfold the evolutionary rate (7.5% per million year, Martin et al., 2002a, b) to estimate the most recent common ancestor (TMRCA) of each genetic group, as this more accurately reflects the rate at which new haplotypes appear (Ho et al., 2011, 2015). This better accounts for the time-dependence of molecular evolution at population level, brings divergence estimates closer to the present and avoids over-estimation of recent divergences when analysing intraspecific lineages.

To estimate divergence times, we utilized BEAST v2.5 with a reduced dataset comprising mitochondrial haplotypes retrieved using DnaSP. We performed a phylogenetic reconstruction combining a phylogenetic mutation rate with a population substitution rate (Ho et al., 2008). We estimated divergence times between the two main clades (Clade A and Clade B) using a phylogenetic rate, as it relates to the substitution rate at which mutation are fixed among clades. This was conducted using an optimized relaxed clock (ORC) model (Douglas et al., 2021) with a substitution rate of 0.75% per million years, a rate proposed for Chironomidae (Martin et al., 2002a, b) and previously used for Parochlus (Allegrucci et al., 2012). For the mutation rate, we used a population rate obtained to estimate the TMRCA of each clade, as it considers the rate at which new haplotypes appear in the clade. Three independent runs of 100 million generations each were conducted, with samples collected every 10,000 iterations. Convergence of results was assessed using Tracer v1.7.2 (Rambaut et al., 2018), and the results were summarized in a single ultrametric tree using TreeAnnotator v2.5.1 (Drummond et al., 2012). The final tree was constructed using the results of both approaches with the Chronos function in ape v5.7-1package (Paradis, 2013) in R, using a sequential secondary age calibration with the time of divergence of the main clades and TMRCA of each lineage (see another example in Pérez-Alvarez et al., 2021).

Within the framework of this study, we sampled from 16 freshwater sites within the species’ reported distribution across Maritime and sub-Antarctic regions and southern South America ( Figure 1A , for more details see Supplementary Material S1 , and Gañán et al., 2021). Two sites were in the MSE (IN = 9, NI = 12 sequences), four in sub-Antarctic South Georgia (BI = 5, TP = 16, STR = 14, HB = 11), and 10 in the South Shetland Islands (KGI = 40, DCI = 18, LVI = 26). The results of Xia’s test showed low substitution saturation (Iss = 0.2; Iss.c = 0.8; p < 0.05) under the assumption of both symmetrical and asymmetrical tree topologies for all OTU subsets. Similarly, PhyloMad results indicated a low risk of substitution saturation for cox1, considering codon positions 1–2 and position 3 independently. Regarding the simple linear regression plots of genetic distances (uncorrected vs. model corrected) identified deviations for cox1 codon position 3 only (R2 = 0.784). Plots for codons 1 and 2 combined indicated strong linear correlations between uncorrected and model-corrected genetic distances (R2 = 1). When codon 3 was combined with codons 1 and 2, the coefficient of determination remained high (R2 = 0.995), supporting the inclusion of codon 3 for estimating evolutionary relationships. For interspecific phylogenetic tree we used a concatenated data considering one fragment of 699 bp length of cox1 (under the Folmer universal primers) and 28S rRNA, resulting in 1347 bp. This concatenated dataset included 264 and 213 segregating sites for cox1, and 28S rRNA, respectively. For the intraspecific analyses, a total of 151 sequences for cox1 and 26 sequences for 28S rRNA were obtained. However, there was no variation in the 28S rRNA sequences within P. steinenii, so no further population analyses were carried out using this locus. We included a final alignment of two cox1 fragments – considering two mitochondrial fragments obtained from Folmer’s universal and UEA set of primers – of a total length of 1470 bp, with 80 variable sites (5.4%), and did not include insertions/deletions or stop codons.

The haplotype tree topology of P. steinenii identified the divergence of two main clades, (A) and (B), with high posterior probability (PP) and bootstrap support (BS) ( Figure 1B ). Clade B showed a well-supported SG clade, separated from remaining but unresolved SSI specimens ( Figure 1B ). The concatenated 28S rRNA and cox1 phylogenetic reconstruction recovered the relationships within the subfamily Podonominae and the monophyly of P. steinenii ( Figure 2 ). At the same time, individuals from the Maritime Antarctic (MA) populations were more closely related to those from sub-Antarctic South Georgia (SG) than to those from the MSE (Allegrucci et al., 2006). Within the Parochlus clade, MA, MSE and SG groups formed a monophyletic group together with previously published sequences of P. steinenii from SG (Cranston et al., 2010). This clade was recovered in an early branching event relative to the other available Parochlus species sequences from northern latitudes of South America (Peru, Chile, Argentina) and New Zealand ( Figure 2 ).

Global mitochondrial haplotype diversity (Hd) was 0.922, with no substantial variation among the three biogeographic regions analyzed ( Table 1 ). The MSE population exhibited the highest level of diversity, in terms of nucleotide diversity (Π = 7.44) and fewer haplotypes (K = 11), despite having a smaller sample size (n = 21). The sub-Antarctic SG and the Maritime Antarctic SSI populations showed similar genetic diversity patterns in terms of nucleotide diversity (Π = 3.04 and 1.2, respectively) and numbers of haplotypes (K = 27 and 28, respectively).

The median-joining haplotype network showed three main groups, without any shared haplotypes ( Figure 3 ). The MSE was the most distinct haplogroup, separated by 10 mutational steps from MA and SG. It also displayed a more expanded genealogy, with several mutational steps between each haplotype, contrasting with those in MA and SG, which displayed a star-like topology and short genealogy.

The Bayesian analysis of the number of genetic groups and the spatial clustering algorithm obtained with GENELAND detected three genetic clusters (K = 3, PP = 0.65, Figure 3 ), in agreement with the three biogeographic regions (i.e., MSE, SG and MA). The posterior probabilities associated with the definition of populations was on average 0.7, corresponding to the boundary between these clusters across the Polar Front (PF), which is located northern of SSI and SG ( Figure 4 ). Similarly, AMOVA of comparison of pairwise differences ( Supplementary Table S5 ), also confirmed the greatest differentiation between MSE and the combination of SG and SSI (70.41%, F_CT = 0.704, p < 0.01), while differentiation within groups only explained 2.72% (F_SC = 0.092, p < 0.01) of the total genetic variance. Pairwise comparison based on two statistics related to genetic (F_ST) and phylogeographic (Φ_ST) differentiation, showed significant structure across the three main biogeographic areas ( Supplementary Table S6 ). In particular, the levels of genetic differentiation within SSI and SG were very low (SSI: F_ST = 0.145, Φ_ST = 0.052; SG: F_ST = 0.106, Φ_ST = 0.199, for more details see Supplementary Table S7 ).

As expected for star-like genealogy, both neutrality tests Tajima’s D and Fu’s Fs were highly negative and significant ( Table 1 ) for MA and SG, but positive and not significant for the MSE. The distribution of pairwise differences analyzed through Mismatch Analyses (MMA), varied between the different demographic units of P. steinenii. The MA and SG showed a typical unimodal curve and smooth distribution ( Figure 5 ), representing a classical population expansion. In case of the SG populations, the curve moved to the right, indicating more mutation accumulation over time. As expected, the Harpending’s raggedness index was not significant in SG (rH = 0.027, p > 0.05) or MA (rH = 0.079, p > 0.05). Likewise, the analysis revealed a non-significant SSD value in SG (SSD = 0.001, p > 0.05) and SSI (SSD = 0.003, p > 0.05), supporting a spatial expansion for both populations In contrast, the MSE population exhibited a multimodal distribution ( Figure 5 ) and its Harpending’s raggedness index was significant (rH = 0.099, p < 0.05). Likewise, the analysis revealed a significant SSD value (SSD = 0.060, p < 0.05), refuting the demographic model of spatial expansion for the South American population and, rather, supporting population equilibrium. Based on Bayesian Inference phylogenetic reconstruction, haplotype network and the genetic differentiation index, we identified the three haplogroups aligning with the three specific geographic areas as separated and independent demographic units. Bayesian Skyline Plot analyses were conducted using the 10× corrected substitution rate ( Figure 5 ), providing historical population dynamic patterns and dates for TMRCA for the three demographic and evolutionary units (MA, SG and MSE). MSE had the oldest TMRCA of ca. 112,000 YBP (95% highest posterior density 95%HPD = 47.8 – 200 ka), followed by SG (ca. 50,000 YBP; 95%HPD: 19.6 – 92.1 ka), and finally MA (ca. 30,000 YBP; 95%HPD: 10.2 – 65.9 ka). The date of the onset of population expansion in SSI was estimated at ca. 10,000 YBP, while SG expanded earlier at ca. 20,000 YBP. For the MSE population, we did not detect a population expansion signal. Values of TMRCA were within the same range values obtained from the τ = 2μt values from the Mismatch Analyses (MSE: τ = 4.005, 36,035 YBP; SG: τ = 3.043, 13,500 YBP; MA: τ = 1.295, 11,655 YBP) ( Figure 5 ).

The haplotype phylogenetic reconstruction based on cox1 shows a clear assortment of haplotypes according to their geographic origin (MSE, SG and MA), as previously observed in the concatenated phylogenetic reconstruction ( Figure 1 ), and the haplotype network ( Figure 2 ). Using sequential calibration points, the divergence times estimated among the main clades of Parochlus suggest the MSE lineage as the first clade-group to diverge (Clade A), around 1.46 Mya (95%HPD): 0.87 – 2.10 Mya), followed by the divergence between MA and SG lineages (Clade B) around 0.84 Mya (95%HPD: 0.50 – 1.25 Mya). The age for each population node was obtained from the TMRCA in BSP analyses ( Figure 6 ).

Given the existence of two clades whose origin is during the Pliocene (1.46 Mya), the presence of private alleles in each population, the deeply diverged haplogroups/clades and distinct reciprocally monophyletic cluster, suggest the presence of separated and multiple refugia within the present distribution of the species. The diversification of P. steineni agrees with the start of the intensification of the glacial cycles during the Pliocene–Pleistocene (Zachos et al., 2001). These results support the long-term persistence of P. steinenii in sub-Antarctic and Maritime Antarctic regions, adding further supporting evidence to the recurring pattern of long-term persistence across a range of Antarctic terrestrial and freshwater invertebrates (see Convey et al., 2020 and citation therein; Verleyen et al., 2021; Maturana et al., 2022; Short et al., 2022). Though phylogeographic approaches can identify one of the outcomes of isolation from ancestral panmixia due to refugia (allopatry among populations), we must also consider factors that maintain regional distinction during interglacial periods, such as physical, environmental or selective mechanisms that may limit gene flow among populations (Maggs et al., 2008).

The representatives of Clade A were collected on Navarino Island in the Cape Horn Biosphere Reserve (CHBR), south of Tierra del Fuego. The lake sampled is located below Bandera Hill at 700 m a.s.l, and represents the only population sampled for genetic analyses within the MSE. The species has also been recorded at Horn Island, slightly further south in the Cape Horn archipelago, from a lake at similar altitude (see Gañan et al., 2020). The MSE population is characterized by high genetic diversity, an extended network topology and older TMRCA, with no signal of past demographic growth. Such genetic and demographic patterns suggest that this population maintained a significant effective population size during the intensification of glacial cycles throughout the Pleistocene, and has not experienced an interglacial population expansion unlike SG and MA. Additionally, the MSE population had a positive neutrality test, meaning a deficit of low frequency haplotypes and high divergence between them. Therefore, we proposed that P. steinenii is a polar/cold-adapted species with effective interglacial refugia at high altitude in MSE (see cryptic southern refugia section, Stewart et al., 2010), as described for other terrestrial examples in the Northern Hemisphere (Hirao et al., 2015). This plausible scenario of P. steinenii as a cold-adapted species allows us also to explain 1) the disjunct distribution and high percentage of divergence from representatives from SG and MA (which represent polar refugia at higher southern latitudes); 2) the persistence in high altitude ecosystems such as small rivers and marginal meltwater channels, similar to the habitats in which they are found today in their current MSE distribution range (Contador et al., 2023); and 3) the abundance and frequency of occurrences within South Shetland Islands, in particular Kitiesh Lake in King George Island where the population density is extremely high (T. Contador personal observations).

Conversely, the contemporary distribution of clade B included a broader allocation (SG, MA). In both populations, we found low levels of nucleotide diversity, unimodal curves in their Mismatch Analyses, and significant demographic expansion. Such a demographic pattern may reflect past contraction-expansion processes in response to climate oscillations, as has been frequently cited in phylogeographic studies of high-latitudes affected by the LGM (Provan and Bennett, 2008; Marko et al., 2010; Maturana et al., 2020). Considering that both SG and SSI are located south of the oceanic Antarctic Polar Front (APF), they would experience greater impact from ice-sheet advances during the LGM (Fraser et al., 2012; Hodgson et al., 2014). In consequence, these populations may have experienced strong bottlenecks due to drastic reductions in their population sizes. Further and ongoing research at small geographic scale, the use of more molecular markers (e.g. mitochondrial genomes), or more contemporary and polymorphic markers (e.g. SNPs), will be required to generate a possible scenarios with greater resolution for P. steinenii.

The archipelagoes of the Scotia Arc form a geological connection between the southern part of South America and the Antarctic Peninsula. The region is also one of Earth’s major ocean gateways and is critical for understanding the behavior of the main oceanographic currents and their influence on Southern Hemisphere biodiversity (Livermore et al., 2004; Scher and Martin, 2006). It is particularly complex in terms of its geology, glacial processes and the spatial distribution of biodiversity (Terauds et al., 2012). Studies based on available distribution records of vegetation (primarily bryophyte and lichen species) and certain groups of terrestrial invertebrates erected 16 distinct Antarctic Conservation Biogeographic Regions (ACBR, Terauds et al., 2012) in the region of Antarctic Treaty governance, three of which cover the western Antarctic Peninsula, South Shetland Islands and South Orkney Islands. This terrestrial biogeographic classification does not extend to the maritime Antarctic South Sandwich Islands, sub-Antarctic South Georgia, or the Magellan sub-Antarctic Ecoregion, contrary to the complex marine biogeographic regions, which have identified patterns and linkages across the Scotia Arc region (De Broyer and Koubbi, 2014; Koubbi et al., 2014; González-Wevar et al., 2022). Understanding the extent to which these biogeographical regions apply to any particular taxon, and how distributions may change given the appearance of new ice-free areas (Lee et al., 2022), are key to informing how biosecurity protocols should be implemented to minimize risks of anthropogenic transfer of native biota between different parts of this region (van Vuuren et al., 2018; Cukier et al., 2023) and in particular in freshwater ecosystems, which have not previously been included in such studies.

Based on our P. steinenii mtDNA sequences obtained from 16 locations across Maritime Antarctica (SSI) and different sub-Antarctic regions (SG, MSE), we detected strong genetic and phylogeographic differentiation combined with ancient population divergence. Our data support the deep divergence between Clade A and B, which may provide evidence of a speciation processes (genealogical lineages). This hypothesis of species differentiation within the overall distribution was also noted and discussed by Allegrucci et al. (2006). However, speciation processes must be considered carefully, since the nDNA did not detect any differentiation between both clades. The lack of any differentiation in the nuclear 28S rRNA gene versus the mtDNA is likely to reflect lower mutation rates in the nuclear gene and should be addressed in future studies by the inclusion of more nuclear genes.

Finally, our phylogeographic and genetic divergence estimates for P. steinenii are comparable to those reported in studies of other sympatric Maritime Antarctic invertebrates (McGaughran et al., 2010; van Vuuren et al., 2018; McGaughran et al., 2019), and the endemic Antarctic midge B. antarctica (Allegrucci et al., 2012; Edgington et al., 2023). Additionally, these divergence time estimates are consistent with those recently reported for the freshwater copepod, Boeckella poppei, based on cox1 sequences obtained across the same three biogeographic regions (Maturana et al., 2022).