Pan-Arctic plankton community structure and its global connectivity

The Arctic Ocean (AO) is going through a rapid change due to anthropogenic warming, marked by the alarming rate of decline in the extent and volume of sea ice. Some models predict that by 2050 the AO will be ice-free in summer (Notz and SIMIP Community, 2020). Planktonic microorganisms are highly sensitive to climatic forcing (Hays et al., 2005); their population dynamics in the AO are tightly linked to sea ice-ocean interactions (Ardyna and Arrigo, 2020). In conjunction with changes in the Arctic icescape, major changes to the oceanography of the region are underway, from processes affecting local mixing and stratification to basin scale circulation (Timmermans and Marshall, 2020). The consequences of these effects on unicellular eukaryotic plankton are complex to delineate. Although future trajectories have been proposed for phytoplankton in general (Ardyna and Arrigo, 2020), ecological knowledge is still lacking at a fine taxonomic level, making predictions and the detection of shifts in these communities difficult.

In addition to environmental variations, plankton communities in the AO are influenced by the influx of Atlantic and Pacific waters that carry their physical and chemical signatures, as well as subpolar and temperate species, into the AO (Wassmann et al., 2015; Ardyna and Arrigo, 2020; Csapó et al., 2021). The transport of plankton by ocean currents (hereafter bioadvection) is now recognized as an important mechanism underlying the observed poleward expansions of some phytoplankton and zooplankton populations in the Barents Sea (Wassmann et al., 2019; Oziel et al., 2020; Greco et al., 2022). Unicellular plankton have high dispersal potential and broad geographic ranges, aided by high abundances, asexual reproduction and—in some cases—the ability to form dormant stages (Finlay, 2002). The success of advected species in establishing populations in the AO, however, depends on their ability to survive and reproduce under the extreme seasonality in day length and sea ice cover, as well as low temperatures and variable salinity. Thus, locally adapted species ought to outcompete or even take advantage of advected organisms from lower latitudes as a food source (Wassmann et al., 2019). As the warming of North Atlantic water continues and the sea ice volume that exits the Fram Strait decreases, the transport of Atlantic water into the AO is expected to intensify (Wang et al., 2020). In the Bering Strait, in situ mooring data already have documented a long-term increase (1990–2015) in the annual mean transport of Pacific water into the Arctic (Woodgate, 2018). Because the volume of Atlantic water that enters the AO through the Fram Strait and Barents Sea openings is about 10 times larger than the volume of Pacific water that enters via the shallow and narrow Bering Strait (Timmermans and Marshall, 2020), the influx of Atlantic taxa is expected to be higher.

How changes in the current balance between bioadvection and local environmental filtering processes within the AO will affect plankton communities of microbial eukaryotes is a relevant and timely question. Efforts to compile published data sets from different regions of the AO have provided initial pan-Arctic inventories of the diversity of this phylogenetically rich group of microorganisms (Lovejoy et al., 2006; Bluhm et al., 2011; Poulin et al., 2011; CAFF, 2017; Rämä et al., 2017; Stoecker and Lavrentyev, 2018), but only a few of these studies have examined Arctic taxa in a global context. The biogeographic structure of these communities beyond the AO remains elusive. Based on microscopy records, the conclusion has been drawn that Arctic dinoflagellates form a coherent Arctic-boreal metacommunity, and that tropical-boreal and Antarctic-tropical-boreal species detected in the AO are allochthonous (Okolodkov and Dodge, 1996). An Arctic-boreal metacommunity was also suggested for the diatom order Thalassiosirales based on a molecular analysis (Luddington et al., 2016). A recent analysis of diatom communities from five Svalbard fjords showed that Arctic genotypes display a variety of biogeographic patterns, ranging from endemic to cosmopolitan distributions (Šupraha et al., 2022). These observations for selected taxa highlight the need for a more extended analysis of biogeographic patterns—on a pan-Arctic scale and across a wider range of Arctic taxa—to better understand the connectivity of AO communities with other oceans and the biogeographic processes that shape them. Such analysis could also help identify target species for monitoring and assessments of biotic responses to climate change.

Here, we have leveraged the globally collected Tara Oceans 18S rRNA gene amplicon sequencing data to study biogeographic patterns of operational taxonomic units (OTUs) from 75 deep-branching lineages that were recovered from the AO. OTUs were grouped into 5 biogeographic categories based on a statistical analysis of their abundance and occupancy patterns within and beyond the AO. These include Arctic/Polar indicators, globally ubiquitous OTUs, non-Polar indicators, and a group of less abundant OTUs that was further divided into a group that was only detected in AO samples (Arctic-exclusives) and a group of background OTUs with a broader range of distributions. Considering these groups, we asked the following questions: (1) What is the relative contribution of each biogeographic category to the contemporary pan-Arctic community and do taxa with high biogeographic specificity to the AO account for a higher proportion of AO biodiversity? (2) Do the 75 deep-branching lineages differ in their biogeographic structure and are there lineages that are more enriched in a specific biogeographic category? (3) On a pan-Arctic scale, are taxa that show geographical preference to low latitudes more likely to be detected in the Arctic gateways compared to other regions of the AO?

Data were collected during the Tara Oceans expedition; overviews of the expedition are provided in Sunagawa et al. (2020) and references therein. The analyses presented here are based on samples collected at 89 globally distributed stations, 18 of them located in the AO (TARA stations 158–209; Figure S1; Table S1). Here we provide additional details on the sampling program in the AO and the specific analyses performed in this study.

In June through October 2013, the Tara Oceans expedition circumnavigated the AO (Figures 1A and S1) using its standardized protocols to study plankton ecosystems (Pesant et al., 2015; Alberti et al., 2017). Contextual hydrographic parameters used here were obtained from daily deployments of a CTD rosette package that provided vertical profiles of temperature, practical salinity, chlorophyll fluorescence, and nitrate (Picheral et al., 2014; Pesant et al., 2015). Surface chlorophyll a concentrations were derived from fluorescence measurements (Wet Labs ECO sensor) that were calibrated against chlorophyll a concentrations obtained from high pressure liquid chromatography analysis. In situ measurements were augmented with remote sensing data of surface chlorophyll a and surface ice concentrations to produce dynamic maps (based on a 250-km latitudinal band) of their variations in the region and time of sampling (Figure 1). Satellite-derived chlorophyll a and sea ice concentration data were obtained, respectively, from the European Space Agency GlobColour project (http://www.globcolour.info) and the National Snow and Ice Data Center (https://nsidc.org). We used bottom depths from the National Oceanic and Atmospheric Administration (NOAA, USA) at 4-min spatial resolution from the R-package “marmap” (version 1.0.5). Additional environmental parameters used in the analysis included the mixed layer depth (MLD) that was calculated from our own potential density profiles following de Boyer Montégut et al. (2004); concentrations of iron, derived from the ECCO2-DARWIN ocean model (Menemenlis et al., 2008); and daylength, calculated from the astronomical applications of the U.S. Naval Observatory (https://aa.usno.navy.mil/data/index) at the sampling date and location. Tara Oceans expedition full contextual data are available at https://doi.org/10.1594/PANGAEA.875582.

Samples for DNA metabarcoding analyses were collected with Niskin bottles from 18 stations that encompass different continental shelf types (inflow, interior and outflow shelves; Carmack and Wassmann, 2006), hydrographic and nutrient regimes, and periods of melt and ice formation (Figures 1 and S2; a list of AO stations indicating sampling dates, locations, station depths and daylength is provided in Table S1). Most samples were collected in open water over Arctic shelves, with some stations coinciding with the marginal ice zone (approximately 30% sea ice cover) and slope waters. The sampling protocol was the same as for the global ocean (Pesant et al., 2015), including surface (5 m) and subsurface chlorophyll maxima (SCM) samples when the latter was detected with a fluorometer mounted on the CTD-rosette. Because of the low number of SCM samples from the AO (8 stations), we chose to focus the main analyses of this paper on surface samples. Seawater samples for DNA analysis were size-fractionated (Pesant et al., 2015). For the purpose of this study we have used the integrative fraction 0.8–2000 µm because it retrieves most functional groups and allows the characterization of the community across a wide range of cell sizes. Filters were immediately deep-frozen in liquid nitrogen for later DNA extraction and sequencing (Alberti et al., 2017).

Our analysis focuses largely on the pan-Arctic assemblage, but to address question #3 (Introduction) stations were grouped into 5 Arctic regions following the framework of the Arctic Marine Biodiversity Monitoring Plan (Gill et al., 2011; CAFF, 2017; Figure 1A): (1) Atlantic inflow shelf, which includes the Greenland and Barents seas (TARA stations 158, 163, 168 and 175); (2) inner Siberian shelves, which include the Kara, Laptev and East Siberian seas (TARA stations 173, 178, 180, 188 and 189); (3) Pacific inflow shelf, which includes the Chukchi Sea (TARA stations 191, 193, 194) and one station on the Beaufort Shelf (TARA station 196); (4) Canadian Archipelago/western Baffin Bay, a region characterized by a net outflow of modified Pacific and Arctic water (TARA stations 201 and 205); and (5) eastern Baffin Bay/Greenland fjord, an area influenced by the northward flowing West Greenland Current with an Atlantic water signature (Münchow et al., 2015; TARA stations 206, 208, 209).

DNA extraction and high-throughput amplicon sequencing protocols are described in de Vargas et al. (2015) and Alberti et al. (2017). For this work we used reads derived from the V9 region of the 18S rRNA gene, and compared our main results to data derived from the V4 region of the same gene (Table S2). For the latter, the utilized primers were TA-Reuk454FWD1 (5′-CCAGCA(G/C)C(C/T)GCGG-TAATTCC-3′, S. cerevisiae position 565–584) and TAReukREV3 (5′-ACTTTCGTTCTTGAT(C/T)(A/G)A-3′, S. cerevisiae position 964–981). In both cases, reads were grouped into OTUs using the following swarm-based pipeline: paired-end reads were merged with vsearch’s --fastq_mergepairs command (version 2.15.1, allowing for staggered reads; Rognes et al., 2016), and trimmed with cutadapt (version 3.0; Martin, 2011), keeping only reads containing both forward and reverse primers. After trimming, the expected error per read was estimated with vsearch’s command --fastq_filter and the option --eeout. Each sample was then de-replicated; that is, strictly identical reads were merged, using vsearch’s command --derep_fulllength, and converted into fasta format. Clustering was performed at the sample level with swarm 3.0 using default parameters (Mahé et al., 2015). Prior to global clustering, individual fasta files (one per sample) were pooled and further dereplicated with vsearch. Files containing per-read expected error values were also dereplicated to retain only the lowest expected error for each unique sequence. Global clustering was performed with swarm (using the fastidious option). Cluster representative sequences were then searched for chimeras with vsearch’s command --uchime_denovo using default parameters (Edgar et al., 2011).

Clustering results, expected error values, taxonomic assignments, and chimera detection results were used to build a “raw” occurrence table. Reads without primers, reads shorter than 32 nucleotides, and reads with uncalled bases (“N”) were discarded. For a “filtered” occurrence table, non-chimeric sequences, sequences with an expected error per nucleotide below 0.0002, and clusters containing at least 2 reads were retained. As primer trimming is not perfect, some sequences can still contain primer fragments or be excessively trimmed. These sub- or super-sequences were identified using vsearch and merged with their closest, most abundant perfectly trimmed sequence. Finally, occurrence patterns throughout our sample collection were used to further refine the occurrence table. Clusters that contain sub-clusters with only a single-nucleotide difference but with different ecological patterns (defined here as uncorrelated abundance values in at least 5% of the samples) were turned into distinct clusters (https://github.com/frederic-mahe/fred-metabarcoding-pipeline). On the other hand, clusters with similar sequences that had correlated abundance values in at least 95% of the samples were merged using a re-implementation of lulu’s method (Frøslev et al., 2017; https://github.com/frederic-mahe/mumu).

To ensure that the V9 and V4 data sets were comparable, we took into account differences in the number of samples and depth of sequencing. The V9 dataset had sequencing data for 18 surface samples and 8 SCM samples from the AO, while for the V4 dataset only 13 surface samples were available. The sequencing depth was approximately 10⁶ reads per sample for the V9 region and approximately 10⁵ reads per sample for the V4 region (Figure S3). These two aspects, together with the inherent differences related to the nucleotide sequence natural variation (e.g., GC content) and length that could affect steps such as the PCR or the subsequent construction of OTUs, motivated us to: (a) center our analyses on the V9 dataset and use the V4 for the validation of emerging patterns, and (b) exclude the very least abundant OTUs by keeping only those with at least one sample exhibiting a relative abundance of 0.001% in surface and SCM samples. We nevertheless recognize that rare OTUs are an important component of microbial communities, and the consequences of their removal were evaluated (Figure S4). From an Arctic perspective, the latter affected mainly the number of OTUs that were detected exclusively in the AO (Figure S4).

Both global datasets were rarefied at the Arctic minimum of 2.6 × 10⁵, which was greater than the global minimum of 1.8 × 10⁵, as a trade-off between reducing the effects of uneven sequencing depths for our ecological analyses and avoiding information loss from the Arctic plankton community. The slightly lower sequencing effort of certain samples from the non-Arctic ocean should be compensated by the higher number of non-Arctic stations (>3x). A sensitivity test was nevertheless performed by also rarefying at the global minimum (1.8 × 10⁵), with no significant changes in the proportion of biogeographic categories.

A taxon name was assigned to each OTU by comparing its representative sequence (the most abundant one) against PR2_V9, a customized PR2 reference database (Guillou et al., 2013; de Vargas et al., 2015; http://doi.org/10.5281/zenodo.3768951), using vsearch --usearch_global command. OTUs classified as Bacteria, Archaea and Metazoa were removed from the analyses, as well as a few OTUs with low representation, pluricellular life-forms (Ulvophyceae, Phaeophyceae, Streptophyta, Fornicata, Glomeromycota, Parabasalia, Phytomyxea, and Rhodophyta, except for genus Rhodella) or nucleomorphs of both chlorarachneans and cryptophytes. OTUs with an identity matching <80% of a reference sequence were excluded as well. Trophic functions (photosynthetic, phagotrophic, parasitic, photosynthetic endosymbiont) were assigned to OTUs for which information is known from literature. The photosynthetic function refers here to the presence of a permanent chloroplast, based on a literature review (de Vargas et al., 2015), which does not exclude the possibility of mixotrophy. Altogether, the V9 region contained more deep-branching lineages than the V4 region. Mismatches to the primer set or the amplified fragment being longer than the sequencing capabilities (2 × 250 bp) explain why Diplonemida and Euglenida from the Excavata supergroup would be absent from the V4 data (Vaulot et al., 2022). Comparisons at the level of deep-branching lineages were performed on the 47 coinciding groups (Figure S5). These groups have consistent relative abundance regardless of the marker used (Spearman’s rho = 0.93, p << 0.001). Raw reads of 18S rRNA gene metabarcoding are deposited at the European Nucleotide Archive (ENA) under projects PRJEB6610 and PRJEB9737.

OTU abundances were normalized by dividing each value by the total sum per sample, followed by a transformation using the “decostand” function. Bray-Curtis pairwise dissimilarities of samples were then calculated and used as input for a principal coordinates analysis. Differences between polar and nonpolar communities were assessed through a Permutational Multivariate Analysis of Variance (PERMANOVA) using the function “adonis.” A selection of contextual variables (temperature, chlorophyll a, iron, mixed layer depth, daylength and station depth) was fit on the ordination with the “envfit” function to explore potential driving factors of community differences. A goodness of fit (R²) and a p-value were retrieved from the function as a result from a permutative procedure (Table S3). This selection of contextual variables was chosen as general habitat descriptor and with the criterion of reducing collinearity (see multiple pairwise correlation analysis between all available habitat descriptors performed in Ibarbalz et al., 2019). Functions used for the above analyses are included in the R-package “vegan” (Oksanen et al., 2019) and were used with default parameters.

OTUs were classified according to broad biogeographic distributions into 5 mutually exclusive categories. We used a species indicator approach to identify OTUs with strong affinity for the AO (De Cáceres and Legendre, 2009). This effort was performed with the function “multipatt” in the R-package “indicspecies” v.1.7.8 with parameter “IndVal.g.” Specifically, an OTU was considered an indicator of the Arctic Ocean (s.Region_AO = 1) and not of any other non-polar oceanic region (s.Region_i = 0), when the association value was >0.7, and the p-value was <0.05, including false discovery rate correction for multiple comparisons. This association value threshold selected OTUs with at least 50% AO occupancy. Because many of these OTUs were also detected in the few samples of the Southern Ocean (n = 6), in the following text we refer to this group as Arctic/Polar indicators. Using the same approach, we defined non-Polar indicators by retrieving those OTUs with high affinity for regions below the Arctic circle (specifically, North Atlantic, South Atlantic, North Pacific or South Pacific Ocean; s.Region_AO = 0, s.Region_i = 1, same parameters). We emphasize that in these two categories the term “indicator” refers to a statistical definition based on their detection patterns throughout the world’s oceans, in contrast with other ecological definitions related to particular functions or sensitivities (Siddig et al., 2016). We did not determine specific Southern Ocean indicators because of the small number of samples (n = 6). A third category was set for globally ubiquitous OTUs that were defined as OTUs having >50% occupancy at Arctic and >50% occupancy at non-Arctic stations, and were not defined previously as indicators. The remaining OTUs that were detected in the AO were found at relatively lower abundances. If they were only detected in the AO, they were considered Arctic-exclusive; otherwise, they were called background. Background OTUs were not differentiated by a degree of cosmopolitanism (i.e., a single detection outside of the AO was sufficient) given their overall low detection in samples.

OTU biogeographic annotations are available in the Appendix (http://dx.doi.org/10.17632/zb96mzvhv5.1). Estimates of the relative contributions of each biogeographic category to the pool of OTUs derived from the V9 variable region of the 18S rRNA gene were compared with estimates for the V4 region (Table S4). The circular panel delineating OTU distributions across taxa, biogeographic categories and Arctic regions (Figure 4) was done with R-package “circlize” (v. 0.4.10; Gu et al., 2014), where each biogeographic category of each lineage was sorted according to a model-based clustering obtained through functions “mclustBIC” and “Mclust” from R-package “mclust” (v. 5.4.6; Scrucca et al., 2016). A full version of Figure 4A is available at https://doi.org/10.6084/m9.figshare.14562606. We acknowledge the challenges and uncertainties of placing drifting (planktonic) organisms into biogeographic categories and use the term “biogeography” in a broad sense.

Gross change in OTU identity and gross change in dominance structure or abundance were calculated for samples collected in the North Atlantic and Arctic oceans following the approach of Hillebrand et al. (2018). Briefly, in a pairwise comparison of samples, gross change in identity accounts for the sum of species that are not shared between samples, relative to the total number of species. This index corresponds to the Jaccard index. Gross change in abundance corresponds to a Simpson-based turnover index. Both indices vary between 0 and 1, with 1 meaning a complete change of community composition. Since gross change in identity can be sensitive to the sampling effort, we rarefied the data at the minimum sum among the subset of samples used prior to calculations (5 × 10⁵ reads per sample; function “rrarefy” from R-package “vegan”). The effect of rare OTUs was evaluated by repeating the analysis starting from a lower OTU percentage cutoff (1 × 10⁻⁴; Figure S4), allowing the number of rare OTUs to increase. This results in an increase in the average of gross change in identity by only 7%, with the overall patterns remaining the same.

Assemblages of eukaryotic plankton were sampled across heterogeneous conditions. The majority of the samples were collected in open waters at varying times since ice melt and varying distances from the seasonal ice edge (Figure 1B). Nitrate and chlorophyll a concentrations also varied greatly, reflecting different stages in the annual cycle of phytoplankton (Figures 1B and S2). Hydrographic profiles of salinity and temperature show the key water masses that were sampled during the Tara Oceans Polar Circle expedition (Figure S2). These included (i) the warm and saline Atlantic water that enters the AO through the Fram Strait and the Barents Sea and spreads across the Eurasian and Canada basins at intermediate depths, (ii) the relatively cold and fresh lens of surface water over the Siberian, Chukchi and Canadian shelves, (iii) the modified Arctic and Pacific waters (cold and fresh) that are advected southward along Baffin Island, and (iv) the warmer and saltier Atlantic water advected along western Greenland. Thus, while each station represented only a snapshot in time of a dynamic environment, their combination and the deep sequencing per site (>10⁶ sequence reads per sample; Figures S3 and S4) captured a significant portion of the diversity of unicellular eukaryotes that occupy Arctic shelves.

An average of 934 (± 264) OTUs per station were recovered in the AO, resulting in a total of 3,082 different OTUs across the AO (for the V9 region; Figure S3; Table S2). These results were in good agreement with the numbers of OTUs derived from the V4 dataset (Table S2; Ibarbalz et al., 2019). The majority of the reads belonged to photosynthetic and phagotrophic organisms (Figure S6A), and likely included mixotrophic groups that can be both autotrophic and heterotrophic. Parasites and putative endo-photosymbionts were detected at lower abundances and mostly in regions associated with northward-flowing Atlantic water (regions 1 and 5; Figure S6A). Relative contributions of major photosynthetic and phagotrophic taxa to the pan-Arctic pool are reported in Figure S6B–D to provide a general overview of the communities that were sampled during the Tara Oceans Polar Circle expedition.

In this study, we focused on the surface layer. Nevertheless, a relatively small number of available samples from the subsurface chlorophyll maximum (SCM, n = 8, V9 region) allowed us to explore how AO photic communities vary with depth, as previous studies have shown that SCM communities often differ from those in surface water (Monier et al., 2013; Ardyna et al., 2017). Our analysis suggests that variability due to geographic location of each station (as well as temporal variability due to the progression of the season during the circumnavigation) surpassed variations due to vertical zonation (Figure S7; PERMANOVA for surface versus SCM: R² = 0.06, not significant). Still, a fraction of low-abundant OTUs (32 ± 9% of OTUs per station, representing 1.9 ± 1.7% of read abundance) were only detected in the SCM, indicating that this layer contributes to AO biodiversity.

When comparing community compositions of all Tara Oceans samples across the world’s oceans, a clear pattern emerges, with AO communities clustering separately from most other communities (Figure 1C–D; see also Sommeria-Klein et al., 2021). Southern Ocean communities and some Atlantic Ocean communities, however, clustered close to AO samples along the first axis (31.9% of total variation). This pattern is consistent with previous reports for other groups of plankton such as bacteria (Ghiglione et al., 2012; Salazar et al., 2019), archaea (Galand et al., 2009a) and viruses (Gregory et al., 2019; Endo et al., 2020), as well as overall distinct plankton co-occurrence networks (Chaffron et al., 2021). Overall, community dissimilarities co-varied strongly with temperature and station depth, and to a lesser extent with other selected environmental parameters (Figure 1D; Table S3). The high proportion of samples collected over continental shelves in the AO introduces a bias in the analysis, as some of the differences from the rest of the predominantly oceanic Tara Oceans samples could be attributed to differences between shelf and oceanic communities. However, the clustering of all AO samples with each other, regardless of station depth (Figure 1D), and with oceanic samples from the Southern Ocean (station depth near 400 m for the two shallowest stations and between 2300–4100 m for the remaining stations; Table S3; PERMANOVA for “Polar” as factor: R² = 0.3, p < 0.001), implies that the more coastal nature of Arctic shelves is not the only driver of the observed dissimilarities.

Throughout the evolution of the modern Arctic, tectonic processes, such as the opening of the major Arctic gateways and glacial controls on sea level, combined with climatic changes have provided opportunities for species from the global ocean to establish populations in the AO. Contrarily, physical barriers for dispersal and isolation in sub-Arctic and Arctic refugia created during glaciation periods allowed genetic divergence in some populations (Hardy et al., 2011). These processes likely gave rise to a distinct mix of organisms, with different biogeographic histories, that is reflected in the contemporary pool of Arctic unicellular eukaryotic plankton. To further understand how the pan-Arctic community differentiates from and is influenced by communities from other oceans, we analyzed the biogeographic structure of the pan-Arctic community (Figure 2A).

On a pan-Arctic scale, OTU abundance and occupancy correlated positively (Figure 2B), a pattern repeatedly observed across taxa in other environments (Gaston et al., 2000). The most abundant and widespread OTUs belonged to the Arctic/Polar indicator and the globally ubiquitous groups. As such, members of these groups likely play pivotal roles in the AO ecosystem, although accounting for only 9–12% and 6% of the pool of OTUs found in the AO, respectively (Figure 2B, Table S4). Winter samplings, in spite of a general decrease in biomass and changes in relative abundances, could determine if these OTUs are (as we think they are) still detected at most of the stations. Approximately half of the Arctic/Polar indicators (55–57%) were detected in our Southern Ocean samples (examples in Figure 2C). This result may indicate that a connectivity between populations of an OTU present in the two polar oceans existed, and perhaps still exists, or a difficulty to distinguish different ecotypes from their V9 or V4 rRNA gene sequence (Šupraha et al., 2022). This connectivity or difficulty may also be the case for the relatively few globally ubiquitous OTUs. Future strain isolations and single-cell sequencing could clarify these indications (Dorrell et al., 2022).

In comparison, Arctic-exclusive OTUs were less abundant, showed a more restricted distribution, and accounted for 10–15% of the OTU richness (significantly affected by the chosen abundance cutoff; see Methods; Figure 2B, Table S4). We cannot preclude the possibility that Arctic-exclusive OTUs might increase in abundance at certain times of the year not covered in this study, display interannual variations in population dynamics, or might simply exist at low abundances (Galand et al., 2009b; Ser-Giacomi et al., 2018). The small proportion of Arctic-exclusive OTUs that were found in this study is consistent with the notion, brought from microscopy observations, that true endemic planktonic species are few (Okolodkov and Dodge, 1996; Bluhm et al., 2011). However, their relatively low abundance, the limitations in the microbial species definition based on rRNA marker genes, and the difficulties in defining geographic boundaries for marine plankton preclude us from reaching firm conclusions about endemicity.

Non-Polar indicators were also present in low abundance and had restricted distribution within the AO (discussed below), yet their overall contribution to the observed richness was not negligible (16–24%; Figure 2B). These OTUs may represent a “newer” phase of the recolonization of the AO (on paleoclimatic time scales), with populations that are not fully established (and may never will). Finally, background OTUs contributed the most to OTU richness (48–54%; Figure 2B; Table S4). In light of their low detection in the global dataset, further investigation will be required to understand their biogeography (Figure 2C).

An intriguing question that follows the previous discussion is how taxonomies and biogeographic categories relate to each other. Of the deep-branching lineages of unicellular eukaryotes represented in the global ocean V9 metabarcoding data set, 88% (75 out of 85) were detected in at least one AO station. A complementary analysis of the V4 data set produced a similar result (Table S2; Figure S5). Ten lineages from the global data set (V9) were not detected in our pan-Arctic sampling and, according to our analysis, their occurrence was found to be confined largely to low latitudes (Appendix). Overall, large variations were found in the representation of the different biogeographic categories among different taxa (Figure 3A; Appendix). Centrohelida, Picomonadida, Mamiellophyceae, and the clades MOCH-1,2 were the lineages with the highest fraction of Arctic/Polar indicator OTUs (accounting for more than 25% of those OTUs). Conversely, non-Polar indicator OTUs made up high proportions (≥50%) of the OTUs associated with Prasinophyceae Clade 7 (Chloropicophyceae), RAD-A, MALV-IV and Collodaria. Highest proportions of Arctic-exclusive OTUs (>40%) were detected among Euglenida (mostly photosynthetic or mixotrophic), Dactylopodida (heterotrophic) and the parasitic lineage Ascetosporea, and fungi in the Chytridiomycota and Oomycota. Globally ubiquitous OTUs were detected in more than a third of the Arctic lineages, rarely exceeding 20% of the OTUs per lineage, while background OTUs were present in all lineages except one (Prasinophyceae Clade 7), with varying contributions per lineage (on average 58 ± 22%). Among the five biogeographic categories, variation in taxonomic and trophic composition was marked (Figure S8). Most notable were the higher proportions of OTUs assigned to MALV-I & II (parasites) and Diplonemida detected as globally ubiquitous, non-Polar indicator and background compared to the two categories with higher affinity to the AO (Arctic/Polar indicator and Arctic-exclusive). For example, the relative contribution of MALV II to the globally ubiquitous OTUs is 16%, to the non-Polar indicator OTUs, 22%, and to the Arctic/Polar indicator OTUs, 2% (Figure S8). The relatively high contribution of OTUs associated with parasitic lifestyles within the non-Polar indicators and background groups points to the potential role of bioadvection in altering trophic dynamics in the AO, if those populations and their hosts become established, and calls for more in-depth investigation of parasitic interactions in the AO.

While the diversity of many of the taxa mentioned here has been described previously in Arctic waters (Poulin et al., 2011; Terrado et al., 2013; Thaler and Lovejoy, 2015; Hassett et al., 2017; Šupraha et al., 2022), it was examined for a smaller subset of taxa and rarely placed in a global context. Generally, global trends show a decline in diversity toward the poles in many groups of organisms, including plankton (Ibarbalz et al., 2019). Impoverished richness (low P_AO/P_G; Figure 3B and C) is seen across all taxa, but is particularly evident among lineages within Alveolata and Rhizaria. Yet, if we consider Arctic/Polar indicator- and Arctic-exclusive OTUs as the specific contribution of the AO to global diversity, our analysis shows that the AO is an important source of biodiversity (Figure 3B and C). This importance is exemplified in lineages such as Euglenida, Mamiellophyceae, Crymonadida, Ebriida and Chytridiomycota, for which Arctic/Polar indicator and Arctic-exclusive OTUs constitute more than a quarter of their global biodiversity (Figure 3C). Given that the AO represents <5% of the global ocean’s surface, this contribution to global plankton biodiversity is notable. While environmental sequencing data should be interpreted with caution (i.e., limited resolution of the DNA marker genes that may underestimate the number of Arctic species, and data that are not strictly quantitative), the high degree of coherence between the V9 and V4 datasets increases confidence in the ecological patterns presented here. The frequently observed low similarity to sequences from known species (Figure 3D) points to important gaps in the reference database and is a reminder that further work on microbial diversity (e.g., cultures or single cell isolations) is highly needed (Šupraha et al., 2022).

Concerns of climate change effects on Arctic biota often focus on loss of biodiversity. The large contribution of OTUs with broad biogeographic signatures (i.e., non-Polar indicators and background) to the pool of OTUs in the AO indicates a role for bioadvection in maintaining biodiversity. Annual exchange of Arctic surface water is on the order of 10% per year (Wefing et al., 2021), implying a continued influx of species from lower latitudes into the AO. Depending on the changing balance between extinction and colonization rates, diversity of unicellular eukaryotes in the AO could maintain or even increase in the near future. Successful colonization, however, depends on the ability of species to adapt to the unique conditions in the AO. While warming could favor low latitude species, they may not be able to flourish under extreme conditions of seasonality in sunlight. We hypothesize that non-Polar indicator OTUs represent species with high dispersal potential, but physiologies that are not well suited for the current AO conditions. By mapping the distributions of all OTUs across the different regions of the AO, we observed that non-Arctic OTUs were more prevalent at the Arctic gateways compared to the inner Siberian shelves (Figure 4), a pattern indicative of the competing roles of advection (dispersal), environmental filtering, and/or biological interactions (Wassmann et al., 2015). This prevalence at the gateways was particularly noticeable among the lineages MALV-II and MALV-I, haptophytes, MAST-3,12, Diplonemida and Dinophyceae. The larger number of non-Polar indicators detected over the Atlantic inflow shelf compared to the Pacific inflow shelf is consistent with the significantly larger volume exchange at the Atlantic sector (Timmermans and Marshall, 2020). In contrast, Arctic-exclusive OTUs were generally less prevalent on shelves influenced by Atlantic water, although distribution patterns varied greatly among lineages (Figure 4). This reduced prevalence may point to a combination of habitat specificity and rarity within this category of OTUs, which could reflect highly specialized species with a narrow environmental niche or stochastic fluctuations in populations (Kunin and Gaston, 2012).

Indirect insights on environmental filtering can be gained from an examination of changes in community structure between the neighboring North Atlantic Ocean (NAO) and AO communities. In the absence of strong environmental gradients across the NAO-AO continuum, we may expect changes in species abundances (e.g., reflecting seasonal changes), but not necessarily a large shift in species identity. However, changes in both species abundances and identity would suggest a strong environmental filtering. We found large, concurrent changes in both species identity and relative abundances within and between the two oceans (Figure 5; see Methods). Although there are large spatial variations in community composition within each ocean basin, changes between AO and NAO communities were significantly greater than those within each basin (Figure 5; Table S5). These changes co-varied with changes in absolute richness (Figure S9), suggesting that the reduced biodiversity of Arctic plankton is contributing to the observed change in species identity. Overall, these results support the notion of environmental filtering at the North Atlantic gateway, a process that is also likely to occur at the Pacific inflow gateway and over outflow shelves (Kalenitchenko et al., 2019).

The use of consistent sampling protocols allowed us to examine assemblages of unicellular eukaryotic plankton in the AO from a global perspective. We have shown that, based on 18S rRNA gene amplicons of both V9 and V4 regions, AO communities consist of a mixture of unique and broadly distributed OTUs with a characteristic community structure. Among the shared OTUs, a considerable fraction (10–16%) have high affinity for lower, warmer latitudes, which could be a sign of change in progress. The importance of advected populations into the AO (Atlantic and Pacific gateways) has been examined primarily for mesozooplankton (Wassmann et al., 2015; Wassmann et al., 2019). Our results further demonstrate the strong connectivity of unicellular eukaryotic plankton communities with the global ocean. Surface velocities are increasing for the poleward-flowing North Atlantic current (Oziel et al., 2020), and the annual mean transport of Pacific water into the AO is also increasing (Woodgate, 2018). Together with other environmental changes over the AO shelves, they will likely alter the current dynamic of community assembly processes, with consequences for the structure and function of the AO ecosystem.

Results of this work are compiled in a resource table in which we provide the taxonomic and biogeographic affiliations for each OTU (Appendix). This resource could become a valuable framework for developing monitoring strategies of unicellular eukaryotic taxa. As an example, regarding the observed higher proportion of non-Polar indicators at the Atlantic inflow shelves, changes to their overall contribution to diversity, or more specifically to one of their deep-branching lineages (e.g., haptophytes; Figure S6; Ardyna and Arrigo, 2020; Egge et al., 2021) and/or corresponding OTUs (e.g., several assigned to haptophyte genera Phaeocystis, Chrysochromulina and the uncultured Clade D; abundance patterns not shown) could be monitored. Of socio-economic relevance, species implicated in harmful algal events are included in the resource table as well (e.g., Dinophysis acuminata and Alexandrium and Pseudo-nitzschia spp. among Arctic/Polar indicators; Dinophysis acuta, Protoceratium reticulatum, and Margalefidinium polykrikoides—previously Cochlodinium polykrikoides—, the haptophyte Chrysochromulina leadbeateri and the pelagophyte Aureococcus anophagefferens among other biogeographic categories). Continued monitoring of target OTUs from this study, as well as an improved understanding of their physiologies and traits, could prove useful for predicting and detecting future changes in species ranges.

This work represents an extensive effort to map occurrence patterns among the 75 different deep-branching lineages that were retrieved from the AO. The majority of the samples in this study were collected from Arctic shelves, which account for 50% of the total area of the AO. Arctic shelf regions are affected most immediately by the decrease in the extent of summer sea ice (Onarheim et al., 2018) and the influx of allochthonous biomass from the neighboring oceans (i.e., the Atlantic and Pacific inflow shelves). In order to further understand the implications on the plankton community and the AO ecosystem, our biogeographic perspective should integrate the winter season and the central sector of the Arctic as well. Under-ice communities of microbial eukaryotes from the central Arctic are diverse and dynamic, with apparent exchanges between sea ice and water column communities (Hardge et al., 2017). The ice-associated taxa and mixotrophic plankton would likely contribute to our Arctic/Polar indicator or Arctic-exclusive categories. As the ice retreats beyond the shelf break, exchanges between shelf and central basin communities are expected to increase due to mixing processes at the shelf break (Charette et al., 2020), which will likely contribute to continued shifting of the AO plankton community.

Raw reads are deposited at the European Nucleotide Archive (ENA) under projects PRJEB6610 and PRJEB9737 where they are publicly available (https://www.ebi.ac.uk/ena/browser/home).

The supplemental files for this article can be found as follows:

Figures S1–S9. Tables S1–S5. PDF

Appendix, pages 1–6. Available at http://dx.doi.org/10.17632/zb96mzvhv5.1.

Tara Oceans (which includes both the Tara Oceans and Tara Oceans Polar Circle expeditions) would not exist without the leadership of the Tara Ocean Foundation and the continuous support of 23 institutes (http://oceans.taraexpeditions.org). We further thank the commitment of the following sponsors: CNRS (in particular Groupement de Recherche GDR3280 and the Research Federation for the study of Global Ocean Systems Ecology and Evolution, FR2022/Tara Oceans-GOSEE), European Molecular Biology Laboratory (EMBL), Genoscope/CEA, the French Ministry of Research, and the French Government “Investissements d’Avenir” programmes OCEANOMICS (ANR-11-BTBR-0008), FRANCE GENOMIQUE (ANR-10-INBS-09-08), MEMO LIFE (ANR-10-LABX-54), and PSL* Research University (ANR-11-IDEX-0001-02). Funding for the collection and processing of the TARA data set was provided by NASA Ocean Biology and Biogeochemistry program under grants NNX11AQ14G, NNX09AU43G, NNX13AE58G and NNX15AC08G to the University of Maine, and Canada Excellence Research Chair on Remote sensing of Canada’s new Arctic frontier (MB). CB additionally acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 835067; Diatomic). FMI is currently a researcher at CONICET, thanks his current PI, Pedro Flombaum, and acknowledges funding from ECOS SUD (AT08ST18), ANPCYT (PICT-2017-3020), UBACYT (20020170100620BA) and MINCYT Pampa Azul (PIDT-A6). We also thank the support and commitment of agnès b. and Etienne Bourgois, the Prince Albert II de Monaco Foundation, the Veolia Foundation, Region Bretagne, Lorient Agglomeration, Serge Ferrari, Worldcourier, and KAUST. The global sampling effort was enabled by countless scientists and crew who sampled aboard the Tara from 2009 to 2013, and we thank MERCATOR CORIOLIS and ACRI-ST for providing daily satellite data during the expeditions. We are also grateful to the countries who graciously granted sampling permissions. Finally, we thank the two anonymous reviewers for their constructive feedback that greatly helped to improve the manuscript. The authors declare that all data reported herein are fully and freely available from the date of publication, with no restrictions, and that all of the analyses, publications, and ownership of data are free from legal entanglement or restriction by the various nations whose waters the Tara Oceans expeditions sampled in. This article is contribution number 144 of Tara Oceans.

Silvia G. Acinas, Marcel Babin, Peer Bork, Emmanuel Boss, Chris Bowler, Guy Cochrane, Colomban de Vargas, Mick Follows, Gabriel Gorsky, Nigel Grimsley, Lionel Guidi, Pascal Hingamp, Daniele Iudicone, Olivier Jaillon, Stefanie Kandels, Lee Karp-Boss, Eric Karsenti, Fabrice Not, Hiroyuki Ogata, Stéphane Pesant, Nicole Poulton, Jeroen Raes, Christian Sardet, Sabrina Speich, Lars Stemmann, Matthew B. Sullivan, Shinichi Sunagawa, Patrick Wincker.

The authors declare no competing interests.

Substantial contributions to conception and design: FMI, LK-B, CB.

Acquisition of data: MA, EB, LK-B, NH, FM, AZ, ES, CdV, PW.

Analysis and interpretation of data: FMI, NH, MA, SM, SC, LZ, CB, LK-B.

Drafting the article or revising it critically for important intellectual content: FMI, LK-B, NH, FM, MA, AZ, ES, CL, FL, OJ, DI, SM, MBS, SC, EK, MB, EB, PW, LZ, CdV, CB.

Final approval of the version to be published: FMI, LK-B, NH, FM, MA, AZ, ES, CL, FL, OJ, DI, SM, MBS, SC, EK, MB, EB, PW, LZ, CdV, CB.

How to cite this article: Ibarbalz, FM, Henry, N, Mahé, F, Ardyna, M, Zingone, A, Scalco, E, Lovejoy, C, Lombard, F, Jaillon, O, Iudicone, D, Malviya, S, Tara Oceans Coordinators, Sullivan, MB, Chaffron, S, Karsenti, E, Babin, M, Boss, E, Wincker, P, Zinger, L, de Vargas, C, Bowler, C, Karp-Boss, L. 2023. Pan-Arctic plankton community structure and its global connectivity. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2022.00060

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Kevin R. Arrigo, Department of Earth System Science, Stanford University, Stanford, CA, USA

Knowledge Domain: Ocean Science