Exploring the Microdiversity Within Marine Bacterial Taxa: Toward an Integrated Biogeography in the Southern Ocean.

Most of the microbial biogeographic patterns in the oceans have been depicted at the whole community level, leaving out finer taxonomic resolution (i.e., microdiversity) that is crucial to conduct intra-population phylogeographic study, as commonly done for macroorganisms. Here, we present a new approach to unravel the bacterial phylogeographic patterns combining community-wide survey by 16S rRNA gene metabarcoding and intra-species resolution through the oligotyping method, allowing robust estimations of genetic and phylogeographic indices, and migration parameters. As a proof-of-concept, we focused on the bacterial genus Spirochaeta across three distant biogeographic provinces of the Southern Ocean; maritime Antarctica, sub-Antarctic Islands, and Patagonia. Each targeted Spirochaeta operational taxonomic units were characterized by a substantial intrapopulation microdiversity, and significant genetic differentiation and phylogeographic structure among the three provinces. Gene flow estimations among Spirochaeta populations support the role of the Antarctic Polar Front as a biogeographic barrier to bacterial dispersal between Antarctic and sub-Antarctic provinces. Conversely, the Antarctic Circumpolar Current appears as the main driver of gene flow, connecting sub-Antarctic Islands with Patagonia and maritime Antarctica. Additionally, historical processes (drift and dispersal limitation) govern up to 86% of the spatial turnover among Spirochaeta populations. Overall, our approach bridges the gap between microbial and macrobial ecology by revealing strong congruency with macroorganisms distribution patterns at the populational level, shaped by the same oceanographic structures and ecological processes.

Introduction

Biogeography has traditionally investigated the geographic distribution of macroorganisms in the Eukaryota domain. Yet, the last two decades have witnessed a growing number of studies focused on the biogeography of microorganisms, and repeatedly reporting non-random community assemblages of various prokaryotic microorganisms (Martiny et al., 2006; Wilkins et al., 2013; Karimi et al., 2018). Unlike contemporary driving factors (i.e., environmental selection) that have been extensively studied (Gilbert et al., 2012; Stegen et al., 2012; van der Gast, 2015), the role of historical processes – past ecological and evolutionary events – onto the present-day distribution patterns of microorganisms remains poorly investigated. Initially, the consensus was that the rapid and widespread dispersal of microbes should erase any signal of past historical events (Martiny et al., 2006). Nevertheless, it is now clear that historical processes, such as the dispersal barriers and geographic distance, might substantially contribute to microbes’ biogeography instead of environmental filtering (Hanson et al., 2012, 2019). Similarly to larger organisms, biogeographic regionalization, isolation, and endemism have been reported for microbes, reflecting the predominant effect of geographic distance over environmental variations (Papke et al., 2003; Whitaker et al., 2003).

To date, most of the microbial biogeographic patterns have been depicted at the whole community level, taking advantage of the extensive survey capacity provided by the next-generation sequencing (NGS) technologies (Almasia et al., 2016; Karimi et al., 2018; Wu et al., 2019). Nevertheless, as observed in various empirical studies, a finer taxonomic scale generally allows better detection of geographic patterns (Hanson et al., 2012; Chase and Martiny, 2018; Bay et al., 2020). Moreover, the ecological processes driving the biogeographic patterns at the community-level would mainly result from the accumulation of micro-evolutive processes, i.e., mechanisms contributing to the genetic composition and diversity within populations, and how they vary in space and time (Hanson et al., 2012; Larkin and Martiny, 2017). Hence the comprehensive description of these micro-evolutive processes requires considering the intra-population diversity, as commonly applied in phylogeographic studies of macroorganisms. In other words, microbial assembly processes need to be investigated at a finer taxonomic resolution than usually done by microbial biogeographic surveys and consider the “microdiversity” within groups (Larkin and Martiny, 2017; Chase and Martiny, 2018).

The oceans have been considered among the most challenging environments to test hypotheses about microbial biogeography, mainly due to the speculated transport of organisms over large distances by marine currents and the absence of perceivable marine barriers impeding potential dispersal events (Fenchel and Finlay, 2004). However, oceanic fronts separating different water masses have been recently identified as major microbial dispersal barriers (Flaviani et al., 2018). The Southern Ocean (SO) is a vast region representing approximately 20% of the world ocean surface. It surrounds Antarctica, and its northern limit is the Subtropical Front (Gruber et al., 2019). Two main oceanographic structures shape the SO biogeography; the Antarctic Polar Front (APF) and the Antarctic Circumpolar Current (ACC). The APF is classically considered a harsh north–south obstacle for dispersing marine organisms due to the brutal change in water temperature and salinity (Griffiths, 2010; Halanych and Mahon, 2018). Phylogenetic reconstruction achieved on various vertebrate and invertebrate taxa clearly supports the role of the APF on their respective diversification processes (Díaz et al., 2011; Poulin et al., 2014; Hüne et al., 2015; González−Wevar et al., 2017; González-Wevar et al., 2019). Accordingly, and based on the described distribution of species, biogeographers have traditionally recognized Antarctica and sub-Antarctica as the two main biogeographic provinces of the SO, even if several provinces have been proposed within each of them (Koubbi et al., 2014). Contrarily, outside the APF, the ACC is generally described as the driver of genetic connection across the sub-Antarctic zone due to its eastward circulation (Nikula et al., 2010; Cumming et al., 2014; Moon et al., 2017; González−Wevar et al., 2018). Intraspecific genetic and phylogeographic studies of macroorganisms have demonstrated the ACC’s role in connecting geographically distant sub-Antarctic provinces (Gérard et al., 2008; Fraser et al., 2009; González−Wevar et al., 2018; Frugone et al., 2019). The marine biota distribution in this region has been synthesized in the Biogeographic Atlas of the SO, providing updated biogeographic information of a wide range of benthic and pelagic taxa from Metazoan, macroalgae, and phytoplankton (De Broyer et al., 2014). Despite being the most abundant and diverse domains on Earth, Bacteria, and Archaea are not included in the SO Atlas (Shade et al., 2018). Studies conducted at the whole community-level support (1) the role of ACC as a likely efficient mechanism of circumpolar microbial transport and dispersal (Murray and Grzymski, 2007; Wilkins et al., 2013) and (2) the role of APF as the main dispersal barrier separating Antarctic and sub-Antarctic microbial assemblages (Wilkins et al., 2013; Flaviani et al., 2018; Raes et al., 2018). However, even when geographic distributions of marine microbial communities have been characterized in the region, the underlying evolutionary processes remain unclear, and their comprehensive understanding may rely on a higher taxonomic perspective, exploring bacterial populations’ microdiversity.

Targeting the intraspecific microdiversity using NGS data requires specific computational methods to discriminate the stochastic noise caused by random sequencing errors from those associated with biologically significant diversity (Hanson et al., 2012; Edgar and Flyvbjerg, 2015; Chase and Martiny, 2018). For this purpose, an algorithm called “Minimum Entropy Decomposition” (MED) relying on the oligotyping method has been proposed by Eren et al. (2015). This algorithm allows to identify true sequence variants (i.e., oligotypes) within the “operational taxonomic units” (OTUs), classically defined at 97% identity of the bacterial 16S rRNA gene. This approach has already been successfully used to unravel fine-grained biogeographic patterns of bacterial microdiversity in Arctic sediments, such as variations in oligotype distribution according to spatial and environmental parameters (Buttigieg and Ramette, 2014). Moreover, focusing on the sulfate-reducing genus Desulfotomaculum in Arctic marine sediments, Hanson et al. (2019) showed clear biogeographic patterns – attributed to historical factors associated with past environments – were only evident at the microdiversity level achieved with the oligotyping method. However, the microevolutionary processes causing the microdiversity were not assayed, and the study did not encompass large-scale distribution among different biogeographic provinces, as it was restricted to the west coast of Spitsbergen, Svalbard in the Arctic Ocean.

In the present proof-of-concept study, we aim to elucidate the evolutionary processes driving microbial biogeography across different provinces of the SO by combining (1) community-wide surveying provided by the high-throughput sequencing of the 16S rRNA gene, (2) intra-species microdiversity resolution obtained through the oligotyping method implemented in the MED pipeline, and (3) phylogeographic analysis as traditionally developed for macroorganisms as models. Considering the SO as an outstanding idoneous frame, we investigated the geographic distribution of genetic diversity of marine bacterial taxa across three biogeographic provinces: maritime Antarctica (King George Island, South Shetland Islands, West Antarctic Peninsula), sub-Antarctic Islands of the Indian Ocean (Kerguelen archipelago), and southern South America (Patagonia). We selected sites separated by the APF, i.e., maritime Antarctica and Patagonia, and others connected through the ACC, i.e., Patagonia and Kerguelen archipelago.

As the contribution of geography to biological diversity patterns (i.e., dispersal limitation) is stronger on habitat-specialists (i.e., taxa found in habitat with high selective strength) (Lindström and Langenheder, 2012; Szekely et al., 2013), and emphasized within homogeneous habitats distributed across large spatial scales (Astorga et al., 2012; Hanson et al., 2012; Langenheder and Lindström, 2019), we focused our study on the bacterial community associated to a specific habitat: the gut of Abatus irregular sea urchins. The Abatus genus is distributed across the SO and gathers various sibling species homologous in ecology and habitat, such as Abatus cavernosus in southern South America, Abatus cordatus in Kerguelen Islands, and Abatus agassizii in maritime Antarctica (Poulin and Feral, 1995; David et al., 2005; Díaz et al., 2012; Guillaumot et al., 2018). Since these species lack specialized respiratory structure, they are restricted to the well-oxygenated coarse sediments found at shallow depth (typically 1–3 m depth) in sheltered bays, protected from the swell (Poulin and Feral, 1995). Within the Abatus hosts, we focused on a specific micro-environment – the gut tissue – previously described to act as a selective filter of the external sediment microbiota, as illustrated by the reduction of bacterial diversity at both taxonomic and functional levels (Schwob et al., 2020). Working on the gut community with supposedly more limited dispersal capacity, through a high sequencing depth, is expected to (1) provide robust coverage of the bacterial diversity, (2) minimize the relevance of environmental filtering between provinces, (3) emphasize the contribution of geographic and oceanographic factors, and therefore (4) enhance the detection of phylogeographic signals across the SO (Hanson et al., 2012). As a model taxon to explore bacterial phylogeography in the SO, we chose the Spirochaeta genus (phylum Spirochaetes). Spirochaeta bacteria are recognized as the most prevalent and abundant genus in the Abatus gut tissue (Schwob et al., 2020). Moreover, spirochaetes are classically reported in marine benthic sediments (Bowman and McCuaig, 2003; Beiruti et al., 2017) and, to a lesser extent, in the water column (Ocean Barcode Atlas[1]). Thus, due to its ease of detection and ubiquity across biogeographic provinces, Spirochaeta represents an illustrative model to validate our methodology and explore marine bacteria’s spatial genetic patterns, from genus to population level. We hypothesized that the strong biogeographic barrier between South America and maritime Antarctica classically observed in the literature for macroorganisms (i.e., vicariance process) also affects the fine-scale genetic structure and the phylogeographic patterns within Spirochaeta OTUs. Besides, the ACC-mediated connectivity among sub-Antarctic provinces should be reflected by a greater genetic homogeneity of Spirochaeta populations between South American sites and the Kerguelen Islands, rather than with maritime Antarctica. Alternatively, the potential high dispersal capacity of Spirochaeta taxa may result in the absence of genetic and phylogeographic structure across the SO.

Materials and Methods

Sampling Collection, DNA Extraction, and 16S rRNA Gene Library Preparation

Adult Abatus individuals were sampled from four localities across the SO, including two sites in Patagonia, southern South America (Possession Bay, PAT1 and Puerto Deseado, PAT2), one site in Kerguelen Islands (Port-aux-Français, KER), and one site in the West Antarctic Peninsula (King George Island, KGI) (Figure 1 and Table 1). Marine surface sediments (0–5 cm, referred here as “external sediment”) were also sampled in each Abatus population’s immediate vicinity as the ingested food source of the sea urchins. Due to logistic constraints, it was not possible to collect external sediment in the PAT2 site. All individuals were dissected under sterile conditions to collect the whole digestive tract minus the caecum (identified as “gut tissue”) following Schwob et al. (2020). Gut tissue samples were gently rinsed with nuclease-free sterile water to remove the content (i.e., in sediment) and were then individually homogenized using mortar and pestle under a laminar-flow cabinet. Genomic DNA was extracted from both external sediments and the totality of the homogenized gut tissue samples using the DNeasy PowerSoil® Kit (Qiagen, Hilden, Germany) following the manufacturer’s recommendations.

FIGURE 1

Sampling localities across the Southern Ocean, encompassing Possession Bay and Puerto Deseado in Atlantic Patagonia (PAT1 and PAT2, respectively), King George Island in Maritime Antarctica (KGI), and Port-aux-Français in Kerguelen Islands (KER). The Antarctic Polar Front (APF) and the Southern Antarctic Circumpolar Front (SACCF) are represented by the limit between the white and the blue areas, and by the dashed line, respectively.

TABLE 1

Experimental design and sequencing data.

Locality	Province	GPS coordinates	Date	Designation	Sample types	N	Nseq. (Relat. Abund.)
King George Island	Maritime Antarctica	62°12′55.3″S 58°56′43.8″W	01–2019	KGI	External sediments	8	255786 (10%)
					Gut tissue	31	563383 (22%)
Bahía Posesión	Patagonia	52°19′52.97″S 69°29′10.50″W	07–2019	PAT1	External sediments	6	271828 (11%)
					Gut tissue	15	447892 (18%)
Puerto Deseado	Patagonia	47°45′07.0″S 65°52′04.0″W	12–2016	PAT2	External sediments	NA	NA
					Gut tissue	10	470087 (18%)
Port-aux-Français	Kerguelen Island	49°21′13.32″S 70°13′8.759″E	11–2017	KER	External sediments	5	92564 (4%)
					Gut tissue	14	440498 (17%)

A metabarcoding approach was used to assess the bacterial community composition in the external sediment and Abatus gut tissue samples. Briefly, extracted genomic DNA was used as the template for PCR amplification using the primers 515F 5′-GTGYCAGCMGCCGCGGTA-3′ and 926R 5′-CCCCGYCAATTCMTTTRAGT-3′ (Parada et al., 2016). The PCR conditions and 16S rRNA gene library preparation were the same as described in Schwob et al. (2020).

Metabarcoding Data Processing

External sediment and gut tissue amplicons were sequenced using the paired-end sequencing technology (2 × 250 bp) on the Illumina MiSeq platform at the UWBC DNA Sequencing Facility (University of Wisconsin–Madison, United States). Reads of 16S rRNA were processed using the open-source software Mothur v1.44.0. Briefly, 3′ and 5′ reads were paired and trimmed according to their length and quality as described in Schwob et al. (2020). Chimeric sequences were removed using Uchime implemented in Mothur (Edgar et al., 2011). Reads were clustered into OTUs at 97% identity threshold, and a filter of relative abundance at >0.0001% was applied, as recommended by Bokulich et al. (2013). Following this, a taxonomic classification was performed with the classify.otu function and the SILVA database v138 implemented in Mothur. An OTU table of Spirochaeta was edited (i.e., all OTUs assigned to the genus Spirochaeta), and converted into a presence/absence matrix to detect the shared or exclusive OTUs across the four selected sites, discarding abundance variations resulting from short-term environmental conditions. Bray–Curtis and Unweighted Unifrac distances were calculated from the OTU presence/absence matrix and used to perform a non-metric multidimensional scaling (NMDS) with the metaMDS function of the ade4 package (Chessel et al., 2009). The NMDS was plotted through a scatter diagram using the s.class function implemented in ade4. The locality’s contribution to the Spirochaeta OTUs composition in gut tissue samples was tested with permutational multivariate analysis of variance (permanova), using adonis and pairwise.adonis functions implemented in vegan and pairwiseAdonis R packages, respectively (Oksanen et al., 2011; Martinez Arbizu, 2017).

A subset of the three most abundant Spirochaeta OTUs present in the four localities was retained for further analysis (Supplementary File 2). All the sequences assigned to the selected Spirochaeta OTUs were retrieved using the bin.seqs command in Mothur. Finally, the resulting fasta files were processed independently through the Minimum Entropy Decomposition pipeline following Eren et al. (2015).

Minimum Entropy Decomposition

Minimum Entropy Decomposition pipeline was used to identify nucleotidic polymorphism at fine-scale resolution (<3% of differences) within the 16S rRNA gene sequences from Spirochaeta OTUs. Briefly, MED employs the Shannon entropy algorithm to discriminate biologically meaningful variations of closely related sequences from the stochastic noise caused by random sequencing errors, focusing on informative-rich variable nucleotide positions (Eren et al., 2013, 2015). The resulting taxonomic units will be referred to as Spirochaeta oligotypes. Unsupervised oligotyping was carried out individually on Spirochaeta OTUs using the default dynamically computed threshold from which entropy is considered as zero (−m). Additionally, each identified oligotype had to have a default minimum relative abundance of 2% in the OTU sequences dataset (Eren et al., 2015). Accumulation curves of oligotypes’ richness were computed for each Spirochaeta OTU at a 95% confidence interval using the package iNext (Hsieh et al., 2016) in r v3.6.0 (R Core Team, 2012). Pie charts of the relative site contributions in the total abundance of the Spirochaeta OTUs oligotypes were performed with the pie function in the package graphics in r v3.6.0.

Genetic Diversity and Structure of Spirochaeta Populations

The number of oligotypes (k), the oligotype diversity (H), the number of discriminant sites (S) and the pairwise difference between sequences (Φ) were estimated individually for each Spirochaeta OTU using the packages pegas (Paradis, 2010) and ape v5.3.0 (Paradis et al., 2015) in R v3.6.0. For comparative purposes among sites with unequal sample sizes, a composite bootstrapping script was written to rarefy the sequence datasets to the minimum number of sequences per site and repeat the rarefaction with 1,000 re-samplings. Confidence intervals at 95% of genetic diversity indices were then calculated using these iteration values. The genetic differentiation (Fst and Φst) among Spirochaeta populations was analyzed using the software arlequin v3.5.2 (Schneider et al., 2000) with 1,000 permutations and a significance threshold at 0.05. Phylogeographic differentiation was also estimated with the nearest-neighbor statistic Snn (Hudson, 2000), and the significance of Snn estimates was tested with a permutation test through DnaSP v5.10.01 (Librado and Rozas, 2009). The reconstruction of oligotype networks was performed using the Median Joining method with the software Populational Analysis with Reticulate Trees v1.7.0 in PopART (Leigh and Bryant, 2015). Because of the differences in Spirochaeta OTUs’ abundances across the four sites, oligotypes frequencies were calculated for each site and then summed per OTU to reconstruct the oligotype networks.

Quantification of Selection, Dispersal, and Drift

The relative contribution of stochastic (i.e., dispersal, drift) and deterministic (i.e., selection) processes, on Spirochaeta oligotype assembly was measured for the selected OTU, following the analytical framework described in Stegen et al. (2015) and illustrated by Feng et al. (2018). Briefly, the approach relies on the comparison of the phylogenetic turnover between communities across samples (β mean nearest-taxon distance, βMNTD) to a null distribution of βMNTD, and denoted as the β-nearest taxon index (βNTI). βNTI values indicate that taxa between two communities are more (i.e., βNTI < −2) or less (i.e., βNTI > +2) phylogenetically related than expected by chance, thus suggesting that communities experience homogenizing or variable selection, respectively (Stegen et al., 2015). βNTI values ranging from −2 to +2 indicate a limited selection effect and point to dispersal limitation and ecological drift out as possible community composition drivers.

Then, the respective effect of dispersal limitation and ecological drift were disentangled using the pairwise Bray–Curtis-based Raup–Crick dissimilarity index (RCBray) among sites (Chase et al., 2011), weighted by oligotype abundance (Stegen et al., 2013). The RCBray values < −0.95 and > +0.95 correspond to communities that have –respectively– more or fewer taxa in common than expected by chance, and therefore indicate that community turnover is driven by homogenizing dispersal (RCBray < −0.95) or dispersal limitation plus drift (RCBray > +0.95). On the contrary, RCBray values > −0.95 and < +0.95 are indicative of ecological drift (Ramoneda et al., 2020).

The βMNTD/βNTI and RCBray matrices, and the respective contributions of the four ecological processes, were calculated using an optimized version of the initial script of Stegen et al. (2013), developed by Richter-Heitmann et al. (2020). The phylogenetic trees required for the βMNTD/βNTI matrix were generated using PhyML v3.0 (Guindon et al., 2010), and the oligotype sequences of Spirochaeta previously aligned with MUSCLE (Edgar, 2004). Comparisons with βNTI and RCBray null models included 999 randomizations.

Testing for Isolation by Distance and Environment

To disentangle the relative effect of geographic distance and abiotic environmental differences on the Spirochaeta oligotype composition between samples, we used the distance-based multiple matrix regression with randomization (MMRR) approach (Wang, 2013). We extracted a set of 9 environmental variables for each of our sampling site from the Bio-ORACLE database (Assis et al., 2018), including pH, the means of nitrate, silicate, and phosphate concentrations, and the means at the mean depth of seawater salinity, dissolved oxygen concentration, seawater temperature, seawater temperature range and chlorophyll concentration. All environmental variables were standardized ((, and were then analyzed through principal components analysis (PCA). As a high percentage of the variation among localities was explained by the first component (PC1, >91%, Supplementary Files 6, 7), we transformed the scores of PC1 into Euclidean distance using the vegdist function in the vegan package in R to use it as the environmental distance matrix further. The longitude and latitude coordinates were converted into kilometers using the earth.dist function implemented in the fossil package in R (Vavrek, 2011). The geographic distances were transformed using the Hellinger method through the decostand function of the vegan package in R. The dissimilarity matrix of Spirochaeta oligotype composition among samples was also created from Bray–Curtis distances using the vegdist function of the R package vegan. Finally, to evaluate the relative weight of environmental and geographic distance matrices, an MMRR was performed using the R package PopGenReport (Adamack and Gruber, 2014), and the correlation coefficients and their significance were estimated based on 10,000 random permutations.

Connectivity Among Spirochaeta Populations

The amount and direction of gene flow among Spirochaeta populations were estimated using the coalescent-based program Lamarc v2.1.10 (Kuhner, 2006). A total of 10 runs was performed for each Spirochaeta OTU, consisting of likelihood searches of 20 initial and 2 final chains, with a minimum of 500 and 10,000 recorded trees, respectively, and sampling every 20 generations after a burn-in of 1,000 genealogies. The effective number of migrants per generation (Nm) among Spirochaeta populations was calculated by multiplying the maximum likelihood estimates (MLE) of the mutation parameter (θ) by the migration parameter (M), both estimated through the Lamarc program. We present the mean and standard deviation of the estimated Nm values obtained from the 10 runs for each Spirochaeta OTU.

Results

Sequencing Performance and OTU-Based Analysis

A total of 4,184,226 raw reads was generated from the 91 samples of external sediments and gut tissues. Once processed, 2,542,038 cleaned sequences distributed into 727 OTUs were obtained from the external sediment and gut tissue samples (details provided in Table 1). Out of this condensed dataset, 425,613 sequences associated with the Spirochaeta genus were retrieved, representing a total of 10 OTUs.

Both Bray–Curtis and Unweighted Unifrac distance methods did not reveal any difference in Spirochaeta OTU composition between the Patagonian sites (PAT1 and PAT2) (Figures 2A,B and Supplementary File 1). Conversely, Kerguelen Islands (KER) and martime Antarctic (KGI) sites were significantly different from each other in terms of Spirochaeta OTU composition and with the Patagonian ones (Supplementary File 1).

FIGURE 2

Non-metric multidimensional scaling scatter diagram of the Spirochaeta OTUs composition in gut tissue samples across the localities. Presence/absence matrix converted in Bray–Curtis distances (A) and Unweighted Unifrac distance (B). Colors are assigned to the locality.

The relative abundance analyses among the 10 Spirochaeta OTUs (Figure 3) showed four of them were shared among all the SO’s sampled provinces. Three OTUs (OTU6, OTU7, and OTU532) were more abundant in maritime Antarctica (KGI), four (OTU40, OTU42, OTU278 and OTU561) were more abundant in the Patagonian locality PAT1, two (OTU23 and OTU221) were predominantly found in the Patagonian locality PAT2, and a single one (OTU349) was predominant in Kerguelen Islands (Figure 3).

FIGURE 3

Clustering of Spirochaeta OTUs based on their relative abundances in each site. Clustering of Bray–Curtis distances matrix using the complete linkage method. The size of the circles indicates the repartition of a given OTU sequences among the four localities. The percentages indicate the OTUs’ proportions in the Spirochaeta dataset. Framed labels and black arrows indicate the selected OTUs that were selected to be process through the Minimum Decomposition Entropy pipeline (MED).

To test the genetic and phylogeographic structure of Spirochaeta at the broadest geographic scale available through our dataset, we selected the three most abundant Spirochaeta OTUs which were detected within the four localities of the dataset (i.e., OTU6, OTU7, and OTU40). These three selected co-distributed OTUs are good targets to constitute a metapopulation, which is a meaningful ecological unit of distinct local populations separated by gaps in habitats and interconnected to some extent via dispersal events of individuals (Halsey et al., 2015). The relative abundance in sample types, the closest sequence retrieved from Blast analysis, and the distribution of these OTUs among the localities are provided in the Supplementary Files 2, 3.

Microdiversity Within Spirochaeta OTUs

A total of 48, 96, and 48 oligotypes were defined for OTU6, OTU7, and OTU40, respectively (Supplementary File 4). Accumulation curves of OTU6, OTU7, and OTU40 oligotypes reached saturation in almost all localities indicating that the overall majority of Spirochaeta microdiversity has been found within the analyzed samples (Supplementary File 5). Diversity indices measured as N, k, S, h, and Π for each OTU in each locality are provided in Table 2. The genetic diversity (H) ranged from 0.0681 (OTU40 in KER) to 0.8036 (OTU7 in PAT1) across localities (Table 2). Patagonian sites exhibited higher oligotype and nucleotide diversity for OTU6 and OTU7 than maritime Antarctica and Kerguelen Islands localities. In contrast, the genetic diversity of the OTU40 oligotypes was higher in the maritime Antarctic site and lower for the Kerguelen population.

TABLE 2

Summary of oligotypes number and genetic indices per OTU and per site for the three most abundant Spirochaeta OTUs found in all sampling localities.

OTU	Site	N	k	S	H	Π
OTU6	KGI	125,112	26 ± 0	24 ± 0	0.4879 ± 0.0004	0.7919 ± 0.0012
	PAT1	6,452	18 ± 0	8 ± 0	0.6701 ± 0.0000	1.5772 ± 0.0000
	PAT2	22,545	22 ± 0	12 ± 0	0.6574 ± 0.0002	1.5426 ± 0.0006
	KER	29,448	33 ± 1	29 ± 0	0.5677 ± 0.0003	0.9687 ± 0.0011
OTU7	KGI	53,645	31 ± 0	32 ± 0	0.5555 ± 0.0003	1.4443 ± 0.0021
	PAT1	7,249	69 ± 0	37 ± 0	0.8036 ± 0.0000	2.0247 ± 0.0003
	PAT2	44,509	60 ± 0	33 ± 0	0.7958 ± 0.0002	1.7336 ± 0.0011
	KER	7,021	43 ± 0	37 ± 0	0.6306 ± 0.0000	1.0261 ± 0.0000
OTU40	KGI	47	4 ± 0	4 ± 0	0.6984 ± 0.0000	1.6606 ± 0.0000
	PAT1	24,612	32 ± 0	14 ± 0	0.3741 ± 0.0007	0.4922 ± 0.0011
	PAT2	2,423*	34 ± 0	11 ± 0	0.3816 ± 0.0000	0.4903 ± 0.0000
	KER	11,017	12 ± 0	8 ± 0	0.0681 ± 0.0004	0.0913 ± 0.0006

Populations Differentiation and Phylogeographic Structure of Spirochaeta Oligotypes

Independently of the OTU considered, the genetic (Fst) and phylogeographic (Φst) structures between the two closest localities (i.e., Patagonian sites PAT1 and PAT2) were extremely to moderately low. In the case of the OTU40, the genetic diversity and frequencies of Spirochaeta oligotypes were fully homogenous between PAT1 and PAT2, as indicated by the non-significant values of Fst and Φst comparisons (Table 3). Contrarily, higher values of Fst and Φst comparisons were recorded among the three provinces considered in this study (Patagonia, PAT1/PAT2; maritime Antarctica, KGI; Kerguelen Islands, KER) (Table 3). Consistently, the distribution of the Spirochaeta oligotypes was geographically discontinuous across the localities, with various province-specific oligotypes (Supplementary File 4). Two exceptions were observed, in the case of maritime Antarctica and Kerguelen Islands (KGI and KER) for OTU6, and in the case of Kerguelen Islands (KER) and Patagonia (PAT1) for OTU7, with relatively lower values of Fst and Φst (Table 3). The Snn test for phylogeographic structure among all sites was significant with statistic values ≥ 0.5 (OTU6; Snn = 0.50, p-value < 0.0001, OTU7; Snn = 0.57, p-value < 0.0001, OTU40; Snn = 0.50, p-value < 0.0001). All in all, these results indicate the existence of both genetically and geographically differentiated Spirochaeta populations across the three biogeographic provinces sampled.

TABLE 3

Genetic (FST) and phylogeographic structure (ΦST) of the Spirochaeta populations among localities.

OTU	Index	Locality	KGI	PAT1	PAT2	KER
OTU6	FST	KGI	–	0	0	0
		PAT1	0.4428	–	0	0
		PAT2	0.4417	0.0026	–	0
		KER	0.0654	0.3856	0.3829	–
	ΦST	KGI	–	0	0	0
		PAT1	0.5371	–	0	0
		PAT2	0.5360	0.0027	–	0
		KER	0.1041	0.5063	0.4933	–
OTU7	FST	KGI	–	0	0	0
		PAT1	0.3436	–	0	0
		PAT2	0.3254	0.0726	–	0
		KER	0.3902	0.0568	0.1734	–
	ΦST	KGI	–	0	0	0
		PAT1	0.4393	–	0	0
		PAT2	0.5361	0.1341	–	0
		KER	0.3967	0.1312	0.3429	–
OTU40	FST	KGI	–	0	0	0
		PAT1	0.5549	–	0.4505	0
		PAT2	0.5445	<0.0001	–	0
		KER	0.8491	0.7358	0.8614	–
	ΦST	KGI	–	0	0	0
		PAT1	0.8326	–	0.2793	0
		PAT2	0.8278	<0.0001	–	0
		KER	0.9354	0.7398	0.8652	–

Within the 48 oligotypes identified in the OTU6, 11 (∼23%) were exclusive to one of the three provinces, and more than half were exclusive to Kerguelen Islands (Supplementary File 4). Maritime Antarctic and Kerguelen Islands (KGI and KER) shared 27 (∼66%) of their oligotypes. The Patagonian localities (PAT1/PAT2) shared 18 of the 23 total oligotypes (∼78%) observed in this province (Figure 4 and Supplementary File 2). A total of five (∼10%) oligotypes were broadly distributed across all localities. Oligotype network of OTU6 showed short genealogies and the presence of at least five dominant oligotypes. The dominant oligotype in Patagonia (PAT1/PAT2), and the dominant oligotype in maritime Antarctica and Kerguelen Islands (KGI/KER), were separated by a single substitution (Figure 4).

FIGURE 4

Median-joining oligotype networks of the three selected OTUs of Spirochaeta. Each circle represents a distinct oligotype. Colors indicate the locality of provenance. Circle size is scaled on the sum of the oligotype frequencies in each locality, within the corresponding OTU dataset. Raw abundances are provided in Supplementary File 2.

In the case of OTU7, the percentage of exclusive oligotypes was almost the same as the OTU6, with 21 (∼22%) oligotypes exclusive to one of the provinces (Supplementary File 4). The dominant oligotype was different between maritime Antarctic and Kerguelen Islands localities (KGI and KER) (Figure 4). While most oligotypes from Kerguelen Islands (KER) were detected in at least one of the Patagonian sites (PAT2 or PAT1) (∼94%), fewer oligotypes from maritime Antarctica (KGI) were observed in Kerguelen Islands (∼68%) and even fewer in the Patagonian localities (∼62%) (Figure 4).

For OTU40, we recorded a predominant group of oligotypes specific to the Patagonian sites representing 65% of the oligotypes identified within the OTU40 (Supplementary File 4 and Figure 4). A clear separation was observed in the oligotype network between the KGI/KER and PAT1/PAT2 localities (Figure 4), with only four shared oligotypes (∼8%) (Supplementary File 2). In maritime Antarctica (KGI), three of the four oligotypes were exclusive, whereas the dominant one from Kerguelen Islands (KER) was shared with maritime Antarctica (KGI) (∼8%) (Figure 4 and Supplementary File 4).

Gene Flow Under a Migration–Drift Equilibrium Model

All the analyzed OTUs showed high genetic similarities between the analyzed Patagonian populations (PAT1 and PAT2). Gene flow analyses identified a bidirectional pattern from PAT2 to PAT1 (effective number of migrants per generation, Nm > 4) and from PAT1 to PAT2 (Nm > 14) (Table 4). The connectivity between the Patagonian and maritime Antarctic localities showed relatively low values of Nm, ranging from 0.001 (OTU40, from PAT2 to KGI) to 0.64 (OTU7, from KGI to PAT1) (Table 4). The OTU6 and OTU7 were both characterized by a substantial gene flow between maritime Antarctica and Kerguelen Islands that was stronger in the direction KGI to KER (OTU6, Nm = 9.81 and OTU7, Nm = 1.39) than in the direction KER to KGI (OTU6, Nm = 2.63 and OTU7, Nm = 0.41) (Table 4). Contrarily, an unidirectional and low gene flow from KER to KGI (Nm = 0.87) was recorded for the OTU40 (Table 4). Finally, the connectivity between Patagonian (PAT1 and PAT2) and Kerguelen Islands (KER) localities was illustrated by three distinct patterns; a low-intensity flow (Nm < 0.5) predominant in the direction PAT1/PAT2 to KER for the OTU6, a substantial flow (Nm values from 2.00 up to 24.41) predominant in the direction KER to PAT1/PAT2 for the OTU7, and an absence of connectivity (Nm < 0.03) in the case of the OTU40 (Table 4). The gene flows are summarized in Figure 5.

TABLE 4

Effective numbers of migrants per generation (Nm) among Spirochaeta populations of the OTU6, OTU7, and OTU40.

OTU	From	θ± SD	To	Nm ± SD
OTU6	KGI	0.003 ± 0.0003	KER	9.81 ± 3.14
			PAT1	0.12 ± 0.06
			PAT2	0.14 ± 0.09
	KER	0.004 ± 0.0004	KGI	2.63 ± 1.19
			PAT1	0.12 ± 0.05
			PAT2	0.07 ± 0.05
	PAT1	0.003 ± 0.0004	KGI	0.14 ± 0.05
			KER	0.43 ± 0.15
			PAT2	24.99 ± 7.10
	PAT2	0.005 ± 0.0005	KGI	0.20 ± 0.07
			KER	0.27 ± 0.16
			PAT1	14.57 ± 4.93
OTU7	KGI	0.002 ± 0.0003	KER	1.39 ± 0.60
			PAT1	0.64 ± 0.23
			PAT2	0.03 ± 0.02
	KER	0.006 ± 0.0006	KGI	0.42 ± 0.13
			PAT1	24.41 ± 8.23
			PAT2	2.00 ± 0.54
	PAT1	0.015 ± 0.0029	KGI	0.16 ± 0.06
			KER	0.41 ± 0.16
			PAT2	4.69 ± 1.82
	PAT2	0.004 ± 0.0005	KGI	0.05 ± 0.03
			KER	0.13 ± 0.09
			PAT1	19.95 ± 9.57
OTU40	KGI	0.002 ± 0.0002	KER	0.09 ± 0.04
			PAT1	0.01 ± 0.01
			PAT2	0.02 ± 0.02
	KER	0.001 ± 0.0002	KGI	0.87 ± 0.21
			PAT1	0.02 ± 0.02
			PAT2	0.02 ± 0.02
	PAT1	0.005 ± 0.0008	KGI	0.01 ± 0.01
			KER	0.03 ± 0.02
			PAT2	28.65 ± 10.15
	PAT2	0.005 ± 0.0008	KGI	0.00 ± 0.00
			KER	0.06 ± 0.02
			PAT1	15.94 ± 10.48

FIGURE 5

Gene flows summary and proposed dispersal routes across the Southern Ocean for each Spirochaeta taxa. Only the gene flows with Nm values > 0.25 are graphically represented. Discontinuous arrows represent Nm values > 0.25 and <1, continuous arrow represent Nm values > 1. Continuous arrows’ width is proportional to Nm values.

Contribution of Contemporary Selection Versus Historical Processes in Shaping the Spirochaeta Microdiversity

For each of the three selected Spirochaeta OTUs, neutral ecological processes were essential in shaping the population composition turnover in Abatus gut membrane. According to the quantitative parsing of ecological processes, the composition of Spirochaeta population was mostly driven by ecological drift, ranging from 50% (OTU40) to 74% (OTU6) of turnover, followed by dispersal limitation, ranging from 12% (OTU6) to 20% (OTU40) of turnover, and homogenizing selection, ranging from 9% (OTU6) to 19% (OTU40) of turnover (Table 5). Overall, deterministic processes (i.e., homogeneous and variable selection) did not account for more than 10% of the populations’ turnover.

TABLE 5

Quantitative parsing of ecological processes driving populations turnover within Spirochaeta OTUs.

	Ecological processes contributions
Spirochaeta OTU	Homogeneous selection (%)	Homogenizing dispersal (%)	Ecological drift (%)	Dispersal limitation (%)	Variable selection (%)
OTU6	2.7	8.8	74.0	12.1	2.3
OTU7	0.4	10.3	63.7	17.6	8.0
OTU40	0.3	19.4	49.6	21.7	9.1

The MMRR approach was used to disentangle the relative effect of geographic distance environmental abiotic differences on the Spirochaeta oligotype compositions between samples. The geographic distance matrix was linearly correlated to the abundance-based similarity matrix of Spirochaeta population composition for OTU7 and OTU40, explaining about 31 and 67% of the observed variation, respectively (Table 6). In contrast, the geographic distance did not significantly impact Spirochaeta oligotype composition for OTU6 (Table 6). Whatever the OTU considered, the environment distance had a significant but slight contribution (<4% of the observed variation) to the Spirochaeta population composition (Table 6), and the global R2 of the model (environmental and geographic distance) did not exceed 24%.

TABLE 6

Multiple matrix regression with randomization (MMRR) to quantify the relative effects of isolation by distance (IBD) and isolation by environment (IBE) on oligotypes assemblage within Spirochaeta OTUs.

Spirochaeta OTU	Model	Coefficient	t Statistic	t p-Value	F statistic	F p-value	R2
OTU6	IBE	0.037	17.82	<0.001	187.32	<0.001	0.208
	IBD	−0.158	−2.59	0.078	NA	NA	NA
OTU7	IBE	0.020	12.10	<0.001	162.43	<0.001	0.165
	IBD	0.309	6.11	<0.001	NA	NA	NA
OTU40	IBE	0.022	6.45	<0.001	102.08	<0.001	0.235
	IBD	0.667	7.97	<0.001	NA	NA	NA

Discussion

In this study, we coupled 16S rRNA metabarcoding and oligotyping algorithm to reveal the microdiversity within three bacterial OTUs affiliated to the Spirochaeta genus, and co-distributed across Patagonia, Kerguelen Islands, and maritime Antarctic provinces of the SO. Through this innovative approach, we identified numerous oligotypes within each of the Spirochaeta OTUs. These oligotypes, corresponding to Spirochaeta sub-taxa, were characterized by contrasting geographic distribution and high levels of 16S rRNA gene similarity (>97%). Taking advantage of the populational taxonomic resolution provided by the oligotype definition, we depicted the Spirochaeta biogeographic patterns across the analyzed provinces in the SO, using various tools adapted from population genetics classically applied in phylogeographic study of macroorganisms’ models. Despite its low substitution rate [approximately 1% in 50 million years (Espejo and Plaza, 2018)], our study demonstrates that the 16S rRNA gene shows value in the evaluation of the microdiversification, as it offers the best compromise between an informative genetic signal, and robust screening of global microbial diversity at intra-OTU level, in a wide range of barely unknown habitats (Chase and Martiny, 2018).

Unlike the studies with macroorganisms, which are usually more demanding in terms of individual sampling effort, we benefit here from the high sequencing depth provided by the metabarcoding of a low-diversity habitat (i.e., the Abatus gut tissue). This allows a robust coverage of the Spirochaeta diversity (up to 180,000 sequences per OTU), and thus high precision of the oligotypes frequencies. Our methodology echoes the “metaphylogeographic” approach recently proposed by Turon et al. (2020) to investigate eukaryotic intraspecies diversity through COI (cytochrome c oxidase subunit I) gene amplicon-sequencing and an oligotyping-like cleaning protocol of the reads based on entropy variation. We propose to expand the concept of “metaphylogeography” to the prokaryotes since it permits phylogeographic inferences of uncultured microbes from a wide range of habitats.

The β-diversity analysis performed at the Spirochaeta genus level revealed that each of the three geographic provinces might host specific Spirochaeta OTUs representing distinct phylogenetic lineages. We also reported a non-random distribution trend with contrasting patterns of Spirochaeta OTU compositions across the localities. Nevertheless, about half of the Spirochaeta OTUs exhibited a broad distribution encompassing Patagonia, maritime Antarctica, and the Kerguelen Islands located more than 7,000 km to the east. This result suggests that despite being mostly detected in Abatus gut, and to a lesser extent in marine benthic sediments, some Spirochaeta representatives would disperse through the SO currents. Concordantly, previous campaigns of high-throughput sequencing of the ocean water column have consistently reported the presence of free-living Spirochaeta OTUs in surface to mesopelagic water, away from the coastlines (Pesant et al., 2015).

For each of the three assessed OTUs, the Spirochaeta populations were expected to be remarkably homogeneous between the two Patagonian sites due to the geographic vicinity and the absence of an evident oceanographic barrier. Consistently with this assumption, and for each of the three OTUs, most of the Spirochaeta oligotypes were shared, the lowest genetic and phylogeographic structures were reported, and high levels of gene flow were recorded between these two sites. Similarly, low or absent differentiation patterns along the Atlantic coast of Patagonia were previously reported for marine Patagonian macroorganisms, including notothenioid fishes (Ceballos et al., 2016), scorched mussels (Trovant et al., 2015), and pulmonate gastropods (Iriarte et al., 2020) presumably due to their high dispersal potential and the ecological continuum of the sampled localities that may conform a same biogeographic province connected through the equator-ward Falkland current (Phillpot, 1985; Arkhipkin et al., 2004). Further phylogeographic studies focusing on microbial taxa of additional sampling sites from Atlantic Patagonia should confirm the microbial biogeographic consistency of this province.

Between Patagonian and maritime Antarctic provinces, Spirochaeta populations exhibited strong genetic and phylogeographic structures, and low levels of gene flow were estimated between these two provinces. These results corroborate our hypothesis that the APF hinders individual dispersion and genetic homogeneity among bacterial populations and suggest that the geographically structured Spirochaeta populations from these two provinces are genetically diverging over time (Nagylaki, 1980; Charlesworth, 2003). Previous studies focusing on diverse macroorganisms taxonomic groups of the SO have evidenced the critical role of the APF on biogeographic patterns, as an open-ocean barrier inducing a genetic break between South America and Antarctica (e.g., ribbon worms (Thornhill et al., 2008); brittle stars (Hunter and Halanych, 2008); notothenioid fishes (Hüne et al., 2015); limpets (González-Wevar et al., 2012, 2019; González−Wevar et al., 2017); sea urchins (Díaz et al., 2011). Regarding the microbial distribution patterns, significant β-diversity differences between prokaryotes assembly from both sides of the APF have been reported in the past, but most of the studies focused on global community in the water column, at high taxonomic resolution [summarized in Flaviani et al. (2018)]. Here, we extend this discontinuity in bacterial diversity to a fine taxonomic resolution (i.e., intra-OTU), revealing province-restricted oligotypes and strong genetic and phylogeographic structure between Patagonian and maritime Antarctic Spirochaeta populations.

Contrarily, and despite the substantial geographic distance separating the sub-Antarctic Kerguelen Islands and the Patagonian and Antarctic sites (>6,500 km), population genetic analyses suggest the existence of some level of connectivity between Kerguelen and the other sites. These findings support a potential dispersion of Spirochaeta taxa from Patagonia and maritime Antarctica to the Kerguelen Islands, and contrariwise, from the Kerguelen Islands to Patagonia and maritime Antarctica. As evidenced by the numerous shared oligotypes, such connectivity would maintain a sufficient gene flow among these provinces to partially counteract the genetic divergence driven by selection, mutation, and genetic drift, inducing oligotypes mixing, and limiting the spatial differentiation of oligotypes assembly (Martiny et al., 2006; Orsini et al., 2013). Moreover, we suggest that this gene flow is not bidirectional, but governed by exclusively eastward oriented dispersion routes (Figure 5), following the major and constant flow of the ACC (Güller et al., 2020). Under this scenario, Spirochaeta individuals from Kerguelen Islands may seed toward Patagonia following the ACC eastward flow around Antarctica. Such ACC-mediated connectivity among sub-Antarctic provinces (Patagonia and Kerguelen Islands) is well known in a wide range of benthic macroorganisms populations, such as buoyant kelps Durvillaea antarctica and Macrocystis pyrifera (Macaya and Zuccarello, 2010), and several kelp-associated macroinvertebrates (Leese et al., 2010; Nikula et al., 2010; Cumming et al., 2014; González−Wevar et al., 2018; Güller et al., 2020). Occasionally, Spirochaeta individuals from Kerguelen Islands may also be able to reach the maritime Antarctic province. Such pattern has been recently reported for the southern bull kelp D. antarctica, a typical sub-Antarctic macroalgae, which is transported by rafting to as far as the West Antarctic Peninsula coasts, pushed by the circumpolar flow of the ACC or by storms leading to the occasional crossing of the APF (Fraser et al., 2020).

Several studies have provided evidence of a high dispersal capacity of marine bacteria by comparing community composition mostly at high taxonomic resolution (e.g., class, genus, or OTU) among various water masses and oceanic regions (Sul et al., 2013; Sunagawa et al., 2015; Milici et al., 2017; Logares et al., 2020). Particularly, the most abundant marine bacteria are supposed to migrate between adjacent regions through passive transport (Sul et al., 2013). An innovative conceptual framework called “Microbial Conveyor Belt” (MCB) has been proposed by Mestre and Höfer (2020), to emphasize that the marine microorganisms’ dispersion would not merely rely on passive and stochastic dispersal, but instead on the adaptation of life-history traits (e.g., dormancy stage). These traits would allow microorganisms to successfully and recurrently disperse in unfavorable habitats through specific dispersion avenues (Locey, 2010). Here, we provided empiric results from Spirochaeta population based on genetic data supporting a partial MCB in the SO driven by the ACC. Unfortunately, details about the benthic Spirochaeta taxa’s ecology are scarce, with a single isolated strain from subseafloor sediment (Miyazaki et al., 2014). Thus, the life-history traits of Spirochaeta, as the sporulation capacity, remain to be investigated to further understand its distribution pattern in the SO. Nevertheless, in order to disperse, we propose that Spirochaeta individuals (enriched in the digestive tract) could be released from the host gut toward the surrounding benthic sediments through fecal pellets. Such enrichment of the digesta with taxa from the host microbiota, as well as the presence of Spirochaeta within the fecal material, have been demonstrated in the sea urchin species Lytechinus variegatus (Hakim, 2015). The released Spirochaeta individuals may be resuspended in the water column through the action of one or several processes such as upwelling, bioturbation by the benthic deposit-feeders, or water column mixing during winter (Meysman et al., 2006; Petro et al., 2017; Mestre and Höfer, 2020). Once in the water column, these Spirochaeta individuals may disperse over large geographic scales, transported through oceanographic features (e.g., currents, punctual meteorological events) (Mestre and Höfer, 2020). The attachment to suspended particulate matter, either biotic [e.g., hitchhiking on zooplankton (Grossart et al., 2010) and seaweed (Serebryakova et al., 2018)] or abiotic [e.g., microplastics (Bowley et al., 2020), known to have a long-distance dispersion potential], may also contribute to the bacterial spreading in the oceans (Milici et al., 2017; Mestre and Höfer, 2020).

The marine prokaryote communities are usually considered widely dispersed and mainly shaped by contemporary ecological processes such as environmental filtering (Sunagawa et al., 2015; Louca et al., 2016). By applying the ecological framework developed by Stegen et al. (2013) to oligotypes data, we found contrarily that ecological drift was the predominant stochastic mechanism shaping intra-populations turnover within Spirochaeta taxa across the SO. Our previous study of the Abatus gut microbiota showed that non-neutral processes drove the bacterial community at the OTU level in the host gut tissue (Schwob et al., 2020). While deterministic processes are usually prevalent in structuring microbial communities’ assembly at a higher taxonomic resolution (Martiny et al., 2011; Wang, 2013; Sintes et al., 2015; Larkin and Martiny, 2017), the stochastic mechanisms tend to have a more significant contribution at finer taxonomic scales (Logares et al., 2020), since niche overlapping and functional redundancy enhance the susceptibility of populations to drift (Zhou and Ning, 2017). Thus, the biogeographic structure observed among Spirochaeta populations might result from stochastic birth, death, disturbance, emigration, and immigration events rather than oligotype-sorting through the biotic and abiotic environmental variations (Martiny et al., 2006; Nemergut et al., 2013; Vellend, 2016). Consistently, the MMRR analysis revealed that the isolation-by-environment (IBE) model might account for a low percentage of the Spirochaeta oligotypes turnover. Altogether, these results tend to validate the strategy applied in our study, that is, to focus on specialist bacterial taxa hosted in sibling sea urchin species with the same habitat preferences, in order to homogenize the environment, to reduce the diversity, to soften the deterministic selection driven by environmental variations, thus leading to the maximization of the detection of neutral micro-evolutive processes associated with biogeography (Langenheder and Lindstrom, 2019).

By analogy with the genetic drift, whereby changes in gene frequencies occur solely by chance in a population (Chase et al., 2011), our result suggests that the microdiversity observed within the Spirochaeta taxa would be mostly generated by genetic drift without any adaptive implications. An earlier study reported that microdiversity observed in the 16S rRNA gene of marine coastal Vibrio splendidus isolates was ecologically neutral (Thompson et al., 2005). Nevertheless, we cannot discard that, while the microdiversity within the V4–V5 of 16S rRNA gene-targeted here is likely to be acquired through neutral processes (Chust et al., 2016), it may also be associated with substantial modifications in niche-defining traits and functional attributes specific of the Spirochaeta strains, driven by deterministic processes, in order to cope with local conditions (Zure et al., 2017; Chase and Martiny, 2018). Further studies will need to focus on other loci (e.g., functional genes), potentially under selection, as they are expected to display a higher degree of differentiation among populations and to provide an insight into the ecology of the Spirochaeta sub-taxa (White et al., 2010).

Given the highly specific targeted habitat (i.e., gut tissue) and the intricate association between Spirochaeta and Abatus, it remains unclear whether or not the host contributes to the Spirochaeta genetic and phylogeographic pattern reported in the gut populations. Indeed, the Abatus guts may have individual colonization histories also contributing to the microbial community assembly process at the population level. However, due to the limited sequencing coverage of Spirochaeta diversity outside the host, we cannot evaluate to what extent the most abundant oligotypes from the gut tissue, also detected in the external sediment samples, would be representative of the composition of free-living populations, thus preventing a comprehensive understanding of the host contribution to Spirochaeta oligotype assembly. Ultimately, conclusion about the predominant effect of stochastic ecological processes onto Spirochaeta oligotype composition need to be nuanced, as potential historical contingencies related to the host can eventually shape the final state of microbiota composition through both stochastic and/or deterministic processes (Obadia et al., 2017). In strong instance, host early-life microbial colonization by early colonizers benefiting from a “priority effect” may influence the order and timing of the historical sequence of species arrival within the local community (Fukami, 2015). Under the hypothesis of a predominantly stochastic gut colonization by random microbes, diverging chronologies of the ensuing microbial species would be expected, leading toward different stable states and a noisy community assembly across Abatus populations (Obadia et al., 2017). Contrastingly, we reported patterns of high dominance and prevalence of Spirochaeta across Abatus individuals, and high homogeneity of oligotypes composition within populations (data not shown), rather suggesting that Spirochaeta may benefit from a priority effect during Abatus early-life, and that the subsequent interactions may deterministically constrain gut colonization success (Fukami, 2015; Higgins et al., 2020). Moreover, while there is no evidence of vertical transmission in Abatus, the transovarial transmission of a bacterial endosymbiont was reported in the regular sea urchin species Heliocidaris erythrogramma (Carrier et al., 2021) and maternal-inherited bacterial communities were detected in unfertilized eggs of Strongylocentrotus purpuratus (Carrier and Reitzel, 2019). The presence of preexisting Spirochaeta strains within the gut may also limit the ensuing horizontal colonization and constrain the microbiota composition toward a homogenous and stabilized state (Renelies-Hamilton et al., 2021). The presence of vertically transmitted Spirochaeta, the compositional transition states of Abatus gut microbiota, and the biogeographic structure of free-living Spirochaeta populations in the external sediment remain to be further investigated.

Notwithstanding the consistency of the global phylogeographic and connectivity patterns depicted across the three tested OTUs, we also reported some differences according to the taxa considered, which might be related to different ecotypes with distinct ecological niches or different dispersal capacity. For instance, various marine bacterial taxa, such as the cyanobacteria Synechococcus or the Vibrio populations, demonstrate fine-tuning of their physiology by accumulating microdiversity in functional genes through duplication events, SNPs, and allelic variants (Everroad and Wood, 2012; Shapiro et al., 2012). Alternatively, these differences may also be related to the Spirochaeta OTU abundance, since the more relatively abundant the Spirochaeta populations were (i.e., higher number of sequences retrieved from the gut tissue through the metabarcoding approach), the more they tend to exhibit cosmopolitan oligotypes (i.e., detected across each of the four localities). It is not unreasonable to infer that a larger population may have more chance to migrate and successfully reach a suitable habitat, while small-size populations may be more likely diluted along the dispersal route with no/too few dispersive particles to establish in the new habitat (Sul et al., 2013).

The diversity units defined by 16S rRNA gene sequences are generally considered as insensitive to diversification resulting from dispersal limitation (Hanson et al., 2012). Contrastingly, we reported that the dispersal limitation was the second most crucial ecological factor driving the turnover of Spirochaeta oligotypes, and by extension, their genetic divergence. Dispersal limitation is classically considered as a historical factor since current oligotypes assemblage results from past dispersal limitations (Martiny et al., 2006). Our result indicates that the potentially suitable habitats are too distant (Heino et al., 2015), or inaccessible due to the existence of oceanic currents (Sul et al., 2013; Logares et al., 2018), hence limiting the homogenization of Spirochaeta oligotypes’ frequencies across populations and allowing the neutral genetic divergence of genomic regions overtime via genetic drift (Nosil et al., 2008; Hanson et al., 2012). Note that our results obtained from distinct methodologies [i.e., the genetic differentiation and phylogeographic structure, the contribution of dispersal limitation from Stegen et al. framework (Stegen et al., 2013), and the contribution of the geographic distance isolation by distance (IBD) from the MMRR analysis] were highly consistent with each other, and across the three selected Spirochaeta OTUs. For instance, the OTU40 that harbored the overall highest value of genetic divergence was also characterized by the highest estimated contribution of geographic distance and dispersal limitation, thus supporting the interrelation between genetic divergence and oligotypes population turnover, and the overall consistency of the approach implemented.

Conclusion

Our study highlights the application of V4–V5 16S rRNA gene metabarcoding and oligotyping approach as rapid, robust, and resolutive enough to unravel marine bacterial phylogeographic patterns and detect genetic connectivity among the SO provinces. Taken together, the three Spirochaeta OTUs analyzed evidence three consistent phylogeographic patterns, classically observed in the studies involving benthic macroinvertebrates across the SO: (1) a high populational and genetic homogeneity within the Patagonia province, (2) a strong barrier to dispersal between Patagonia and maritime Antarctica due to the APF, resulting in a high differentiation of Spirochaeta populations, and (3) the existence of connectivity between sub-Antarctic provinces of the Kerguelen Islands and Patagonia, and from Kerguelen Islands to the maritime Antarctic, due to the ACC-mediated connectivity. Nevertheless, as connected as these provinces are, the gene flow does not seem to be strong enough to prevent the ongoing intraspecific microdiversification of the Spirochaeta taxa. The microdiversity of Spirochaeta, underlying these biogeographic patterns, is essentially driven by historical processes, such as ecological and genetic drift, and dispersal limitation related to the SO’s oceanographic features. In the future, extending this framework to other localities and taxonomic groups will contribute to the comprehensive understanding of the SO microbiota.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/ PRJNA658980, https://www.ncbi.nlm.nih.gov/ PRJNA590493, and https://www.ncbi.nlm.nih.gov/ PRJNA659050.

Author Contributions

EP, JO, and LC designed the study. EP and JO organized the sampling missions. EP, JO, LC, and GS collected samples. GS extracted the DNA, organized sequencing and managed data mining and analyses. NS contributed to the gene-flow and the MMRR analyses and designed the illustrative maps. GS, JO, NS, CG-W, and EP interpreted the results. GS wrote the manuscript. All authors contributed substantially to manuscript revisions and read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.