Chaotic Genetic Patchiness in the Highly Valued Atlantic Stalked Barnacle Pollicipes pollicipes From the Iberian Peninsula: Implications for Fisheries Management

The stalked barnacle Pollicipes pollicipes inhabits rocky shores from the Atlantic coasts of Brittany (France) to Senegal. Because of the culinary traditions of southern Europe, stalked barnacles represent an important target species for local fisheries on the Iberian Peninsula. To manage this fishery sustainably, it is therefore important to assess the dynamics of local populations over the Iberian coast, and how they are interconnected at a wider scale using finely tuned genetic markers. In this work, a new enriched library of GT microsatellites for P. pollicipes was prepared and sequenced using Ion Torrent™ Next Gen-Sequencing Technology. 1,423 adults and juveniles were sampled in 15 localities of three geographic regions: southern Portugal, Galicia and Asturias (both in northern Spain). Twenty polymorphic loci arranged in five multiplex PCRs were then tested and validated as new molecular tools to address the spatial and temporal genetic patterns of P. pollicipes. Our results revealed high genetic diversity among adults. However, juveniles were genetically more structured than their adult counterparts, which alternatively displayed much more connectivity among the three studied regions. The lack of spatial genetic heterogeneity in adults may be due to the overlapping of several generations of settlers coming from different geographic origins, which mainly depends on the orientation of residual currents along the coast during reproduction. The genetic differentiation of juveniles may indeed be congruent with Iberian Peninsula hydrodynamics, which can produce chaotic genetic patchiness (CGP) at small temporal scales due to sweepstake reproductive success, collective dispersal and/or self-recruitment. Remarkably, most of the genetic heterogeneity of juveniles found in this work was located in Galicia, which could represent an admixture between distinct metapopulations or an old refuge for the most northern populations. To conclude, high genetic variation in P. pollicipes can lead to the false impression of population panmixia at the Iberian scale by masking more restricted and current-driven larval exchanges between regions. This possibility should be taken into consideration for further specific management and conservation plans for the species over the Iberian Peninsula.

The stalked barnacle Pollicipes pollicipes inhabits rocky shores from the Atlantic coasts of Brittany (France) to Senegal. Because of the culinary traditions of southern Europe, stalked barnacles represent an important target species for local fisheries on the Iberian Peninsula. To manage this fishery sustainably, it is therefore important to assess the dynamics of local populations over the Iberian coast, and how they are interconnected at a wider scale using finely tuned genetic markers. In this work, a new enriched library of GT microsatellites for P. pollicipes was prepared and sequenced using Ion Torrent TM Next Gen-Sequencing Technology. 1,423 adults and juveniles were sampled in 15 localities of three geographic regions: southern Portugal, Galicia and Asturias (both in northern Spain). Twenty polymorphic loci arranged in five multiplex PCRs were then tested and validated as new molecular tools to address the spatial and temporal genetic patterns of P. pollicipes. Our results revealed high genetic diversity among adults. However, juveniles were genetically more structured than their adult counterparts, which alternatively displayed much more connectivity among the three studied regions. The lack of spatial genetic heterogeneity in adults may be due to the overlapping of several generations of settlers coming from different geographic origins, which mainly depends on the orientation of residual currents along the coast during reproduction. The genetic differentiation of juveniles may indeed be congruent with Iberian Peninsula hydrodynamics, which can produce chaotic genetic patchiness (CGP) at small temporal scales due to sweepstake reproductive success, collective

INTRODUCTION
The percentage of stocks exploited at biologically unsustainable levels increased from 10% in 197410% in to 34.2% in 201710% in (FAO, 2020 after decades of management strategies based on catchrate limitations (i.e., the EU Common Fisheries Policy). As an alternative or complementary approach, management practices are increasingly incorporating the spatial allocation of fishing intensity through marine protected areas, marine zoning, or spatial user rights, particularly for sessile or low-motility species (Lorenzen et al., 2010;Rassweiler et al., 2012). Optimization of these processes depends on the accurate estimation of the connectivity among management units, mediated by the dispersal of the planktonic larval stages (Silva et al., 2019). In this regard, a fundamental issue concerns whether the dispersal scales are consistent with the management scales (Ouréns et al., 2015). Although advection by ocean currents should lead to long dispersal distances exceeding the scale of management, there is increasing evidence that long-distance dispersal may be rare on ecological time scales (Palumbi, 2003;Selkoe et al., 2010;D'Aloia et al., 2015). This phenomenon can be explained by a combination of seascape characteristics such as eddies, gyres or upwellings of deep water bodies and specific larval behavior that would favor local retention and reduced dispersal (Morgan et al., 2009(Morgan et al., , 2018Barshis et al., 2011;Kough et al., 2013). An additional line of evidence reveals surprising patterns of spatial and temporal genetic structure observed in some marine species at a scale where genetic variation should be efficiently homogenized by gene flow via larval dispersal, collectively coined chaotic genetic patchiness (CGP) (Johnson and Black, 1982;Hedgecock and Pudovkin, 2011;Eldon et al., 2016).
The stalked barnacle (Pollicipes pollicipes) is a pollicipedomorph cirriped (Chan et al., 2021) inhabiting rocky coasts that are highly exposed to waves in the northeast Atlantic. Its range extends from southwestern England through the coasts of Brittany (France), Spain, Portugal, and West Africa to Dakar (Senegal) (Barnes, 1996(Barnes, , 2008Southward, 2008;Fernandes et al., 2010). In the Iberian Peninsula, it represents a highly valued resource that reaches very high market prices due to an old gastronomic tradition (Molares and Freire, 2003;Jacinto et al., 2011;Rivera et al., 2014). Remains of its consumption Abbreviations: CGP, chaotic genetic patchiness; TURF, territorial use rights for fishing.
have been found in early Holocene archeological sites, mainly associated with Mesolithic and Neolithic shell-middens on both the Atlantic and Mediterranean coasts (Álvarez-Fernández et al., 2010(Álvarez-Fernández et al., , 2013Álvarez-Fernández, 2011). Between 2013 and 2016, the European stalked barnacle fisheries have an annual economic value of EUR 10 million, involving approximately 500 t of landings and 2,100 professional fishers (Aguión et al., 2021). At some localities, the pressure exerted by poachers can be extremely high (more than 60% of their catches), especially in banned areas or periods (Jacinto et al., 2010;Rivera et al., 2014;Ruiz-Díaz et al., 2020).
Management of the stalked barnacle fishery in the Iberian Peninsula is highly heterogeneous (Aguión et al., 2021). In Galicia (NW Spain) since 1992, the regional government has developed a co-management system between fishers' guilds ("cofradías") and the fisheries authority through territorial user rights for fishing (TURFs) (Molares and Freire, 2003;Macho et al., 2013), where exclusive right of access are granted to fishing communities (Costello et al., 2010;Rivera et al., 2014). Similarly, in the West coast of Asturias (N Spain), the barnacle fishery has been managed through a co-management system with TURFs since 1994 (Rivera et al., 2014(Rivera et al., , 2017. Both Galicia and western Asturias present adaptive spatial management with nested scales at regional, local and patch/rock levels; recognized to promote fisheries sustainability (Aguión et al., 2021). However, on the eastern coast of Cape Peñes (eastern Asturias) and Portugal, the fishery is managed at a regional scale through general regulations without management plans (Aguión et al., 2021). In Portugal, however, there are two protected areas subjected to specific regulations for harvesting P. pollicipes: the Reserva Natural das Berlengas (RNB) and the Parque Natural do Sudoeste Alentejano e Costa Vicentina (PNSACV) (Sousa et al., 2013;Cruz et al., 2015;Carvalho et al., 2017). The first one (RNB) is subjected to local management, resembling a TURF in many aspects (Aguión et al., 2021). Currently, there is interest and potential for developing co-management systems similar to the one in Galicia and western Asturias in both Portuguese protected areas Sousa et al., 2020). Among different management approaches, TURFs represent the best option for the sustainable management of small-scale sessile fisheries (Gutiérrez et al., 2011;Rivera et al., 2017;Aguión et al., 2021). However, the design of management areas mandates a good understanding of population renewals for which estimates of connectivity are crucial (Aceves-Bueno et al., 2017;Silva et al., 2019). Dispersal, settlement, and subsequent recruitment are decisive processes in the population dynamics of marine invertebrates with planktonic larval stages, allowing the connection between remote populations and leading to meta-populations that are globally viable (Cowen and Sponaugle, 2009).
Pollicipes pollicipes larvae go through six planktotrophic nauplius stages before turning into a lecithotrophic stage, called cypris. According to Molares et al. (1994) and Franco et al. (2016Franco et al. ( , 2017, the pelagic larval development is finalized after 15 days to 1 month under optimal conditions in the laboratory, whereas in the natural environment, the total pelagic larval duration is estimated to last 2 months (Cruz, 2000;Macho, 2006). The presence of stalked barnacles on the shore might favor the settlement of cyprids, because recruitment is intense on conspecifics (e.g., Cruz et al., 2010;Fernandes et al., 2021). The dynamics of ocean circulation are recognized as important aspects in shaping connectivity patterns among marine populations (Treml et al., 2008). In this situation, significant effort is required to study population dynamics locally to adequately manage the resource (Molares and Freire, 2003). For P. pollicipes, a minimum potential passive migration distance of 600 km during the planktonic stage has been suggested (Quinteiro et al., 2007); nevertheless, reanalysis of genetic data and basic biophysical modeling point to modest dispersal distances in the range of tens of kilometers in the Asturian region (Rivera et al., 2013). At a large spatial scale, it has been suggested that P. pollicipes displays a metapopulation structure, where disconnected adult populations share a common larval pool (the n-islands model hypothesis) (Molares and Freire, 2003). However, the metapopulation structure has not yet been addressed. Alternatively, species with long larval dispersal potential, such as P. pollicipes, may exhibit surprising patterns of spatial and temporal genetic structure. CGP (Johnson and Black, 1982;Hedgecock and Pudovkin, 2011;Eldon et al., 2016) has been consistently reported in marine species that broadcast larvae at a scale where genetic variation should be efficiently homogenized by gene flow via larval dispersal. Eldon et al. (2016) reviewed and discussed how selection, sweepstake reproductive success, collective dispersal, and temporal shifts in local population dynamics may play a crucial role in generating such unexpected patterns. Moreover, Pineda et al. (2006) reported the existence of "recruitment windows" in a close barnacle species (Semibalanus balanoides), in which after a recruitment period of approximately 3 months, only recruits able to settle in just a couple of weeks survive after settlement and mature into adults. In spite of its interest for the management of this species, the processes that shape the genetic structure of P. pollicipes in the Atlantic Ocean have not been studied.
Genetic markers are a powerful tool for fisheries management because they present an array of very useful applications: they can address the correct identification of species, delimit distinct fish stocks (Borrell et al., 2012;Papa et al., 2020), assess relatedness levels within populations (Veliz et al., 2006;, expose population connectivity (Pascual et al., 2017;Muñoz-Ramírez et al., 2020), estimate larval dispersal (Van Wyngaarden et al., 2017) and larval diversity (Chen et al., 2013;Wong et al., 2014;Alshari et al., 2021) or the source-sink dynamics within the population structure (Pineda et al., 2007;Brault et al., 2013;Lindegren et al., 2014). Genetic data, however, integrate information on the past demographic history of populations and are not always easily applicable for the present-days management for marine species with high fecundity and dispersal capabilities (Gagnaire et al., 2015). Estimating some of the population parameters that are crucial for stock management imposes the need to develop numerous highly polymorphic markers. These will help to discriminate between past and present-day processes that shape populations of species with highly effective population sizes (e.g., Hongjamrassilp et al., 2020). Despite the economic relevance of the P. pollicipes fishery, only a few articles have been published on the genetics of the stalked barnacles, most of which are based on mitochondrial markers (Quinteiro et al., 2007;Campo et al., 2010;Rivera et al., 2013). According to Quinteiro et al. (2007), the panmixia hypothesis is rejected, and 5 population groups are established: (1) Brittany; (2) Asturias-East; (3) Galicia, Portugal and Morocco; (4) Canary Islands; and (5) Cape Verde Islands, with the latter being extremely divergent. The Cape Verde population was later considered a new species  and described as Pollicipes caboverdensis . Campo et al. (2010) revealed genetic differences among populations between Brittany (France) and the rest of the species distribution range, while Rivera et al. (2013) described small-scale, asymmetric connectivity in gooseneck barnacle populations, when reanalyzing data from Campo et al. (2010) for the Cantabrian coast. Microsatellites usually display high levels of genetic variation and can detect subtle genetic differentiation among populations separated by only a few hundred kilometers (Borrell et al., 2012). Moreover, they seem to be very useful to detect parentage/familial structures, when assessing the origin of recruits (St-Onge et al., 2015;Couvray and Coupé, 2018;Dubé et al., 2020). Microsatellite markers have been previously developed and, in some cases, used to infer the population genetic structure for several closely related acorn barnacles, such as S. balanoides (Dufresne et al., 1999;Flight et al., 2012); Chthamalus montagui (Pannacciulli et al., 2005;Fontani, 2009); Tetraclita spp. (Dawson et al., 2010;Chen et al., 2015); Megabalanus coccopoma (Reigel et al., 2015); Chelonibia testudinaria (Ewers-Saucedo et al., 2016; Notochthamalus scabrosus (Barahona et al., 2019) and two stalked barnacle species: Pollicipes elegans  and P. pollicipes (Seoane-Miraz et al., 2015;Fernandes et al., in prep.). The latter appear to have shown positive results with specific crossamplifications in the congeners P. elegans, Pollicipes polymerus, and P. caboverdensis.
The aim of the present study was to revisit and test the previously described genetic homogeneity of P. pollicipes at the scale of the Iberian Peninsula with highly polymorphic microsatellite markers. The final goal is to provide support for the design of adequate and sustainable fishery management plans, using an in-depth analysis of genetic patterns inferred from a hierarchical geographic sampling of the barnacle populations along the Iberian coastline. However, preliminary tests using published microsatellite markers have provided inconsistent and non-reproducible PCR results in two different and independent genetic labs, necessitating the development of new highly variable genetic markers for the species P. pollicipes (this study).

Study Area and Sampling
A total of 1,423 individuals from 15 different localities belonging to three Atlantic regions of the Iberian Peninsula covering the most important spots of the barnacle fishery were sampled. These three regions were SW Portugal, Galicia (NW Spain), and W Asturias (N Spain). Thus, a hierarchical sampling of populations was performed in which five distinct localities were sampled within each region (Figure 1). The five localities belonging to Portugal are Aljezur (AL), Azenha do Mar (AZ), Cabo Sardão (CO), Malhão (MA) and Sines (CS). The five localities belonging to Galicia are Baiona (BA), Cabo Home (CH), Aguiño (AG), Camelle (CA) and A Coruña (AC). The five localities belonging to Asturias are La Cruz (PC), El Cuerno (CU), Las Llanas (LM), La Erbosa (ER), and El Corviru (EC) (Figure 1). Samples were transferred to the laboratory and frozen on the same day of collection until further individualization and labeling.
Within each of the targeted localities, one hundred individuals collected in September and October 2017 were randomly sampled within two distinct developmental cohorts according to their rostrum-carinal (RC) length [see Figure 4 in Parada et al. (2013)] (50 adults of commercial size greater than >18 mm; 50 juveniles between 2 and 4 mm). As barnacles are usually found in groups of sessile individuals, fixed on primary rocky substrates with small juveniles attached to adult peduncles , juveniles were first removed from the adults, avoiding the collection of more than one juvenile by adult and then treated secondarily. Each barnacle was put individually in a tube previously labeled and preserved in absolute ethanol at room temperature. In the laboratory, a small portion of the peduncle muscle was dissected from each individual for genomic DNA extraction. In the case of adults, special care was taken to dissect the tissue from the inner part of the peduncle to avoid possible contamination by attached post-larvae (cyprids) and juveniles.
To characterize the typical upwelling circulation during the stalked barnacle larval season in summer/autumn 2017 along the coasts of northern and western Iberia, sea surface temperature (SST) along with modeled sea surface currents datasets were retrieved during the peak of meridional Ekman transport at central Portugal on 11-08-2017. Daily 4 km SST data were obtained from the all-satellites combined Copernicus' product 1 . The 5-days averaged meridional and zonal components of the surface currents were obtained from the OSCAR model with a spatial resolution 2 of 0.33 • .

Microsatellite Markers and Multiplex PCR Development
Genomic DNA from five adult individuals was extracted using the EZNA R Mollusk Kit (Omega Bio-Tek Inc., Norcross, GA, United States). An enriched biotin-labeled CT/GT library for dinucleotides was obtained using the methodology described by Bloor et al. (2001) and Sotelo et al. (2007), where DNA was digested with HaeIII (NEB). Digestions were run in 1.5% agarose gels stained with ethidium bromide. Fragments between 400 and 800 bp were excised from gels and purified using a QIAquick Gel Extraction Kit (Qiagen). Fragments were ligated to a doublestranded adaptor using ligase (NEB) and enriched by PCR using oligoA. Purified PCR products were denatured and incubated with 200 pmol of 5 biotinylated (CT) 12 and (GT) 12 probes (Invitrogen) attached to streptavidin-coated magnetic beads (Streptavidin MagneSphere Paramagnetic Particles, Promega). Hybridization was carried out in 6 SSC for 30 min at 60 • C in a thermocycler. Specific fragments were recovered after washing the bead suspension with solutions progressively desalted at 60 • C, and subsequently amplified using Oligo A. A DNA library was prepared using an Ion Plus Fragment Library Kit (Thermo Fisher Scientific, Austin, TX, United States) according to the manufacturer's protocol. Next-generation Ion Torrent sequencing of the library was conducted using the Ion Torrent platform on an Ion PGM System (Life Technologies, Foster City, CA, United States) using Ion PGM 400 sequencing reagents and Ion 318v2 chips following the manufacturer's instructions at the University of Vigo Central Services (CACTI). Quality control procedures and filtering of the resulting reads were afforded using PRINSEQ software (Schmieder and Edwards, 2011). Tag Sequence Check and Sequence Duplication routines were used to trim adapters and eliminate duplicates. Sequences shorter than 100 bp with a mean quality Phred score lower than 20 were removed. Tandem Repeats Finder (Benson, 1999) was used with all the parameters by default for locating and displaying tandem repetitions in DNA sequences. Forward and reverse primers were designed for effective microsatellite amplifications using FastPCR 6.5 software following Kalendar et al.'s (2009) recommendations. Finally, primers were proposed and tested by individual PCR on 30 individuals of P. pollicipes from 3 distinct geographic populations: 10 individuals from Baiona (Galicia, Spain), 10 individuals from Los Xatos (Asturias, Spain) and 10 individuals from Toulbroc'h (Brittany, France) (Figure 1). PCR tests were equally subdivided between different research laboratories with 41 primers tested per laboratory at the University of Vigo, the University of Oviedo and the Roscoff Marine Station. In this way, microsatellite markers were calibrated between geographic regions and the three institutes involved in the project.
Individual PCRs were conducted in a 20 µL total volume with Green GoTaq R Flexi Buffer (1×) (Promega Corporation, Madison, WI, United States), MgCl 2 (2.5 mM), dNTPs (0.5 mM), 0.2 µM of each primer, 0.1 U of GoTaq R G2 Flexi Polymerase (Promega Corporation, Madison, WI, United States), and 50 ng of DNA in sterile distilled water. The PCR program included an initial 5 min denaturation step at 95 • C, 35 cycles of denaturation at 95 • C for 30 s, annealing at 57 • C for 30 s and elongation at 72 • C for 30 s. PCR products were visualized using electrophoresis on a 2% agarose gel stained with SimplySafe TM (EURx, Gdañsk, Poland). Primer pairs without amplification, leading to a multiband pattern or a band size differing from its expected size, were discarded. Twelve microsatellite loci were amplified reliably and arranged in three multiplex PCRs (M1, M2, and M3 with four microsatellite markers per multiplex each) using Multiplex Manager 1.2 software (Holleley and Geerts, 2009) according to the dye colors and expected amplicon sizes. In addition to the twelve microsatellite markers retained with this screening, eight microsatellite markers previously developed in a parallel study (Fernandes et al., in prep) were tested, calibrated, and added in two supplementary multiplexes (M4 and M5) following the previously detailed methodology. This process resulted in a total of 5 multiplex PCRs. Forward primers were labeled using fluorescent dyes: 6-FAM TM , NED TM , VIC TM , and PET R (Applied Biosystems, Foster City, CA, United States) ( Table 1). PCR products were sequenced at the Genomer platform of the Roscoff Marine Station and at Servicios Científico-Técnicos of the University of Oviedo. Allele sizes were manually scored using GeneMapper v.4.0 (Applied Biosystems, Foster City, CA, United States).

Multiplex PCR and Microsatellite Genotyping
As explained above, all adult DNA was extracted with the EZNA R Mollusk Kit (Omega Bio-Tek Inc., Norcross, GA, United States). Juvenile DNA was extracted using the Chelex R 100 (Bio-Rad Laboratories Inc., Hercules, CA, United States) method (Estoup et al., 1996). PCRs were carried out following a unidirectional workflow that started in a pre-PCR room to prepare PCR plates. Amplification by PCR and processing of the subsequent PCR products always took place in a post-PCR area to avoid any possible contamination. M1, M2, and M3 multiplex PCRs were conducted using the QIAGEN Multiplex PCR Kit (QIAGEN Inc., Venlo, Netherlands) in a final reaction volume of 13 µL with the following components: 1× QIAGEN Multiplex PCR Master Mix, 1× Q-Solution, 50 ng of DNA template and 0.2-0.5 µM of each primer ( Table 1). PCR conditions consisted of an initial denaturation step at 95 • C for 15 min, followed by 40 cycles at 94 • C for 30 s, an annealing temperature of 60 • C (M1 and M2) or 64 • C (M3) for 1:30 min and 72 • C for 1 min, with a final extension at 60 • C for 30 min. M4 and M5 multiplex PCRs were incorporated and tested in a later stage and they were conducted   using the TouchDown PCR technique (Hecker and Roux, 1996 , and 5-10 ng of DNA in distilled water. The samples were initially heated at 95 • C for 5 min, followed by 10 cycles consisting of 95 • C for 30 s, 60 • C (decreasing incrementally by 0.5 • C per cycle) for 40 s, and 72 • C for 40 s, followed by 25 cycles at 95 • C for 30 s, 55 • C for 40 s, and 72 • C for 40 s, culminating in a final cycle at 72 • C for 10 min. PCR results were checked on a 2% agarose gel. For each multiplex amplification, 2 µL of reaction product (diluted 1/40 with Milli-Q water for M1, M2, and M3) was mixed with 9.5 µL of Hi-Di formamide (Applied Biosystems, Foster City, CA, United States) and 0.50 µL of SM594 molecular weight marker (Mauger et al., 2012). The mixture was heated at 94 • C for 5 min, immediately chilled on ice for 2 min, loaded in an ABI Prism R 3130XL automatic sequencer (Applied Biosystems, Foster City, CA, United States) of 16 capillaries using POP-7 polymer and run at 60 • C, 15 kV, 1,200 s using the sequencing platform Plateforme Genomer (Station Biologique de Roscoff). To ensure that the allele spread calibration held between the set of samples analyzed, controls were included in each plate to be genotyped as reference genotypes. Each genotype was then scored after analyzing the amplification products with Genemapper 4.0 (Applied Biosystems, Foster City, CA, United States).

Population Genetic Analyses
The allele frequencies, number of alleles per locus (k), observed heterozygosity (H O ) and unbiased expected heterozygosity (H E ) were calculated with GENETIX 4.05 (Belkhir et al., 2004). Moreover, possible genotyping errors and null allele frequency estimation were conducted using MICRO-CHECKER 2.2.3 (Van Oosterhout et al., 2004) and FreeNa (Chapuis and Estoup, 2007) with a number of replicates fixed to 10,000. Moreover, to explore the influence of null alleles on data we assessed F IS and F ST correlation, F IS and the number of missing data (putative null homozygotes) correlation and estimated the StrdErrF IS and StrdErrF ST values following the De Meeûs (2018) and Manangwa et al. (2019) [but see Waples (2018)].
The significance of correlations was tested with a unilateral (ρ > 0) Spearman's rank correlation test with Rcmdr package (Fox, 2005(Fox, , 2007 for R. Furthermore, for each population, the number of private alleles was calculated with GENALEX 6.5.03 (Peakall and Smouse, 2012).
Possible deviations from expected proportions in Hardy Weinberg's equilibrium and linkage disequilibrium for each locus and population were assessed using FSTAT 2.94 software (Goudet, 1995). FSTAT 2.94 software (Goudet, 1995) was used to calculate the allelic richness (A R ) and to determine the fixation indices (F-statistics) within and across populations using the method described by Weir and Cockerham (1984). Significance levels of F IS were estimated by permutating alleles between genotypes within samples 2,000 times and adjusted following Bonferroni correction (Rice, 1989) from all tested juvenile and adult samples. To test self-recruitment, the relatedness between individuals (R XY ) was estimated with the "related" package in R (Pew et al., 2015). The relative performance of seven different relatedness estimators was examined (dyadml, lynchli, lynchrd, quellergt, ritland, trioml, and wang) through comparison of the observed values to expected values generated from a simulated sample set of 400 individuals of known relatedness [with one hundred individuals from 4 categories: parent-offspring (R XY = 0.500), full-sib (R XY = 0.500), half-sib (R XY = 0.250), and unrelated pairs (R XY = 0.000)]. The results showed that the dyadic likelihood relatedness estimator (dyadml) provided the most consistent estimates through all possible levels of kinship; therefore, it was performed with 500 iterations. The bottleneck hypothesis was tested using the software BOTTLENECK 1.2.02 (Piry et al., 1999) under the two-phased model of mutation (TPM), taking into account 90% single stepwise mutations with a variance of 12.
Comparisons between regions and between cohorts (adults and juveniles) were conducted using a two-sided statistical analysis included in the FSTAT software for several statistics [A R , H O , H E , F IS , F ST , relatedness (R), and corrected relatedness]. In addition, F ST values were estimated using FreeNa, which estimates unbiased F ST following the ENA method (Chapuis and Estoup, 2007). The F ST values and associated p-values between cohorts and within and between regions were also calculated using FSTAT 2.94 (Goudet, 1995) to test for the regional and local structure. To assess the significance levels of F ST , multilocus genotypes were permutated 2,000 times between pairs of samples, and the significance threshold was obtained by applying a false discovery rate (FDR) over samples (Benjamini and Hochberg, 1995). Partial Mantel tests to estimate the correlation between genetic and geographical distance were performed with FSTAT 2.94 (Goudet, 1995) using the INA correction method for the chord distance (Cavalli-Sforza and Edwards, 1967) (D CSE ) provided by FreeNA (Chapuis and Estoup, 2007) and combining a ln transformation of Haversine geographic distances following Séré et al. (2017) and Rousset's θ/(1-θ) and a log transformation of Haversine geographic distances with 10,000 permutations (Rousset, 1997).
The software BayeScan v2.1. (Foll and Gaggiotti, 2008) was used to identify candidate loci deviating from neutral expectations from genetic data using differences in allele frequencies between populations. Twenty pilot runs of 5,000 iterations each, followed by an additional burn-in of 50,000 iterations and then 5,000 samplings with a thinning interval of 10, were conducted. To correct for multiple testing, the program computes q-values based on the posterior probability for each locus. Loci with α-values significantly > 0 and q-values < 0.05 were defined as "outliers" -, i.e., loci putatively under directional selection. Loci with α-values significantly < 0 were considered putatively under balancing selection. The remaining loci were classified as neutral.
An analysis of molecular variance (AMOVA) implemented in Arlequin 3.5.1.3 (Excoffier et al., 2005) to partition genetic variation across nested levels, regions and sites within regions was used. For the AMOVA, the number of different alleles was used as a measure of genetic variation (F ST -like option in Arlequin), and 10 000 permutations were used to test for statistical significance. Moreover, the "adegenet" package in R was used to estimate the genetic differentiation and visualize individual clustering with principal component analysis (DAPC, Jombart, 2008;Jolliffe, 2011) among adults and juveniles from each of the three regions separately and both among adults and among juveniles for all three regions pooled together. A neighbor-joining (NJ) tree based on the pairwise Nei's genetic distance D A (Nei and Takezaki, 1983) for all microsatellites and localities (15 localities; adults and juveniles grouped together) and then adding temporal cohorts as independent samples (i.e., 15 localities and 2 cohorts, 30 samples) was constructed with the software POPTREEW (Takezaki et al., 2014) using 10 000 bootstraps and visualized in The Interactive Tree of Life 3 (Letunic and Bork, 2019). Finally, STRUCTURE 2.3.4 (Pritchard et al., 2000) was also run to explore the population structure with Bayesian clustering. STRUCTURE was run using the 15 localities and 30 samples using admixture (Gilbert et al., 2012;Novembre, 2016) and also using adults and juveniles taken separately from the three regions (Portugal, Galicia and Asturias) in the same conditions to explore putative genetic units. The settings used were an admixture model from K = 1 to K = 30 in 20 runs following Evanno et al. (2005) and (Gilbert et al., 2012). Assignment clusters were made with burn-in periods of 20,000 and 200,000 Markov chain Monte Carlo repetitions. The most likely value of K was chosen using the delta K statistic (Evanno et al., 2005) using STRUCTURE HARVESTER software (Earl and VonHoldt, 2012), and visualization and grouping of the individual STRUCTURE runs was performed using CLUMPAK (Kopelman et al., 2015).

RESULTS
The typical upwelling circulation during the stalked barnacle larval season in summer/autumn 2017 along the coasts of northern and western Iberia, sea surface temperature (SST) along with modeled sea surface currents datasets revealed that the SST patterns showed strong onshore advection of cold waters (13-15 • C) on the Galician and Portuguese shelves with upwelling filaments extending further offshore especially at the upwelling centers of Fisterra, A Guarda, and Cape da Roca (Figure 1). Slightly onshore cooling indicative of upwelling was also observed along the western Cantabrian coast. Westward and southward currents in the order of few cm/s off the Cantabrian and Atlantic shores, respectively, clearly pointed to upwelling circulation (Figure 1). These flows are weaker close to the coast probably due to friction with the coastal boundary layer. Off southern/central Portugal in between Cape da Roca and Cape San Vincente, an anticlockwise cyclonic eddy was apparent with strong southward currents (>10 cm/s) along its western side. The dynamic structure of this feature matched SST patterns remarkably well, with a warm core (18 • C) surrounded by colder upwelled waters (14 • C).
The microsatellites markers development process produced libraries with a total amount of 42,860 reads showing a mean sequence length of 91.61 ± 103.29 bp (minimum length: 25 bp -maximum length: 517 bp) and a mean GC content of 63.66 ± 18.90%. A total of 10 781 sequences with a mean sequence length of 244.48 ± 97.61 bp, a length range of 418 bp and a mean GC content of 50.30 ± 5.61%, resulted after quality control procedures. A total of 1,140 sequences containing di, tri, tetra, and pentanucleotides were selected after locating and displaying tandem repetitions in DNA sequences. Finally, 123 pairs of primers were proposed and tested in three different research laboratories (University of Vigo, University of Oviedo, and the Roscoff Marine Station). A new set of twelve microsatellite loci currently arranged into three multiplex PCRs (M1, M2, and M3) was developed for the stalked barnacle P. pollicipes in this work (Genbank accession numbers: MW443103-MW443114). Moreover, eight previously developed microsatellite loci by Fernandes et al. (in prep) were also tested and included in another two multiplexes (M4 and M5, Genbank accession numbers: MZ576446-MZ576456). This procedure resulted in a total of 5 multiplex PCRs (Table 1) leading to scorable and reproducible genotypes for all 20 microsatellite loci. None of these loci showed evidence of linkage disequilibrium between alleles (p > 0.05). These loci were highly polymorphic and exhibited approximately 15% private alleles (n = 87) only present at one locality (Tables 1, 2).
The comparative analysis for levels of genetic variation between regions [Portugal (PTL), Galicia (GAL), and Asturias (AST)] revealed no significant differences for expected heterozygosities (H S ) ( Table 2). Slight differences in genetic diversities were, however, observed depending on the population parameter estimated. Galicia showed the highest values for allelic richness and observed and expected heterozygosity in adults and juvenile populations ( Table 2). Portugal showed the highest number of private alleles (mean A P PTL = 3.3), which was mainly attributable to adults (mean A P = 4.4) ( Table 2). Significant differences in allelic richness were also observed for juveniles (A R PTL: 11.569; A R GAL: 12.418; A R AST: 11.900; p < 0.01) and in observed heterozygosity for adults  Table 2). This phenomenon was especially obvious in Portuguese samples, where the average number of private alleles decreased by 50% (A P AD: 4.4 to A P JV: 2.2) ( Table 2). In this later region, significant differences were found between adults and juveniles both in terms of allele richness (A R AD: 12.318; A R JV: 11.569; p < 0.01) and observed heterozygosity (H O AD: 0.640; H O JV: 0.605; p < 0.05) or expected heterozygosity (H S AD: 0.773; H S JV: 0.757; p < 0.01) ( Table 2). Globally, juveniles were also more related to each other (R XY = 0.067) than their adult (R XY = 0.058) counterparts, as indicated by relatedness analyses. Juveniles from Portugal (R XY value = 0.073, p < 0.002) and Asturias (R XY = 0.070, p < 0.002) were much more related than expected from panmixia. Bottleneck software showed that none of the 30 samples tested (15 localities × 2 cohorts) exhibited a significant excess of predicted heterozygotes under the TPM model and could not be considered to have experienced a recent genetic bottleneck ( Table 2). When the bottleneck hypothesis was tested with all juveniles and adults together (15 samples) and at the regional scale (15 samples grouped in 3 regions, for 2 cohorts), the statistics remained non-significant ( Table 2).
According to the overall F ST , there was no significant genetic differentiation of adults between and within regions (Figure 2A). Only 20 out of the 75 possible pairwise F ST values between adult samples from different regions (25.3%) showed p-values lower than the 0.05 cutoff value, and these critical values were more often encountered between Galicia and Asturias (12/25 = 48%) (Figure 2A). However, no p-values remained significant after FDR correction (Figure 2A). In contrast, the overall F ST statistics estimated for the juveniles between and within regions indicated notable regional and spatial structuring ( Figure 2B). Pairwise F ST estimated between juvenile samples from Galicia and Portugal (12/25 = 48% before and 6/25 = 24% after FDR) and between Galicia and Asturias (13/25 = 52% before and 6/25 = 24% after FDR) confirmed this trend and showed clear regional structuring (Figure 2B). Asturias and Portugal were, however, less differentiated from each other, with fewer significant pairwise F ST values (3/25 = 12% before and 2/25 = 8% after FDR) ( Figure 2B). Some spatial structuring within regions was detected for juveniles using pairwise F ST analyses but only in the case of Portugal (2/15 = 13% before and 1/15 = 6% after FDR) ( Figure 2B). The pairwise F ST analyses between adults and juveniles within regions revealed that only 8 out of the 75 possible comparisons (11%) had p-values lower than the 0.05 cutoff threshold for Portugal and Asturias (but not in Galicia), which, however, did not remain significant after FDR correction (Figures 2C-E). The analysis conducted with BayeScan v2.1 for outlier detection resulted in no loci under selection or biased by species admixture and hybridization which have the same expectations in terms of outliers; the twenty loci showed signatures of balanced or purifying selection with negative alpha values. The results of the partial Mantel tests indicated no correlation between genetic and geographic distances, with R 2 = 1.61 and p-value = 0.1916 for adults and R 2 = 3.16 and p-value = 0.0698 for juveniles using the INA correction method for D CSE , and R 2 = 0.06 and p-value = 0.8002 for adults and R 2 = 0.22 and p-value = 0.6241 for juveniles using the Rousset method. The population structure was therefore closer to an n-island model than a stepping stone model, and the pairwise F ST between adjacent sites often exceeded those obtained between geographically distant locations.

DISCUSSION
The analyses using twenty new microsatellite loci aimed to define, more accurately, the temporal and spatial evolution of the genetic structure of stalked barnacle P. pollicipes. This species is highly appreciated in the Spanish and Portuguese markets, and its management must be based on reliable scientific data. Previous studies have suggested that larval dispersal driven by ocean currents, in particular, the Iberian Poleward Current have played a crucial role in determining the population structure, and two distinct regional configurations FIGURE 2 | F ST heatmaps [based on Weir and Cockerham (1984)] following genetic analyses of P. pollicipes using microsatellites along the Iberian Peninsula. The darker the color, the higher the F ST value. Asterisks (*) indicate significant p-values (p < 0.05) while significant values after a FDR correction (Benjamini and Hochberg, 1995)  have been established using mitochondrial DNA for P. pollicipes within its distribution range along the northeastern Atlantic. Quinteiro et al. (2007) suggested that P. pollicipes is structured into four genetically differentiated groups: French populations, eastern Asturian populations, Galician-Portuguese populations, and Canarian populations. Conversely, Campo et al. (2010) suggested the presence of only two groups, among which French populations were highlighted as a peculiar and differentiated genetic entity, as a result of a past population fragmentation during Pleistocene glacial/interglacial periods. Regardless, later studies based on estimates of population migration rates have suggested that barnacle population connectivity occurred on a small scale and in an asymmetric manner in the Cantabrian coast (Rivera et al., 2013). Information based on highly variable nuclear molecular markers can provide crucial information on both population connectivity and stock renewal for this species within the Iberian Peninsula. This information is needed for the delimitation of conservation/management units in this fishery and the improvement of the management plans and the performance of TURFs.
Genetic diversity contributes to the ability of a species to respond to environmental changes, and highly fecund species that release high numbers of small eggs into the environment (the so−called r−strategists) are much more polymorphic than species that produce a small number of relatively large offspring and provide parental care (called K−strategists) (Ellegren and Galtier, 2016). Recent studies in S. balanoides have confirmed that barnacles harbor high levels of genome-wide genetic variation (Nunez et al., 2021). The level of genetic diversity of P. pollicipes found in this work was particularly high. We observed higher levels of genetic variation in P. pollicipes than in other barnacles of the same genus, such as P. elegans . Our results showed that Galicia exhibited the highest values for allelic richness and observed and expected heterozygosity in adult and juvenile populations. Conversely, newly settled cohorts (juvenile) had a lower genetic diversity than adults across all the studied regions, particularly when examining both allelic richness and private alleles.
The main principal assumption of the Hardy-Weinberg principle is that the sample comes from a single, randomly mating population where perturbing forces (such as selection, genetic drift, mutation, and migration) are absent or balanced (Waples, 2014). All loci and populations showed significant deviations from Hardy-Weinberg equilibrium in this work due to, more or less pronounced, heterozygote deficiencies. This phenomenon could be the consequence of local admixtures of genetically differentiated populations (Wahlund effect), assortative mating, inbreeding, selection (Palumbi, 2003) and finally null alleles. The presence of null alleles has been reported in the vast majority of previous microsatellite studies in barnacles (Dufresne et al., 1999;Pannacciulli et al., 2005;Reigel et al., 2015;Abreu et al., 2016;Ewers-Saucedo et al., 2016) as well as in other marine invertebrate species such as clams (Borrell et al., 2014;Chiesa et al., 2016;Rico et al., 2017), octopus (Greatorex et al., 2000;De Luca et al., 2016), sea urchins (Mccartney et al., 2004;Carlon and Lippé, 2007), jellyfish (Aglieri et al., 2014), and polychetes (Jolly et al., 2003(Jolly et al., , 2009(Jolly et al., , 2014. The presence of null alleles is an inherent trait of microsatellite loci and is caused by mutations in the primer sequences, leading to the lack of amplification and the dropout of alleles (Selkoe and Toonen, 2006). In addition, an increase of the null allele frequency would be expected with the increase of alleles per locus and previous studies have indicated that the presence of null alleles seems to be particularly common in populations with high effective population sizes (Chapuis and Estoup, 2007). Although the presence of null alleles leads to an overestimation of both F ST and genetic distances in cases of significant population differentiation (Chapuis and Estoup, 2007), our results showed no differences worth considering for both the F ST or F ST ENA values. It has been argued that the conservative approach of discarding loci deviating from Hardy-Weinberg equilibrium expectations could rob us of our most informative markers, weakening our ability to interpret biological phenomena (Dharmarajan et al., 2013). Moreover, De Meeûs (2018) stated that in case of null alleles, F IS and F ST are augmented and a positive correlation is expected between F IS and F ST as is expected a positive correlation between F IS and the number of missing data (putative null homozygotes), and StrdErrF IS being at least twice StrdErrF ST . If FIGURE 4 | Neighbor Joining trees using DA distance (Nei and Takezaki, 1983) of P. pollicipes populations using microsatellites along the Iberian Peninsula.   Waples (2018) had also argued about this and simulated 10% of null alleles suggesting that caution in interpreting F IS × F ST correlations under conditions where null alleles might be common it is indeed necessary and more efforts will be needed for a comprehensive evaluation of this complex topic. In this work, panmixia is rarely met for any locus (Table 1), we found positive F IS × F ST correlations, StrdErrF IS > StrdErrF ST and we did not find positive correlations between F IS and the number of missing data (putative null homozygotes) pointing out to the fact that, even when null alleles are present, other biological factors also play a fundamental role to explain significant heterozygote deficits in our data.
Heterozygote deficiencies can as well be the result of local admixtures of genetically differentiated cohorts in populations, or due to sweepstake reproductive effort (Waples, 1998;Hedgecock and Pudovkin, 2011). Growth of individuals in P. pollicipes populations is highly heterogeneous Jacinto et al., 2015), so that individuals of similar size may differ greatly in age. Our adult samples likely contained a mixture of cohorts from different reproductive and dispersal events, potentially leading to significant departures from Hardy-Weinberg equilibrium, locally. Genetic heterogeneity of cohorts can potentially blur the genetic signal in adults and may decrease the genetic differences over time, given that the geographic origin of migrants might change throughout the breeding/dispersal seasons depending on prevailing local hydrodynamics during these periods. However, it should be noted that a special care was taken in this work to sample only one cohort of juveniles with a specific size (2-4 mm RC). If the deficiencies of heterozygotes were due the superimposition of cohorts, juveniles should not show such deficiencies. This was clearly not the case here as our results demonstrated that juvenile mean F IS values were higher than those for adults in all the three regions ( Table 2). It has been stated that the surf zone and its surrounding nearshore waters are known to act as selective barriers to the onshore transport of many larval invertebrates on the local scale (Porri et al., 2006;Rilov et al., 2008). The permeability of such barrier is modulated by small scale topography that generates retentive oceanographic features like coastal fronts (Pineda, 1999;Shanks et al., 2003). In fact, the larvae of P. pollicipes and other barnacles have been shown to accumulate in great numbers at internal waves and river plume fronts off the Cantabrian coast only at some specific locations (Weidberg et al., 2014;Höfer et al., 2017). In this topic the available evidence are indeed scarce, however, genetic data seems to confirm it.
Pollicipes pollicipes has asynchronous broods during the reproductive season which usually occurs from March to September (e.g., Cardoso and Yule, 1995;Cruz and Hawkins, 1998;Pavón, 2003;Macho, 2006), where several batches of larvae are produced, and potentially lead to the co-occurrence of different settlement events. Juveniles sampled in this study might, however, come from one to few settlement events. Despite the possibility of several discrete settlement events, post larval mortality might favor one specific batch of survivors, and in the end, the 2-4 mm RC juveniles might become more related than what would have been expected from the mixing of several reproductive events. Pineda et al. (2006) found that recruitment to the reproductive stage of acorn barnacles (S. balanoides) was composed of survivors that settled in a recruitment window. The recruitment window (to reproduction in the case of the Pineda study, to 2-4 mm in our study) might be narrower than the recruitment season. If by some reason these survivors correspond to larvae that are genetically more related, then a pattern of genetic differentiation could occur among recruits. The concept of a "recruitment window" proposed by Pineda et al. (2006) matches quite well with Hedgecock's "sweepstakes-chance matching hypothesis" also known as "sweepstakes reproductive success hypothesis, " which is based in part on the observation of reduced genetic variability in young-of-the-year populations relative to adult populations. This reduced genetic variability among recruits suggests that the surviving young of the year are the products of spawning by only a small fraction of the adult population, which, according to Hedgecock's hypothesis, happened to produce their offspring at a place and time that was suitable for survival (Hedgecock, 1994). Moreover, barnacles rear embryos in bags before hatching and there is also the possibility that the larval release is only efficient for a small proportion of the reproductive adults depending on the local hydrodynamics. In this work, we found evidence indicative of reproductive sweepstakes in adult and juvenile samples. Although globally, the relatedness coefficients estimated for P. pollicipes were in the same range as those from other studies previously conducted with barnacles (Veliz et al., 2006;, they were significantly slightly greater in juveniles (i.e., Asturias and Portugal) compared with adults. Juveniles were significantly more related to each other than expected from random mixing despite their larval entrainment in the water column during the planktonic phase.
Unexpected genetic differentiation in marine invertebrates can occur due to three neutral processes: sweepstake reproductive success (Hedgecock, 1994), collective dispersal (Johnson et al., 1993;Li and Hedgecock, 1998) and asynchronous population dynamics (Eldon et al., 2016), but also selective processes during the settlement process. According to the Hedgecock's "sweepstakes-chance matching hypothesis" or selective sweepstakes (Hedgecock, 1994), only a fortunate combination (hence sweepstakes) of reproductive traits and oceanographic conditions would allow an individual to complete the long mobile phase from spawning and fertilization through larval survival to recruitment back to the adult habitat. In a highly fecund species and a locally heterogeneous oceanographic setting, this would involve strong selection favoring just a handful of genotypes at each locality, leading to a local-scale genetic mosaic but a relatively large-scale uniformity. Postlarval settlement selection under different environmental conditions has been argued to create CGP in coastal areas of temperate regions over a mosaic of contrasting habitats able to impose a strong differential selective sieve or a target for habitat choice in larvae (Eldon et al., 2016). We detected significant genetic differentiation for juveniles among and within regions (but not for adults), together with significant genetic heterogeneity between P. pollicipes generations. However, we did not find evidence of such selective processes for the assayed microsatellites. There seemed to be a genome-wide pattern that was more parsimoniously explained by neutral processes such as sweepstake reproductive success, which may greatly reduce the genetic diversity of a given cohort while provoking unexpected heterozygote deficiencies, as seen previously, by mimicking local bottlenecks (genetic diversity drawn from a small subset of parents). In addition to this phenomenon, genetic differentiation may persist in recruits when dispersal is limited in space, when larvae from different cohorts do not mix completely during dispersal (collective dispersal), or when local conditions may promote self-recruitment (Eldon et al., 2016).
The genetic data obtained in this work, after applying dissimilar approaches (F and φ ST statistics, Discriminant Principal Component and Bayesian analyses), pointed all out to the existence of significant genetic heterogeneity in the Iberian coasts rejecting previous findings using mitochondrial DNA. The results herein highlighted Galicia as a peculiar genetic entity possibly representing a superimposition of two distinct metapopulations or potentially an old refuge for the most northern populations from France (not sampled in this study). Among Galician northernmost populations, Camelle (CA) and A Coruña (AC) are also the most differentiated from Portugal and Asturias and may have a specific demographic history. The sampled P. pollicipes populations are located along the Atlantic Iberian coast, whose hydrodynamic patterns have been well studied. The western peninsular coast (SW Portugal and Galicia) is characterized by a complex current system subject to strong seasonality and mesoscale variability, showing inverse patterns between summer and winter in the upper layers of the shelf and slope. During spring and summer (coinciding with P. pollicipes breeding season), northerly winds along the coast are dominant, causing coastal upwelling and producing a southward current on the surface and a northward undercurrent on the slope. In the Cantabrian Sea (Asturias) the surface currents flow generally eastward in winter and early spring and shift westward in late spring and summer following the wind force with intermittent summer upwelling events west of Cape Peñes (ICES, 2021). Different aspects of the oceanographic circulation in Iberia were reviewed by Relvas et al. (2007). Casabella et al. (2014) divided the upwelling affecting the coasts of Galicia into three regions: Rías Baixas, Fisterra-Bares and Cantabrian. These two locations (CA and AC) would be found in the Fisterra-Bares region, which is the region with a greater intensity of upwelling, although the period favorable for upwelling is longer in the region of Rias Baixas (sampled here i.e., Baiona). Galician juveniles showed clear genetical differences from those of Portugal and Asturias. The main explanation for this distinction is that the Biscay Bay Current, characterized by a wide gyre, can trap larvae, and thus should favor selfrecruitment and perhaps local importations from the French and Cantabrian populations. This ocean circulation could also be responsible for the differentiation between juveniles from Asturian and Galicia. Previous studies on adult barnacles have found significant differences between the Asturian and Galician localities (Quinteiro et al., 2007). However, it should be noted that most of the Asturian sites sampled in this study are located to the West of Cape Peñes, while the site sampled by Quinteiro et al. (2007) was located to the East of the same cape, which has been described as a biogeographic barrier (Anadón and Niell, 1980). Rivera et al. (2013) showed that during a year of high upwelling activity (2009), the theoretical P. pollicipes recruitment success was 94%, with a recruitment peak predicted 56 km west of the emission point. Consistently, migration rates derived from genetic analyses showed that westward dispersal was much more likely along the Cantabrian coast, which matches the upwelling driven circulation typical of the stalked barnacle larval season in summer/autumn (Figure 1). Thus, the recurrence of upwelling may not only define the spatial scale and direction of the dispersal process but also the genetic structure of the barnacle metapopulation.
The Western Iberian upwelling system represents an important crossroad between Lusitanian and boreal temperate species (Jolly et al., 2006;Maggs et al., 2008). Upwelling/downwelling wind-driven circulation and tides are recurrent physical processes along the Atlantic Iberian coastlines and are among the most energetic phenomena that can affect near-shore circulation during the spring and summer periods when reproduction occurs and, during the summer and beginning of autumn in the case of recruitment (Queiroga et al., 2007). However, when studying a strong upwelling region in the northeastern Pacific coast, Morgan et al. (2009) observed that the larvae of most invertebrate species remain close to the shore even during strong upwelling, where high local retention and limited connectivity have been evidenced in populations of several species, such as Petrolisthes cinctipes (Hameed et al., 2016) or in the red rock lobster Panulirus interruptus (Iacchei et al., 2013). Despite this phenomenon, upwelling areas have been pointed out as probable climate change refuges for the distribution of Fucus guiryi, other barnacles such as S. balanoides and other sessile marine species (Gómez and Lunt, 2007 for a review ;Hoarau et al., 2007;Provan and Bennett, 2008;Lourenço et al., 2016;Herrera, 2019). In addition, Campo et al. (2010) suggested the existence of a Pleistocene refuge area off the coast of North Africa and two additional northern glacial refuges for P. pollicipes, in the English Channel/Brittany region and in the northwestern Iberian Peninsula.
Previous studies have mentioned that the southern region of Portugal also represents a well-known upwelling area (Lourenço et al., 2016) with a high level of barnacle larval settlement (Queiroga et al., 2007) and recruitment (Aguión et al., in prep.). Remarkably, the number of private alleles was significantly higher in adults there when compared with those from Galicia and in Asturias. Portuguese juveniles were, however, significantly less genetically diversified and more related to each other than expected based on random mating. Moreover, Nolasco et al. (in prep.) show that connectivity matrices integrated over the period of the observations (July 2017 to July 2019) indicate high levels of larval retention. Such retention is probably caused by the recurrent eddies driven by upwelling circulation observed off southern-central Portugal in between Cape Roca and Cape San Vicente (Figure 1; Haynes et al., 1993;Batteen et al., 2000;Sánchez and Relvas, 2003;Peliz et al., 2004). These findings suggest that Portuguese populations are likely to export more migrants than they receive. As Queiroga et al. (2007) hypothesized, regular exchanges of larvae over the distance separating the southern and northern parts of Portugal are unlikely. Conversely, the Portugal Current, which shows a northor southward direction, depending on the season, could be an important factor in promoting gene flow between our sampling locations in southern Portugal and other, unsampled, P. pollicipes southernmost areas such as the Canarian and North African coasts. Nevertheless, microsatellite markers have recently shown a genetic differentiation between European and African P. pollicipes populations (Fernandes et al., in prep).
The correct management of marine ecosystems relies on understanding the scale and magnitude of connectivity among populations through the identification of adaptive genetic differences (Almany et al., 2009;Aceves-Bueno et al., 2017), because locally adapted populations should be considered poorlyconnected, separate management units (Waples, 1998). Our results suggested that P. pollicipes populations in the Iberian Peninsula possibly exhibit a "CGP" structure, which extends from a few kilometers apart to as much as hundreds of kilometers apart. This phenomenon has clear consequences for the sustainable management of resources. Currently, an increasing number of small-scale fisheries have successfully implemented co-managed TURFs; a governance arrangement that enables the collaboration across diverse stakeholders, develops new knowledge and increases the capacity of the system to deal with new drivers (Rivera et al., 2014). However, the design of TURFs does not usually account for the spatial configuration of resources (Aceves-Bueno et al., 2017) due to the multi-species nature of fisheries. This mismatch between management and biological scales can compromise the sustainability of sessile stocks (Ouréns et al., 2015), like barnacles. However, a better understanding of the spatial structure and larval dynamics of the population, permits the redefinition of management units according to population boundaries. In addition to these management measures, it would be interesting to implement networks of protected areas at detailed scales to ensure that propagules are available when and where conditions are favorable for their survival (Larson and Julian, 1999;Ouréns et al., 2015).

CONCLUSION
New molecular markers have been developed in the highly valued species P. pollicipes and offer useful tools to provide a better finetuning assessment of its population dynamics along the Iberian Peninsula. P. pollicipes displays high genetic diversity, which is attributable to large effective population sizes representing a well-connected network of local populations. However, temporal and spatial genetic differentiation of populations over regional scales, on one hand, and a significant reduction in genetic diversity in juveniles, on the other hand, clearly indicate that patterns of exchanges together with seasonal wind-induced upwelling may induce genetic differences between settlers throughout generations. Such patterns of CGP are likely due to sweepstake reproductive success with possible collective dispersal or episodic self-recruitment events. Therefore, our P. pollicipes genetic dataset suggests that recruitment may be stochastic and highly dependent on climatic conditions with multiple sources of emissions. These phenomena may have strong implications in terms of management plans over the whole Iberian Peninsula with the need to protect a series of putative sources within each region. Future research should combine genetic information at broader spatiotemporal scales with larval dispersal models based on ecological and biological characteristics of P. pollicipes. This means, among others, mapping the complete species distribution and tracking the genetic structure of age groups over time and space. It also means applying new sequencing technologies to fully understand the dynamics of larval exchanges and the post-larval settlement of the stalked barnacle but also to better apprehend how environmental variations shape genomic variation in this species.

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 in the article/ supplementary material.

AUTHOR CONTRIBUTIONS
MP: methodology, formal analysis, investigation, writingoriginal draft, and writing-review and editing. PM: conceptualization, supervision, methodology, and writingreview and editing. MB and NW: methodology and writingreview and editing. JLA: conceptualization, supervision, project administration, funding acquisition, and writing-review and editing. AA, JA, JNF, LG-F, and KJG: writing-review and editing. JC: methodology. TC and GM: conceptualization and writing-review and editing. EG-V: conceptualization, funding acquisition, and writing-review and editing. ET and DJ: conceptualization, supervision, and writing-review and editing. YJB: conceptualization, supervision, funding acquisition, and writing-review and editing. All authors contributed to manuscript revision, read, and approved the submitted version. This work also had the support of Fundação para a Ciência e Tecnologia (FCT), through the strategic project UIDB/04292/2020-MARE granted to MARE and Xunta de Galicia (GRC, ED431C 2020-05). AA was supported by a FPU fellowship (Ministerio de Ciencia, Innovacion y Universidades de España, Grant No. FPU2016-04258) and KJG was supported through the Severo Ochoa Ph.D. program (Principado de Asturias, PA-18-PF-BP17-184). GM was supported by postdoctoral contracts from project PERCEBES (PCIN-2016-063) and MARISCO (CTM2014-51935-R, Ministerio de Economía y Competitividad). NW was funded by NASA grant 80NSSC20K0074 during data processing and manuscript elaboration. This is a contribution of the Marine Observatory of Asturias (OMA) and the Biotechnology Institute of Asturias (IUBA).