Trophic strategies of intertidal foraminifera explored with single‐cell microbiome metabarcoding and morphological methods: What is on the menu?

Abstract In mudflats, interactions and transfers of nutrients and secondary metabolites may drive ecosystems and biodiversity. Foraminifera have complex trophic strategies as they often rely on bacteria and eukaryotes or on potential symbionts for carbon and nitrogen resources. The capacity of these protists to use a wide range of adaptive mechanisms requires clarifying the relationships between them and their microbial associates. Here, we investigate the interactions of three foraminiferal species with nearby organisms in situ, by coupling molecular (cloning/Sanger and high‐throughput sequencing) and direct counting and morphological identification with microscopy. This coupling allows the identification of the organisms found in or around three foraminiferal species through molecular tools combined with a direct counting of foraminifera and diatoms present in situ through microscopy methods. Depending on foraminiferal species, and in addition to diatom biomass, diatom frustule shape, size and species are key factors driving the abundance and diversity of foraminifera in mudflat habitats. Three different trophic strategies were deduced for the foraminifera investigated in this study: Ammonia sp. T6 has an opportunistic strategy and is feeding on bacteria, nematoda, fungi, and diatoms when abundant; Elphidium oceanense is feeding mainly on diatoms, mixed with other preys when they are less abundant; and Haynesina germanica is feeding almost solely on medium‐large pennate diatoms. Although there are limitations due to the lack of species coverage in DNA sequence databases and to the difficulty to compare morphological and molecular data, this study highlights the relevance of combining molecular with morphological tools to study trophic interactions and microbiome communities of protists at the single‐cell scale.


| INTRODUC TI ON
Intertidal mudflats host abundant and diverse microbial communities that play major roles in primary production, food web, biogeochemical cycles, and sediment stabilization (Cesbron et al., 2016;Lebreton et al., 2019;Lubarsky et al., 2010;MacIntyre et al., 1996;Miller et al., 1996). Among these communities, the microphytobenthos (MPB) is composed of an assemblage of benthic photosynthetic microalgae and cyanobacteria often dominated by diatoms Méléder et al., 2007). It is a major contributor of mudflats primary production and a food source for heterotrophs (Blanchard et al., 2001;Miller et al., 1996;Underwood & Kromkamp, 1999). Furthermore, microbial eukaryotes/prokaryotes interactions are essential to marine ecosystems as they facilitate nutrient recycling, photosynthetic activity, and secondary metabolite production (Amin et al., 2012).
Nevertheless, it has not been proven yet that the kleptoplasts are photosynthetically active in E. oceanense and E. selseyense. Haynesina germanica (Ehrenberg, 1840) has been shown to feed on large diatoms (Austin et al., 2005) and to be a photosynthetically active kleptoplastic species (Jauffrais et al., 2016;Jesus et al., 2021;LeKieffre et al., 2018;Lopez, 1979). A recent study using a metabarcoding approach confirmed the omnivorous diet of Ammonia and the kleptoplastic activity of E. selseyense and H. germanica .
Investigating these mechanisms requires clarifying the relationships between foraminifera and their microbial associates. Improving our knowledge on foraminiferal trophic interactions would allow to better understand this understudied group and its role in the ecosystem functioning and in the biogeochemical cycles.
In the present study, we combine molecular (cloning/Sanger sequencing and high-throughput sequencing or HTS) with morphological (granulometric measurements and optical microscopy observations) approaches to investigate the identity of organisms interacting with foraminifera in situ. Both approaches have pros and cons. Microscopy allows to count specimens, but is more limited for species identification, whereas eDNA is more precise for species identification, but has only semiquantitative resolution. Here, we define the microbiome as the nonforaminiferal DNA sequenced from foraminifera, which can originate from symbionts, commensals, parasites, decomposers, or preys. Three to five specimens of three foraminiferal species are collected from three sites in the Bourgneuf mudflat (France). Single foraminifer extractions are used to sequence bacteria and chloroplasts with the 16S rDNA marker and eukaryotes with the 18S rDNA marker to investigate organisms interacted with the foraminifera.
Samples are sequenced using HTS to get an overview of the diversity of taxa found associated with foraminifera. In addition, we used a cloning/Sanger sequencing on one sample/species/site to gather the most abundant and representative sequences present in the different foraminiferal microbiomes. We did so because this approach allows obtaining longer sequences than HTS, and therefore a taxonomic assignment of much higher quality and resolution. In parallel, optical microscopy observations are performed to count foraminifera and diatoms from fixed volumes of sediment to get an estimation of the densities of both groups in situ at the time of sampling. We expect that combining these methodologies, which is new for foraminiferal studies, will improve our results by getting advantages of eDNA for species identification, and microscopy for density estimation. For the granulometric analysis, the superficial sediment layer (~ first centimeter) was scrapped with a spoon in the three stations and brought back to the laboratory in a cooling box and frozen on arrival.

| Sampling site and granulometry
For each station, 1 g of sediment was prepared and analyzed through liquid dispersion with a laser diffraction particle size analyzer (Malvern Mastersizer 3000E, Malvern Instruments) at the UMR 6112 LPG (University of Angers). This analysis allowed to define sediment grain size by the relative abundance (% volume) of silt (Ø < 63 μm) and sand (63 < Ø < 2000 μm) according to the Udden-Wentworth's scale.

| NDVI and microphytobenthic assemblages in the first millimeters
To retrieve the MPB biomass for each station, a SPOT image was analyzed following Méléder et al. (2003) and Echappé et al. (2018).
This image was acquired from the sampled area by the SPOT 7 satellite, with 6 m of spatial resolution, on the 12th of September 2015 at 10:54 GMT, that is, 1:28 after the maximum low tide, which was 1.30 m (Echappé et al., 2018). Reflectance data from each pixel of the image were translated into NDVI (Normalized Difference Vegetation Index) values (Figure 1c), used as a proxy of the MPB biomass due to the chlorophyll a absorption (Benyoucef et al., 2014;Brito et al., 2013;Méléder et al., 2003), and averaged over an area of ~600 m 2 (i.e., 20 pixels) around each sampling station and are compared (ANOVA and Tukey-test).
F I G U R E 1 Localization of the three sampling stations in the Bay of Bourgneuf: (a) France with a star indicating the studied region, (b) Bay of Bourgneuf with a star indicating the sampling area, (c) the stations H17, H18, and H19 near to oyster reefs (indicated by arrows). In (b) light gray surfaces indicate the lower level of spring tide and dark gray surfaces the rocky areas, including oyster reefs. The intertidal zone in (c) is covered with an NDVI map retrieved from a Spot image (15/09/12, 10:54 GMT).

(a) (b) (c)
The first millimeters of sediment with biofilm (~10 ml) were scraped using a clean spoon for the three stations. Samples were kept in cooling boxes during the few hours of transportation and were stored at −20°C back in the laboratory until processing. The organic matter, including microorganisms, was isolated from the sediment following a method adapted from Blanchard et al. (1988) by Méléder et al. (2007) using Ludox® HS-40 colloidal silica (see Méléder et al., 2007 for details)  and following previous reference works (Ribeiro, 2010;Witkowski et al., 2000). In addition, samples of cleaned frustules were mounted on cover slips fixed to metallic supports and coated with platinum In parallel, a biometric analysis was done on few specimens (n > 3) of the more abundant morphospecies to estimate their lengths and widths.

| Foraminiferal assemblages in the first centimeter
In each of the three stations, three replicates were cored at one meter  (Dunn, 1964) was applied to identify which density differs from the others.  (Pawlowski, 2000). Five specimens of each species were sampled in each of the three stations. Additional superficial sediment was scraped directly with 50 ml Falcon™ tubes from sites H17, H18, and H19 and immediately placed in the cooling box for the journey back to the laboratory. There, it was stored at −20°C until DNA was extracted from 10 g of sediment from each site with the DNeasy PowerSoil Kit (Qiagen) according to the manufacturer's instructions.

| High-throughput sequencing
The primers 515f and 806r were used to amplify a ~250 bp fragment in the V4 region of 16S rDNA (Caporaso et al., 2011)

| Sanger sequencing
To better identify the taxa found in the foraminifera, longer DNA fragments from a subset of foraminiferal specimens (one per species per site) were amplified, cloned, and sequenced with the Sanger method. Extractions were amplified with three different sets of primers to amplify fragments of the SSU rDNA gene for different groups: foraminifera, prokaryotes/chloroplasts (16S), and eukaryotes (18S). For foraminifera, taxon-specific primers s14F3-J2 and s14F1-N6 (Darling et al., 2016;Pawlowski, 2000) were used with two rows of PCR following the protocol described in Darling et al. (2016).
The amplified region (~500 bp) is situated at the 3′ end of the SSU rDNA, in the 18S V9 region, and is used for foraminiferal barcoding (Pawlowski & Holzmann, 2014). 16S rDNA and 18S were also amplified from the same DNA extractions following the protocol de- In addition, the diatomaceous 18S rDNA sequences were placed with a representative selection of sequences belonging to diatoms taken from GenBank and aligned with MUSCLE (Edgar, 2004) implemented in Seaview v.4 (Gouy et al., 2010). Four subsets were prepared with the GenBank sequences most closely related to the studied sequences. Molecular phylogenetic trees were built with the PHYML program (Guindon & Gascuel, 2003) implemented in Seaview v.4, choosing the GTR (General Time Reversible) evolutionary model (Tavaré, 1986) and the approximate likelihood ratio test (aLRT) for branch support estimation (Anisimova & Gascuel, 2006).
To further analyze the Sanger sequences retrieved from the nine foraminiferal specimens, a microbiome interspecies comparison was done using the Gephi software (http://gephi.github.io/; Bastian et al., 2009) on 16S and 18S rDNA data. Gephi is a software often used in biology allowing the visualization of network (Jacomy et al., 2014;Serive et al., 2017). The DNA network analysis associated rDNA data (16S and 18S) extracted from the three studied foraminiferal species; that is, the software grouped the DNA data in communities sharing a common foraminiferal species. Circles with a high diameter represent foraminiferal species, while circles with a smaller diameter represent diatoms and other taxa found in their cytoplasm and identified with 16S and 18S rDNA. We used the ForceAtlas2 algorithm of the software to carry out the network analysis; it is a force-directed layout where nodes repulse each other while edges attract their nodes. The final network helped to interpret intuitively the data through a community-based analysis network (Jacomy et al., 2014). For clarity, 16S and 18S rDNA data mainly found in either Ammonia sp., Haynesina germanica or Elphidium oceanense were presented using colors identical to the ones of the foraminifera in which they were mainly encountered.

| Granulometry of the sampled stations
Stations H17 and H18 have a similar distribution of grain sizes, whereas station H19 is slightly coarser (Figure 2). Stations H17 and H18 contained, respectively, 90.6% and 91.2% of mud and 9.4% and 8.8% of sand, while station H19 contained 85.6% of mud and 14.4% of sand.

| NDVI and microphytobenthic assemblages in the first millimeter
The biomass map ( Figure 1c) showed NDVI values on the intertidal mudflat ranked as expected from 0.1 (no biomass) to 0.3 (maximum of MPB biomass) (for comparison, see Echappé et al., 2018;Méléder et al., 2003). NDVI mean values, calculated for H17, H18, and H19 were, respectively, 0.27 ± 0.01 (n = 18), 0.26 ± 0.009 (n = 10), and 0.15 ± 0.009 (n = 25), indicating that H19 was the station with the lowest biomass whereas both H17 and H18 were similar (n = 3; ANOVA: F = 39.14, p ≤ .001; Tukey-test: t = 1.02, p = .76 between H17 and H19, t = 10.3 and t = 11.31, both p ≤ .001 between H19 and, respectively, H18 and H17). In the three sediment samples analyzed by microscopic observations, the biomass at the surface was mainly due to MPB assemblages only composed of diatoms.   Table 2), was the most abundant morphospecies in H17 and was only found at this station. Navicula spartinetensis, a smaller pennate ( Figure 3d, Table 2) dominated H18 and P. delicatulum ( Figure 3h, Table 2) H19, respectively, although they were found in lower numbers in the other stations (Table 1). Centric diatoms were the second major group for the three stations with Thalassiosira and Odontella spp. (Figure 3i, Table 2). Cymatosira belgica was the third most abundant morphospecies in H17 and H18 and the sixth in H19 (Tables 1 and 2; Figure 4). This species is a small colonial diatom ( Figure 3g,  Figure 5). However, in H17 and H18, there were also several diatoms larger than 100 μm (e.g., P. seriata and G. wansbeckii) and even close to 400 μm in some cases (Nitzschia sigma and Gyrosigma balticum), but it was not the case in H19 with only small (e.g., P. delicatulum, C. belgica, and Cocconeis) and medium (N. spartinetensis and N. phyllepta) morphospecies (Table 2; Figure 5).

F I G U R E 2
Granulometric characterization of the sediment in H17, H18, and H19 with comparison of particle diameter (μm) against class weight (%).
The means of foraminiferal total densities in the stations H17, H18, and H19 are 841 ± 513, 1457 ± 159, and 1257 ± 490 per 50 cm 3 , respectively (Table 3). Although the standard deviation gave an indication of spatial heterogeneity at each station, the densities of foraminifera were not statistically different between the three sta-    (Table 6). The most numerous sequences belonged to diatoms (191 sequences), followed by fungi (69 sequences) and animalia (two sequences of nematodes). The main sequenced taxa of diatoms were Thalassiosira (76 sequences) and Gyrosigma (66 sequences), representing >75% of the diatom sequences.
The percentages of diatoms, fungi, and animalia varied between foraminiferal species; either the fungal or the diatomaceous sequences dominated ( Diatomaceous nuclear sequences belonging to several phylotypes of the same genera were retained for further phylogenetic analyses; four different alignments were made for these subsets and phylogenetic trees were built with these data sets (Figures 9-12).
For Entomoneis (Figure 9), three different phylotypes were recognized. For the Naviculales (Figure 10), four phylotypes were recognized for Gyrosigma, one for Navicula, and two for Pleurosigma.
Five phylotypes of Nitzschia have been identified (Figure 11). The last phylogenetic tree concerned Thalassiosira (four phylotypes) and Odontella (two phylotypes) (Figure 12). None of these phylotypes were 100% identical to sequences identified at the species level.

| Foraminiferal microbiome network analysis
To further investigate the molecular data, a microbiome interspecies comparison was done with DNA network analyses associating rDNA data (16S and 18S) extracted from their respective foraminiferal species and sequenced with the Sanger method ( Figure 13).
The community analysis with 16S rDNA data highlights microbiome differences for the three studied species (Figure 13a The 18S rDNA community analysis based on Sanger sequences also highlights species-specific microbiomes between the three foraminiferal genera with differences and similarities compared with the 16S rDNA community analysis. Ammonia sp. T6 holds diatoms, fungi and nematods. Elphidium oceanense contains fungi and diatoms, and its diatom 18S rDNA sequences belong mainly to taxa often encountered in mudflats (e.g., Nitzschia, Entomoneis, and Navicula) and ubiquitous ones such as Odontella and Thalassiosira. Haynesina germanica exclusively retains large pennate diatoms (Gyrosigma and Pleurosigma) and thus shows a lower taxonomic diversity than the 16S rDNA data.

| DISCUSS ION
In tidal mudflats, the species diversity of foraminifera is rather low compared with other environments such as the top of the continental margin (e.g., Bignot, 1985;Fontanier et al., 2002;Mojtahid et al., 2010). In the Bay of Bourgneuf, seven species have been rec-

| Densities of MPB, diatoms and foraminifera in stations H17, H18, and H19
NDVI values calculated from satellite data (Figure 1) et al., 1998;Méléder et al., 2007;Paterson & Hagerthey, 2001;Ribeiro et al., 2013). H17 and H18 are muddy stations (mud >90%) and can be viewed as similar. H19 is considered as sandy mud with less mud than the other stations (mud = 85.6%). The strong presence of Planothidium delicatulum in H19 (Figure 4) can be explained by a higher proportion of sand in this station, as this morphospecies is epipsammic. However, the granulometry of this site is not coarse enough to expect an important change in MPB communities (Méléder et al., 2007). Therefore, in the present case, MPB density changes between stations H17-H18 and H19 can probably not be attributed to grain sizes changes between the three stations. Other factors influencing these site differences could be the distance to the oysters or different currents or interactions with other organisms. For example, the positive feedback of oyster dejections on the microphytobenthos was shown by Méléder et al. (2007) and Echappé et al. (2018).
As MPB is almost only composed of diatoms, it is rather logical that the direct counts of diatoms follow a similar trend as NDVI values, with the highest density in H17 and a decrease in H18 and H19. This is also observed with sediment eDNA data where the percentage of diatoms decreases from H17 to H19 (Figure 8, Sediment 18S). In the three stations, most of the diatoms have a size around 100 μm or lower, but H19 is the only station with no diatom bigger than 100 μm ( Figure 5).
When combining diatoms and foraminifera data (Tables 1-3, 5 and 6; Figures 5 and 8), we can see that Ammonia sp. T6 and E. oceanense contain all sizes of diatoms from small to large, whereas H. germanica harbors medium to large diatoms from two genera ( Figure 14). Moreover, H. germanica and E. oceanense hold diatoms in stations H17, where diatom density is higher, and H18 (Figures 8   and 14). In station H19, with the lowest density of diatoms and an absence of large diatoms (Table 2, Figure 5), H. germanica continues to harbor diatoms with a lower percentage and two specimens of E. oceanense still have a majority of diatom DNA (Figures 8 and 14).
Station H19 has the lowest diatom density and the highest density of H. germanica, which could be explained by a top-down control, that is, when populations of organisms from lower trophic levels (diatoms here) are controlled by the organisms of higher trophic levels (forams here).

| Comparison between morphological and molecular identifications of diatoms
The comparison between morphological (Table 1), 16S (Table 5) and 18S rDNA ( were not recognized with DNA datasets (Sanger and HTS). This could either be explained by the absence of these genera in the foraminiferal microbiomes (at least in the most numerous taxa, as only 50 clones were selected for each foraminifer), the taxa selectivity of primers during amplification (primer bias) or possibly by a discrepancy between morphological and molecular taxonomies. Conversely, some of the genera identified with 16S rDNA (Asterionellopsis, Haslea, and Lithodesmium) and 18S rDNA (Ditylum and Entomoneis) were not recognized with the other datasets. As 16S rDNA has a lower phylogenetic resolution than 18S rDNA (Pillet et al., 2011), different diatom species or genera may share a common 16S rDNA sequence for their chloroplasts. Therefore, previously unsequenced diatoms from Bourgneuf could have the same sequences as Asterionellopsis, Haslea, and Lithodesmium, which could explain why these genera were not retrieved from 18S rDNA and morphological analyses. For 18S rDNA, Entomoneis is a diatomaceous genus present in Bourgneuf, but as its frustule is very fragile, it usually disappears during the processes used to prepare the material for morphological observation (see Section 2.2). This fragility could explain the absence of Entomoneis from the list of common morphospecies, whereas the genus was recognized with DNA. The absence of Ditylum in the morphospecies list could be explained by a discrepancy between morphological and molecular taxonomies or identification problems (e.g., Amato et al., 2007;Kaczmarska et al., 2007).
These results agree well with what is known on Ammonia in the literature. This genus is thought to be omnivorous, feeding on organic detritus, bacteria, microalgae, and meiofauna such as nematods (Dupuy et al., 2010;Mojtahid et al., 2011;Pascal et al., 2009;Wukovits et al., 2018). A study using a metabarcoding approach with Several ultrastructural studies have shown the total ingestion of diatom frustules by Ammonia sp. LeKieffre et al., 2017). In addition, diatom chloroplasts quickly become nonfunctional in this taxon as they are digested, demonstrating that this foraminifer is not kleptoplast (Jauffrais et al., 2016 Nomaki et al., 2014;Salonen et al., 2019).  (Bird et al., 2020;Hayward et al., 2004), which can now be distinguished morphologically (Richirt et al., 2019). The present study and the ones of Chronopoulou et al. (2019), Jauffrais et al. (2016Jauffrais et al. ( , 2018 and LeKieffre et al. (2017)

| Microbiome and possible trophic strategies of Elphidium oceanense
With 16S data, the specimens of E. oceanense from H17 had the highest percentage of chloroplastic sequences with 22%-66%, whereas specimens from H18 and H19 had more than 99% of bacterial DNA (Figure 8, 16S). It was more contrasted for 18S data, ten specimens had a majority of diatomaceous sequences (44%-92% of the total) and the five remaining specimens had a majority of fungal sequences (80%-98% of the total) (Figure 8, 18S). H19-21 had 100% of bacterial sequences and 95.2% of fungal sequences ( Figure 8). As no foraminiferal DNA could be amplified from this replicate, it may have been dead at the time of collection (see Schweizer, 2015) and its microbiome would be the reflection of the decay mechanisms happening after its death with bacteria and fungi acting as decomposers.
The same may be true for specimens H17-22, H18-25, and H19-25, which also contained a majority of fungal DNA (no data for 16S).
The densities of E. oceanense are similar in the three stations, and  Figure 7). The diatoms caught by this foraminifer have small to large sizes and a higher taxonomic diversity than in Ammonia and Haynesina (Figure 14). The microbiome network analysis showed that E. oceanense contained bacteria, diatoms, and fungi, some similar to Ammonia sp. T6 ones, and others (bacteria Herminiimonas and Delftia, diatoms Nitzschia, Entomoneis, Navicula, and Odontella) not shared by other foraminifera (Figure 13). Herminiimonas was isolated only from H19-21, and as this specimen could have been dead at the time of collection, the bacterium could be linked to decay processes.
It is very difficult to find information on E. oceanense in the literature, as this species was often mixed with other ones in the E. excavatum morphospecies. Comparisons in this group are difficult. For example, Elphidium selseyense, which was also included in the E. excavatum F I G U R E 9 Partial 18S rDNA phylogeny of Entomoneis inferred using the ML method with the GTR model and the aLRT SH-like branch support. Sequences coming from this study are indicated in bold; other sequences come from GenBank. Amphora sequences were used as out-group. 804 out of 866 sites were used and 81.4% of these sites had no polymorphism.

Entomoneis sp. Bn1
group, is a kleptoplastic species Jauffrais et al., 2018). However, E. oceanense is different morphologically, genetically, and physiologically from E. selseyense (Darling et al., 2016;Jauffrais et al., 2018). There is only one ultrastructural study with a clear identification of E. oceanense where the ingested chloroplasts are in degraded states , showing that they are eaten, but probably not used for photosynthesis. Unpublished data would be needed to investigate their roles.
F I G U R E 1 0 Partial 18S rDNA phylogeny of Gyrosigma, Pleurosigma, and Navicula (Naviculales) inferred using the ML method with the GTR model and the aLRT SH-like branch support. Sequences coming from this study are indicated in bold; other sequences come from GenBank. Nitzschia sequences were used as out-group. 842 out of 877 sites were used, and 75.4% of these sites had no polymorphism. Pleurosigma sp. Bn2

HQ121419-Nitzschia
Therefore, we hypothesize that E. oceanense is a probable heterotrophic herbivorous foraminifer feeding mainly on diatoms (Figures 8 and 13). Where diatoms are less abundant (e.g., H19), E. oceanense possibly mixes diatoms with other preys (bacteria and fungi, unless they have a different role; Figures 8 and 13).

| Microbiome and possible trophic strategies of Haynesina germanica
The replicates of H. germanica had a majority (57.9%-92.1%) of chloroplastic sequences in all stations, contrary to other foraminiferal species (Figure 8, 16S). All specimens from H17 and four specimens from H18 had more than 90% of diatom nuclear sequences, whereas in the remaining specimens, diatom sequences were either still in majority (H19-02, H19-03, and H19-04) or equal with other eukaryotes (H18-09), or fungal sequences were prominent (H19-01 and H19-10) ( Figure 8, 18S). Although diatom densities are decreasing from H17 to H19 (Table 1), H. germanica is equally present in the three stations (Table 3, Figure 7) and is relatively more abundant than Ammonia sp.
T6 in H19, the station with the lowest diatom density. As nuclear DNA is not supposed to be retained by kleptoplastic species, the 18S rDNA may reflect the recent feeding activity of H. germanica contrary to chloroplastic DNA. The discrepancy with 16S rDNA sequence identification could also be explained by the low phylogenetic resolution of 16S data, either short or long amplicon (Figure 13a vs. b) Figure 13a). Our data show a clear preference of H. germanica for diatoms among other preys (Figure 8, 18S).
The diatoms caught by this foraminifer have medium to large sizes ( Figure 14), even when these taxa/sizes are scarce or virtually absent in the environment (Tables 1 and 2, Figure 5).
Haynesina germanica has been shown to crack frustules of large diatoms and suck out the cell content while keeping the frustule outside their shell (Austin et al., 2005;Jesus et al., 2021). Moreover, it is known to be kleptoplast-bearing, which means able to steal chloroplasts from its diatom preys and use them to perform photosynthesis (Jauffrais et al., 2016Jesus et al., 2021;LeKieffre et al., 2018;Lopez, 1979 et al., 2005). This shows that foraminifera living in the same habitats may use different ways to feed on the same kind of preys (diatoms).
There are discussions about the presence/absence of diatom nuclei in kleptoplastic foraminifera (e.g., Jauffrais, LeKieffre, Pillet et al., 2011), but as most of the photosynthesis genes are in the nucleus instead of the chloroplast (Eberhard et al., 2008), the diatom nucleus may be needed in the foraminiferal cytoplasm to keep the kleptoplasts running. As H. germanica is probably continually feeding on diatoms in its natural environment, diatom nuclei could always be present in its cytoplasm. Our 18S rDNA data and the ones of

| CON CLUS ION
Our results provide new information about foraminiferal ecology with an original combination of molecular and morphological data of foraminifera and diatoms, revealing the complex interactions between these protists through symbiosis or trophic relationships. To summarize, three different trophic strategies can be deduced for the foraminiferal species investigated here. Ammonia sp. T6 is a heterotrophic omnivorous foraminifer with different trophic strategies depending on resources availability and is feeding on diatoms only when they are very abundant. Elphidium oceanense is a probable heterotrophic herbivorous foraminifer, preferably feeding on diatoms when they are abundant. Haynesina germanica is a mixotrophic herbivorous foraminifer, specialized in medium to large pennate diatom preys and performing kleptoplasty. We can conclude that Ammonia sp. T6 is probably more opportunistic than E. oceanense and H. germanica (Figures 8, 13 and 14), with a wider diet as it can even prey on nematods (Dupuy et al., 2010) and other metazoa .
This study, together with other recent ones on foraminifera Pillet et al., 2011;Prazeres et al., 2017;Schmidt et al., 2016;Tsuchiya et al., 2015), highlights the importance of molecular tools to study trophic interactions and microbiome communities of protists at the single-cell scale. In addition, this study shows the relevance of combining molecular and morphological tools for studying trophic interactions and relationships between protists and their microbial associates using single-cell analysis and morphological counting methods to assess densities. Nevertheless, as mentioned earlier, some limitations linked to the lack of data in DNA databases and to the difficulty to compare morphological and molecular data may be noticed. These limitations will require further dedicated studies to be able to tackle the holobiont of single-cell eukaryotes with a higher accuracy. Moreover, additional studies focusing more on metazoa, fungi, and bacteria in the environment and in the foraminifera, using labeled organic matter and/or investigating the ultrastructure of foraminifera, are needed to go further.

F I G U R E 1 4
Diagram summarizing the diatomaceous genera and other taxa identified in the microbiome of Ammonia sp. T6, Elphidium oceanense, and Haynesina germanica (Sanger sequencing, 16S and 18S data merged) at the three stations H17, H18, and H19 with a size estimation of the diatoms based on data from formal analysis (equal); writing -original draft (lead); writing -review and editing (lead). Thierry Jauffrais: Conceptualization (equal); formal analysis (equal); writing -original draft (equal); writing -review and editing (supporting). Constance Choquel: Conceptualization (equal); data curation (equal); formal analysis (equal); writing -original draft (supporting); writing -review and editing (supporting).

ACK N OWLED G M ENTS
Fieldwork was done during the project COSELMAR funded by the We thank Bruno Jesus (ISOMer, University of Nantes) for discussions on an early version of the manuscript and his help to produce Figures 3 and 5 and Eric Bénéteau (LPG, University of Angers) for sediment grain size measurements. We also thank four anonymous reviewers for fruitful comments.

CO N FLI C T O F I NTE R E S T
We have no conflict of interest.