Samples were collected at the active hydrothermal vents located on the Hatiba Mons volcano, in the Red Sea Rift. The active Hatiba hydrothermal fields are situated at ≈1,000 m water depth and are characterized by widespread diffuse venting of low-temperature fluids (up to 40°C). Forty-five individual sites are mainly built of iron-oxyhydroxide mounds and chimneys, surrounded by warm Red Sea bottom water with a temperature of 21.7°C and high salinity (41 psu), flourishing microbial mats, and the absence of vent-specific macrofauna. The Hatiba vent fields are considered the largest active low-T vent field area observed so far (van der Zwan et al., 2023).

A total of 38 subsamples from crusts, sediments, and microbial mats from five hydrothermal sites within Hatiba Mons (Figure 1) were collected aboard the KRSE Aegaeo RV cruise in May 2022 (Table 1). Detailed information of the sampling sites in the Hatiba Mons vent fields is available in van der Zwan et al. (2023). At the Field 8 hydrothermal vent site, three hydrothermal crust samples (KRSE4-2 ROV HYD) were collected using a remotely operated vehicle (ROV). Each crust was subsampled into three non-homogenized parts, resulting in nine subsamples. Push cores (4 cm in diameter and 60 cm in length) were collected from the Kabir Field. Each core was divided into two depth intervals for subsampling: a microbial mat layer (0–5 cm, flocculant to loosely consolidated) and an underlying precipitate layer (5–10 cm, more consolidated). Each interval was subsampled into three non-homogenized parts, resulting in six subsamples per core and 18 across the three cores. A hydrothermal chimney crust sample (KRSE4-4 ROV HYD X) was collected and subsampled into three parts. A microbial mat sample (KRSE4-6 ROV NET) from the Farwah Safaa Ridges was retrieved using an ROV sampler and divided into three non-homogenized zones; replicate aliquots from these zones yielded eight subsamples in total.

Subsampling and transfers were conducted on board the research vessel as rapidly as possible in a designated clean workspace, adhering to aseptic technique throughout the process. The workspace consisted of a dedicated bench area that was cleaned with 70% ethanol before and during handling, with sterile tools and consumables used throughout. During all sampling procedures, personnel wore gloves and face masks, and work surfaces and gloves were routinely disinfected with ethanol prior to and during processing. All samples were aseptically transferred into sterile Whirl-Pak bags and sterile Falcon tubes by using sterile spatula. Whenever possible, subsampling focused on the interior portions of mats, cores, and crusts, with the exteriors avoided to minimize potential contamination. Herein, field blanks during subsampled were not taken. Following collection, all subsamples were promptly sealed and immediately frozen at −20°C onboard, then subsequently transferred to −80°C in the laboratory for long-term storage.

The DNA was extracted from each of the 38 samples (~10 g per sample) using a DNeasy PowerMax Soil Kit (Qiagen, Germany), following the manufacturer’s protocol. DNA concentration and integrity were checked using a Qubit dsDNA HS Assay and Qubit fluorometer 4.0 (ThermoFisher Scientific, Waltham, USA) and agarose gel (1%) electrophoresis containing Sybr Safe (ThermoFisher Scientific, Waltham, USA), respectively.

The extracted DNA was subjected to internal transcribed spacer (ITS) rRNA gene sequencing at the MR DNA (Molecular Research LP, Shallowater, USA). The DNA was amplified using a commonly used ITS primer set for fungal studies, ITS1 (5’-CTTGGTCATTTAGAGGAAGTAA-3’) and ITS2 (5’-GCTGCGTTCTTCATCGATGC-3’) (White et al., 1990). Amplification was performed using the HotStarTaq Plus Master Mix Kit (Qiagen, USA), as follows: initial denaturation at 95°C for 5 min, followed by 30 cycles at 95°C for 30 s, 53°C for 40 s, and 72°C for 1 min, and a final extension step at 72°C for 10 min. PCR products were then run in 2% agarose gel, and subsequently used to prepare the barcoded library. Barcoded libraries were quantified using a bioanalyzer and pooled at an equimolar ratio based on their molecular weight and ITS concentrations and later purified using AMPure XP beads. For the ITS sequencing (2 × 300 bp), MiSeq reagent kit v3 (Illumina, USA) was employed, in accordance with the manufacturer’s instructions.

Extraction controls were included throughout the DNA extraction process and gene amplification. These controls were assessed using the Qubit dsDNA HS Assay and Qubit fluorometer 4.0 (ThermoFisher Scientific, Waltham, USA), as well as agarose gel (1%) electrophoresis containing Sybr Safe (ThermoFisher Scientific, Waltham, USA). All extraction and PCR controls consistently showed negative results, with no detectable DNA or PCR product bands on the gel and concentrations below the detection limit. We emphasize here that the blank controls from DNA extraction, which were negative by gel electrophoresis and high-sensitivity assays for quantification, were nevertheless sequenced and sequences were deposited alongside the dataset (Supplementary Table 1). However, reads derived from blank were not included in further data processing steps and were not extracted from the sequences of the biological samples.

The ITS raw reads were processed using QIIME2 (v2021.11) (QIIME 2 Development Team, 2021). Raw ITS sequencing reads were initially inspected for quality using FastQC (Andrews, 2010) and the QIIME 2 demux summarize visualization. Reads were then processed and denoised using the DADA2 plugin in QIIME 2 (version 2021.11; Bolyen et al., 2019). As part of the DADA2 workflow, forward and reverse reads were truncated to 220 bases, and reads containing more than two expected errors were discarded. Reads were further truncated at the first base with a Phred quality score below 2. Only paired reads with a minimum overlap of 12 bp were merged, and chimeric sequences were identified and removed using the consensus method (Supplementary Figure 1). This process generated a feature table of Amplicon Sequence Variants (ASVs) and their representative sequences. Taxonomic classification of ASVs was performed using a pre-trained Naive Bayes classifier implemented in QIIME 2, with the UNITE database for fungi (version 7, clustered at 97% similarity, released December 2017; Kõljalg et al., 2013) as the reference database for taxonomic assignment. QIIME2 output files, including the feature table, taxonomy assignments, phylogenetic tree, and sample metadata, were imported in R (v4.4.1) using ‘qiime2R’ (v0.99.6) and ‘phyloseq’ (v1.50.0) packages.

Before all downstream analyses, non-fungal taxa were excluded by subsetting for sequences assigned to the Kingdom Fungi. To evaluate sequencing depth and sampling completeness, rarefaction curves were generated using the rarecurve()function in the vegan R package, based on the ASV abundance table. To visualize the relative abundances of fungal communities at different taxonomic levels, the data were aggregated at the phylum, class, order, genus, and species levels. The relative abundances were calculated by normalizing the read counts within each sample. The most abundant taxa were plotted for each taxonomic level, while low-abundance taxa were grouped under ‘Other’ to simplify interpretation. Stacked bar plots were generated to display the relative abundances across sample groups, using ‘ggplot2’ v3.5. 1. To generate species heatmaps across sample groups, the ASV table was processed using the ‘phyloseq’ v1.50.0 package. The table was collapsed to the species level with the tax_glom() function, and sample counts were normalized to relative abundances within each sample. For visualization, the dataset was aggregated by species and sample groups, summing relative abundances across samples in each group. The top 15 most abundant fungal species, based on cumulative relative abundance, were retained to improve readability. The resulting matrix, with species as rows and sample groups as columns, was visualized using the ‘pheatmap’ v1.0.12 package. Row-wise Z-score scaling highlighted relative differences in species composition. A continuous white-to-dark red color gradient represented low to high abundance, and hierarchical clustering was applied to both rows and columns to reveal distribution patterns.

Alpha-diversity metrics, including Observed richness, Shannon and Simpson diversity index, were calculated using the estimate_richness() function from the ‘phyloseq’ package. The calculated diversity indices were merged with the sample metadata for downstream analysis. Beta-diversity was assessed using the Bray-Curtis dissimilarity index to quantify compositional differences among fungal communities. Principal Coordinates Analysis (PCoA) was conducted to visualize variation in community structure. To statistically evaluate differences in beta-diversity, Permutational Multivariate Analysis of Variance (PERMANOVA) was conducted using the adonis2 function in vegan, with 999 permutations, and the coefficient of determination (R²) was reported to indicate the proportion of variance explained by groupings. To confirm homogeneity of variance, Permutational Analysis of Multivariate Dispersions (PERMDISP) was performed using the betadisper function. If PERMDISP was non-significant (p > 0.05), the assumption of equal variance was met, validating PERMANOVA results. A PCoA plot was generated using ‘ggplot2’ v3.5.1, with sample points colored by group, and 95% confidence ellipses included only if PERMANOVA results were significant.

To characterize the ecological functions of the identified fungal taxa, we utilized the FUNGuild tool v1.1 (Nguyen et al., 2016; available at https://github.com/UMNFuN/FUNGuild). Taxonomic classifications of fungal ASVs were first obtained by exporting taxonomy assignments from QIIME2 (Bolyen et al., 2019). Afterward, the FUNGuild was used to assign ecological functions to the ASV based on its taxonomic information, providing details such as trophic modes, specific functional guilds, and associated confidence rankings. In parallel, the fungal abundance table containing ASV counts across all samples was also exported from QIIME2. Finally, the abundance data were merged with the FUNGuild output, to better understand the fungal communities’ taxonomic composition and functional roles. For visualization of the patterns in fungal functional diversity, two metrics were employed to describe guild composition, in which ASV richness, representing the proportion of unique ASVs assigned to each guild, and sequence richness, representing the proportion of total sequence abundance, were attributed to each guild. These metrics were calculated separately for each sample type and hydrothermal vent site. The results were summarized as percentage values and visualized using ‘ggplot2’ v3.5. 1. It is important to acknowledge that precise community-wide conclusions for ecological function remain challenging due to several factors, such as multiple trophic strategies of the fungal taxa, guild data are still limited for numerous fungal groups (Nilsson et al., 2019), and the ecological roles of marine fungal taxa are underexplored (Velez et al., 2022). Nevertheless, FUNGuild can provide a valuable foundation toward advancing our understanding of deep-sea and hydrothermal vent fungal communities.

Raw sequencing yielded a total of 8,271,048 reads (4,135,524 paired-end read pairs) across all samples. After quality control, 228,043 reads were retained from hydrothermal sediment, microbial mats and chimney crust samples collected in low-temperature hydrothermal vent fields in the Hatiba Mons, with 104,274 for microbial mat samples, 58,547 for crust samples and 65,222 for hydrothermal sediments (Supplementary Table 2). A total of 34,940 reads corresponded to KRSE4-2 ROV HYD, 33,137 for KRSE4-3 ROV SED 1, 48,032 for KRSE4-3 ROV SED 2, 46,991 for KRSE4-3 ROV SED 4, 23,607 for KRSE4-4 ROV HYD X, and 41,336 for KRSE4-6 ROV NET.

The retained sequences were classified into 550 ASVs out of which 534 represented fungal taxa (214,707 reads). All of the rarefaction and extrapolation curves of the fungal assemblages in the samples from the different sites reached the asymptote (Supplementary Figure 2), indicating that the data generated from the 38 samples provided a good description of the fungal diversity. However, it is important to note that technical factors, such as DNA extraction method, primer specificity, or PCR bias, may have influenced the recovery and detection of certain fungal groups (Tedersoo et al., 2015; Portela et al., 2025).

Out of 228,043 reads, a total of 21,4707 were taxonomically assigned to Fungi. In all the studied sites, the fungal community was dominated by the Basidiomycota, followed by Ascomycota. The crust sample from Field 8 exhibited more diversity at the phylum, also harboring Mucuromycota and Zoopagomycota, while the crust from Farwah Safraa Ridges presented only the most abundant fungal phyla and higher abundance of unidentified ASVs. In the crust of Field 8 and sediment samples of Kabir Field, Mortierellomycota was also observed (Figure 2A). At the order level, dominance patterns varied across the hydrothermal fields. For instance, sequences of Poliporales were vastly observed in the microbial crust and mat sample of Farwah Safraa Ridges with ASVs showing 96.0–98.3% identity to Resinoporia piceata, a known terrestrial member of Polyporales, while Filobasidiales and Sacharomycetales were highly represented in the hydrothermal sediments of Kabir Field. High relative sequence abundance of Agaricales were found in the crust samples (Field 8 and Farwah Safraa Ridges). Higher sequence abundance of the other observed orders was found in the Field 8 and Kabir Field (Figure 2B). Farwah Safraa Ridges samples (crust and microbial mat) showed higher sequence abundance of the class Agaricomycetes, followed by Saccharomycetes, while for the samples from Field 8 and Kabir Field, the abundance was more equally distributed between the different classes (Figure 2C). Prevalence of Lipomyces was higher in the samples of Field 8, Farwah Safraa Ridges, and one sample group of Kabir Field, and for the other samples of Kabir, Filobasidium were in higher sequence abundance (Figure 2D). The genus Coprinus was only observed in the crust samples, especially in the Field 8, with ASVs showing 96% identity to known Coprinus sequences, while Kluyveromyces was found only in one sample of the hydrothermal sediment of Kabir Field, and Rigidoporus was uniquely found in Farwah Safraa Ridges crust (Figure 2E). Both Coprinus and Rigidoporus are examples of taxa known from the terrestrial environment and while these, as well as other close relatives of terrestrial taxa may be able to survive marine conditions, some may also persist in the deep sea as spores. Overall, at the species level, nine fungal species were assigned to the ASVs (Figure 2F), whereas each species was found in higher abundance in different samples. For instance, Coprinus cordisporus was abundant in the crusts of Field 8, Filobasidium wieringae, Trichoderma longibrachiatum, Debaryomyces udenni, Malassezia globosa, and Neoascochyta paspali were dominant species in Kabir Field, and Rigidoporus ulmarius were prevalent in the microbial mat of Farwah Safraa Ridges.

Measures of fungal alpha-diversity (number of observed ASV, Shannon and Simpson indexes) were significantly higher in the crust of the Field 8 across the other Red Sea hydrothermal vent type samples and locations. Microbial mat from Farwah Safraa Ridges and samples from Kabir Field showed similar richness and evenness, without significant difference between their fungal communities (Figure 3A). The same pattern was observed when analyzing the fungal beta-diversity explored by Bray–Curtis and weighted Unifrac distances. The PcoA analysis revealed a clear distinction between the crust and the rest of the samples and locations that formed a single cluster (Figure 3B).

Prediction of the fungal communities’ functions were performed using FUNGuild, and characterized in trophic mode in the sample’s location, and also across the sample type. In all of the samples’ locations and types, the dominance was for unassigned mode, whereas around 50% of the sequences were not assigned to any function. For the assigned reads, it was observed that putative saprotrophs were dominant in all the sites and sample types, except for microbial mat of Kabir Field (sample KRSE4-3 ROV SED 2), whereas it was observed putative pathotroph-saprotroph-symbiotroph was highest enriched with relative to the other trophic modes (Figure 4A). Crust samples, in comparison with hydrothermal sediment and microbial mat, presented the highest abundance for putative symbiotrophs (Figure 4B).