Protists Exhibit Stronger and More Recent Oceanic Genetic Structure than Archaeplastida and Metazoa

The global distribution of biodiversity is shaped by a complex interplay of evolutionary history and ecological processes. While biogeographic patterns are well defined for animals and plants, the global distributions of protists remain unclear. A key question is whether protists follow the same broad biogeographic principles as macroscopic life. To address this, we compiled 88 marine COI metabarcoding studies performing population‑genetic analyses across ocean basins. Our results reveal that most protist phyla exhibit pronounced genetic structure among oceans, a pattern exceeding that reported for Archaeplastida and Metazoa. This likely reflects recent, and potentially human-mediated, introductions , influencing protist dispersal and contemporary community assembly. By demonstrating that protist distributions are not historically cosmopolitan, our study supports the existence of common eukaryotic biogeographic patterns that transcend organismal size. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Main Text: Unicellular eukaryotes (protists), which represent most of the diversity of eukaryotes from which multicellular groups evolved, have classically been viewed as fundamentally different from multicellular plants and animals with regard to their biogeographic patterns (1). Early theories postulated that protists have cosmopolitan distributions and virtually unlimited dispersal capacity, giving rise to the famous concept “everything is everywhere, but the environment selects” (2). This view, however, was based largely on limited morphological data. The advent of large -scale molecular techniques, such as metabarcoding, now permits better characterization of protist distributions (3). Some recent marine metabarcoding studies have begun to reveal that protists exhibit biogeographic patterns more similar to those of multicellular eukaryotes than previously assumed (4–6). Nevertheless, these investigations have predominantly focused on presence-absence matrices, beta diversity patterns, and relative abundances inferred from read depth. A critical knowledge gap persists regarding population-level genetic diversity and the historical biogeographic processes shaping these distributions (7). This gap hinders our ability to determine whether the biogeographic principles observed in multicellular animals and plants are truly universal across all eukaryotes , or if smaller taxa follow distinct rules in marine environments. To resolve this metabarcoding offers a promising approach (11), as the large amount of molecular and geographic data from different organisms enables the integration of geographical and population genetic theories within a unified biogeographic framework (8–10). However, its effectiveness depends on meeting certain conditions: i) the use of molecular markers with intraspecific resolution (for example , the mitochondrial cytochrome oxidase subunit I gene, COI, which has proven invaluable in animal biogeography (12)); ii) the availability of a comprehensive global database for such markers that includes data from both protists and multicellular eukaryotes (13); iii) the use of standardized taxonomic units that mitigate the inherent biases of metabarcoding (14); and iv) the integration of population genetic

and theory into metabarcoding analyses. Construction of a database with informative OTUs to test biogeographic patterns To investigate global eukaryotic biogeographic patterns at the population level, a critical first step is the creation of a comprehensive metabarcoding database amenable to population-genetic analyses. So, we here compiled data from 88 independent metabarcoding studies into an environmental -COI (eKOI metabarcoding) database, comprising 976,865 amplicon -sequence variants (ASVs) with no abundance filter ( Fig. S1) [see supplementary data (15)]; and 302,809,791 reads from 4,102 samples collected at 957 unique sites that span all oceans (Fig. 1a). We then taxonomically annotated them at phylum level, clustering them into 51,558 operational taxonomic units (OTUs). To ensure robust population -genetic analyses and minimize potential artifacts, we stringently filtered the data to retain only informative OTUs, defined as those i) containing at least four sequences from a minimum of two distinct locations , and ii) containing at least one pair of sequences that differ by  ≥ 1% (uncorrected p‑distance). This rigorous filtering resulted in a final dataset of 5,970 informative OTUs, comprising 276,303 ASVs and 101,898,078 reads distributed across 41 phyla (Fig. 1b, 2a). Crucially, inter-phylum comparisons of nucleotide diversity (π) and Tajima’s D revealed minimal significant differences (0.001% and 6.9% of pairwise comparisons, respectively; Fig. S2; Supplementary Data), confirming broadly comparable molecular-evolutionary pressures across phyla. Consequently, this constructed eKOI metabarcoding resource provides .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint standardized taxonomic units, thereby establishing a solid foundation for rigorous, population-level investigations of marine biogeographic hypotheses. Protists have larger geographic ranges than animals and plants A long -standing debate in eukaryote biogeography concerns whether protists possess larger geographic ranges than macro -organisms (1). Testing this hypothesis has been hindered by the sheer volume of molecular and geographic species-level data required at a global scale (16). Our eKOI database allows for a direct assessment of this question. We analyzed the distributional patterns of all eukaryotic phyla by comparing the range size of every informative OTU within each eukaryotic phylum, quantifying i) the mean pairwise geographic distance among sampling sites , and ii) the maximum distance between any two sites. The results show that while Metazoa and Archaeplastida show mean inter-sample distances around 1,000 km, most protist phyla average over 5,000 km, with some individual OTUs exceeding 10,000 km ( Fig. S3). This demonstrates that protists, at the OTU level defined by COI, possess considerably larger geographic ranges than plants and animals. Molecular and geographic expansions in protists To understand why protists exhibit these distinctive geographic distances patterns, additional data, particularly molecular data, are essential, since they enable the detection of within-OTU genetic patterns linked to geographic range (9). Therefore, we compared genetic versus geographic distances within each informative OTU to test whether the observed spatial patterns were non -random and consistent across different phyla. Our analyses using both linear and nonlinear regressions revealed that most OTUs displayed significant non -random spatial patterns. Specifically, when compared to null models, 81.5% of OTUs in the linear models and 92.6% in the nonlinear models had R² values that were statistically significant (p < 0.05; Supplementary Data). Furthermore, linear and nonlinear models (Fig. S4), as well as Mantel tests (Fig. S5), consistently showed positive relationships between genetic and geographic distances across all eukaryotic phyla examined. In other words, genetic divergence tended to increase with greater geographic separation. The rate of genetic divergence per unit distance, quantified by the slope of the linear regression and by the parameter “a” of the nonlinear (exponential) model, was largely uniform across phyla. Significant differences among phyla in these rates were rare, occurring in only 5.7% of pairwise comparisons for the linear slopes and 3.6% for the nonlinear “a” parameters. Moreover, the geographic extent of an OTU (mean inter - sample distance) was only very weakly co rrelated with its genetic divergence rate, explaining merely 3.4% of the variance in slope and 4.0% in the “a” parameter. Together, these results indicate that geographic distance alone does not strongly limit genetic differentiation in protists, even across vast oceanic distances. We propose that the primary explanation for this weak isolation-by-distance pattern is the recent introduction and rapid geographic expansion of many protist species (17, 18). Such rapid expansion over short evolutionary timescales would limit the accumulation of substantial genetic divergence, thus maintaining high genetic homogeneity across geographically distant populations. Nevertheless, this interpretation should be treated with caution, as multiple introductions can sometimes reduce genetic founder effects (19), yet molecular diversity often remains lower compared to that observed in the centers of origin. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Metazoa exhibit higher richness of haplotypes among eukaryotes To follow up on those rapid geographic expansions, we then inquired how and when current species distributions were established. To achieve this, we performed explicit population‑genetic analyses to determine how molecular diversity is structured within each informative OTU. Our approach involved examining the genetic structure of these OTUs across different oceanic regions, treating each region as a distinct population. This allows for inferences regarding both contemporary and historical processes shaping current species distributions. To ensure unbiased analysis, population -level units within each OTU were defined irrespective of sequencing depth variations. We observed that haplotype richness per informative OTU had only a weak correlation with sequencin g depth (R² = 18.2%), therefore being a robust metric for inferring population -genetic patterns in metabarcoding studies. Using the haplotypes identified within each OTU, we constructed haplotype networks for each informative OTU, without imposing a strict ly bifurcating tree topology (20). This network -based approach allowed us to infer intra - OTU genealogies while avoiding the assumption of purely bifurcating relationships inherent to traditional phylogenies, which may not accurately capture genetic structure at the population level (21). As a first step, we characterized haplotype diversity distributions across major phyla using three complementary metrics: haplotype diversity (Hd) (22), branch diversity (Bd), and a combined haplotype -branch diversity index (H bd) (23). These metrics were calculated based on the number of ASVs (unique haplotypes), the total read abundance per haplotype, and log -transformed read counts. B d and H bd values showed relative consistency across metazoan phyla, with few significant inter -phylum differences ( Fig. S6). In contrast, H d exhibited greater variation across phyla, primarily driven by comparisons involving Metazoa, which consistently displayed near -maximal Hd values. This indicates a significantly higher richness of unique haplotypes and, consequently, greater intra-OTU genetic diversity within metazoan OTUs compared to those of other eukaryotic phyla. Genetic Structuring Among Oceanic Regions suggest recent OTU introduction events in protist phyla Following the characterization of within -phylum haplotype diversity, we next assessed genetic structuring across oceanic regions. For each informative OTU, we analysed haplotype frequency and distribution patterns across sampled oceans. We quantified population-level differentiation using four standard metrics of population genetic structure: total genetic diversity (H t), within -ocean genetic diversity (H s), the fixation index (Fst) (a measure of among-ocean differentiation) (24), and the effective number of migrants (Nm) (an indirect estimate of gene flow) (25). These metrics revealed significant differences among phyla in their degree of genetic structuring, with 20% of inter -phyla comparisons showing a significant difference in H t, 9.7% in H s, 21.9% in F st, and 9.3% in Nm (Fig. S7). The most pronounced disparities were observed in comparisons involving Metazoa and Archaeplastida versus other phyla. Specifically, Metazoa and Archaeplastida showed higher H s values, reflecting greater average haplotypic diversity within individual oceans. This inter -oceanic divers ity was coupled with very low F st values in both clades (often approaching zero and never exceeding 0.5), signifying minimal genetic differentiation among oceanic regions. So, Metazoa and Archaeplastida haplotypes tend to be geographically restricted. In contrast, protist phyla exhibited much stronger genetic structuring, with consistently lower Hs and high F st values (≥ 0.5). This .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint pattern indicates that although lineages in those phyla occur across multiple oceans, their populations remain genetically isolated between regions. We propose that this strong structure in protists may reflect recent geographic expansions or introductions rather than long-term historical isolation. To address the limitations of H s, in capturing haplotype distribution across oceans, we developed the Population-level Haplotype Diversity (PHd) metric, which integrates intra- oceanic diversity and the evenness of haplotype distribution across oceans [see supplementary materials and methods (15)]. PHd penalizes uneven distributions, and the difference between PHd and Hs (PHd - Hs) indicates the degree to which an OTUs genetic diversity is structured between oceans. PH d analysis revealed significant inter -phylum differences, 20.6% of comparisons, largely mirroring H s and F st trends. Metazoa, Archaeplastida, and Ochrophyta showed near -zero PH d values and negative PH d - Hs differences, suggesting concentrated haplotypic diversity in limited oceanic regions (Fig 2a and b). However, some metazoan group s such as Arthropoda (PH d - Hs = 0), Gastrotricha (-0.17), Placozoa (0.09), and Rotifera (-0.21), exhibited more even haplotype distributions. Protist phyla generally showed slightly higher PH d values and PH d - Hs differences close to zero, indicating moderate population genetic structuration with traceable dominant oceans. Overall, PH d analysis supports the hypothesis of recent geographic expansions or introductions shaping genetic patterns in these marine eukaryotes. Further analysis using the PH d metric to identify dominant source oceans, harboring the majority of that OTU’s haplotypic diversity, and introduced (recipient) oceans revealed clear differences among taxonomic groups. Metazoa, Archaeplastida, and Ochrophyta showed minimal evidence of multi-ocean distributions (Fig. 2C, S8). In contrast, certain metazoan and protist phyla frequently displayed one or more introduced oceans per OTU, indicating recent inter -ocean colonization. For example, Arthropoda, Gastrotri cha, Porifera, Rotifera, and Xenacoelomorpha, often had broader distributions, averaging around one introduced ocean per OTU (Fig. S8). Similarly, protists frequently displayed one or more introduced oceans per OTU, indicating a strong signature of recent inter - ocean introduction, with low haplotypic diversity, in these phyla. However, caution is required when characterizing species as native or introduced, as the evolutionary history of some species remains unclear with the current data (i.e, cryptogenic species) (26). To minimize misclassifications, we only considered OTUs with clear patterns of intra -OTU diversity across oceans [see supplementary materials and methods (15)]. To validate the capacity of our approach to identify genuine introduction events, we examined specific OTUs within the phylum Chordata, where introduction histories are often well -documented (27). Our analysis successfully highlighted several OTUs corresponding to species known for human -mediated introductions or widespread dispersal, such as Acipenser gueldenstaedtii, Botryllus schlosseri, Salmo trutta, Siganus rivulatus, and Styela plicata [see supplementary data (15)]. This reinforced the reliability of our approach in identifying recent introductions from metabarcoding data. Furthermore, the inferred geographic scope of these colonization events differed markedly between groups. In Metazoa and Archaeplastida, most inferred source -sink connections occurred between geographically adjacent oceans and within the same hemisphere, suggesting predominantly regional spread (Fig. S9). In contrast, protist phyla OTUs often showed long -distance dispersal between non -adjacent oceans and different hemispheres, implying that recent introductions of these microorganisms frequently span distant oceanic regions, far exceeding the regional spread seen in larger eukaryotes (28). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Finally, by aggregating dominant and introduced ocean assignments across all phyla, we uncovered broad biogeographic trends that align with known global patterns of species introductions [see supplementary materials and methods (15)]. The highest numbers of introduced informative OTUs were detected in Northern Hemisphere oceans, particularly the Atlantic and North Pacific (Fig. 3), consistent with these regions’ historically high rates of species introductions (27, 29, 30). In contrast, the main source regions (dominant oceans) across taxa were the Mediterranean Sea and Indian Ocean, both noted for their high levels of endemic diversity (31, 32). Overall, these findings underscore that recent colonization, most likely facilitated by anthropogenic activities, have left a distinct imprint on global marine genetic diversity, with identifiable source areas and introduction hotspots explaining the distribution patterns we observe.

Our global COI analysis, leveraging informative OTUs and the innovative PHd metric for assessing population -level genetic diversity, definitively reveals fundamental distinctions in the contemporary biogeography of marine eukaryotes. Metazoa and Archaeplastida generally exhibit significantly restricted ranges and weaker ocean -basin genetic structure compared to protist phyla, underscoring divergent molecular population dynamics across major lineages. These patterns likely reflect the profound impact of contemporary events, notably anthropogenic dispersal (33). While global shipping is a known vector for metazoan introductions (34), our data robustly indicate that protists and smaller metazoans (Arthropoda, Rotifera, Gastrotricha, Placozoa) show heightened susceptibility to long - distance dispersal, evidenced by elevated PHd values with respect other metazoan phyla. Moreover, the lower PH d in these smaller metazoans compared to protist s phyla likely stems from protists rapid generation times (35), or multiple introductions (19), could facilitate the genetic homogenization in newly colonized oceans , reducing the founder effect. Conversely, restricted distributions and weak structure in some small and unicellular Archaeplastida suggest lower transoceanic transport survival, potentially due to light and resource limitations, corroborated by the differential recovery of taxa in ballast water studies. Consistent with this explanation, studies of ballast water frequently recover resilient photosynthetic taxa such as Bacillariophyta and dinoflagellates, which are capable of producing resistant forms (36), whereas groups like Chlorophyta appear underrepresented (37, 38). This suggests that ballast water may be a less effective vector for long-distance dispersal in certain Archaeplastida lineages. Ultimately, our global study establishes a unified framework for interpreting broad biogeographic and population genetic patterns across diverse eukaryotic phyla. The congruent patterns of marine geographic distribution and genetic structure within each group, strongly suggest shared historical processes, such as recent introductions, shaping their contemporary ranges. These findings challenge the traditional view of protists as ubiquitously dispersed entities (1), at the population genetic level, and reveal potentially universal principles governing biogeographic dynamics and community assembly , that transcend organismal size and complexity, particularly in the context of the Anthropocene.

eKOI metabarcoding database .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint To generate the eKOI metabarcoding database we selected 88 marine metabarcoding studies, with publicly accessible raw sequencing data, that had been generated on an Illumina paired-end platform using the cytochrome oxidase subunit I (COI) molecular marker, and employed the following primers: HCO (TAAACTTCAGGGTGACCAAAAAATCA) (39), LCO (GGTCAACAAATCATAAAGATATTGG) (39), jgHCO2198 (TAIACYTCIGGRTGICCRAARAAYCA) (40), mlCOIintF (GGWACWGGWTGAACWGTWTAYCCYCC) (41), and mlCOIintF -XT (GGWACWRGWTGRACWITITAYCCYCC) (42). After obtaining the raw data for each study, we processed each dataset separately according to the protocol outlined in (43). As a first step, primer trimming and demultiplexing were performed with Cutadapt (version 2.8) (44), when necessary. Next, the resulting reads for each sample in each metabarcoding study were processed independently using the DADA2 R package (45), with minor modifications applied depending on the specific metabarcoding study. Chimeric sequences were subsequently removed using the DADA2 removeBimeraDenovo function. After inferring the amplicon sequence variants (ASVs), we stored them in a FASTA file using the following header format: >sequence_ID merged_sample={ 'LOCALITY_ID1': reads, 'LOCALITY_ID2': reads, ... }. A separate metadata file contained information associated with each locality ID (see metadata.csv in the supplementary data). For each sample in the metadata, we then assigned the corresponding ocean or sea. To achieve this, we downloaded shapefile datasets from the “Global Oceans and Seas” database (version 1, accessed 05/11/2024) and “Natural Earth” (accessed 05/11/2024). Using each locality’s latitude and longitude coordinates, we determined and recorded the respective ocean or sea for that sample. For subsequent biogeographic analyses, we utilised the broad ocean regions defined in the “Global Oceans and Seas” database, which groups the world into 11 major seas and oceans. This approach was chosen to avoid excessive subdivision of oceanic regions. All FASTA files generated by each metabarcoding study, in the above format, were then merged into a single file ( eKOI_metabarcoding_database.fasta in the supplementary data) using a custom Python script (1_fastas_combine.py). This script utilises the Biopython library (46) to identify identical ASVs across studies and combine their locality information. Each unique ASV was assigned an identifier in the format eKOIX (where X is a unique number for each ASV). Once the final merged FASTA file was assembled, we checked it for c himeric sequences using the VSEARCH uchime_denovo command (version 2.14.1) (47). Taxonomic assignment of the ASVs was performed using the eKOI taxonomy database (13). For this purpose, we used the script (5_taxonomic_assignation.py) develop in that study. This script created a separate folder for each FASTA file in the working directory. Within each folder, an Excel file was generated containing the taxonomic assignme nt information for each ASV (eKOI_metabarcoding_database_taxonomic_annotations.xlsx in the supplementary data), obtained via the VSEARCH usearch_global command. ASVs with less than 84% similarity to any reference sequence were not considered further. Finally, we generated separate FASTA files for each taxonomic level of interest; in this case, we chose the phylum level. The resulting ASV sequences for ea ch phylum can be downloaded from the “Phyla_data” folder in the supplementary data as eKOI_database.fasta file per phylum. Sequence Data Processing and Analysis .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint OTU Clustering Once FASTA files were generated for each phylum, they were analysed independently. First, all metadata for each ASV was extracted using the script 2_generate_abundances_unique.py. Biopython was used to extract each ASV’s ID and match it with the metadata f ile, producing a DataFrame containing the sequence ID, geographic coordinates, reads and ocean for each record (abundances_unique.csv in the supplementary data per phylum). Because some ASVs were found in multiple localities, duplicate entries were created for each locality where an ASV occurred by appending a suffix “_dupN” (where N is the duplicate number) to the ASV’s ID, yielding a unique identifier for each locality occurrence of that ASV. To avoid biases from multiple sequencing replicates per locality (common in metabarcoding studies), only one ASV entry per locality was retained (duplicate entries with identical ASV ID and coordi nates were removed). After this filtering, a comprehensive FASTA file was compiled, including all ASV sequences (with each locality -specific duplicate as a separate entry). This file was aligned using MAFFT version 7.490 (48), with optimised parameters (--auto, --ep, - -op, --maxiterate, --large) (aligned_sequences_mafft.fasta in the supplementary data per phylum). Operational Taxonomic Units (OTUs) were then delineated from the alignment using the script 3_generate_OTU.py. In particular, VSEARCH clustering at 97% identity was employed to group sequences into OTUs. We set the clustering threshold at 97 % identity (i.e., a 3 % divergence cut -off) because COI barcoding studies of Amoebozoa (43), Cercozoa (49) and diverse animal taxa report similar values (12). Although each lineage can display its own diversification rate, published barcoding gaps generally fall between 2 % and 3 %, so we applied a uniform 3 % threshold across all phyla. Representative ASV (centroids) for each OTU were saved, and a mapping of each ASV to its OTU was generated as a .uc file ( otus.uc in the supplementary data per phylum). The .uc file was subsequently processed to produce a tabular mapping file (otus_mapping.txt in the supplementary data per phylum) that records the assignment of each ASV ID to its corresponding OTU. Molecular and geographic distance Matrix Construction and OTU Filtering After obtaining the OTUs, the relationship between genetic and geographical distances among ASVs was analysed using the R script 4_1_generate_molecular_and_geographic_distances.R. We used the Biostrings R package to convert the phylum-level ASVs alignments into DNAbin objects, suitable for distance calculations. Genetic distance matrices were computed under the K80 model (Kimura 2 -parameter) (50), with gaps (insertions/deletions) excluded using the pairwise.deletion method. In parallel, a matrix of geographical distances between sampling localities was calculated from the coordinate data using the distVincentyEllipsoid function of the geosphere R package. The genetic and geographical distance matrices were both converted to long format and merged into a single data frame, so that each row contained a pairwise genetic distance, the corresponding geographical distance, an d the associated OTU (and phylum). To ensure robust analysis, we first filtered the data to include only OTUs represented by at least four independent ASVs. Next, we further filtered these OTUs, called the informative OTUs, to retain only those in which at least one pair of sequences had a genetic distance ≥ 0.01 and at least one pair had a geographical distance ≥ 1 m. The ratio of the number of total OTUs and filtered informative OTUs per phylum can be downloaded from ratio_total_OTUs_informative_OTUs_by_phylum.csv in the supplementary data, and the ID of the informative OTUs in the file informative_OTUs.txt of each phylum. This .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint two-step filtering removed OTUs with very few sequences or with negligible variation, which could otherwise introduce noise. Metabarcoding studies that uses general primers lack targeted taxonomic resolution and often retrieve fragmented sequences or sequences from multiple species, potentially biasing interpretations (14). Filtering out OTUs with uniformly low molecular divergence could mitigates these issues. While this approach might also exclude some genuine low genetic diversity cases (for example, populations that have undergone bottlenecks or are in decline), it pri marily reduces spurious OTUs arising from methodological biases, improving the reliability of biogeographic patterns interpretations. The results of each analyses of each informative OTUs were compiled into the file informative_OTUs_results.csv, in the supplementary data. In this file, each row corresponds to a informative OTU and each column the results of the analysis or metrics, such as: Phylum (the taxonomic phylum of that OTU), OTU (the OTU identifier, unique within each phylum), Max_Distance (the maximum geographical distance, in metres, between any two localities), Mean_Distance (the mean geographical distance among all locality pairs), Std_Dev_Distance (standard deviation of geographical distances among localities) and Count (number of pairwise ASV comparisons, genetic/geographical distance pairs). Within-OTU Diversity and Tajima’s D To quantify genetic diversity within each informative OTU, we calculated the nucleotide diversity (π). This was done using the previously generated alignments per phylum: the nuc.div function (51) was applied with pairwise.deletion = TRUE to compute π for each OTU using the script 9_calculate_nucleotide_diversity_by_otu.R. The resulting values were recorded in the Nucleotide_Diversity_Pi column of informative_OTUs_results.csv. We also calculated Tajima’s D for each informative OTU to investigate potential signals of selection or demographic history (e.g. recent expansion, bottleneck) using the script 10_calculate_tajimasD_by_otu.R. The tajima.test function (51) was applied to the alignment data for each informative OTU. In theory, a negative Tajima’s D indicates an excess of rare variants consistent with recent population expansion or purifying selection, whereas a positive Tajima’s D suggests a deficiency of rare variants, which may indicate a population bottleneck or balancing selection (52). Given our filtering criteria (which ensured multiple sequences per OTU and some level of genetic divergence), we expected Tajima’s D values to trend negative. Each OTU’s Tajima’s D value was recorded in a Tajima_s_D column, with its associated p -values recorded in P_Value_Normal (normal approximation) and P_Value_Beta (beta distribution approximation) columns. Subsequently, we examined whether nucleotide diversity and Tajima’s D varied significantly among phyla. (For Tajima’s D, this comparison was restricted to OTUs with a Beta p -value < 0.05) We first checked the distribution of each variable: the Shapiro -Wilk test indicated that the data were not normally distributed. Therefore, a non -parametric Kruskal -Wallis test (using the kruskal.test function in R) was employed to test for differences among phyla. A significant Kruskal-Wallis p-value was taken as evidence of variation among phyla for the given metric. In cases of significance, post-hoc pairwise comparisons were performed using Dunn’s test (dunnTest function in R) with a Bonferroni adjustment for multiple testing (53). The results of the pairwise comparisons for each variable were saved to significant_group_comparisons_Dunn_test.csv (available in the Supplementary Data). The file contains the columns Z (Z-Statistic), P_unadj (Uncorrected p-value) and P_adj (Adjusted p-value) for every variable in the phyla pairwise comparison. Phyla with an .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint insufficient number of informative OTUs were excluded from these comparative analyses (Filasterea, Hemichordata, Nibbleridia, Prasinodermophyta and Ustilaginomycotina). Additionally, for infer differences between and within groups, phyla were grouped into three broad categories (Archaeplastida, Metazoa, and all other phyla) to compare the prevalence of significant differences between and within these groups (significant_group_comparations_Dunn_test.csv in the supplementary data). The same statistical approach was applied in later analyses for other variables and metrics of interest, as described below. Read Count Filtering Impact It is common in metabarcoding studies to discard ASVs with very low read counts (below a certain threshold) to avoid spurious ASVs (54). To evaluate how such filtering might affect the recovery of low -abundance haplotypes and OTUs, we characterised the occurrence of rare haplotypes within each OTU. First, for each informative OTU, the total number of ASVs and the total number of reads was recorded in the Num_Sequences and Total_Reads columns, respectively in the informative_OTUs_results.csv file. In addition, separate columns captured the number of reads contributed by each ocean named _reads. We then evaluate the potential loss of data and in reads counts filtering by counting the number of haplotypes with fewer than 10 reads for each informative OTU, recorded in the haplotypes_lt10 column. We also determined how many of these rare haplotypes were supported by at least two independent ASVs (each with <10 reads); this number was recorded in haplotypes_recovered_multiple column. By requiring multiple independent ASVs, we can identify rare haplotypes that are consistently observed, thereby reducing the likelihood of counting sequencing artifacts as true haplotypes. We also assessed whether the low -read filter could inadver tently affect more abundant haplotypes. For each OTU, we counted the number of haplotypes with total reads > 10 that nevertheless contained at least one constituent ASV with <10 reads. This count was recorded in haplotypes_gt10_with_lt10_seq column. Finally, to examine the biogeographic distribution of low -abundance ASVs, we counted how many oceans had a total read count < 10 for each OTU. This value (the number of distinct oceans in which the OTU is represented by fewer than 10 reads) was re corded in the global_lt10 column. Given that metabarcoding studies are typically constrained to specific regions/oceans, any potential “tag jumping” across oceans (a form of cross -sample contamination) can be ruled out at this stage. Overall, these metrics allowed us to understand the extent to which stringent read filtering might exclude genuine but low - frequency haplotypes and the veracity of ASVs with low number of reads. Geographic and genetic distances Linear regression and non -linear analyses were conducted for each informative OTU to evaluate the relationship, strength and statistical significance of the correlation between genetic and geographic distances using the R scripts 5_calculate_linear_models_by_otu.R and 6_calculate_nls_models_by_otu.R, respectively. The pa rameters of the linear model (lm) were calculated for each informative OTU using the dplyr R package, and the results were saved in the file informative_OTUs_results.csv, with the Columns: Intercept_lm (intercept), Slope_lm (slope), R2_lm (coefficient of determination), P_value _lm (p-value of the slope), AIC_lm (Akaike Information Criterion), and BIC_lm (Bayesian Information Criterion). Additionally, a non -linear regression analysis (nls) was performed using the model: .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint “Genetic Distance” = a * log(“Geographical Distance” + 1) + b. Parameters a and b were estimated from predefined initial values. The output included the a_nls (slope parameter), b_nls (vertical offset parameter), RSS_nls (residual sum of squares), R2_nls (coefficient of determination), AIC_nls (Akaike Information Criterion), and BIC_nls (Bayesian Information Criterion). To interpret the differences in AIC and BIC between these two models, categories were established based on the magnitude of the differences (Δ): (1) no evidence to prefer either model (|Δ| ≤ 2), indicating similar fit; (2) moderate evidence (2 7), indicating strong support for one model ( models_comparison_summary.csv file in supplementary data). Based on these results we continued with the non-linear models for the graphical representations, as it better explained the genetic and molecular distance correlations per phylum. Next, the linear and non-linear models were compared with null models to assess whether the relationship between genetic and geographic distances was consistent with a stochastic pattern. The null models were generated by simulating random distributions of genetic and geographic distances using the Python script 7_calculate_null_models_by_otu.py, fitting linear models with the scikit -learn library and optimising non-linear models with SciPy’s curve_fit function. For each OTU, 1,000 simulations were performed. In each simulation, the number of data points equalled the number of pairwise comparisons (count column in informative_OTUs_results.csv file) for that OTU. Geographical distances were sampled uniformly between 0 and the maximum recorded distance ( max_distance column), and genetic distances were sampled uniformly between 0 and 0.03 (the specified OTU cutoff). From these simulations, the mean, median and standard deviation of R^2, AIC and BIC were calculated for comparison with the observed models. Results were saved in informative_OTUs_results.csv file with the following columns: Mean_R2_lm_null, Median_R2_lm_null, Std_R2_lm_null, Mean_AIC_lm_null, Median_AIC_lm_null, Std_AIC_lm_null, Mean_BIC_lm_null, Median_BIC_lm_null, and Std_BIC_lm_null; and for non -linear models: Mean_R2_nls_null, Median_R2_nls_null, Std_R2_nls_null, Mean_AIC_nls_null, Median_AIC_nls_null, Std_AIC_nls_null, Mean_BIC_nls_null, Median_BIC_nls_null, and Std_BIC_nls_null, Mean_a_nls_null, Median_a_nls_null, Std_a_nls_null, Mean_b_nls_null, Median_b_nls_null, and Std_b_nls_null. P -values were calculated to determine whether the observed relationships differed significantly from random expectations, based on the coefficient of determination, and were saved in the columns P_R2_lm_null (linear models) and P_R2_nls_null (non-linear models). To further examine whether a correlation exists between genetic and geographic distances and to assess its consistency across phyla, Mantel tests were performed using the R script 8_calculate_mantel_test_by_otu.R. For each OTU, genetic and geographic dista nce matrices were generated, and the Mantel test was conducted using the mantel function from the vegan package (version 2.6) (55) with 999 permutations, the results were recovered in the columns Mantel_R and P_Value_Mantel. However, the Mantel test can be problematic due to inherent dependencies in the data (56). In this context, the permutations used to calculate the p -value may not adequately reflect the expected null distribution under true independence; therefore, this results should be interpreted with caution. After completing these analyses, only OTUs with a Mantel test p -value < 0.05 were retained, and comparisons between phyla were conducted as described above. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Biogeography Based on Haplotype Networks Haplotype network construction and diversity metrics Unique haplotypes were first identified from the alignment for each informative OTU using the haplotype function (51). The number of haplotypes per informative OTU was then recorded in a column named num_haplotypes. Next, using the haploNet function (51), haplotype networks based on Minimum Spanning Networks (MSNs) were constructed for each informative OTU and used in subsequent analyses. We then calculate the Haplotype diversity (Hd), branch diversity (Bd) and their combined metric (Hbd) (23) using R script 11_calculate_haplotype_diversity_by_otu.R. We saved the results in columns named Hd, Bd and HBd, respectively, in the informative_OTUs_results.csv file. Additionally, the maximum value of the squared haplotype frequencies per OTU was recorded in a column Efh2_Hd. The above metrics treat each ASV equally, regardless of its abundance (i.e. they are based only on ASV presence/absence within each haplotype). To account for the relative abundance of each ASV (as a proxy for its prevalence in the sample), we incorporated the associated read counts. Specifically, we used two metrics (1) the total reads count for each ASV and (2) a logarithmic transformation, (log(reads+1)), using the R script 12_calculate_weighted_haplotype_diversity_by_otu.R. By weighting sequences according to their total read counts and a logarithmic transformation, haplotypes with higher abundance have a greater influence on the calculated diversity metrics. Therefore, if a sequence-based diversity metric exceeds its read-weighted counterpart, this indicates that rarer haplotypes disproportionately influence the haplotypic diversity (i.e. diversity is overestimated in the unweighted case). This approach is particularly useful for communities with high variability in abundances, as it allows patterns of d ominance or inequality to emerge. However, it may underestimate diversity in OTUs that contain many low-abundance haplotypes. The results of these abundance -weighted calculations were recorded in the following columns: Hd_reads and Hd_log_reads (haplotype diversity); Bd_reads and Bd_log_reads (branch diversity); Hbd_reads and Hbd_log_reads (combined haplotype and branch diversity); and Efh2_reads and Efh2_log_reads (maximum squared haplotype frequencies). Oceanic genetic structure After calculating haplotype diversity metrics, we evaluated genetic structuring between and within populations (each corresponding to a different ocean). Analysis of Molecular Variance (AMOVA) was performed for each OTU using R script 13_calculate_amova_by_otu.R. OTUs were selected for AMOVA based on the following criteria: (i) Presence in at least two different oceans; (ii) At least two populations; and (iii) At least two distinct haplotypes within those populations. We used the popr.amova function to perform the AMOVA. Variance components and sums of squares were calculated for between-population and within-population comparisons. The results were saved in the informative_OTUs_results.csv file, in columns named Between_samples_SumSq, Within_samples_SumSq, Total_SumSq, Between_samples_Sigma, Within_samples_Sigma and Total_Sigma. These metrics were then compared across phyla as described above. Because the sequence length is relatively short (approximately 313 bp), the AMOVA results should be interpreted with caution. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint To mitigate potential biases introduced by short sequences, we evaluated genetic structuring among oceans using the haplotype frequencies and distributions for each informative OTU. The haplotypic diversity metrics total diversity (Ht), within-population diversity (Hs), genetic differentiation coefficient (Fst) and gene flow (Nm) were calculated for each OTU using script 14_calculate_Hs_fst_nm_by_otu.R. These four metrics were stored in columns Ht, Hs, Fst, and Nm; and comparisons between phyla were conducted as described above. A modified haplotype diversity metric (PH d) was created to infer the distribution of haplotypes among populations (here defined by different oceans) using script 15_calculate_PHd_and_global_ocean_dominance_by_otu.R. This new metric incorporates the relative distribution of haplotypes across populations. By using weighted relative haplotype frequencies, PH d quantifies the contribution of population -level differences to the genetic diversity of each OTU. The classical haplotype diversity formula was adapted to account for the distribution of haplotypes across populations. In this adaptation, the relative frequency of each haplotype is redefined as its frequency in a given population divided by its total frequency in the OTU. The modified PHd is therefore calculated as: 𝑃𝐻𝑑 = 1 − ∑ 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛(𝑜𝑐𝑒𝑎𝑛𝑠) (𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 ℎ𝑎𝑝𝑙𝑜𝑡𝑦𝑝𝑒𝑠 𝑖𝑛 𝑡ℎ𝑒 𝑝𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑇𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 ℎ𝑎𝑝𝑙𝑜𝑡𝑦𝑝𝑒𝑠 𝑖𝑛 𝑡ℎ𝑒 𝑂𝑇𝑈 ) ^2 The PHd values were saved in a column named PHd. The relative haplotype frequencies for each ocean were stored in columns named _sqfreq. These PHd metrics were then compared across phyla using the same procedure as for the other diversity metrics. Furthemore we compare PHd with Hs metric to helps to identify whether the populations of an informative OTU are genetically structured or relatively homogeneous across oceans. If PHd exceeds Hs (a positive difference), this indicates that the populations have distinct genetic compositions and are partially isolated from one another, reflecting significant genetic structuring and no single population dominating the haplotype diversity. Conversely , if PH d is less than H s (a negative difference), one or a few populations account for most of the global haplotypic diversity in the OTU, suggesting that these populations dominate the haplotypic diversity. For each informative OTU, we identified the ocean with the highest haplotype relative contribution (the maximum sqfreq value) as a proxy for the putative population of source and vice versa for recent introduction populations. However, this analysis is considered a preliminary proxy to guide further investigation. To distinguish likely source populations from recent introductions, we calculated ratios of relative haplotype frequencies between oceans using script 20_introductions_vs_dominant_oceans_byOTU.R. Specifically, within each OTU we used the ocean with the highest relative frequency as the reference (ratio = 1), and we computed for each other ocean the ratio of the reference frequency to that ocean’s frequency. We then examined the distribution of thes e ratios across all informative OTUs to detect shifts toward lower values. A notable change in the distribution occurred at a ratio of 0.25. Consequently, for each OTU we classified any ocean with a ratio below 0.25 as a “recently introduced ocean/populati on,” since such oceans have very low relative haplotype frequency compared to the highest haplotypic frequency ocean. Furthermore, to identify a clear source ocean for each informative OTU, we applied additional filtering criteria. We selected OTUs in whic h exactly one ocean had a ratio of 1.0 (i.e. a single ocean had the maximum relative frequency) and excluded informative OTUs where multiple oceans shared this maximum frequency, since that .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint would preclude a clear source. We also excluded informative OTUs in which any other ocean had a ratio exceeding 0.5, to avoid ambiguity. Although these thresholds are somewhat arbitrary, they provide a more restrictive criterion that allows a more reliable inference of the putative source ocean for the remaining informative OTUs. We create a pie chart graph of the proportion of dominant oceans per phylum ancestral_oceans.pdf. Furthermore, we combine the number of OTUs dominant and introduced oceans in the file dominant_introduced_oceans.csv. To validate this proxy for source (dominant) and introduced ocean populations, we used phylum Chordata taxonomically annotated informative OTUs as a positive control. We included only those informative OTUs that had at least one ocean classified as “recent ly introduced”. For each selected informative OTU, the first amplicon sequence variant (ASV) was compared to entries in the NCBI database using BLAST (57). We retained those OTUs for which the BLAST hit showed ≥ 97% sequence identity to a species in the database, ensuring that species -level identification was feasible. Taxonomic annotation

for these OTUs are recorded in file introduced_chordata_otus.csv. For each informative OTU, we inferred connectivity between different oceans from its haplotype network to evaluate genetic dispersal and exchange. This analysis was performed using the script 16_generate_ocean_connectivity_by_otu.R. In each haplotype network, nodes represent haplotypes and edges represent connections between them. We annotated each edge with the ocean of origin for each haplotype in the pair. Only connections linking haplotypes from different oceans were retained; intra -oceanic connections were discarded because they do not inform inter-ocean connectivity. We then aggregated the data to count the number of connections between each pair of oceans. For each informative OTU, we recorded the total number of oceanic haplotypic connections in the column Total_Connections and the number of connections per ocean in columns named _connections. Then, we classified oceanic pairs as “nearby” if they are geographically adjacent and “distant” if they are geographically separated using the script 17_calculate_ocean_conectivity.R, the number of connections were recovered in the columns nearby_oceans_connections and distant_oceans_connections. OTUs in which all haplotypes occur in the same ocean therefore have only “nearby” connectivity and no “distant” connectivity. We also characterised intra -hemispheric and inter - hemispheric connectivity, focusing on connections between oceans in the Northe rn and Southern Hemispheres, recovered in the columns Inter-Hemisferica and Intra- Hemisferica. For each phylum, we generated network graphs illustrating connections among oceans using the script 18_oceanic_network_by_phyla.py and 19_global_oceanic_network.py. In these graphs, each node represents an ocean and each edge represents the total number of connections between two oceans, weighted by relative strength. We summed all connections across OTUs for each phylum and applied a log1p transformation (the natural logarithm of one plus this sum) to normalise the data. We then applied the Louvain community detection algorithm to identify clusters (communities) within each network, colouring nodes by their community membership. The networks were constructed using the networkx library, with edge weights reflecting connection strength the network per phylum were recovered as network_.pdf and .jpg. Finally, we analysed correlations among the different calculated variables for each informative OTU. We computed correlation matrices and associated p -values using Pearson’s correlation coefficient, generating a table listing only those correlations that were statistically significant (p < 0.05) significant_correlations.csv file. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint

and Notes 1. J. Whitfield, Is Everything Everywhere? Science 310, 960–961 (2005). 2. L. G. M. B. Becking, Geobiologie of inleiding tot de milieukunde (W.P. Van Stockum & Zoon, 1934). 3. E. Pinseel, K. Sabbe, E. Verleyen, W. Vyverman, A New Dawn for Protist Biogeography. Global Ecology and Biogeography 33, e13925 (2024). 4. C. de Vargas, S. Audic, N. Henry, J. Decelle, F. Mahé, R. Logares, E. Lara, C. Berney, N. Le Bescot, I. Probert, M. Carmichael, J. Poulain, S. Romac, S. Colin, J. -M. Aury, L. Bittner, S. Chaffron, M. Dunthorn, S. Engelen, O. Flegontova, L. Guidi, A. Ho rák, O. Jaillon, G. Lima -Mendez, J. Lukeš, S. Malviya, R. Morard, M. Mulot, E. Scalco, R. Siano, F. Vincent, A. Zingone, C. Dimier, M. Picheral, S. Searson, S. Kandels-Lewis, Tara Oceans Coordinators, S. G. Acinas, P. Bork, C. Bowler, G. Gorsky, N. Grimsley, P. Hingamp, D. Iudicone, F. Not, H. Ogata, S. Pesant, J. Raes, M. E. Sieracki, S. Speich, L. Stemmann, S. Sunagawa, J. Weissenbach, P. Wincker, E. Karsenti, Eukaryotic plankton diversity in the sunlit ocean. Science 348, 1261605 (2015). 5. G. Sommeria-Klein, R. Watteaux, F. M. Ibarbalz, J. J. Pierella Karlusich, D. Iudicone, C. Bowler, H. Morlon, Global drivers of eukaryotic plankton biogeography in the sunlit ocean. Science 374, 594– 599 (2021). 6. P. C. Junger, H. Sarmento, C. R. Giner, M. Mestre, M. Sebastián, X. A. G. Morán, J. Arístegui, S. Agustí, C. M. Duarte, S. G. Acinas, R. Massana, J. M. Gasol, R. Logares, Global biogeography of the smallest plankton across ocean depths. Sci. Adv. 9, eadg9763 (2023). 7. G. S. Bradburd, P. L. Ralph, Spatial Population Genetics: It’s About Time. Annual Review of Ecology, Evolution, and Systematics 50, 427–449 (2019). 8. L. J. Rissler, Union of phylogeography and landscape genetics. Proceedings of the National Academy of Sciences 113, 8079–8086 (2016). 9. S. Manel, M. K. Schwartz, G. Luikart, P. Taberlet, Landscape genetics: combining landscape ecology and population genetics. Trends in Ecology & Evolution 18, 189–197 (2003). 10. J. C. Avise, J. Arnold, R. M. Ball, E. Bermingham, T. Lamb, J. E. Neigel, C. A. Reeb, N. C. Saunders, INTRASPECIFIC PHYLOGEOGRAPHY: The Mitochondrial DNA Bridge Between Population Genetics and Systematics. Annual Review of Ecology, Evolution, and Systematics 18, 489–522 (1987). 11. P. Taberlet, E. Coissac, M. Hajibabaei, L. H. Rieseberg, Environmental DNA: ENVIRONMENTAL DNA. Molecular Ecology 21, 1789–1793 (2012). 12. P. D. N. Hebert, S. Ratnasingham, J. R. de Waard, Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species. Proc. R. Soc. Lond. B 270 (2003). 13. R. González-Miguéns, A. Gàlvez-Morante, M. Skamnelou, M. Antó, E. Casacuberta, D. J. Richter, D. Vaulot, J. del Campo, I. Ruiz -Trillo, A Novel Taxonomic Database for eukaryotic Mitochondrial Cytochrome Oxidase subunit I Gene (eKOI): Enhancing taxonomi c resolution at community-level in metabarcoding analyses. bioRxiv [Preprint] (2024). https://doi.org/10.1101/2024.12.05.626972. 14. M. Couton, F. Viard, F. Altermatt, Opportunities and inherent limits of using environmental DNA for population genetics. Environmental DNA 5, 1048–1064 (2023). 15. Materials and methods and additional supplementary text sections are available as supplementary materials. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint 16. K. J. Gaston, K. J. Gaston, The Structure and Dynamics of Geographic Ranges (Oxford University Press, Oxford, New York, 2003)Oxford Series in Ecology and Evolution. 17. N. D. Tsutsui, A. V. Suarez, D. A. Holway, T. J. Case, Reduced genetic variation and the success of an invasive species. Proceedings of the National Academy of Sciences 97, 5948–5953 (2000). 18. C. E. Lee, Evolutionary genetics of invasive species. Trends in Ecology & Evolution 17, 386–391 (2002). 19. J. Roman, J. A. Darling, Paradox lost: genetic diversity and the success of aquatic invasions. Trends in Ecology & Evolution 22, 454–464 (2007). 20. D. Posada, K. A. Crandall, Intraspecific gene genealogies: trees grafting into networks. Trends in Ecology & Evolution 16, 37–45 (2001). 21. D. H. Huson, D. Bryant, Application of Phylogenetic Networks in Evolutionary Studies. Molecular Biology and Evolution 23, 254–267 (2006). 22. M. Nei, Molecular Evolutionary Genetics (Columbia University Press, New York, 1987). 23. E. Garcia, D. Wright, R. Gatins, M. B. Roberts, H. T. Pinheiro, E. Salas, J. -Y. Chen, J. R. Winnikoff, G. Bernardi, Haplotype network branch diversity, a new metric combining genetic and topological diversity to compare the complexity of haplotype networks. PLOS ONE 16, e0251878 (2021). 24. B. S. Weir, C. C. Cockerham, Estimating F -Statistics for the Analysis of Population Structure. Evolution 38, 1358–1370 (1984). 25. M. Slatkin, Gene Flow in Natural Populations. Annual Review of Ecology, Evolution, and Systematics 16, 393–430 (1985). 26. J. T. Carlton, Biological Invasions and Cryptogenic Species. Ecology 77, 1653–1655 (1996). 27. J. L. Molnar, R. L. Gamboa, C. Revenga, M. D. Spalding, Assessing the global threat of invasive species to marine biodiversity. Frontiers in Ecology and the Environment 6, 485–492 (2008). 28. B. W. Bowen, M. R. Gaither, J. D. DiBattista, M. Iacchei, K. R. Andrews, W. S. Grant, R. J. Toonen, J. C. Briggs, Comparative phylogeography of the ocean planet. Proceedings of the National Academy of Sciences 113, 7962–7969 (2016). 29. S. Katsanevakis, I. Wallentinus, A. Zenetos, E. Leppäkoski, M. E. Çinar, B. Oztürk, M. Grabowski, D. Golani, A. C. Cardoso, Impacts of invasive alien marine species on ecosystem services and biodiversity: a pan-European review. AI 9, 391–423 (2014). 30. N. Streftaris, A. Z. & E. Papathanassiou, “Globalisation in Marine Ecosystems: The Story of Non - Indigenous Marine Species Across European Seas” in Oceanography and Marine Biology (CRC Press, 2005). 31. J. D. DiBattista, J. Howard Choat, M. R. Gaither, J. -P. A. Hobbs, D. F. Lozano -Cortés, R. F. Myers, G. Paulay, L. A. Rocha, R. J. Toonen, M. W. Westneat, M. L. Berumen, On the origin of endemic species in the Red Sea. Journal of Biogeography 43, 13–30 (2016). 32. M. Coll, C. Piroddi, J. Steenbeek, K. Kaschner, F. B. R. Lasram, J. Aguzzi, E. Ballesteros, C. N. Bianchi, J. Corbera, T. Dailianis, R. Danovaro, M. Estrada, C. Froglia, B. S. Galil, J. M. Gasol, R. Gertwagen, J. Gil, F. Guilhaumon, K. Kesner -Reyes, M.-S. Kitsos, A. Koukouras, N. Lampadariou, E. Laxamana, C. M. L.-F. de la Cuadra, H. K. Lotze, D. Martin, D. Mouillot, D. Oro, S. Raicevich, J. Rius-Barile, J. I. Saiz -Salinas, C. S. Vicente, S. Somot, J. Templado, X. Turon, D. Vafidis, R. Villanueva, E. V oultsiadou, The Biodiversity of the Mediterranean Sea: Estimates, Patterns, and Threats. PLOS ONE 5, e11842 (2010). .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint 33. C. Capinha, F. Essl, H. Seebens, D. Moser, H. M. Pereira, The dispersal of alien species redefines biogeography in the Anthropocene. Science 348, 1248–1251 (2015). 34. H. Seebens, N. Schwartz, P. J. Schupp, B. Blasius, Predicting the spread of marine species introduced by global shipping. Proceedings of the National Academy of Sciences 113, 5646–5651 (2016). 35. S. P. Lawler, P. J. Morin, Food Web Architecture and Population Dynamics in Laboratory Microcosms of Protists. The American Naturalist 141, 675–686 (1993). 36. A. McMinn, A. Martin, Dark survival in a warming world. Proceedings of the Royal Society B: Biological Sciences 280, 20122909 (2013). 37. H. Wu, C. Shen, Q. Wang, R. B. Aronson, C. Chen, J. Xue, Survivorship characteristics and adaptive mechanisms of phytoplankton assemblages in ballast water. J. Ocean. Limnol. 37, 580–588 (2019). 38. A. Ardura, J. Rick, J. L. Martinez, A. Zaiko, E. Garcia -Vazquez, Stress resistance for unraveling potential biopollutants. Insights from ballast water community analysis through DNA. Marine Pollution Bulletin 163, 111935 (2021). 39. O. Folmer, M. Black, W. Hoeh, R. Lutz, R. Vrijenhoek, DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol 3, 294–299 (1994). 40. J. Geller, C. Meyer, M. Parker, H. Hawk, Redesign of PCR primers for mitochondrial cytochrome c oxidase subunit I for marine invertebrates and application in all-taxa biotic surveys. Molecular Ecology Resources 13, 851–861 (2013). 41. M. Leray, J. Y. Yang, C. P. Meyer, S. C. Mills, N. Agudelo, V. Ranwez, J. T. Boehm, R. J. Machida, A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for characterizing coral reef fish gut contents. Frontiers in Zoology 10, 34 (2013). 42. O. S. Wangensteen, C. Palacín, M. Guardiola, X. Turon, DNA metabarcoding of littoral hard -bottom communities: high diversity and database gaps revealed by two molecular markers. PeerJ 6, e4705 (2018). 43. R. González-Miguéns, E. Cano, A. Guillén-Oterino, A. Quesada, D. J. G. Lahr, D. Tenorio-Rodríguez, D. de Salvador-Velasco, D. Velázquez, M. I. Carrasco -Braganza, R. T. Patterson, E. Lara, D. Singer, A needle in a haystack: A new metabarcoding approach to survey diversity at the species level of Arcellinida (Amoebozoa: Tubulinea). Molecular Ecology Resources n/a (2023). 44. M. Martin, Cutadapt removes adapter sequences from high -throughput sequencing reads. EMBnet.journal 17, 10–12 (2011). 45. B. J. Callahan, P. J. McMurdie, M. J. Rosen, A. W. Han, A. J. A. Johnson, S. P. Holmes, DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13, 581–583 (2016). 46. P. J. A. Cock, T. Antao, J. T. Chang, B. A. Chapman, C. J. Cox, A. Dalke, I. Friedberg, T. Hamelryck, F. Kauff, B. Wilczynski, M. J. L. de Hoon, Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics 25, 1422–1423 (2009). 47. T. Rognes, T. Flouri, B. Nichols, C. Quince, F. Mahé, VSEARCH: a versatile open source tool for metagenomics. PeerJ 4, e2584 (2016). 48. K. Katoh, K. Misawa, K. Kuma, T. Miyata, MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Research 30, 3059–3066 (2002). 49. R. González-Miguéns, C. Soler-Zamora, F. Useros, S. Nogal- Prata, C. Berney, A. Blanco-Rotea, M. I. Carrasco -Braganza, D. Salvador -Velasco, A. Guillén - Oterino, D. Tenorio -Rodríguez, D. Velázquez, T. J. Heger, I. Sanmartín, E. Lara, Cyphoderia ampulla (Cyp hoderiidae: Rhizaria), a tale .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint of freshwater sailors: The causes and consequences of ecological transitions through the salinity barrier in a family of benthic protists. Molecular Ecology (2022). 50. M. Kimura, A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol 16, 111–120 (1980). 51. E. Paradis, pegas: an R package for population genetics with an integrated –modular approach. Bioinformatics 26, 419–420 (2010). 52. F. Tajima, Statistical Method for Testing the Neutral Mutation Hypothesis by DNA Polymorphism. Genetics 123, 585–595 (1989). 53. O. J. Dunn, Multiple Comparisons Among Means. Journal of the American Statistical Association 56, 52–64 (1961). 54. N. A. Bokulich, S. Subramanian, J. J. Faith, D. Gevers, J. I. Gordon, R. Knight, D. A. Mills, J. G. Caporaso, Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. Nat Methods 10, 57–59 (2013). 55. P. Dixon, VEGAN, a package of R functions for community ecology. Journal of Vegetation Science 14, 927–930 (2003). 56. P. Legendre, M.-J. Fortin, D. Borcard, Should the Mantel test be used in spatial analysis? Methods in Ecology and Evolution 6, 1239–1247 (2015). 57. S. F. Altschul, W. Gish, W. Miller, E. W. Myers, D. J. Lipman, Basic local alignment search tool. Journal of Molecular Biology 215, 403–410 (1990). Acknowledgments: we express our gratitude to A. Berlinches, E. Lara, G. Bercedo and M. Villar-DePablo and the MCG lab group for helpful discussions. Work funded by the European Union (ERC, MISSINGRELATIVES, 101097659 ). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them . We also acknowledge support to Departament de Recerca i Universitats de la Generalitat de Catalunya (exp. 2021 SGR 00751) and support by PIE - 202120E047- Conexiones-Life, and MICIU/AEI /10.13039/501100011033 - FSE+ (JDC2023-050439-I). Author contributions: Conceptualization: RGM, IRT Methodology: RGM Investigation: RGM Visualization: RGM Funding acquisition: RGM, IRT Project administration: RGM Supervision: EC, IRT Database maintenance: AGM, GT, CB .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Writing – original draft: RGM, IRT Writing – review & editing: RGM, AGM, GT, CB, EC, IRT Competing interests: Authors declare that they have no competing interests. Data and materials availability: All data, code, results, including eKOI metabarcoding database, are available via Figshare (10.6084/m9.figshare.29144645). All other data are available in the manuscript or the supplementary materials. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. 1. eKOI metabarcoding database . ( A) Map showing the locations of samples included in the database. The size and color of the points represent the log(total reads) for each sample. To the left of the map, a scatter plot displays the relationship between log(ASVs) and log(reads) for the sampled localities, colored according to log(reads). (B) rose plots representing the proportion of taxonomic annotated ASV as log(ASV) on the inner layer and log(reads) on the outer layer in the different eukaryotic phyla used in the analysis. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. 2. Population structure between eukaryotic phyla. (A) Boxplots representing the PHd values for each OTU between phyla, where values close to 0 indicate low structuring in terms of intra-OTU haplotype diversity. Red circles with crosses represent the median Hs values for each phyla. The number near each phyla name represents the number of OTU recovered after filtration. The numbers at the top indicate the difference between PHd and H s. ( B) PHd values grouped between different phyla. ( C) Number of recent introduced oceans per OTU grouped between different phyla. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. 3. Number of recently introduced and dominant informative OTUs per ocean. (A) Boxplots showing the number of informative OTUs considered invasive in each ocean and (B) the number of informative OTUs per ocean identified as their dominant/sink. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S1. Rare haplotypes (< 10 reads). Boxplot represents the proportion of haplotypes with less than 10 total reads. The median is indicated by black diamonds, and the box colors correspond to different taxonomic groups. Blue circles represent the median proportion of haplotypes with more than 10 reads but composed of at least one ASV with less than 10 reads. Red squares represent the median proportion of oceans composed of less than 10 reads. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S2. Nucleotide diversity metrics. (A) Boxplot representing the values of Tajima’s D for OTUs grouped by phyla, with p -values less than 0.05. ( B) Boxplot of nucleotide diversity values for OTUs grouped by phyla. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S3. Geographic distances . The boxplots represents the maximum distribution between two localities per informative OTU. The red circles represents the median of the mean geographic distances between all the localities for each informative OTU. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S4. Parameter “a” of the nonlinear (exponential) model . Boxplots representing the values of parameter “a” from the non -linear models for each informative OTU, relating genetic and geographic distances. Orange points represent null models generated for each phyla, with squares marking the median values. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S5. Significant Mantel test . Boxplots representing the R values from the Mantel test, including only informative OTUs with a p-value less than 0.05. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S6. Haplotypic diversity metrics . (A) Boxplot representing haplotype diversity values, (B) branch diversity, and (C) haplotype-branch diversity, all based on the number of sequences. In all three cases, the median considering the number of sequences is represented by diamonds, the median based on total reads is represented by squares, and the median using log(reads) is represented by circles. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S7. Population structure. (A) Boxplot representing the Fst values, and (B) boxplot of Hs values per informative OTU grouped by phyla. Circles represent the median of Nm values, and squares represent the medians for Ht values. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S8. Characterization dominant and introduced ocean per OTU. (A) Distribution of haplotype relative frequency ratios between oceans of all informative OTUs. ( B) number of introduced oceans per informative OTU grouped by phyla. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint Fig. S9. Ratios of haplotype oceanic connections. Bar plot representing the composition of the ratios for the different types of ocean connections, based on the haplotypic networks, for each phylum. The left side of the bar plot illustrates the proportion of connections between inter- versus intra-hemispheric oceans, while the right side shows the proportion of connections between separated versus adjacent oceans. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted June 2, 2025. ; https://doi.org/10.1101/2025.06.02.657477doi: bioRxiv preprint