Changes in the Microbial Community Associated with the Large Benthic Foraminifera Cycloclypeus Carpenteri, Along a Depth Gradient

Large benthic foraminifera (LBF) are marine protists which harbor photosynthetic microalgae as endosymbionts that provide energy to the host (Prazeres & Renema, 2019). These algae include diatoms (Holzmann et al., 2006) that likely show strict host-specificity (Leutenegger, 1984; Pawlowski et al., 2001). As a result, LBF are light dependent and restricted to photic environments, yet some species, such as Cycloclypeus carpenteri Brady, 1881, occur in the deepest part of the photic zone where light is extremely sparse (Hohenegger et al., 2000; Renema, 2018). These deep occurrences have been hypothesized to be possible because of nutrient poor conditions in which bacteria feed on the photosynthates of the eukaryotic photosymbionts and the foraminiferal hosts feed on the bacteria, possibly in addition to food gathered from the environment (Pomar et al., 2017). However, little data are available on either the identity of the photosymbionts or the endobiotic microbial communities of deep-living LBF.

While eukaryotic endosymbionts of LBF have been studied intensely in the past (e.g., reviewed in Prazeres & Renema, 2019), fewer studies have addressed the prokaryotic community composition associated with LBF. Most studies characterized the prokaryotic microbiome of foraminifera collected at shallow depths (Bourne et al., 2013; Martin et al., 2019; Macher et al., 2021), while other studies analyzed changes in the microbiome as an effect of temperatures and nutrient concentrations (Prazeres et al., 2017a), or assessed the depth-dependent response of LBF-associated microbial communities to bleaching (Prazeres, 2018). However, few data are available on the microbial community associated with foraminifera in the mesophotic zone.

Understanding the association of LBF with microbial taxa is important, as light intensity and nutrient availability in the deep photic zone are likely to be impacted by global and local drivers of environmental change (Baker et al., 2008). This might impact light-dependent, deep-living LBF such as C. carpenteri, which are an important part of biodiversity in mesophotic coral ecosystems (MCEs; 30 to >150 m; Renema, 2019; Rowley et al., 2019). Recent research has revealed that Indo-Pacific MCEs are thermally dynamic environments yet harbor incredibly biodiverse, vibrant reefs (Rowley et al., 2019). While mobile species might be able to trace their ecological optimum by altering their depth range, sessile benthic organisms such as foraminifera are less flexible. Available data suggest that the endosymbiotic microalgae-foraminiferal host relationship is conserved (Prazeres et al., 2017b), which begs the question how LBF populations accommodate changes in environmental conditions along a depth gradient and whether their microbial community plays a role in this.

Large benthic foraminifera have a sexual macrospheric (A-form) generation and an asexual microspheric (B-form) generation, which can be most easily separated based on size. A-forms are small, and thin shelled with a larger proloculus (hence their name: macrospheres). Microspheres can grow to large diameters, yet their proloculus is smaller in diameter; they are heavily calcified and surrounded by fragile outer whorls. The macrospheric tests of C. carpenteri are rarely larger than 1.5 cm (Briguglio et al., 2015; Krüger et al., 1996), while the microspheric tests can reach a diameter of up to 10 cm, but usually do not grow beyond 6 cm (Renema, 2018).

In this study we analyzed the eukaryotic and prokaryotic microbial community associated with the diatom-bearing LBF C. carpenteri along a depth-gradient, using molecular metabarcoding. Cycloclypeus carpenteri is the largest species of the lamellar perforate family Nummulitidae, which occurs at mid-lower mesophotic depths (Hohenegger et al., 2000; Renema, 2018) in the tropical Pacific ocean from Japan to New Caledonia (Renema, 2018). The species has been shown to occur in high abundance between 80 and 100 m depth and to be restricted to the 20°C isotherm at that depth (Koba, 1978). Based on previous studies on shallow-living LBF, we hypothesize that the prokaryotic community associated with C. carpenteri is variable and changes with depth (Prazeres & Renema, 2019), while the eukaryotic symbiont-host relationship is obligatory and stable.

Field collections were conducted during July 2019 at two study sites in the Federated States of Micronesia (FSM). The first study site was a steep outer-reef slope on the northwest side of Pohnpei (6°59′30.08″N, 158°8′13.17″E). The second study site was an outer reef wall on the northwest side of Ant Atoll (6°45′21.19″N, 157°55′9.58″E; Fig. 1). Both sites are characterized by low attenuation gradients in light and wave surge. The mean diffuse attenuation coefficient was K_o 0.043±0.003 m⁻¹ and K_o 0.038±0.002 m⁻¹ at Pohnpei and Ant Atoll, respectively (Rowley et al., 2019). Ant Atoll was less turbid with a 1% subsurface irradiance (the compensation depth, or the lower limit of the euphotic zone where the rate of photosynthesis is equal to cellular respiration) at 110 m water depth, (23 μE m⁻² s⁻¹) compared with 95 m water depth (29±3 SE μE m⁻² s⁻¹) at Pohnpei (see Rowley et al., 2019). Significant thermal variability on a diurnal, seasonal, and annual basis led to mean temperature values at 50 m depth of 29.1°C (max/min: 30.6–23°C ±0.8 SD), at 90 m of 24.2°C (max/min: 31–13.2°C ±4.4 SD), and at 130 m of 19.5°C (max/min: 30–11.8°C ±3 SD; Rowley unpublished data; from August 2018–July 2019).

Live Cycloclypeus specimens were collected by hand using closed circuit rebreather diving technology (Liberty CCR, Divesoft LLC under the auspices of the British Sub-Aqua Club [BSAC], permit no: RD-0017). Collection was targeted to sample at each 10 m depth increment; at missing depths it was not possible to find specimens within the time constraints of the dive plan. The C. carpenteri individuals were placed in pre-labelled polyethylene bags along a depth gradient as follows (Table 1): Pohnpei (n = 29) (50, 70–130 m) and Ant Atoll (n = 30) (60, 90–120 m). Collected specimens were preserved in 98% ethanol and stored at −20°C prior to analysis.

Prior to DNA extraction, specimens were removed from ethanol and inspected for adhesions. No obvious adhesions were observed, but specimens were rinsed with lab grade sterile water (MilliQ). Then specimens were air dried under a sterile fume hood. The dried specimens were crushed to a fine powder with a TissueLyser II (Qiagen, Hilden, Germany). All individuals were homogenized by 1-minute of rigid shaking at 30 Hz. Microspheric specimens were crushed in a stainless-steel jar (Qiagen, Hilden, Germany) to which a 20-mm stainless steel grinding ball (Qiagen, Hilden, Germany) was added. To minimize the loss of material, macrospheric foraminifera were crushed in their original 2-mL sampling tubes, to which four 3.2 mm Stainless Steel Beads (Next Advance, Troy, United States) were added. The resulting fine powder was stored at −20°C. Genomic DNA extraction was conducted using the PowerSoil DNA Isolation Kit (Qiagen, Venlo, Netherlands), following the manufacturer’s instructions. Cell lysis was improved by implementing an extra 10-minute incubation at 70°C after addition of the proprietary aqueous lysis solution (C1). The 10-minute vortex step was executed using the TissueLyser II for 3 x 3 minutes at 30 Hz.

The quantity and quality of the extracted DNA was examined on a DropSense96 spectrophotometer (Trinean NV/SA, Gentbrugge, Belgium) corrected for turbidity.

The 18S rRNA V9 region (for eukaryotes) and the 16S rRNA V3-V4 region (for prokaryotes) were used for metabarcoding. The extracted DNA and controls were amplified using a two-step PCR protocol. The first PCR was performed using primer pairs with a Nextera overhang nucleotide sequence. The Euk1391f (5′-GTACACACCGCCCGTC-3′) and EukBr (5′-TGATCCTTCTGCAGGTTCACCTAC-3′) primer pair was used for amplification of the 18S rRNA V9 region (Earth Microbiome Project, 2019). The 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTA-3′) primer pair (see Appendix Table S1 for primers plus index sequences) was used for amplification of the 16S rRNA V3-V4 region (Klindworth et al., 2013). The PCR reaction master mix consisted of 3 µL Milli-Q H₂O, 10 µL TaqMan^TM Environmental Master Mix 2.0 (Thermo Fisher Scientific, Waltham, United States), 1 µL 10 µM forward primer, 1 µL 10 µM reverse primer and 5 µL DNA template/control.

Escherichia coli genomic DNA was used as a positive control for the 16S rRNA V3-V4 region. Milli-Q H₂O was used as a negative control throughout the entire extraction and PCR process. The reaction mixes were kept on ice prior to thermocycling to prevent degradation. Amplification was performed with C1000 Touch Thermal Cyclers (Bio-Rad, Hercules, United States) under the following conditions. For 18S rRNA: 95°C for 10 min, then 35 cycles of 95°C for 45 sec, 57°C for 60 sec, and 72°C for 90 sec, followed by a final elongation step of 5 min at 72°C (Earth Microbiome Project, 2019). For 16S rRNA: 95°C for 10 min, then 30 cycles of 95°C for 15 sec, 55°C for 30 sec, and 72°C for 40 sec, followed by a 5 min at 72°C final elongation step (adapted from Prazeres, 2018). To verify that the fragment size matched the expected ∼260 bp (for 18S rRNA) and ∼530 bp (for 16S rRNA), the PCR products were examined using gel electrophoresis (1% Agarose gel in 0.5x TBE buffer with the addition of cybersafe). PCR product purification was conducted using magnetic beads, whereby a volume of 15 µL PCR product was mixed with 13.5 µL NucleoMag® NGS-Beads (Macherey Nagel, Düren, Germany) according to the manufacturers protocol and using the VP 407AM-N 96 Pin Magnetic Bead Extractor stamp (V&P Scientific, San Diego, United States).

The purified PCR product was transferred to a 96-well 0.2-mL PCR plate and stored at −20°C pending further processing. The second PCR added unique Illumina Nextera XT indexes (Illumina, San Diego, United States) to the fragments obtained from the first PCR (Appendix Table S1). This PCR was conducted under the same conditions as the first PCR, but with 15 PCR cycles. Volumes of 20 µL 1:10 diluted PCR product were loaded on the QIAxcel platform (Qiagen, Hilden, Germany) to verify the size and the concentration of the PCR product, using the library prep process profile with the QX Alignment Marker 15 bp - 3 kb and the QX Size Marker 50 bp - 800 bp (10 µL, 1:3 diluted).

Sequencing libraries were prepared by equimolar pooling of 18S and 16S samples, respectively. The QIAgility platform (Qiagen, Hilden, Germany) was used for equimolar pooling of the 16S DNA samples. The negative controls were added with 4.8 µL for each to the sub pools, the same volume as the PCR product with the lowest concentration. All samples were amplified in duplicate (i.e., with one technical replicate). The 18S DNA samples were equimolar pooled, and 10 µl of negative control was added to the final library. Pools were cleaned with NucleoMag® NGS-Beads as described above. Concentration and amplicon size of the final pools were checked using a Bioanalyzer High sensitivity DNA chip (Agilent Technologies, Santa Clara, United States). The 18S library was sequenced on the HiSeq X Ten platform at Macrogen Europe (2 × 150 bp length; Amsterdam, Netherlands). The 16S library was sequenced on the Illumina MiSeq platform at Baseclear (2 × 250 bp length; Leiden, Netherlands).

We also amplified a 310 bp-long mitochondrial cytochrome c oxidase I (COI) fragment for a subset of 25 Cycloclypeus specimens to assess potential molecular differences between specimens sampled at different depths and between specimens from Ant Atoll and Pohnpei Island (Table S2). Specimens for COI barcoding were selected based on availability of DNA left after metabarcoding, sampling sites, and depth. Fourteen specimens from Ant Atoll and 11 specimens from Pohnpei were sequenced. The sampling depths of these specimens ranged from 90 m to 120 m (Ant) and 50 m to 130 m (Pohnpei). PCR was conducted with the primers Foraminifera_COI_fwd1/Foraminifera_COI_rev (Macher et al. 2021). The PCR reaction mix consisted of 11.3 µL Milli-Q H₂O, 2 µL Qiagen PCR buffer CL, 0.4 µL Qiagen MgCl₂, 0.8 µL Bovine Serum Albumine (BSA), 0.2 µL Qiagen Taq, 0.4 µL dNTP, 1 µL 10 µM forward primer, 1 µL 10 µM reverse primer, and 2.5 µL DNA template. Amplification was performed with C1000 Touch Thermal Cyclers under the following conditions: 96°C for 3 min, then 35 cycles of 96°C for 15 sec, 50°C for 30 sec, and 72°C for 40 sec, followed by a final elongation step of 5 min at 72°C. The PCR products were sent for bidirectional Sanger sequencing at Baseclear (Leiden, The Netherlands).

The obtained 18S and 16S reads were quality trimmed, processed, and taxonomically annotated using the Galaxy platform (Afgan et al., 2018) according to the workflow and parameters described in Appendix Table S3. Following pre-processing, read counts per sample were normalized and all Operational Taxonomic Units (OTUs) with relative abundance lower than a single read in the sample with the lowest total read count (abundance cut-off of 9.3E-6 for 18S and 6.2E-5 for 16S). This approach is commonly used in metabarcoding studies (Elbrecht & Leese, 2015; Weigand & Macher, 2018; Pereira-da-Conceicoa et al., 2019). The OTUs present in the negative controls were removed from the final OTU table if the average read count in the negative controls was <500 times smaller than the average read count in the actual samples. To further increase reliability of results, a conservative approach was followed, and only OTUs that were present in both technical replicates per sample were retained. Reads from the two replicates per specimen were pooled to build the final dataset. Following bioinformatic processing, the dataset was analysed for OTUs that were consistently part of the foraminifera’s microbiome (“core” OTUs). Core OTUs were defined as OTUs present in at least 90% of the analysed specimens. Reads identified as foraminiferal were further analysed by comparing them with the single available Cycloclypeus carpenteri reference sequence in GenBank (NCBI accession AJ879133.1 from Japan) spanning the analysed 18S fragment. The foraminiferal reads were aligned using the Clustal Omega (EMBL-EBI; Madeira et al., 2019) multiple sequence aligner ClustalW, resulting in a percent identity matrix. The generated sequence alignment was visualized using Jalview (Waterhouse et al., 2009).

The 18S and 16S datasets were separately analysed for the sampling locations Pohnpei and Ant Atoll, as macrospheric and microspheric organisms were not equally distributed between the two sampling locations (Table 1). The effect of foraminiferal generation (macrospheric or microspheric) on the relative abundance of the dominant diatom and two main foraminiferal OTUs was visualized in a Multi two-group Cumming plot created using the DABEST tool of Ho et al. (2019), using a bootstrap 95% confidence interval.

Differences in prokaryotic and eukaryotic community composition between locations and the pooled samples above and below 1% surface irradiance (50–90 vs 100–130 m at Pohnpei and 60–100 vs 110–120 m at Ant; Table 1), as well as between the two morphological types (macrospheric or microspheric) by a permutational multivariate analysis of variance (PERMANOVA) test on log-transformed read-counts derived Bray-Curtis dissimilarity matrices using the ‘adonis’ function from the R package vegan (Oksanen et al., 2019). In Pohnpei, enough specimens were available to use the same test for eight water depths along the depth gradient (Table 1).

For analysis of the dissimilarity between eukaryotic communities between locations, depths, and morphological type, the two main foraminiferal OTUs and the dominant diatom symbiont OTU were removed from the 18S rRNA dataset to focus on the remaining taxa, for which it is unknown whether they are symbionts. For the following analyses, the 16S and 18S rRNA dataset were log-transformed, as the data were highly skewed. Bray-Curtis distances among the specimens were calculated and visualized using non-metric multi-dimensional scaling (NMDS) using the MetaMDS function (iterations = 1000, dimensions = 2) from the R package vegan (Oksanen et al., 2019).

The COI sequences of C. carpenteri were analyzed in Geneious Prime (v2020.2). Read quality was manually checked, primer sequences were trimmed, forward and reverse reads were assembled, and final reads were aligned using the Geneious aligner. Genetic distances between specimens were assessed with the Geneious distance matrix.

The 59 collected C. carpenteri specimens were from two life-stage associated morphological types: macrospheric and microspheric (39 and 20, respectively). In total, 49 specimens were used for DNA extraction and analysis (Table 1). The sequenced COI barcoding fragment, analyzed for a subset of 25 specimens, showed no differences among the analyzed specimens. The 18S amplification of specimens collected from 60 and 90 m water depths at Ant Atoll, and one microspheric specimen collected at 70 m water depth from Pohnpei failed and the specimens were not further processed. The final 18S OTU-table consisted of 3262 eukaryotic OTUs. Out of these, 2263 OTUs were assigned to a reference in NCBI Genbank using blastn, but of those, 842 OTUs were identified as ‘unknown class’ (i.e., have reference sequences that could not be assigned to a taxonomic unit). The five most OTU-rich eukaryotic classes were Bacillariophyceae (179 OTUs), Dinophyceae (151 OTUs), Agaricomycetes (124 OTUs), Polychaeta (97 OTUs), and Anthozoa (91 OTUs). Fifteen OTUs were assigned to foraminifera [order Rotaliida (n = 10), order Astrorhizida (n = 4), and order Textulariida (n = 1)]. The 16S OTU-table consisted of 3607 prokaryotic OTUs. Out of these, 3399 OTUs could be assigned to a reference, but 2688 of these were assigned to ‘unknown class’. The five most common prokaryotic classes were Gammaproteobacteria (422 OTUs), Actinobacteria (131 OTUs), Planctomycetia (65 OTUs), Bacilli (37 OTUs), and Deltaproteobacteria (32 OTUs).

An OTU was perceived to be part of the core community if it was present in >90% of the analysed C. carpenteri specimens. Analysis of the eukaryotic OTUs revealed three to be part of the core community. Two of these were foraminiferal OTUs, which matched to NCBI accession AJ879131 (Heterostegina depressa) with 98.4% identity, and NCBI accession DQ440526 (Operculina ammonoides) with 97.6% identity, respectively. The third was a diatom [NCBI accession AM497721 (Fragilaria delicatissima), with 93.6% identity]. The two main foraminiferal core OTUs made up 68.5% and 29.3%, respectively, of all foraminiferal reads. The sequences of the two most abundant foraminiferal OTUs had a difference of one base pair (Appendix Table S4, Appendix Fig. S1) and were 97.6% and 96.8% similar to the available reference of Cycloclypeus carpenteri (NCBI accession AJ879133.1). The other thirteen foraminiferal OTUs were low in abundance and were only sporadically present (Appendix Figs. S2, S3). The diatom (class Bacillariophyceae) was identified as the dominant symbiont, with an average abundance of 60.1% of all eukaryotic reads per sample. We identified two core OTUs in the prokaryotic dataset, one cyanobacterium (NCBI accession KC246081, 97.8% identity) and an unidentified bacterium (NCBI accession JQ654955, 98.8% identity), which is similar (97.8% identity) to the core cyanobacteria.

The macrospheric and microspheric specimens differed in the relative abundance of foraminiferal reads (Fig. 2). The relative abundance of diatom reads did not significantly differ between microspheric and macrospheric specimens (two-sided permutation t-test, P = 0.257; Fig. 2). In contrast, the average fraction of foraminiferal reads was much lower in the macrospheres (type A) than in the microspheres (type B; 0.052 and 0.258, respectively; two-sided permutation t-test, P = 7.6·E⁻⁶; Fig. 2).

The differences in community composition between location, morphological type, and depth were analyzed with PERMANOVA. Reads from foraminifera and the dominant, putative symbiont diatom were removed to focus on the remaining eukaryotic community (22.8% of total reads remaining). The sampling location had a significant influence on the prokaryotic (PERMANOVA pseudo-F = 4.574; P<0.001) and eukaryotic communities (PERMANOVA pseudo-F = 5.404; P<0.001; see overview of statistics in Table 2). However, these differences were reduced (not significant in eukaryotes and just significant in prokaryotes (F= 1.963; p=0.004), when including the interaction between site and communities above and below the 1% irradiance cut-off.

Due to the differences in community composition between the two locations, and the sampling not being evenly distributed [only macrospheric specimens from Ant Atoll, and both micro- and macrospheric specimens from Pohnpei (Table 1)], the following analyses focus on samples from Pohnpei only.

The eukaryotic community in specimens from Pohnpei was significantly different between depths (PERMANOVA pseudo-F = 1.916, P<0.001). Furthermore, the eukaryotic community did not differ significantly between microspheric and macrospheric C. carpenteri (PERMANOVA pseudo-F = 1.399, P = 0.075). Measures of dissimilarity indicated that prokaryotic communities associated with Cycloclypeus were different for the different collection depths (PERMANOVA pseudo-F = 2.137, P<0.001), but not for the morphological types (PERMANOVA pseudo-F = 1.942, P = 0.023).

The eukaryotic and prokaryotic community compositions were visualized in a two-dimensional NMDS plot based on the Bray-Curtis dissimilarity (Fig. 3). Both NMDS plots show a division in the samples, with specimens collected at 100 m and deeper on the left side of the plot and specimens collected at 90 m and shallower depths concentrated on the right side (Figs. 3a, b). This indicates that the community composition changes between 90 and 100 m depth at Pohnpei, corresponding to the 1% subsurface irradiance mark for this site.

The eukaryotic and prokaryotic OTUs were grouped by taxonomic class, and the abundance of each class was analyzed over the depth, and thus irradiance gradient. The read abundance of the six most abundant eukaryotic classes, excluding the ubiquitous diatom symbiont, was visualized in a stacked bar plot (Fig. 4). The Polychaeta were more abundant in the shallower depths (50–90 m), and the Bacillariophyceae were more abundant at deeper depths (100–130 m). The same difference in the associated community composition was observed between microspheric and macrospheric specimens collected at the same depth (Appendix Fig. S4). Polychaeta were less abundant in macrospheric specimens (Appendix Fig. S5). The Bacillariophyceae were more abundant in macrospheres compared to microspheres (Appendix Fig. S6).

Analysis of the five most abundant prokaryotic classes (Fig. 5) shows that Actinobacteria and Bacilli were more abundant at shallower depths (50–90 m) in Pohnpei, while Gammaproteobacteria, Deltaproteobacteria, and Planctomycetia were more abundant at deeper depths (100–130 m). The same difference in community composition was observed between microspheric and macrospheric specimens collected at the same depth (Appendix Fig. S6). Deltaproteobacteria and Planctomycetia were more abundant in macrospheric specimens collected from deeper depths compared to the microspheres (Appendix Fig. S7). Gammaproteobacteria were more abundant in specimens collected at deeper depths, while Actinobacteria and Bacilli were more abundant in the C. carpenteri specimens collected at shallower depths (Fig. 5, Appendix Figs. S7A, B).

We studied the eukaryotic and prokaryotic microbial community associated with the large benthic foraminifera Cycloclypeus carpenteri at two sites in the Federated States of Micronesia and tested the hypothesis that C. carpenteri hosts a diverse eukaryotic and prokaryotic community that changes with depth. We found that C. carpenteri specimens are host to a single highly abundant diatom species, which is likely the main photosymbiont, and an associated community of both eukaryotes and prokaryotes that changed with depth.

Variation within C. carpenteri specimens in 18S, possibly due to the presence of multiple nuclei per specimen, has been reported for foraminifera and can be as high as 5% (Weber & Pawlowski, 2014; Girard et al., 2022). The molecular data of the two most abundant and core foraminiferal OTUs did not distinguish between Heterostegina and the only Cycloclypeus sequence deposited in the NCBI GenBank reference database. Since the 18S rRNA sequences of Foraminifera are often highly variable compared to other eukaryotes (Holzmann et al., 1996; Pawlowski & Lecroq, 2010) and the fragment used in this study is short and targets a wide range of eukaryotic taxa, we conclude that a longer fragment is needed for reliable molecular species identification of the foraminiferal host. A few base pairs difference in the analyzed 125 bp fragment led to relatively large genetic differences. This can explain why several foraminiferal OTUs were found. Intra-specific and intra-specimen variability (possibly due to the presence of multiple-nuclei specimen) has been reported for foraminifera and can be as high as 5% (Weber & Pawlowski, 2014). Within the most frequent foraminiferal OTUs, we did not observe a pattern of genetic differentiation with depth or between microspheric and macrospheric specimens. We therefore conclude that the 15 foraminiferal 18S OTUs observed in our study likely represent the molecular diversity in this population of C. carpenteri. This is consistent with the absence of variability between specimens from different sampling locations and between microspheric and macrospheric specimens in the analysed CO1 fragment, and we conclude that we sampled a single, genetically homogeneous population of C. carpenteri. However, we point out that future studies could use Amplicon Sequence Variants (ASVs), which can potentially detect intraspecific variability in greater detail (Callahan et al., 2017).

The C. carpenteri specimens hosted a variety of core OTUs that are consistently part of their associated microbiome. The main eukaryotic core organism was a diatom (Class Bacillariophyceae) that, on average, constituted 60.1±20% of the reads per specimen. The relatively short 18S rRNA read length and 93.6% identity to the closest reference sequence does not allow for adequate identification at lower taxonomic levels. The reported diatom symbiont of Cycloclypeus is a Thalassionema sp. (Holzmann et al., 2006), for which no reference is available for the 18S rRNA V9 region. Due to a lack of references, we refrain from drawing conclusions with regard to the identification of the diatom symbiont at a lower taxonomic level. For species identification of the main diatom symbiont we suggest analyzing a genetic marker commonly used for diatom identification, such as rbcl (Rimet et al., 2019). The abundance and diversity of polychaete OTUs was remarkable, especially in the microspheric specimens. Although visual inspection of specimens did not show the presence of polychaetes, they may have been present in the biofilm on top of the thick, heavily calcified center of the foraminiferal tests, or be the result of microbioerosion. Large, microspheric tests can provide a suitable substrate for colonization by benthic organisms such as polychaetes and bryozoa (Berning et al., 2009).

Only two prokaryotic OTUs were observed in >90% of the C. carpenteri specimens and were therefore considered core OTUs. However, based on their read count these two OTUs were rare, with only 0.8% (cyanobacteria) and 0.6% (unidentified bacteria) of all reads, respectively. The consistent presence of these OTUs might indicate that they fulfil a symbiotic function, but we cannot exclude the possibility that these OTUs were present in the environment at all sampling sites and served as food for the foraminiferal host or are part of the biofilm on the foraminiferal test.

Microspheric and macrospheric specimens differ in the relative abundance of foraminifera DNA. We observed a lower relative abundance of foraminiferal reads and a higher relative abundance of non-diatom eukaryote reads in macrospheres compared to microspheres (Fig. 2). However, the average depth at which our macrospheres were sampled was deeper (>100 m) than for the microspheres (<100 m). Our finding that the deep-living macrospheric specimens had low read counts of foraminiferal DNA and high read counts of diatom DNA might indicate that deep-living C. carpenteri need more diatom symbionts to efficiently use the little light that is available at that depth. Alternatively, it might be that the thick-walled center of microspheres contains only few photosymbiotic symbionts, which are concentrated in the thin-walled outer region of the test, resulting in a higher foraminifera-to-symbiont ratio.

The prokaryotic community associated with C. carpenteri changes along the depth gradient. We observed only two prokaryotic OTUs that were present in most specimens, suggesting that most of the prokaryotic OTUs are not obligatory associates of C. carpenteri. Heterostegina, a LBF closely related to Cycloclypeus (Holzman et al., 2006), had the most diverse prokaryotic community compared to other reef-dwelling taxa hosting photosymbionts in the Great Barrier Reef, which was suggested to be due to strong environmental control (Bourne et al., 2013). In the C. carpenteri specimens analyzed in our study, the largest shift in community composition was observed at 100 m water depth at Pohnpei, which matches the 95-m compensation depth in Pohnpei, the depth at which the rate of photosynthesis is equal to cellular respiration (Rowley et al., 2019). The prokaryotic classes Gammaproteobacteria, Planctomycetia, and Deltaproteobacteria were more abundant in the deeper collected specimens, while Actinobacteria and Bacilli were more abundant in the C. carpenteri specimens collected at shallower depths. This is in line with the previously described distribution of prokaryotic classes in seawater, where Gammaproteobacteria, Deltaproteobacteria, and Planctomycetes have been found in higher abundance in deeper waters (below 70 m; Ganesh et al., 2014; Sunagawa et al., 2015; Tseng et al., 2015). This is a further argument that the prokaryotic community found in association with C. carpenteri is mostly shaped by the environment and that the highly diverse community consisted largely of either food organisms or biofilm.

Due to the very challenging conditions during sampling of C. carpenteri at great depths, we did not collect sediment and water samples from the sampling area, which would allow comparing the community of prokaryotes and eukaryotes present in these samples to the community found in and on C. carpenteri specimens. If technically possible, future studies should aim to include such samples as environmental controls.

Furthermore, we point out that while we identified breaks in community composition at about 100 m water depth, most specimens collected at and below this depth were macrospheres, while most specimens collected from shallower depths were microspheres. It is therefore not possible to exclude that the morphological type of the C. carpenteri influences the associated community. Future studies based on a larger number of sampling sites and including analyses of environmental factors, such as light quality and nutrient availability, should address these issues by sampling more specimens of the same type from more sites and depths.

Our results show that molecular methods can help us understand diversity and ecology of Foraminifera by providing insight into the diversity of host and associated microbes. However, it also highlights that intraspecific genetic variability in Foraminifera is still poorly understood. We observed: 1) 15 different foraminiferal 18S RNA OTUs in 49 specimens, which might correspond to the intraspecific variability of Cycloclypeus carpenteri, and a single COI variant; this intraspecific variability was homogeneous throughout the entire population; 2) a single dominant diatom OTU in all C. carpenteri specimens, which likely corresponds to the described diatom symbiont of Cycloclypeus; and 3) a difference in prokaryotic community composition along the depth and, hence, irradiance gradient. The largest change observed in the prokaryotic community composition was observed at ∼1% surface irradiance.

We conclude that there is a strong relationship between the foraminiferal host and a diatom photosymbiont, but that the prokaryotic community is largely driven by environmental conditions. The role of the prokaryote community is not clear. The absence of a core community, and the congruence of the shift in prokaryote community from Actinobacteria and Bacilli in shallow samples, and Gammaproteobacteria, Deltaproteobacteria, and Planctomycetes in deeper water are indications that the prokaryotic community is mostly shaped by the species pool in the ambient environment. Our data highlight the importance of understanding the molecular diversity of the foraminiferal holobiont to address questions on the ecology and distribution of Foraminifera and their symbionts.

We sincerely thank Walter Wilbur and family at the Nihco Marine Park, Pohnpei, the Hawley family at the Pohnpei LP Gas Distributing Company, and Ant Atoll, and the Conservation Society of Pohnpei (CSP), FSM for all their assistance and logistical support. Thank you to the permitting support and guidance of S. Malakai, B. Elias, and S. Lihpai, at the Department of Resources and Development, Division of Natural Resource Management, Pohnpei State Government, FSM. Gratitude is also extended to Divesoft LLC who generously provided technical rebreather equipment support. We thank Elza Duijm, Marcel Eurlings, and Roland Butôt from the Naturalis laboratory team and the Marine Biodiversity group for support and helpful discussions. All raw reads are deposited to NCBI Sequence Read Archive, Bioproject PRJNA775074. All COI barcodes of Cycloclypeus carpenteri specimens are available in NCBI GenBank, Accession numbers OK646000–OK646024. Appendix Figures S1–S7 can be found linked to the online version of this article.

Figure S1. Sequence alignment and consensus of the 10 foraminiferal OTU sequences order Rotaliida. Generated using Jalview.

Figure S2. Stacked bar graph showing the abundance of different foraminiferal OTUs in the individual specimens collected over a depth gradient of collected C. carpenteri specimens in Ant Atoll.

Figure S3. Stacked bar graph showing the abundance of different foraminiferal OTUs in the individual specimens collected over a depth gradient of collected C. carpenteri specimens in Pohnpei.

Figure S4. Stacked bar graph showing the abundance of different eukaryotic classes in the individual specimens collected over a depth gradient of collected C. carpenteri specimens in Pohnpei, excluding the main diatom and two main foraminiferal OTUs. Type A: macrospheric and type B: microspheric specimens.

Figure S5. Stacked bar graphs of the most abundant eukaryotic classes in the analyzed foraminifera. The plot visualizes abundance of eukaryotic classes per depth category based on the sum of reads. Per depth, the average of the sum of reads per class is plotted, excluding the main diatom and two main foraminiferal OTUs. (A) macrospheric (type A) Cycloclypeus carpenteri collected at Pohnpei; (B) microspheric (type B) C. carpenteri collected at Pohnpei; and (C) macrospheric (type A) C. carpenteri collected at Ant Atoll.

Figure S6. Stacked bar graph showing the abundance of different prokaryotic classes in the individual specimens collected over a depth gradient of collected C. carpenteri specimens in Pohnpei. Type A: macrospheric and type B: microspheric specimens.

Figure S7. Stacked bar graphs of the most abundant prokaryotic classes in collected foraminifera. The plot visualizes abundance of different prokaryotic classes at different depths based on the sum of reads. Per depth, the average of the sum of reads per class is plot. (A) macrospheric (type A) Cycloclypeus carpenteri collected at Pohnpei; (B) microspheric (type B) Cycloclypeus carpenteri collected at Pohnpei and; (C) macrospheric (type A) C. carpenteri collected at Ant Atoll.

Table S1. Primers used in this study. Containing Nextera overhang sequences (underlined) and linker sequences, overlapping the indexing primers and the amplicon primers (shown in italics).

Table S2. Workflow and parameters used for sequencing data pre-processing.

Table S3. Sample location, depth and morphotype for the 25 specimens that were used for COI barcoding.

Table S4. Percentage Identity matrix of the fifteen OTUs identified as foraminifera created with the Multiple Sequence alignment of Clustal Omega (EMBL-EBI). The OTUs 0002 and 0004 are the most abundant and make up 68.5% respectively 29.3% percent of all foraminiferal reads.