Structure de l’écosystème planctonique : apport des données à haut débit de séquençage et d’imagerie

Structure globale de l’écosystème planctonique

Préambule Dans l’introduction, la complexité des écosystèmes a été présentée. Ces écosystèmes sont composés de très nombreuses entités en interaction, caractérisables à de multiples échelles spatio-temporelles et biologiques, et sur lesquelles différents processus évolutifs et écologiques s’exercent. Ainsi, comprendre le fonctionnement et les mécanismes qui structurent les écosystèmes est un réel défi, qui pourrait bénéficier des approches intégratives (Michener and Jones, 2012, Garcia, 2019, Zipkin et al., 2021). Comme souligné dans l’introduction, ces dernières décennies ont été marquées par une augmentation très rapide de la quantité et de la diversité de données acquises. Ces données massives pourraient permettre de mieux comprendre les écosystèmes planctoniques. Dans les deux précédents chapitres, différentes sources de données ont été utilisées. Des structures de l’écosystème différentes ont été identifiées, mettant en évidence que les multiples sources de données contiennent des informations variées mais potentiellement complémentaires. Intégrer l’ensemble de ces données dans une mˆeme analyse pourrait fournir un signal plus riche et représentatif de la structure réelle de l’écosystème. Cette étude a été menée dans la région subtropicale de l’Océan Atlantique Nord. L’océan Atlantique Nord abrite le système Atlantic Meridional Overturning Circulation (AMOC), profondément impliqué dans la régulation du climat et dont les modifications  peuvent avoir un impact sur les écosystèmes océaniques et continentaux à l’échelle globale (Caesar et al., 2018; Liu et al., 2020). Dans sa partie subtropicale, l’océan Atlantique Nord présente un large gyre anticyclonique oligotrophe supporté par plusieurs courants (i.e., Gulf Stream, front subpolaire, courant des Canaries, courant nord-équatorial). Cette région présente des gradients environnementaux et biologiques (Aiken et al., 2000) et des régions stables aux bordures fluctuant spatialement au cours du temps y ont été délimitées (Longhurst, 2010; Reygondeau et al., 2013). De nombreuses études ont été réalisées afin de mieux comprendre cet écosystème. Par exemple, la zone a été échantillonnée pendant plus de 30 ans grˆace au Continuous Plankton Recorder (CPR), ce qui a permis un suivi temporel de plusieurs groupes planctoniques (Batten et al., 2019). Ces études ont été le socle d’une augmentation du savoir sur la diversité dans la zone et de ses réactions au changement climatique (e.g., Beaugrand et al., 2009; Barton et al., 2016; Vannier et al., 2016; Brun et al., 2019). Cependant, l’échantillonnage CPR présente un biais spécifique aux espèces en raison du comportement d’évitement des organismes, ce qui peut entraˆıner des problèmes dans la quantification de la composition de la communauté (Batten et al., 2019). De plus, l’échantillonnage CPR peut endommager les organismes (en particulier les gélatineux, Batten et al., 2019), ne fournissant pas d’informations morphologiques précises pour tous les groupes. Les identifications à fine résolution taxonomique peuvent également ˆetre difficiles et longues (Batten et al., 2019) et certains organismes écologiquement importants, comme le picoplancton, ne sont pas échantillonnés. Ainsi, le CPR ne semble pas idéal pour étudier les communautés dans leur ensemble. La combinaison des données de séquen¸cage et d’imagerie pourrait permettre de surmonter certaines de ces limites. En les intégrant avec les données environnementales, de multiples niveaux de structuration de l’écosystème planctonique pourraient ˆetre explorés simultanément. L’objectif de ce chapitre est de réaliser l’analyse simultanée des différents types de jeu de données, afin d’identifier la structure globale de l’écosystème. Les sous-parties suivantes correspondent au manuscrit intitulé Multi-scale exploration of the planktonic ecosystem structuring in the North Atlantic Ocean. 94 5.2. Multi-scale exploration of the planktonic ecosystem in the North Atlantic Ocean 5.2 Multi-scale exploration of the planktonic ecosystem in the North Atlantic Ocean Authors Ophélie Da Silva1 , Alberto Baudena1 , Guillaume Bourdin2 , Paul Frémont3 , Rachel L. Kelly4 , Louis Caray–Counil1 , Laetitia Jalabert1 , Nicolas Henry5,6, Frédéric Mahé7,8, Julie Poulain5 , Sarah Romac5,6, Sakina-Dorothée Ayata1 , Chris Bowler6,9, Seth John4 , Samuel Chaffron6,10, Colomban de Vargas5,6, Damien Eveillard ´ 6,10, Olivier Jaillon3,6 , Matthew B. Sullivan6,11, Patrick Wincker3,6, Fabien Lombard1,6,12, Lionel Guidi1,6 1 – Sorbonne Universités, CNRS, Laboratoire d’Oceanographie de Villefranche, LOV, F-06230 Villefranche-sur-Mer, France. 2 – School of Marine Sciences, University of Maine, Orono, ME, United States 3 – Génomique Métabolique, Genoscope, Institut Fran¸cois Jacob, CEA, CNRS, Univ Evry, Université Paris-Saclay, Evry, France 4 – Department of Earth Sciences, University of Southern California, Los Angeles, CA, USA 5 – Sorbonne Université, CNRS, Station Biologique de Roscoff, UMR 7144, ECOMAP, 29680 Roscoff, France 6 – Research Federation for the study of Global Ocean Systems Ecology and Evolution, FR2022/Tara Oceans GOSEE, 3 rue Michel-Ange, 75016 Paris, France 7 – CIRAD, UMR PHIM, F-34398 Montpellier, France 8 – PHIM, Univ Montpellier, CIRAD, INRAE, Institut Agro, Montpellier, France 9 – Institut de Biologie de l’ENS (IBENS), Departement de biologie, Ecole normale supérieure, CNRS, INSERM, Universite PSL, 75005 Paris, France 10 – Université de Nantes, Centrale Nantes, CNRS, LS2N, F-44000 Nantes, France 11 – Department of Microbiology, The Ohio State University, Columbus, OH, USA 12 – Institut Universitaire de France (IUF), Paris, France Abstract Deciphering the structure of ecosystems is essential for predicting their response to their drivers’ variations. Yet because ecosystem structure results from the interweaving processes co-occurring at several spatio-temporal scales, its study is challenging. Planktonic ecosystems embrace different levels of heterogeneous information. They include biotic components that can be characterized by taxonomic or morphological approaches, and  components including biogeochemical or hydrodynamic properties of water masses. The proper integration of all components is a necessary step to understand planktonic ecosystem better. We propose an integrative methodology to a North Atlantic transect based on the high-resolution Tara Pacific imaging and 18S V9 metabarcoding datasets. We have highlighted the existence of two main scales of ecosystem structuring and identified associations between components. Taxonomic composition, morphology of the multicellular-enriched part of communities and biogeochemistry are similarly structured at the basin-scale. However, finer-scale structuring was detected for morphology of the protist-enriched part of the community and for hydrodynamics. Such differences in scales of the components could lead to a highly complex response of the ecosystem and its potential restructuring under climate change. By highlighting the multi-scale structuring of planktonic ecosystems, this study confirms the importance of considering processes occurring over multiple scales. In order to better understand the impacts of these coexisting scales on ecosystem structure, the couplings between the components and the mechanisms by which they interact must be further explored.In natural ecosystems, diversity and abundance patterns result from interactions within and between local communities exchanging individuals (Metacommunity concept, Leibold et al. 2004). For planktonic organisms, dispersal mostly relies on physical transport (advection and diffusion) by ocean currents. In addition, physico-chemical environmental characteristics also shape diversity patterns since they determine survival and fitness of organisms under the environmental filtering process (Kraft et al., 2015). Biotic interactions, on the other hand, are also crucial to understand the structuring of communities because competitive phenomena can lead to the exclusion or dominance of some species in an a priori favorable abiotic environment (i.e., exclusive competition, Kraft et al. 2015). However this process is criticized within the plankton framework since, even in resource limiting conditions, the diversity of plankton tends to remain high in the continuously evolving oceanic environment (Plankton paradox ; Hutchinson 1961; Roy, Chattopadhyay 2007). Patterns of diversity, abundance and interactions of eukaryotic plankton communities in the ocean have been extensively studied thanks to oceanographic campaigns, and the emergence of new technologies such as environmental sequencing and imaging. Recent global surveys have revealed that a large part of eukaryotic taxa and genes are unknown (De Vargas et al., 2015; Carradec et al., 2018). Interactions between taxa have also Multi-scale exploration of the planktonic ecosystem in the North Atlantic Ocean been inferred, leading to the identification of key groups and the importance of biological interactions in the structuring of the community (Lima-Mendez et al., 2015). These different levels of biotic organization have been related to the environment in which organisms evolve. Among others, temperature (Ibarbalz et al., 2019; Salazar et al., 2019) or iron concentrations (Caputi et al., 2019) have been shown to be major drivers of the structuring of plankton communities. At the functional level, diversity declines towards the poles which is critical in a context of climate change (Ibarbalz et al., 2019; Salazar et al., 2019). Finally, the morphology of the organisms seems to condition their dispersal, which impacts their distribution patterns (Villarino et al., 2018; Richter et al., 2020; Sommeria-Klein et al., 2020). All of the above examples illustrate advances of knowledge about planktonic community structuring at global-scale. However, these data have rarely been used for finer-scale community studies, even though the importance of local heterogeneity in their structuring has already been highlighted (d’Ovidio et al., 2010; McGillicuddy, 2016; Lévy et al., 2018). Moreover, most of the studies dealt separately with community abundance, diversity or morphology, and only rarely have attempts been made to link them with the environment or processes of transport and mixing. To our knowledge, to date no study has explored simultaneously all levels of ecosystem structuring. Their interplay thus remains unknown. Yet, ecosystems are complex systems whose emergent properties cannot be studied at the level of the entities composing it (Karsenti, 2008; Ladyman et al., 2013). Thus, such systems could benefit from integrative approaches (Karsenti et al., 2011; Michener, Jones, 2012; Sunagawa et al., 2020; Zipkin et al., 2021), as it is the case in other fields of living sciences (Leonelli, 2013). The Tara Pacific expedition (2016 to 2018), during which the global ocean was sampled at a near mesoscale spatial resolution for omics and morphological data, was conceived to overcome previous limitations (Gorsky et al., 2019). The different observations collected will allow integrating taxonomy, morphology and environment into a unified framework, leading to a better understanding of the interplay between these different components of the ecosystems. Here we focused on a fine resolution North Atlantic transect realized during this expedition. The Atlantic Ocean is a crucial oceanic basin as it houses the Atlantic Meridional Overturning Circulation system, deeply involved in climate regulation and whose variation can impact surrounding continents as well as global-scale ecosystems (Lenton et al., 2008; Liu et al., 2020). Therefore, the North Atlantic Ocean has been intensively sampled in order to better understand this ecosystem. In particular, numerous Continuous Plankton Recorder (CPR) surveys have been carried out over time, allowing temporal monitoring of several planktonic groups (Reid et al., 2003; Batten et al., 2019). Among others, important results include evidence that Calanoid copepods are migrating northward in response to temperature variations (Beaugrand et al., 2009). The impacts of copepod community structure on carbon fluxes have also been investigated, and indicate that climate change, by modifying large-bodied copepod distribution, has modified activity of the biological carbon pump in the North Atlantic Ocean (Brun et al., 2019). Only a few studies have targeted protists but some studies suggest the northward migration of diatoms and dinoflagellates (Barton et al., 2016). While CPR data can provide information on fine spatial and temporal variations of some planktonic groups or species, it does not allow the study of the entire plankton community composition because of the sampling method limitations (Batten et al., 2019). In order to fill this gap, here we analysed environmental, metabarcoding and imaging data generated during the Tara Pacific expedition to simultaneously address community composition and morphological structuring of the planktonic ecosystem in the subtropical North Atlantic Ocean. 5.2.2 Material and Methods Atlantic transect of the Tara Pacific expedition During the Tara Pacific expedition, which took place from May 2016 to October 2018, biological samples were collected and physicochemical characteristics of the water masses were measured. Data were collected following the Tara Pacific Expedition protocols detailed in Gorsky et al. (2019). In our study, we focused on the outward North Atlantic Ocean transect, comprising 27 stations from [46.667°N ; 6.279°W] sampled on May 29th, 2016, to [24.180°N ; 81.406°W] sampled on July 6th, 2016 (see maps Figure 5.1). Environmental context : biogeochemistry and hydrodynamics We used the Tara Pacific biogeochemical context (version 1, DOI :10.5281/zenodo.4068292), including continuous and discrete measurements enriched by satellite and modeling data (Gorsky et al., 2019). Variables with more than 10% of missing values were discarded. For the variables recorded by several instruments, only one representative variable was retained, based on measurement quality (from high to low quality : in situ water sampling, inline system, satellite with a resolution of 1/12° and 1/24°, Copernicus assimilative model). The resulting dataset is composed of 24 environmental variables (D.1). Stations 1, 7 and 8 contained, respectively, 3, 10 and 8 missing values (out of 24) and were discarded from the biogeochemical analysis. Salinity measured with a thermosalinograph was missing at station 9. A linear model was adjusted between this variable and salinity provided by Copernicus, and the missing value was replaced by the predicted one. 98 5.2. Multi-scale exploration of the planktonic ecosystem in the North Atlantic Ocean Particulate Inorganic Carbon (PIC) variable also showed one missing value which was replaced by the average of PIC values at the previous and at the next station. Data were then normalized with a box-cox transformation and standardized. In order to characterize the water mixing and transport properties of the area, satellitederived velocity fields with geostrophic and Ekman components were analysed. From their processing, eight Eulerian and Lagrangian diagnostics were computed, listed as follows (further details are provided in Supp. Info. D.1). Finite Time Lyapunov Exponents (FTLE ; Shadden et al. 2005), a Lagrangian quantity associated with mixing and horizontal stirring, were calculated to identify areas of intense frontal or turbulent activity. The Lagrangian betweenness (Ser-Giacomi et al., 2015) was employed to detect the structures which are crucial for the connectivity of the region : high betweenness values identify bottleneck-like regions in which water particles arrive from various origins and are then spread toward several different destinations. The horizontal divergence was used to identify areas of possible upwelling or downwelling. Indeed, this diagnostic detects water convergence and can be used as a proxy of vertical movements in the upper part of the water column. The horizontal divergence was computed both locally (Eulerian divergence) and along water parcel trajectories (Lagrangian divergence ; Hernandez-Carrasco et al. 2018). The presence and the direction of rotation of vortexes were evaluated from the Okubo-Weiss metric and from the vorticity, respectively. Both these diagnostics, analogously to the horizontal divergence, were computed in Eulerian and Lagrangian ways. In order to provide an estimation of the time spent by the water parcels inside an eddy, in the months prior to the observations, their retention time was calculated (d’Ovidio et al., 2013). Finally, absolute speed and kinetic energy, providing information about the intensity of the horizontal currents, were estimated. Each diagnostic was computed using different sets of parameters. Every set was tested and, in total, 9 metrics synthesizing the information were kept for the analyses (more details in Supp. Info. D.1). Each of these metrics was log modulus transformed (L(x) = sign(x) ∗ log(|x| + 1)) and standardized. Characterization of plankton communities Biological samples collected during the Tara Pacific expedition span from viruses to metazoans. We investigated the metabarcoding datasets including 97 samples (data available on zenodo, doi :10.5281/zenodo.4073035) from 27 surface stations and four consecutive size fractions, ranging from 0.2 µm to 2000 µm, in which eukaryotes were collected (i.e.,0.2-3 µm, 3-20 µm, 20-2000 µm and 300-3000 µm size-fractions). For each sample, the V9-18S rDNA metabarcodes were sequenced, filtered and clustered into Operatio99 Chapitre 5. Structure globale de l’écosystème planctonique nal Taxonomic Units (OTUs). Taxonomic assignments of OTUs were determined by global alignment against the 18S database EukRibo D.1. Since the primer used for amplification can amplified bacteria, a filtering step based on taxonomic assignments was performed to remove bacteria potentially present in the datasets. The resulting datasets corresponded to raw abundance matrices of OTUs in each sample. As raw abundances were not directly comparable, they underwent a total sum scaling transformation. The resulting relative abundances were represented by histograms to visualize differences in composition between samples based on the most abundant OTUs. In addition, to take into account the compositional nature of the data, a clr-transformation (Aitchison, 1999) was performed on raw abundances, and the dissimilarity between each pair of communities was estimated by the Aitchison distance between samples (i.e., Euclidean distance based on clr-transformed abundances). Imaging data, generated during the Tara Pacific expedition with Flowcam (20-200 µm ; Gorsky et al. 2019; Sieracki et al. 1998) and Zooscan (300-2000 µm ; Gorsky et al. 2010, 2019) imaging technologies, were also investigated. Due to heterogeneous processing of samples, Flowcam results from stations 1 to 10 were not considered. We analyzed 40 imaging samples distributed into two datasets : one derived from 17 samples analyzed with Flowcam, the other derived from the 23 samples analyzed with Zooscan. As the Flowcam datasets were voluminous (5,000-15,000 objects per sample), the images composing them have been treated differently according to the size of the objects (i.e., organisms). All objects larger than a threshold of 500 pixel area were analyzed, only 30% of the objects smaller than 500 pixel area were taken into account. Taxonomic assignment of each organism was predicted with an automatic recognition algorithm and manually validated or corrected by taxonomic specialists using EcoTaxa (Picheral et al., 2017). The images and their taxonomic assignments are available at EcoTaxa (under project numbers 1700, 1737 and 1744). We discarded every image corresponding to non-living objects (including bad focus, detritus, artifacts, fibers and bubbles) or seaweed. For each object, 12 morphological features automatically estimated in Ecotaxa were used (Table 5.11). Taxonomic assignments from imaging technologies are coarser than metabarcoding ones. However, morphology of individuals are accessible by imaging while it is not by metabarcoding. Thus, these two datatype are complementary and only morphological information were exploited from imaging. The morphological features were normalized by a Yeo Johnson transformation (Yeo, Johnson, 2000) and community weighted means (CWM ; D´ıaz, Cabido 2001) were computed.