Fungal Communities in Sediments Along a Depth Gradient in the Eastern Tropical Pacific

Deep waters represent the largest biome on Earth and the largest ecosystem of Costa Rica. Fungi play a fundamental role in global biogeochemical cycling in marine sediments, yet, they remain little explored. We studied fungal diversity and community composition in several marine sediments from 16 locations sampled along a bathymetric gradient (from a depth of 380 to 3,474 m) in two transects of about 1,500 km length in the Eastern Tropical Pacific (ETP) of Costa Rica. Sequence analysis of the V7-V8 region of the 18S rRNA gene obtained from sediment cores revealed the presence of 787 fungal amplicon sequence variants (ASVs). On average, we detected a richness of 75 fungal ASVs per sample. Ascomycota represented the most abundant phylum with Saccharomycetes constituting the dominant class. Three ASVs accounted for ca. 63% of all fungal sequences: the yeast Metschnikowia (49.4%), Rhizophydium (6.9%), and Cladosporium (6.7%). We distinguished a cluster composed mainly by yeasts, and a second cluster by filamentous fungi, but we were unable to detect a strong effect of depth and the overlying water temperature, salinity, dissolved oxygen (DO), and pH on the composition of fungal communities. We highlight the need to understand further the ecological role of fungi in deep-sea ecosystems.


INTRODUCTION
Fungi inhabited the oceans, including the deep-sea ecosystem, long before they conquered terrestrial environments. In addition, considering that the deep sea represents the largest biome on Earth, there is a paucity of studies on the diversity and ecology of fungi in this ecosystem compared to the rest of the ocean. Furthermore, what is known about the microbial ecology in deepsea sediments is mainly about bacteria and archaea (Edgcomb et al., 2011;Nagano and Nagahama, 2012;Dekas et al., 2016;Xu et al., 2018). Therefore, detailed knowledge of deep-sea fungi is required to understand better the overall fungal contribution to marine food webs and biogeochemical cycles at the global scale (Manohar and Raghukumar, 2013;Barone et al., 2018;Drake and Ivarsson, 2018;Grossart et al., 2019;Román et al., 2019;Hassett et al., 2020).
Fungal communities have been studied in only a small part of the great variety of habitats that exist in deep waters. Some of these habitats include sediments of hydrothermal vents, methane-cold seeps, oxygen-minimum zones, and associated with other macro-organisms (Nagahama et al., 2011;Zhang et al., 2016;Batista-García et al., 2017). In addition, some studies have shown that the subseafloor represents a vast ecosystem where micro-aerobic respiration occurs and where microbial life subsist, even hundreds of meters below the seafloor (D'Hondt, 2002;Roy et al., 2012;D'Hondt et al., 2015;Ivarsson et al., 2016a;Nagano et al., 2016).
In recent years, there has been a growing interest in studying fungal communities in deep-sea environments using culturedependent and, to an increasing extent, culture-independent methods. Abundant fungal populations have been observed in a variety of deep-sea locations such as asphalt seeps in São Paulo Plateau (Nagano et al., 2017), methane seeps in the Kuroshima Knoll (Takishita et al., 2006), hydrothermal vents in the Mid-Atlantic Ridge (Le Calvez et al., 2009;Xu et al., 2017), sediments of the Peru Trench (Edgcomb et al., 2011), the East Indian Ocean (Zhang et al., 2014), the High Arctic (Zhang et al., 2015), the Mariana Trench (Xu et al., , 2018, the Yellow Sea (Li et al., 2016), the Mediterranean Sea (Barone et al., 2018), the Yap Trench , and subsurface sediments in Suruga-Bay (Nagano et al., 2016).
In deep waters, fungi must be adapted to the total absence of light, low temperatures, and high hydrostatic pressure. Fungi in the deep-sea sediments may survive on marine snow, which consists of organic matter derived from photosynthesis that takes place in the photic layer (Bochdansky et al., 2017). In addition to performing aerobic respiration, fungi could be capable of carrying out processes such as fermentation, sulfate reduction, methanogenesis (D'Hondt, 2002;Lenhart et al., 2012), and possibly lithoautotrophy (López-García et al., 2003;Nealson et al., 2005;Ivarsson et al., 2016b). Transcriptomic analyses also confirm fungi as active members of deep-sea sediments, performing activities related to complex carbon and fatty acid metabolism (Pachiadaki et al., 2016). These metabolic processes may be more critical for fungi in deep waters since it has been observed that as depth increases, fungal populations exhibit a more multitrophic lifestyle .
Considering the enormous area to be explored for fungal diversity and function in deep-sea sediments, the existing studies are minimal and often lack an adequate spatial and temporal resolution (Grossart and Rojas-Jimenez, 2016;Grossart et al., 2019;Morales et al., 2019). Therefore, there is still a large number of geographical locations that have not yet been studied, including the Eastern Tropical Pacific (ETP). The deep-sea waters of the ETP constitute a particularly important ecosystem in Costa Rica since they represent about 90% of the whole territory (Cortés, 2016(Cortés, , 2019. The Costa Rican ETP comprises a chain of mountains and submarine volcanoes across the subduction zone of the Cocos and Caribbean tectonic plates. Here, there is a high diversity of microhabitats (Lizano, 2001;Protti et al., 2012;Rojas and Alvarado, 2012) including methane seeps (Sahling et al., 2008;Levin et al., 2012Levin et al., , 2015. Previous studies have shown high endemism and diversity of macro-and microorganisms in this region (Rusch et al., 2007;Cortés, 2008Cortés, , 2019Rojas-Jiménez, 2018). Also, the Costa Rican ETP is part of a marine corridor that extends through Isla del Coco to the Galapagos Islands in Ecuador, which represents an essential site for the conservation and regeneration of marine species throughout the ETP (Cortés, 2012).
In this work, we have explored the diversity and composition of fungal communities in deep-sea sediments of the Costa Rican ETP. Two expeditions were carried out with transects of approximately 1,500 km length each, and sediments were sampled at 16 locations at depths between 380 m and 3,474 m. We extracted DNA from subsamples of each sediment core, sequenced the 18S rRNA gene, and performed a subsequent bioinformatic analysis. This work confirms the high abundance and diversity of fungi in sediments of the ETP region. We expect that our results will support current efforts to conserve this region by providing a baseline of the high diversity of fungal species and microhabitats found in its deep-sea waters.

MATERIALS AND METHODS
We analyzed the fungal community composition in sediments along a depth gradient in the ETP of Costa Rica, across two transects of ca. 1,500 km length each (Figure 1). All samples were collected with the permission of the Ministry of Environment and Energy of Costa Rica (SINAC-CUSBSE-PI-R-032-2018; R-070-2018-OT-CONAGEBIO). The RV Atlantis surveyed the Pacific continental margin of Costa Rica from October 24th to November 5th, 2018, from the continental slope to the offshore seamounts across a subduction zone. In this region, several methane-rich seeps have been detected (Sahling et al., 2008;Levin et al., 2012Levin et al., , 2015. All sediment cores were collected by the human-occupied vehicle (HOV) Alvin equipped with mechanical, maneuverable arms. We analyzed eight sediment-cores from this expedition. The following year, the RV Falkor surveyed the seamounts extending from the mainland to the Isla del Coco National Park between January 6th and 21st, 2019. This region comprises several seamounts and natural gas seeps and provides an important corridor for highly specialized biological communities occupying the area. The sediment cores were collected employing the remotely operated vehicle (ROV) SuBastian, which is also equipped with mechanical, maneuverable arms. We analyzed another eight sediment cores from this expedition. The cores consist of an acrylic sleeve (6.7 cm diameter by 25.4 cm long) with a PVC cap and a rubber flap on the top to allow for water to escape while inserting the core while sealing as the core is removed from the sediment. The cores were kept in a "quiver" which is a PVC sleeve with a stopper at the bottom. They are sealed to the outside and are not contaminated by seawater on the way to the surface. Because the cores traveled from a higher pressure to a lower pressure, we rule out seawater intrusion. The transit time of the ROV on the longest recovery (>3,200 m depth) was approximately 2 h. Further details of the sampling sites, dates, depth, temperature, salinity, dissolved oxygen (DO), and pH are shown in Table 1.
We used the top 15 cm of the cores. Nearly one gram of the upper (1-2 cm), middle (6-7 cm), and lower (13-14 cm) parts of each core was deposited into a 1.5 ml tube, stored at −20 • C on board the vessel and at −80 • C in the laboratory. The sediment DNA was extracted with a DNA isolation kit (PowerSoil R , Qiagen, Carlsbad, CA, United States) following the manufacturer's instructions. From some subsamples, unfortunately, it was not possible to obtain enough DNA for subsequent analyzes, so in total, we retrieved 40 DNA samples (out of the 48 possible) from the 16 cores sampled in both transects. The V7 and V8 regions of the 18S rRNA gene were amplified with primers FF390/FR1 (Vainio and Hantula, 2000), using the HotStarTaq Plus Master Mix Kit (Qiagen, Carlsbad, CA, United States). The PCR conditions consisted of 95 • C for 3 min initial denaturation followed by 35 cycles at 95 • C for 45 s, 53 • C for 1 min, 72 • C for 1 min, and a final extension at 72 • C for 5 min. Multiple samples are pooled together in equal proportions based on their molecular weight and DNA concentrations. Pooled samples were purified using calibrated Ampure XP beads. The pooled and purified PCR product of nearly 350 bp were used to prepare illumina DNA library. Sequencing was performed at MR DNA 1 (Shallowater, TX, United States) on a MiSeq sequencer with v3 2 × 250 nt chemistry (Illumina, San Diego, CA, United States).
We used the DADA2 pipeline version 1.16 to process the Illumina-sequenced paired-end fastq files and to generate a table of ASVs, which are higher-resolution analogs of the traditional OTUs (Callahan et al., 2016). Briefly, we removed primers and adapters, inspected the quality profiles of the reads, filtered and trimmed sequences with a quality score <30, estimated error rates, modeled and corrected amplicon errors and inferred the sequence variants. Then, we merged the forward and reverse reads to obtain the full denoised sequences, removed chimeras, and constructed the ASV table. To assign taxonomy 1 www.mrdnalab.com to the ASVs we used the function assignTaxonomy, which is an implementation of a naive Bayesian classifier method using as input the set of sequences to be classified and a training set of reference sequences with known taxonomy, which in this case was Silva SSURef NR 132 2 (Quast et al., 2013). The assignments were verified and further curated using the BLAST tool of NCBI Genbank. All ASVs that appeared only once in the dataset were discarded. The sequence data were deposited into the NCBI Sequence Read Archive under BioProject PRJNA632873 and BioSample accessions: SAMN14924417-SAMN14924456 3 . Statistical analyses and their visualization were performed with the R statistical program (R-Core-Team, 2019) and the RStudio interface. Package Vegan v2.5-6 (Oksanen et al., 2020) was used to calculate alpha diversity estimators and, non-metric multidimensional scaling analyses (NMDS). Data tables with the ASV abundances were normalized into relative abundances and then converted into a Bray-Curtis similarity matrix. To determine if there were significant differences between the fungal community composition according to factors such as depth or transect, we used the non-parametric multivariate analysis of variance (PERMANOVA) and pairwise PERMANOVA (adonis2 function with 999 permutations). For the network analysis, we selected the 10 most abundant fungal ASVs, which corresponded to 82% of the total number of fungal sequences. We considered a valid co-occurrence event if the Spearman's correlation coefficient was >0.5 (Junker, 2008). The resulting correlation matrix was converted into an undirected matrix. We used the R package igraph v1.2.4.2 to generate the network based on the Kamada-Kawai layout algorithm (Csardi and Nepusz, 2006).
The environmental data was collected from measurements performed in the water column overlying the sediment cores and for which various instruments and sensors were used ( Table 1). Temperature and salinity data were obtained from the conductivity-temperature-depth (CTD) sensors on the HOV Alvin (CTD SeaBird SBE49) and ROV SuBastian (CTD Seabird FastCAT SBE49), which were also equipped with Niskin bottles for water sampling. There was a DO optode on the ROV SuBastian (Aanderaa 3841 O2 Optode) as well as the autonomous underwater vehicle (AUV) Sentry which was deployed over some of the sites during the 2018 Atlantis expedition. Niskin rosettes with attached CTDs were also deployed from the Atlantis and Falkor over the sites, and the Falkor CTD had a DO optode as well. DO data were compiled from a combination of these sources. DO data for the samples from the 2018 Alvin dives were derived from either the Sentry data (if available from the site) or calculated from a curve fitted from the closest CTD cast, typically, from the same site. DO data for the 2019 SuBastian push core samples was determined from SuBastian optode. The pH data were exclusively from the water samples obtained by the rosette deployed from the ship or the niskin bottles on the submersibles. Water samples were brought to room temperature and the pHT (total scale) was measured using an Orion 5 Star pH meter and glass electrode (ROSS Ultra pH/ATC Triode 8107BNUMD,

RESULTS AND DISCUSSION
We determined the presence of 787 fungal ASV in marine sediments of the Eastern Tropical Pacific of Costa Rica, obtained from 16 locations (40 subsamples) along a bathymetric gradient from 380 to 3,474 m. Fungi represented 59.72% of the 2,746,436 sequences obtained from the specific primers used for the amplification of the V7-V8 region of the 18S rRNA gene. Ascomycota was the most abundant phylum, which represented 43% of all fungal sequences and 71% of the ASVs. The second most abundant fungal group was Basidiomycota, representing nearly 3% of the sequences but 22% of the ASVs. Most of the ASVs within Basidiomycota were assigned to the order Agaricales. Chytridiomycota represented the third most abundant fungal group, with 3.5% of the sequences and 2.79% of the ASVs. Other less frequent fungal groups observed in this ecosystem were, Blastocladiomycota, LKM11, LKM15, Mucoromycota, and Zoopagomycota (Figure 2).
When analyzing the relative abundances at the class level, we detected a total of 32 classes in the deep-sea sediments, where Saccharomycetes was the most prominent in the majority of the samples. In samples where Saccharomycetes was dominant, they were typically accompanied by the presence of Chytridiomycetes. There was a second group of samples with high abundances of Eurotiomycetes, Dothideomycetes, and Agaricomycetes, but  where Saccharomycetes were practically absent (Figure 3). This observation was consistent with positive correlations within each group. For example, the correlation calculated with the Spearman method between Saccharomycetes and Chytridiomycetes was 0.80, which implies that they were present in almost all the same samples and that they presented high values of their relative abundances. On the other hand, it was also determined that the correlations between the group dominated by Saccharomyces and the other dominated by filamentous fungi were negative (Supplementary Figure 1).
These results are consistent with those obtained, at the phylum level, in deep-sea sediments from places such as the Western and Central Pacific, the Mediterranean Sea or the São Paulo Plateau, which show Ascomycota as the most abundant group, together with the presence of Basidiomycota and Chytridiomycota in lower proportions (Li et al., 2016Xu et al., 2016Xu et al., , 2019Zhang et al., 2016;Nagano et al., 2017;Barone et al., 2018;Wang et al., 2019). However, this is the first work that shows, to our knowledge, the fungal class Saccharomycetes as the most abundant, and also highly correlated with Chytridiomycetes, in deep-sea sediments.
We also observed high variability in the fungal composition within the horizons of some samples. In this sense, the homogeneity or heterogeneity of the horizons could be related to the specific conditions of the sampled site, which include the geochemical characteristics of the region, the sedimentation time, as well as the microbiological activity. A limitation of this study is the lack of geochemical data on sediment cores, since the data on the environmental variables of the overlying water column are not sufficient to explain what is happening in the vertical gradient of sediments. Some studies have shown large variations in physicochemical conditions in the profile of deep-sea sediments (Roy et al., 2012;D'Hondt et al., 2015;Román et al., 2019). Therefore, it will be necessary to continue exploring in more detail the variations in fungal communities that occur along the vertical gradient of sediment profiles.
The samples of the deep-sea environment studied, characterized by high hydrostatic pressure, low temperatures, and the absence of light, presented an average richness of 75 fungal ASVs per sample (range 13-147), while the average value of the Shannon index was 1.77 (range 0.84-2.68). As with the community analyses, there were no significant differences in the alpha diversity estimations between depths and expeditions (Kruskal-Wallis, p > 0.05). The average value of the Pielou's evenness was 0.42 (range 0.21-0.71), indicating a certain uniformity in the abundances of most of the observed ASVs (Supplementary Figure 2).
The genus Metschnikowia was the most abundant within the class Saccharomycetes and also the most abundant in the majority of the sediments analyzed. The genus Metschnikowia comprises single-celled budding yeasts known for its participation in fermentation processes and wine production, reported mainly in terrestrial environments (Kang et al., 2017;Wang et al., 2017;Pawlikowska et al., 2019). There are few references to the presence of this genus in deep waters, although its presence had been previously reported in subtropical Chinese seas, including the southern and northern Yellow Sea and the Bohai Sea (Li et al., 2016), but with lower abundances than those reported in this study. Also, we showed that this fungal genus was present in a wide depth gradient, from 380 to 3,474 m, indicating that it can be highly tolerant to gradients in temperature, DO, food supply, and the hydrostatic pressure associated with this change in depth. However, in six of the studied sediment cores Metschnikowia was almost absent, pointing more to microhabitat variability.
The most abundant genus within Chytridiomycetes was Rhizophidium which can function as parasite and decomposer (Letcher et al., 2006;Kagami et al., 2007;Frenken et al., 2017), while the most abundant genera of Eurotiomycetes and Dothideomycetes were Aspergillus and Cladosporium, respectively. Previous studies have shown that Aspergillus and Penicillium are common inhabitants of deep-sea sediments; likewise, the presence of yeasts in this ecosystem has been frequently detected, but mainly related to genera such as Pichia, Cryptococcus, Malassezia, and Rhodotorula (Takishita et al., 2006;Zhang et al., 2015;Nagano et al., 2016Nagano et al., , 2017Grossart et al., 2019). Within Agaricomycetes, the most abundant ASV had a percentage of identity of 98.73% with Armillaria, a saprophytic genus of wood that was particularly abundant in Coco Canyon (F7) at a depth of 950 m, which is also the site furthest from the coast. With the available information it is difficult to determine if this fungus, which is known to occur in terrestrial ecosystems, is active in these sediments.
Statistical analyzes, at the ASV level, did not show significant differences (PERMANOVA, p > 0.05) in the structure of fungal communities by depth, expeditions or between filtration/non-filtration areas. Neither according to variables of the overlying water columns of sediments such as pH, salinity and DO (Supplementary Table 1 and Supplementary Figure 3). For example, we showed that depth (and, consequently, hydrostatic pressure) does not have an apparent effect on the composition of communities, given the wide distribution range of species. In addition, the temperature, salinity, DO and pH values of the water column overlying the habitats of the two fungal clusters identified were similar ( Table 2). Therefore, it seems that the conditions of the deep waters are not limiting for the growth of the fungi and that other factors, likely more related to the geochemistry of the sediments, can be influencing the composition of the communities.
As an empirical observation note, samples that contained a higher proportion of mud were the ones that exhibited  a higher abundance of Saccharomycetes. In contrast, sandy samples showed higher abundances of Eurotiomycetes and Dothideomycetes, which are filamentous fungi. This observation suggests a possible relationship between fungal morphology and its ability to colonize substrates of different textures. For example, yeasts may directly depend on the type and concentrations of organic matter found in the habitat, but could also perform fermentation processes in muddy sediments (Takishita et al., 2006;Kutty and Philip, 2008;Zhang et al., 2015;Taube et al., 2018). We used network analysis to further explore possible relationships between the fungal groups that coexist in deep marine sediments of Costa Rica (Figure 4). This technique allowed us to visualize positive associations between the most abundant ASVs (representing 82% of the total sequences). We report a single co-occurrence and positive correlation between Metschnikowia and Rhizophydium. The association between these two taxa occurred regardless of the depth, location or conditions of the overlying water column. We have not found previous reports of the strong association between these two genera. We also report another group of co-occurring taxa that includes Cladosporium (Dothideomycetes), Aspergillus and Exophiala (Eurotiomycetes), and Armillaria (Agaricomycetes). The co-occurrence and high abundance of Cladosporium and Aspergillus is relatively common in deep-sea sediments Wang et al., 2019;Xu et al., 2019). However, with the available information it is difficult to determine whether the positive co-occurrence can be coincidental or can indicate a true positive interaction. Based on the results of the network analysis, as a hypothesis generating tool, we hypothesized that the fungi in both clusters can be carrying out mainly heterotrophic activities, but probably in sediments with different physicochemical conditions. The nature of the interactions within clusters should be further explored. Finally, we highlight the high prevalence of fungi in deep-sea sediments of the ETP of Costa Rica. To our knowledge, this is the first work showing a high abundance of Metschnikowia in deep-sea ecosystems. The high abundance of this type of yeasts should be further studied using cultivation-dependent methods to provide better insights into the physiology, genomic makeup, and their contributions to global biogeochemical processes. Since it was difficult to distinguish the association of specific environmental variables with variations in the composition of fungal communities, particularly in the two clusters identified, further research will be necessary to determine how fungal communities in deep-sea waters are structured as well as to determine their ecological role in the largest biome on the planet.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The sequence data were deposited into the NCBI Sequence Read Archive under BioProject PRJNA632873 and BioSample accessions: SAMN14924417-SAMN14924456 (https: //www.ncbi.nlm.nih.gov/bioproject/PRJNA632873).

AUTHOR CONTRIBUTIONS
KR-J, H-PG, EC, and JC designed the study and performed the analysis. JC and EC collected the samples. KR-J wrote the manuscript. All authors helped to revise the manuscript.

FUNDING
This project was funded by Universidad de Costa Rica, DFG project GR1540/33-1 given to H-PG, and NSF OCE 1635219 awarded to EC.

ACKNOWLEDGMENTS
We would like to thank all of the officers and crew of the RV Atlantis and RV Falkor for their assistance, and the pilots of the HOV Alvin and ROV SuBastian for their efforts in the collection of samples. We thank the Ministry of Environment and Energy of Costa Rica, Instituto Costarricense de Pesca y Acuicultura (permit INCOPESCA-CPI-003-12-2018) and Comisión Nacional para la Gestión de la Biodiversidad (permit R-070-2018-OT-CONAGEBIO). Sample processing was carried out with assistance from Odalisca Breedy and environmental data were compiled by Steve Auscavitch, Jay Lunden, and April Stabbins at Temple University. We also thank Lisa Levin and Greg Rouse at Scripps Oceanography, UC San Diego, La Jolla, CA, and Victoria Orphan at the California Institute of Technology, Pasadena, CA for their support during the development of this project. We acknowledge the support of the Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of University of Potsdam. We thank the reviewers for their critical comments on the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020. 575207/full#supplementary-material Supplementary Table 1 | Statistical analysis of the fungal community composition related to different variables. The PERMANOVA tests were performed using function adonis2 and implemented in Vegan package. Data were normalized by converting the ASV counts into relative abundances. Binning of continuous variables Depth, Temperature, Dissolved Oxygen and pH was performed with package Hmisc.