Time or Space? Relative Importance of Geographic Distribution and Interannual Variation in Three Lineages of the Ascidian Pyura chilensis in the Southeast Pacific Coast

Spatial and temporal variation of environmental parameters can affect dispersal, recruitment and population persistence of marine benthic species. Studies including inter-annual comparisons of genetic structure often indicate high/moderate temporal heterogeneity in marine invertebrates, which may be a prevailing pattern. This suggests that temporal studies are necessary to understand the dynamics of marine metapopulations. In this study, we analyzed the spatio-temporal genetic structure of the ascidian Pyura chilensis, a low dispersal sessile marine species endemic from the Southeast Pacific coast and highly demanded for human consumption. We sequenced a fragment of the mitochondrial gene Cytochrome Oxidase I (COI) from 1,005 individuals of six locations (30–40 individuals per site and year) spanning a wide latitudinal range (24°–42°S) and sampled over 5 years (2012, 2014, 2015, 2016, and 2017). The genetic structure of COI indicates the presence of three monophyletic lineages (haplogroups 1–3) previously described for the species, being one of them highly divergent and geographically restricted (~39°S, Los Molinos). Considering the whole dataset, a picture of strong spatial differentiation but temporal stability emerged in Pyura chilensis. However, detailed studies of the two main lineages revealed important differences in the extent of spatio-temporal variation. Analyses using haplotype frequencies sorted by site and year showed that, for haplogroup 1, genetic variation was explained mainly by differences between sites, while for haplogroup 2 differences between years were prevailing. Haplogroup 3 was restricted to the most southern sites, and also showed inter-annual variability in its frequency. These results point to disparate patterns of genetic differentiation, which may reflect different adaptive scope or variation in reproductive and dispersal features and could be a response to extreme events such as El Niño (2015–2016). This work calls for caution when obtaining general trends in species clearly differentiated in lineages, and prompts instead for separate analyses of sub-specific genetic lineages whenever possible.


INTRODUCTION
Marine benthic coastal species inhabit an environment that is not only variable geographically but also temporarily (Menge et al., 2003). Heterogeneous environmental factors exert a strong influence on the dispersal, recruitment and population persistence of marine benthic species (Siegel et al., 2008;Marshall et al., 2010). Yet, how much does genetic structure vary in space vs. time in marine invertebrates is a topic rarely investigated (e.g., Calderón et al., 2012;Pineda et al., 2016a;Couvray and Coupé, 2018).
Phylogeographic studies of marine benthic invertebrates have traditionally focused on the spatial genetic structure using one or more molecular markers, revealing that most benthic marine species show genetic structure in their range of distribution, related to different abiotic and biotic features (Weersing and Toonen, 2009;Kelly and Palumbi, 2010;Schiebelhut and Dawson, 2018). Studies that include inter-annual comparisons indicate that temporal variation in the genetic structure of populations is a prevailing pattern in the marine environment (e.g., Virgilio and Abbiati, 2006;Barshis et al., 2011;Quintero-Galvis et al., 2020).
Temporal studies are necessary to better understand the dynamics of marine metapopulations and how genetic diversity is affected by reproductive patterns of local populations (Crowder et al., 2000;Cobben et al., 2011;Ludford et al., 2012). Many factors can promote or suppress local larval retention, thus the degree of reproductive connectivity between populations is variable and stochastic in nature (Johnson, 2004;Siegel et al., 2008), ranging from highly connected populations, to populations that show variable degrees of isolation. Differential settlement and survival of arriving larvae has been widely documented in marine organisms (e.g., Berger et al., 2006;Penney et al., 2006;Vigliola et al., 2007) and can also be strongly affected by oceanographic factors such as coastal circulation and upwelling (Roughgarden et al., 1988;Shanks and Brink, 2005;Barshis et al., 2011).
Studies on temporal genetic structure in ascidians are scarce, and they show contrasting results, ranging from yearly or seasonal fluctuations in genetic composition (Paz et al., 2003;Reem et al., 2012;Caputi et al., 2019) to instances of temporal stability over seasons and years (Pineda et al., 2016a,b). Here, we studied the degree of interannual genetic variation of six local populations of the ascidian Pyura chilensis in the Chilean coast between 24 and 42 • S, covering the southern portion of its geographic range of distribution (from 10 to 44 • S), with 5 years of sampling. This species is a conspicuous and dominant competitor of hard substrates in the Humboldt Current System (Ambler and Cañete, 1991;Valdivia et al., 2005) and a bioengineer species (Castilla et al., 2004) that is widely distributed in the Southeast Pacific coast (Lancellotti and Vásquez, 2000). It can be therefore considered a keystone species in the shallow sublittoral habitats of these areas. Besides its ecological importance, P. chilensis has been intensively exploited for human consumption (Manríquez and Castilla, 2005).
The spatial genetic structure of P. chilensis was studied by Haye and Muñoz-Herrera (2013) along the same geographic area, and three distinct mitochondrial lineages or haplogroups were found (hereafter Hg1-Hg3); two of them widely distributed (Hg1 and Hg2) and one geographically restricted to the southern area (Hg3; 39 • -42 • S). The frequency of each haplogroup showed strong geographic variation. On the northernmost location, only the most widespread haplogroup was present (Hg1). In intermediate localities (between 26 • S and 36 • S) two haplogroups were present simpatrically (Hg1 and Hg2). Finally, in the southern localities these two clades together with the third haplogroup (Hg3) were present. This southern haplogroup (Hg3) was particularly abundant in Los Molinos (39 • S). The nuclear gene Elongation Factor 1 alpha (EF1α) did not reveal the same structure in distinct groups as COI, and evidence was found of admixture between the two main haplogroups. Only some of the individuals of the more restricted COI haplogroup (Hg3) were divergent with EF1α. These results indicated that the three haplogroups are not reproductively isolated. The authors suggested that these haplogroups were formed during Pleistocene glacial cycles as a consequence of isolation during glaciations when sealevel fell. In the subsequent interglacial periods, populations could interbreed, resulting in the present day make-up of three divergent mitochondrial, yet not nuclear, lineages. In a follow up study, Segovia et al. (2017) analyzed P. chilensis populations on this same area using genome-wide SNP markers, detecting a strong geographical structure and a divergent nuclear lineage at Los Molinos (likely corresponding to haplogroup Hg3). However, in that work the COI haplotype composition of the studied individuals was not taken into account.
The present study tries to tackle the following main questions: (i) What is the relative importance of the temporal vs. the spatial components of genetic variation in P. chilensis, and (ii) Do the different lineages of P. chilensis have the same spatio-temporal pattern of genetic structure? To answer these questions, we obtained sequences of the COI gene from six local populations of P. chilensis at five time points (2012, 2014, 2015, 2016, and 2017). The data were used to assess the spatial-temporal patterns of genetic diversity in this ecologically and commercially important species.  Figure 1). Samples of yearlings were taken from the same area of each locality in each year. Sampling locations were selected based on results of Haye and Muñoz-Herrera (2013), and included populations with a priori different proportions of the three mitochondrial haplogroups of Pyura chilensis. Collections were performed in 2014, 2015, 2016, and 2017 which, together with data available from samples collected in 2012 (see Haye and Muñoz-Herrera, 2013), led to five time periods analyzed for temporal comparison of genetic structure. 1 | Genetic diversity measures for all COI data and for haplogroups 1 2, and 3 (Hg1, Hg2, and Hg3, respectively) of Pyura chilensis from five temporal samples (years 2012, 2014, 2015, 2016, and 2017).   Sampled individuals had a maximum size of 5 cm in order to avoid overlapping generations between samples of different years (Svane and Lundälv, 1982;Valdivia et al., 2005) and all individuals were obtained from clumps separated by at least 2 m. DNA extractions were performed from 25 mg of mantle tissue using the standard phenol-chloroform procedure. Partial sequences of mitochondrial COI gene were obtained using the universal primers HCO-LCO described in Folmer et al. (1994). For both primers, PCRs consisted of 1X PCR buffer, 1.3 Mm of MgCl 2 , 0.6 µM of each primer, 0.2 mM of each dNTP, 0.03 U µL −1 of Taq polymerase (Fermentas), 1.5 mg mL −1 Bovine Serum Albumin, and 20 ng of DNA template. Cycling conditions consisted of an initial soak of 10 min at 94 • C followed by 35 cycles, each of 1 min at 94 • C, 1 min at 40 • C and 2 min at 72 • C, and a final extension of 72 • C for 13 min. Resulting sequences were aligned and inspected in the software Geneious R8 1 .

Data Analyses
We determined for each year the frequency of haplogroups per location. Genetic diversity indices (haplotype and 1 www.geneious.com nucleotide diversity) were obtained with DnaSP 5.10 (Librado and Rozas, 2009).
To compare and examine the effect of both, spatial and temporal variation in genetic diversity, we used two approximations. First, AMOVA analyses were performed in Arlequin 3.5 (Excoffier and Lischer, 2010) using Years and Sites as a priori grouping, respectively, in order to evaluate the influence of both factors under a hierarchical structure (e.g., variation among sites within years and variation among years within sites). AMOVA was run using the F ST estimator of Weir and Cockerham (1984) based on allele frequencies. Second, we carried out a permutation-based multivariate analysis (PERMANOVA) using the Adonis function in the vegan package (Oksanen et al., 2020) in R (R Development Core Team, 2020). With this function, we decomposed the variance for distance matrices using years and sites as crossed factors. We used a Bray-Curtis distance matrix based on haplotype frequencies calculated using vegdist function and transformed using Hellinger Method with the decostand function, in vegan. Significant values were estimated after 10,000 permutations. Both analyses were carried out using the complete dataset and separately for haplogroups Hg1 and Hg2. As haplogroup 3 (Hg3) was spatially restricted to two FIGURE 1 | Temporal variation of COI haplogroups of Pyura chilensis. Geographic distribution and relative frequency of COI haplogroups for temporal samples (years 2012, 2014, 2015, 2016, and 2017). Pie charts represent the total number of sequences analyzed, colors represent COI haplogroups (red: Hg1, green: Hg2, blue: Hg3), and numbers by the pies indicate number of individuals analyzed per site per year.
localities (with most individuals found in Los Molinos), it was not considered for separate analyses.
Population-pairwise (within years) and temporal-pairwise (within sites) differentiation was assessed using F ST and Arlequin 3.5. The significance of the values obtained was tested with 10,000 permutations and a correction for multiple comparisons was applied following the false discovery rate (FDR) method described in Benjamini and Yekutieli (2001). As before, these analyses were performed for the complete dataset, and separately for the two main haplogroups (Hg1 and Hg2). To asses temporal stability of genetic differentiation, a congruence among distance matrices (CAMD) test (Legendre and Lapointe, 2004) was performed based on pairwise F ST distance matrices between localities for each year. For this analysis, we used the CADM.global and CADM.post functions implemented in the ape package (Paradis and Schliep, 2019) in R, which conducted a permutation test (10,000 permutations) of the congruence of each distance matrix with respect to all the other matrices and then an extended Mantel test with a post hoc correction of Holm were conducted. Mantel values go from 0 to 1, with higher values meaning higher concordance between matrices and, thus, temporal stability. The null hypothesis was complete incongruence among matrices. The alternative hypothesis is that the matrix (year) tested is congruent with at least one other matrix of another year.
Finally, haplotype networks were constructed with the software PopArt v1.7 2 using a Median-Joining algorithm for each sampling site including all temporal samples.

RESULTS
COI sequences were obtained for a total of 1,005 individuals of Pyura chilensis from six localities (Figure 1) and five time periods (2012, 2014, 2015, 2016, and 2017). The three COI haplogroups that had been described in a previous study of P. chilensis (Haye and Muñoz-Herrera, 2013) were found to have a temporally stable geographic distribution (Figure 1). Haplogroup 1 (Hg1) was found on all sampled locations in all years, and had the widest distributional range of all the haplogroups. Haplogroup 2 (Hg2) was found from Caleta Pajonales (CP) to the south, and was absent only in the northern location Pan de Azúcar (PA, 26 • S) in all sampled years with the exception of a single individual in 2016 (Figure 1). As expected, Hg3 was the most restricted haplogroup, present only in the two southern locations, mainly in Los Molinos (LM) and at a much lower frequency in Ancud (AC). Haplogroup frequencies per location showed slight variation between years, with the exception of Los Molinos that showed strong variation between some years (LM, Figure 1). Including all years, a total of 386, 111, and 77 haplotypes were found for Hg1, Hg2, and Hg3, respectively. High haplotype diversity was found at all localities and haplogroups, while nucleotide diversity was low when considering each haplogroup separately ( Table 1).
The relative genetic variation explained by sites and years was analyzed with AMOVAs (Table 2 and Figure 2). When comparing the results of years or sites as a priori groupings, 2 http://popart.otago.ac.nz  Frontiers in Marine Science | www.frontiersin.org we found for the complete dataset that the main factor "among sites" explained a significant amount of the variance (F CT = 0.0030, p = 0.012), while sampling year (among years) did not (F CT = 0.0021, p = 0.116). The same pattern was found for the most frequent haplogroup (Hg1), where 2.52% of the total variance was explained by differences among sites (F CT = 0.0252, p < 0.001) while the temporal component had a negative value (F CT = −0.0022, p = 0.555). Interestingly, this pattern was reversed for Hg2: among year differences explained 4.31% of the total variance and was significant (F CT = 0.0431, p = 0.004), while sites as a priori grouping had a negative percentage of explained variance and non-significant (p = 0.567). When sites or years were nested within each other, more subtle patterns emerged ( Table 2). Even if year is not a significant main factor in the complete dataset, the variation among years within sites explains a low but significant amount of variation (F SC = 0.0082, p < 0.0001), indicating that, at least for some localities, temporal variation is significant. This same pattern is found for Hg1. Similarly, for Hg2, while sites as a main factor was not significant, variation between sites within years is (F SC = 0.0310, p = 0.004), indicating a significant spatial component at least in some years. In all cases, most genetic variance was found in the within population component (all percentages of variation explained >92%). In the same page, permutation-based analysis of variance (PERMANOVA) using sites and years as crossed factors showed that sampling site was the most important factor for the complete dataset (Fvalue = 1.6996, p < 0.001; Table 3) and for the most abundant haplogroup (Hg1) (F-value = 1.5507, p < 0.001). Differences between years, however, were also significant in both cases. For Hg2, on the other hand, this analysis showed that only differences between years were significant (F-value = 2.0873, p < 0.001). Pairwise F ST matrices and heatmaps using haplotype frequencies sorted by site and year showed that, for the complete dataset and Hg1, genetic variation was explained mainly by differences between sites (Figures 2A-D and Table 4) and, particularly for Hg1, most of the significant outcomes involved the comparison between southern locations (TH, LM, and AC). Conversely, for Hg2 non-significant values were found between sites ( Table 4) and most variation was explained by differences between years (Figures 2E,F and Table 5), with significant values generally involving comparisons of 2015 and 2016 with other years. These same years showed also significant comparisons for the global dataset and for Hg1, but the F ST values were an order of magnitude higher for Hg2 (Table 5). Likewise, in the congruence among distance matrices analysis (CADM) for the complete dataset, the CADM.global function showed a high and significant Kendall coefficient of concordance (W = 0.76, χ 2 = 53.47, p = 0.001), indicating overall congruence between distance matrices and, thus, temporal stability. Additionally, a posteriori Mantel correlation values (CADM.post) between each year and the rest of years ranged from 0.63-0.73 being significant for all years ( Table 6). For Hg1, the global test also showed a high and significant Kendall coefficient (W = 0.71, χ 2 = 50.23, p = 0.001). The a posteriori Mantel correlation values were all >0.6, but were only significant for 2012, 2017 and marginally significant for 2014 ( Table 6), indicating that these were the years with higher temporal stability. Interestingly for Hg2, there was a low and non-significant Kendall coefficient (W = 0.19, χ 2 = 8.57, p = 0.48), and Mantel correlation mean values ranged from −0.37 to 0.17, being none of them significant ( Table 6). This means that, for Hg2, there is a complete incongruence of    F ST pairwise distance matrices among years, indicating a high temporal variability for this haplogroup. Haplotype networks for each sampling site allowed the evaluation of the relationship between haplotypes and the temporal stability within localities (Figure 3). The networks showed that haplotypes sampled in 2012, 2014, 2017 were similar while those from 2015 and 2016 were often different (Figure 3), forming in general clusters that were separated from central haplotypes of each haplogroup. Similarly, haplotype sharing was common between the years 2012, 2014, 2017, while 2015 and 2016 haplotypes were more commonly shared only among these 2 years (Figure 3).

DISCUSSION
The genetic structure of P. chilensis consisted of three highly differentiated COI haplogroups. Haplogroup 1 (Hg1), the most frequent and present throughout the analyzed geographic area, Haplogroup 2 (Hg2), only missing on the northernmost site (Pan de Azúcar) and Haplogroup 3 (Hg3), the most differentiated and restricted to the two southern sites (mainly in Los Molinos), as shown previously by Haye and Muñoz-Herrera (2013). The proportion of haplogroups of P. chilensis did not vary between the 5 years sampled, with the exception of Los Molinos, where the frequencies of each Haplogroup, particularly Hg3, were more variable, especially in years 2015-2016.
A noteworthy result is the very high haplotype diversity found (in many cases each individual in a sample had a different haplotype). This is in agreement with Haye and Muñoz-Herrera (2013) results and is attributed by these authors to large effective population size and fast mutational rate. It also suggests that the populations are experimenting a demographic expansion in the sampling area. High haplotype diversity in COI, but less than in the present work, has also been found in other solitary ascidians, such as Microcosmus squamiger (Rius et al., 2009; overall Hd = 0.794), or Styela plicata (Pineda et al., 2011;overall Hd = 0.810).
When all the data were considered (all haplogroups together), the overall variability in genetic structure of Pyura chilensis showed a greater spatial than temporal component over the time frame studied (5 years). Despite this general trend, a striking result of this study is that the relative importance of spatial and temporal variation varies between the mitochondrial haplogroups of P. chilensis. The most frequent and spatially extended Hg1, showed a higher spatial than temporal variation between the 5 analyzed years. As in Hg1 most of the variance was explained by differences between sites and since this is the most abundant haplogroup, the trend of Hg1 may bias the general analysis, obscuring the temporal component.
Conversely, for Hg2, which was absent in the northernmost location (PA), most of the variance was explained by differences between years. A temporal genetic variation was also detected in the Hg3 if we consider its frequency in the site LM, which harbors most individuals of this geographically restricted COI haplogroup: the relative frequency of this group ranged from ca 0.2 to 0.9 in the years surveyed. We could not analyze, however, changes in haplotype frequencies due to limited data for this haplogroup. Albeit the spatial or temporal components are the dominant, respectively, in Hg1 and Hg2, there is nevertheless evidence for significant differences in haplotype composition in particular years in Hg1, as evidenced by PERMANOVA and pairwise F ST results. For Hg2, AMOVA showed a significant site within years component, suggesting that spatial differentiation also exists in at least some years.
The existence of lineage specific differences in temporal and spatial components of variation was unexpected. Particularly because, using nuclear sequences of the Elongation Factor 1 alpha (EF1α) gene, Haye and Muñoz-Herrera (2013) showed that there is admixture between Hg1 and Hg2. Other solitary ascidians also show a lack of relationship between mitochondrial and nuclear genetic structure. For instance, in Styela plicata there are two groups based on COI haplotypes and two groups based on nuclear ANT sequences, but these groups are not congruent FIGURE 3 | Haplotype networks per location of COI haplotypes of Pyura chilensis. Each circle represents a sampled haplotype, and the size is proportional to the number of individuals with the haplotype. Colors represent temporal samples (years 2012, 2014, 2015, 2016, and 2017). Lines between haplotypes represent a mutational step, and small crossed lines correspond to additional steps. Hypothetic/non-sampled haplotypes are represented with a black square. The three COI haplogroups (Hg1, Hg2, and Hg3) are shown in a single network per site and numbers represent the number of mutational steps between haplogroups. (Pineda et al., 2011), which can be explained by the different inheritance mechanisms of both genes. Nuclear genes can be shuffled by sexual reproduction leading to a different distribution than mitochondrial lineages. In addition, Pineda et al. (2011) found an excess of individuals bearing one allele of each nuclear group, indicating that selection may have favored admixture in this species. In the solitary ascidian Microcosmus squamiger there are also two main mitochondrial groups (Rius et al., 2009) and two nuclear clusters as defined by microsatellites (Rius and Teske, 2011). However, these two groups are also not congruent (Ordóñez et al., 2013). Interestingly, in this species it was found that nuclear lineages tend to remain separated even in sympatry in long-term introduced populations (Ordóñez et al., 2013), thus in this case interbreeding seems to be selected against.
We cannot explain at present what are the causes of the different spatio-temporal dynamics of the two haplogroups of P. chilensis. One option is that mitochondrial genes may determine reproductive or biological parameters that affect the dispersal capabilities of the species and, thus, the temporal variation. Alternatively, the differentiation between mitochondrial lineages can correspond to some degree to nuclear differentiation as well, and this may not be evident when evaluating a single gene (EF1α in Haye and Muñoz-Herrera, 2013). We reanalyzed the genome-wide SNP markers genotyped in a subset of the individuals here studied (data from Segovia et al., 2017Segovia et al., , 2020, taking into account the mitochondrial information of the individuals. Our results (Supplementary File 1) revealed that the admixture between Hg1 and Hg2 is not perfect. Thus, there is some degree of nuclear differentiation between them, related both to neutral and to potentially selective processes. The processes responsible for the maintenance of imperfect admixture (either pre-zygotic or post-zygotic mechanisms) between haplogroups cannot be assessed at present and would require an in-depth knowledge of the biological features of each haplogroup and experimental crosses between them. Nevertheless, this nuclear differentiation opens a much wider scope of genetic differentiation than if only mitochondrial DNA could account for the differences observed. Selective or non-selective processes with implications in dispersal abilities can therefore translate into differences in the nuclear genetic makeup of the haplogroups that can explain their different spatio-temporal behavior. One such process may be the onset of the reproductive season (for instance, one lineage may reproduce at a time when currents are stronger). Another potential process could be differing thermal tolerance, that is relevant for larval survival (Perez-Valdes et al., 2017).
The general trend of greater spatial than temporal variation, with an overall low to moderate temporal variation, was also found in Styela plicata, an ascidian of similar biological features (Pineda et al., 2016b). This pattern has also been described in several other marine invertebrates such as polychaetes (Kesäniemi et al., 2014); molluscs (Kenchington et al., 2006;Cassista and Hart, 2007;Lee and Boulding, 2009;Ng et al., 2010), crustaceans (Ball and Chapman, 2003;Pascual et al., 2016), and urchins (Flowers et al., 2002;Calderón et al., 2009). On the contrary, a greater temporal than spatial variation similar to what we detected in Hg2 (and, although limited sampling allowed only to evaluate changes in relative frequencies, in Hg3), has been found in some marine invertebrates such as the barnacle Balanus glandula (Barshis et al., 2011), the oyster Ostrea edulis (Hedgecock et al., 2007), periwinkles Littorina spp. with planktonic larval development (Lee and Boulding, 2009), and the anomuran crab Petrolisthes cinctipes (Toonen and Grosberg, 2011).
A strong degree of spatial genetic differentiation is to be expected in marine invertebrates with low natural dispersal potential. This is a common result in population genetic analyses of ascidians (e.g., López-Legentil and Turon, 2006;Pérez-Portela and Turon, 2008), as they have lecitotrophic, short-lived larvae. In addition, temporal genetic variation may also be linked to the duration of the planktonic larval life stages. For example, species of the genus Littorina vary in their degree of temporal variation in accordance with their dispersal potential; species with planktonic larval dispersal showed high temporal variation while species with direct development showed low temporal variation (Lee and Boulding, 2009). In this sense, the finding of a strong temporal signature in Hg2 is somewhat unexpected and suggests strong temporal variation in response to environmental parameters.
Temporal genetic variation can be established by different processes (Calderón et al., 2009). In organisms with high fecundity and high larval mortality many factors affect reproductive success, including biotic variables (adult densities, sperm availability, larval traits) and abiotic parameters (current patterns, temperature changes, upwellings). These factors can determine selective forces with a variable effect on larvae (Watts et al., 1990). Alternatively, the sweepstakes reproductive success hypothesis (Hedgecock, 1994;Hedgecock and Pudovkin, 2011) establishes that stochastic changes in the environment can result in high variance in reproductive success, thereby different adults contribute to the pool of successful recruits in each generation. High variance in reproductive success can explain temporal variability in genetic diversity and structure (Turner et al., 2002;Pujolar et al., 2007;Calderón et al., 2009;Kesäniemi et al., 2014), this reproductive variance potentially leading to temporal chaotic genetic patchiness (Colgan, 1981;Johnson and Black, 1984;Hedgecock, 1994).
In our case, Hg2 had a significant overall temporal variation, but it is particularly strong when comparing the years 2015-2016 with the rest. For Hg1, the significant F ST comparisons were also related to these 2 years (albeit with values one order of magnitude lower than in Hg2). Consistently, haplotype networks evidenced that in 2015/2016 a different suite of haplotypes was encountered in all haplogroups. Environmental conditions along the southeast Pacific coast are highly unpredictable as a result of processes that have temporal variation such as seasonal winds and upwelling, and the strong influence of the El Niño Southern Oscillation (ENSO), among others (Thiel et al., 2007). In this context, significant changes in haplotype frequencies in 2015/16 are coincident with one of the strongest recorded El Niño events (Blunden and Arndt, 2016) comparable only to the events of 1982/1983 and 1997/1998 in terms of intensity and anomalies in temperature (McPhaden and Zhang, 2009;Santoso et al., 2013) and salinity (Corbett et al., 2017). A recent study suggests that P. chilensis is sensitive to temperature changes, being one of the environmental variables that affects its adaptive genetic structure (Segovia et al., 2020). Sea Surface Temperature anomalies during 2015/2016 ENSO event reached˜3.5 • C (Paek et al., 2017), which suggests that temporal stability in haplotype frequencies was significantly affected by these anomalies. Changes in upwelling intensity (Aravena et al., 2014;Pinochet et al., 2019) may also influence temporal patterns of genetic structure in P. chilensis, as reported for other species such as Balanus glandula in the northwest Atlantic coast (Barshis et al., 2011).
In short, our results suggest different mitochondrial lineage-specific patterns of spatio-temporal genetic variation in P. chilensis. The temporal dynamic that has been found is in accordance with the high temporal variability in environmental conditions in the Southeast Pacific. We also confirmed previous reports of spatial genetic differentiation of P. chilensis along its latitudinal range of distribution, reflected in different frequencies of haplogroups and, at least within the main haplogroup, a strong spatial signature with temporal stability. Our results call for caution when general trends are applied to species with highly differentiated lineages. We prompt instead for separate analyses of sub-specific lineages whenever possible, as they can prove a rich source of information otherwise untapped. This has been noted for the analyses of spatial genetic structure (i.e., Carreras et al., 2020), but a lineage-specific pattern of temporal variation has not, to our knowledge, been previously reported for any marine invertebrate.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Genbank accession numbers of DNA sequences of haplotypes per year: MW785940-MW786618.

AUTHOR CONTRIBUTIONS
PH and NS conceived the idea and contribuited to sample collection, laboratory work, and acquisition of data. PH, XT, and NS contributed reagents, materials, analysis tools, conceived and designed the analysis, and wrote the manuscript.

FUNDING
This study was funded by the Chilean National Fund for Scientific and Technological Development through Grants FONDECYT 1140862, FONDECYT 3190482, Millennium Science Initiative Program-ICN2019_015, and partly funded by project PopCOmics CTM2017-88080 (MCIU/AEI/FEDER) of the Spanish Government.