Comparative Metagenomic Analysis of Fungal Communities and Genomic Analysis of Yeasts from the Chilean Altiplano and Central Mountains

1. Introduction

Fungi are ecologically versatile and physiologically flexible organisms capable of persisting across a broad range of environmental conditions, including cold, dry, nutrient-limited, and otherwise stressful habitats. They are essential components of soil ecosystems, where they contribute to organic matter decomposition, nutrient cycling, and plant–microbe interactions. Fungi inhabiting persistently cold environments or experiencing periodic freezing have evolved physiological and biochemical adaptations that enable them to grow and survive at low temperatures while coping with associated stresses, including limited nutrient availability, low water activity, oxidative stress, and freeze–thaw fluctuations [1,2,3]. Cold environments encompass extensive areas of the Earth’s surface, including polar regions, high-mountain ecosystems, and the deep sea [4,5]. Despite their ecological importance, fungal communities inhabiting many of these environments remain comparatively understudied.

Among high-altitude environments, the Tibetan Plateau has received increasing attention, with recent studies based on ITS amplicon sequencing reporting diverse fungal communities dominated primarily by Ascomycota and Basidiomycota [6,7,8,9,10,11,12,13,14]. These studies highlight the considerable fungal diversity in high-altitude ecosystems and provide evidence that fungal communities can persist under the combination of environmental restrictions for life characteristic of these regions. The Altiplano, the second-largest high plateau in the world, extends across the central Andes of South America. The Chilean Altiplano has a mean elevation of approximately 3,700 m above sea level and is characterized by low precipitation, intense solar radiation, low atmospheric pressure and humidity, and pronounced diurnal and interannual temperature variability [15,16]. The simultaneous occurrence of high elevation, aridity, intense radiation, and substantial thermal variability creates a highly demanding environment for microbial life and provides a natural laboratory for investigating microbial communities exposed to multiple environmental stresses. Although microbial research in the Altiplano has historically focused more strongly on bacteria [17,18,19,20], recent studies have begun to characterize the diversity and composition of its fungal communities. In the Bolivian Altiplano, ITS-based analyses identified Ascomycota as the dominant phylum, with Fusarium and Didymella among the most abundant genera [21]. Similarly, ITS2 sequencing of soils from the Salar de Huasco in the Chilean Altiplano revealed diverse fungal communities dominated by Ascomycota and Basidiomycota [22]. More recently, analyses of soil fungal communities from the northern Chilean Altiplano also found a predominance of Ascomycota and Basidiomycota, while unclassified families and genera were prevalent at several sites, highlighting the substantial fraction of fungal diversity that remains taxonomically unresolved [23].

Among fungi, yeasts represent a biotechnologically important microbial group, particularly those adapted to cold environments. Cold-adapted yeasts have evolved a range of physiological and biochemical mechanisms, including changes in fatty acid composition to maintain membrane function, accumulation of cryoprotective compounds, and secretion of cold-active extracellular enzymes that facilitate the utilization of available carbon sources [24,25,26,27]. Recent next-generation sequencing studies have revealed substantial genomic diversity among cold-adapted yeasts, with genome sizes ranging from approximately 11 to 40 Mb and variation in GC content, reflecting diverse genomic and evolutionary characteristics [28,29,30]. Comparative genomic analyses of Glaciozyma antarctica identified lineage-specific genes, including those encoding cold-shock domain proteins [31], while analyses of Naganishia vishniacii identified genes associated with cold adaptation, including those involved in fatty-acid desaturation and trehalose biosynthesis [32] . In Mrakia psychrophila, the presence of numerous Major Facilitator Superfamily (MFS) transporter genes has been proposed to facilitate nutrient uptake under the nutrient-limited conditions characteristic of the Qinghai–Tibet Plateau [33,34]. Similarly, genomic and transcriptomic studies of eight Antarctic and four Alpine yeasts identified numerous orthologs associated with responses to cold and other environmental stresses, including genes involved in membrane remodeling, compatible-solute biosynthesis, protein folding, and responses to oxidative and osmotic stress [29,30,35].

A largely unexplored aspect of environmental yeasts is their mitochondrial genome (mitogenome), which encodes functions essential for cellular metabolism and energy production. The mitogenome is a mosaic of prokaryotic ancestry and eukaryotic innovation, exhibiting considerable evolutionary plasticity [36]. Although a conserved core of oxidative phosphorylation genes, including cytochrome components and ATP synthase subunits, is generally retained, fungal mitogenomes show extensive variation in size, gene order, intron content, and accessory elements. This diversity has been associated with processes such as intron proliferation, the accumulation of mobile genetic elements, horizontal gene transfer, and expansion of intergenic regions [37,38]. Fungal mitogenomes are predominantly circular, although linear forms also occur, and their sizes can range from approximately 12 to 236 kb. Their organization may vary substantially even among closely related or conspecific isolates, with frequent gene rearrangements likely mediated by intrachromosomal recombination, highlighting the dynamic nature of mitochondrial genome evolution [39]. In pathogenic fungi such as Cryptococcus neoformans and Candida albicans, mitogenomic variation has been associated with metabolic regulation, virulence, and antifungal resistance [40]. Similarly, mitogenomes of Trichoderma and Fusarium exhibit lineage-specific gene rearrangements, differences in tRNA gene content, and horizontally acquired or positively selected ORFs, further illustrating the evolutionary dynamism of fungal mitochondrial genomes [38].

In this work, we investigated fungal diversity in soils from 10 localities spanning high-altitude environments in Chile using two complementary approaches. First, we characterized fungal community composition through ITS amplicon sequencing to assess patterns of diversity and taxonomic distribution across these environments. Second, we conducted comparative genomic analyses of 11 yeast isolates recovered from different sites, including three isolates putatively assigned to the same species and eight isolates putatively belonging to the same genus. These analyses aimed to characterize genomic variation and identify genomic features potentially associated with the occurrence of these yeasts across contrasting environmental conditions, including variation in both nuclear and mitochondrial genomes. Together, these approaches provide complementary community- and genome-level perspectives on fungal diversity and genomic variation in high-altitude Chilean soils.

2. Results

2.1. Taxonomic Composition of Fungal Communities Across Sites

The relative abundance of fungal taxa at the phylum, family, and genus levels was calculated for each site. The ten sites showed variation in fungal community composition, with differences becoming more pronounced at lower taxonomic levels (Figure 1; supplementary Table S1). At the phylum level, Ascomycota was predominant at most sites, except at Paso Sico 3 (PS3), where most ASVs were unclassified, and at Salar Maricunga 2 (SM2), where Basidiomycota was slightly more abundant than Ascomycota.

At the family and genus levels, a large proportion of ASVs remained unclassified, with the highest proportions observed at Paso Sico 3 (PS3), Bofedal Maricunga 3 (BM3), and Laguna del Inca (LI1). Among the five most abundant taxa at each site, Herpotrichiellaceae and Hydnodontaceae were the most frequently represented families, while Subulicystidium, Cladosporium, Neophaeococcomyces, and Sporobolomyces were the most frequently represented genera. At the individual-site level, Hydnodontaceae was particularly abundant at Salar Maricunga 2 (SM2), Aspergillaceae at Ojos del Salado 1 (OS1), and Sporormiaceae at Bofedal Maricunga (BM1). Similarly, Subulicystidium was particularly abundant at SM2, Penicillium at OS1, and Preussia at BM1.

The structural organization and co-occurrence patterns of the fungal communities were evaluated through co-occurrence network analyses conducted at the family and genus levels. The resulting networks showed a densely connected topology, with a high density of associations among taxonomic groups at both levels (Figure 2). At the family level, central hub nodes, identified by high node degree (up to 20), were predominantly represented by Hydnodontaceae, Venturiaceae, and Pyrenotrichaceae, whereas at the genus level, the main hub nodes were Neophaeococcomyces, Subulicystidium, Penicillium, and Solicoccozyma.

2.2. Fungal Diversity Metrics and Comparative Analysis

Shannon diversity (H′) varied across sampling localities at both the family and genus levels, ranging from 0.6 at PS3 to 2.1 at Cristo Redentor 5 (CR5) for families and from 0.5 at PS3 to 1.7 at OS1 for genera. These differences were associated with contrasting combinations of taxon richness (S) and Pielou’s evenness (J′) among localities (supplementary Table S2). The site with the lowest family H′ (PS3) was characterized by low richness and evenness, whereas the site with the highest H′ (CR5) showed high family richness (S = 44) and relatively high evenness (J′ = 0.6). The remaining sites exhibited contrasting combinations of these components: BM3 had relatively high richness but lower evenness and H′; Bofedal Maricunga 1 (BM1) had low richness but intermediate H′ and moderate evenness; Salar Maricunga 3 (SM3) combined relatively high richness with moderately high H′; and OS1 showed low richness but relatively high evenness and H′. Similar patterns were observed at the genus level, where OS1 exhibited the highest H′ despite relatively low richness (S = 21) and high evenness (J′ = 0.6), whereas PS3 showed the lowest evenness and H′. When sampling localities were ordered by increasing Shannon diversity, the trajectory of Pielou’s evenness broadly paralleled that of H′, whereas taxon richness showed no consistent increase and fluctuated substantially among localities (Figure 3). This pattern was particularly evident in localities such as BM1 and OS1, which exhibited relatively high H′ despite lower richness than several sites with lower H′, coinciding with comparatively higher evenness.

Sites were compared based on the presence and abundance of fungal families and genera, with pairwise correlations and hierarchical clustering summarized in Figure 4. Positive correlations were observed among sites at both taxonomic levels. At the family level, the highest correlations (r = 0.6) occurred between PS2 and PS3, SM2 and SM3, and OS1 and PS3, whereas the lowest correlations (r = 0.2) were observed for five site pairs, including BM3–CR5 and LI1–PS2. At the genus level, the highest correlation was observed between SM2 and SM3 (r = 0.6), followed by correlations of r = 0.5 among CR6, LI1, PS2, and PS3. The lowest genus-level correlations (r = 0.1) occurred between CR5–OS1, BM3–CR5, and BM1–CR5. Hierarchical clustering revealed groupings consistent with the geographic distribution of the sampling sites, particularly at the genus level. Sites from the Altiplano formed distinct clusters comprising BM1–BM3, SM2–SM3, and PS2–PS3–OS1, whereas sites from the Central Mountain formed a cluster comprising CR5–CR6–LI1.

2.3. Yeast Genomic Characterization and Comparative Analysis

Yeasts isolated from Altiplano and Central Mountain sites in this study and in our previous study [23] were preliminarily assigned based on ITS region sequencing and comparison of the resulting sequences against the NCBI nr_clust database using BLASTn. Isolates showing high-confidence species-level matches to the same species but originating from different sampling sites were selected for genome sequencing (supplementary Table S3). Draft genomes were generated for 11 yeast isolates, all showing >90% complete BUSCOs (C), with genome sizes ranging from 20 to 39 Mbp and GC contents ranging from 51% to 61% (Figure 5A). The draft genomes were subjected to additional taxonomic analyses. Single-copy BUSCO sequences were extracted from each genome and compared against the NCBI nr_clust database using DIAMOND. The resulting hits were tabulated, and the taxonomic identities of the matching organisms were counted to determine the predominant taxonomic affiliations of each isolate. For most isolates, the predominant affiliations were consistent with those obtained from the ITS-based BLASTn analysis. The exception was isolate LI1-17, for which the ITS-based analysis yielded predominantly matches to Sporobolomyces salmonicolor, whereas the BUSCO-based analysis yielded predominantly matches to Naganishia albida and Filobasidium floriforme (supplementary Table S4). The reference genomes for the taxa identified in the previous analyses were downloaded and compared with the draft genomes using average nucleotide identity (ANI). The ANI comparisons yielded the highest similarities to reference genomes of Naganishia for eight isolates: N. antarctica for LC2-15 (83%), N. diffluens for CR5-9 and LC4-24 (79%), N. kalamii for BP6-25 and LI1-17 (89%) and CR2-8 (80%), and N. vaughanmartiniae for BP6-7 (84%) and LC2-21 (93%) (supplementary Table S4). The remaining three isolates showed the highest similarities to S. salmonicolor: BP6-19 (99%), LC2-23 (100%), and NP14-27 (100%). Considering the combined results from ITS-based identification, BUSCO-derived taxonomic affiliations, and ANI comparisons, eight isolates were putatively assigned to Naganishia spp., whereas three were putatively assigned to S. salmonicolor.

As a first comparative assessment, synteny analysis was performed across all 11 draft genomes. As shown in Figure 5B, the genomes formed distinct groups, with high-tier syntenic blocks shared predominantly among isolates assigned to the same genus. This pattern was particularly evident among the S. salmonicolor isolates collected from different Altiplano localities: LC2-23 and NP14-27 shared 319 high-tier syntenic blocks, BP6-19 and LC2-23 shared 263, and BP6-19 and NP14-27 shared 250. Isolates assigned to Naganishia sp. showed variable synteny patterns, with a higher number of shared blocks between the Altiplano isolates BP6-7 (Bofedal Parinacota) and LC2-21 (Lago Chungara), which shared 352 blocks. Notably, high numbers of shared blocks were also observed between isolates from the Altiplano and Central Mountain: BP6-25 (Bofedal Parinacota) and LI1-17 (Laguna del Inca) shared 517 blocks, whereas CR5-9 (Cristo Redentor) and LC4-24 (Lago Chungara) shared 399 blocks. In contrast, CR2-8 (Cristo Redentor) and LC2-15 (Lago Chungara) shared a relatively low number of syntenic blocks (99). The three isolate pairs showing higher levels of synteny block sharing also showed their highest ANI values against the same reference genomes: BP6-25 and LI1-17 against Naganishia kalamii (GCA_022813745.1), CR5-9 and LC4-24 against Naganishia diffluens (GCA_037832395.1), and BP6-7 and LC2-21 against Naganishia vaughanmartiniae (GCA_030039245.1). These findings indicate that isolates collected from different localities and preliminarily assigned to the same genus share greater genomic similarity with the same Naganishia reference genomes. The concordance between synteny patterns and ANI similarities further supports their assignment to Naganishia and suggests that these isolates may represent genetically differentiated lineages within the genus, a possibility that could be investigated through additional taxonomic analyses.

2.4. Yeast Genomic Functional Characterization and Comparative Analysis

The number of predicted genes ranged from 6,624 in the genome of the isolate NP14-27 to 15,921 in BP6-25 (Figure 5A). Predicted genes were functionally classified into COG categories, and the number of genes assigned to each category was normalized across isolates (Z-scores) and compared across the 11 yeast isolates. Hierarchical clustering based on these profiles revealed patterns of functional similarity among the genomes (Figure 6A). The resulting clustermap separated the isolates into two main branches. The larger branch comprised six Naganishia isolates with broadly similar overall COG profiles. Within this branch, LI1-17 and CR2-8 (both from Central Mountain), LC4-24 and CR5-9 (from Altiplano and Central Mountain, respectively), and LC2-21 and LC2-15 (both from Altiplano) clustered more closely. The smaller branch comprised the three Altiplano S. salmonicolor isolates and the putative Naganishia isolate BP6-25.

To further compare functional genomic profiles among the yeast isolates, hierarchical clustering was performed based on the percentage completeness of predicted KEGG pathways and functional modules. Both analyses produced similar clustering patterns, separating the three S. salmonicolor isolates (BP6-19, LC2-23, and NP14-27) from the putative Naganishia isolates (Figure 6B). Overall, functional profiles were broadly similar among isolates; however, several KEGG pathways showed notable differences in mean completeness between the two groups (Figure 6C). GPI-anchor biosynthesis (map00563), basal transcription factors (map03022), other glycan degradation (map00511), and valine, leucine, and isoleucine degradation (map00280) showed, on average, 10–13% higher completeness in S. salmonicolor than in the putative Naganishia isolates. In contrast, steroid biosynthesis (map00100) and mitophagy—yeast (map04139) showed, on average, 10–11% higher completeness in the putative Naganishia isolates than in S. salmonicolor. Regarding module completeness, overall similarity was also observed among isolates; however, larger differences in completeness were observed for some modules between the S. salmonicolor and putative Naganishia isolates (Figure 6 D and E). The largest difference was observed for the UDP-Gal biosynthesis, Gal → UDP-Gal module (M00554), which showed 0% completeness across all putative Naganishia isolates, compared with 50–100% completeness among the S. salmonicolor isolates. In addition, the Galactose degradation, Leloir pathway, galactose → α-D-glucose-1P module (M00632) showed, on average, 42% higher completeness in the S. salmonicolor isolates than in the putative Naganishia isolates.

The predicted protein sequences from each genome were compared against the NCBI nr_clust database using DIAMOND BLASTp. The taxonomic distribution of the resulting hits broadly reflected the patterns described above. The isolates BP6-7, BP6-25, CR2-8, CR5-9, LC2-15, LC2-21, LC4-24, and LI1-17 yielded the highest numbers of hits associated with Naganishia albida (1,188–1,374), followed by hits associated with Filobasidium floriforme (811–865) (Figure 7A). In contrast, BP6-19, LC2-23, and NP14-27 showed the highest numbers of hits associated with Sporobolomyces salmonicolor (2,596–3,192). The second most frequent taxonomic association was with Cystobasidium minutum MCA 4210 in BP6-19 and with Rhodotorula toruloides in LC2-23 and NP14-27. Principal component analysis (PCA) of the protein similarity profiles revealed clear separation among the isolates along the first two principal components (Figure 7B). PC1 accounted for 56.1% of the total variance and separated the isolates into two main groups. BP6-25, LI1-17, CR2-8, LC4-24, LC2-15, CR5-9, BP6-7, and LC2-21 clustered together, whereas LC2-23, NP14-27, and BP6-19 formed a distinct group. PC2 explained an additional 10.67% of the variance and provided further differentiation within these groups. BP6-25 and LI1-17 were positioned toward positive PC2 values, whereas LC2-21 and BP6-7 clustered near the negative extreme of this axis. Among the isolates associated with S. salmonicolor, BP6-19 showed a slight displacement along PC2 relative to the closely positioned LC2-23 and NP14-27. Overall, the PCA pattern was consistent with the separation observed in the protein-level taxonomic profiles and further supported the distinction of these isolates into two broad groups.

Among the approximately 7,300 protein hits detected, we focused on those associated with stress responses and sterol homeostasis. The distribution of hits associated with different stress types was relatively homogeneous across the yeast isolates (Figure 7C). The highest numbers of hits were associated with the general stress, proteotoxic, and starvation response categories across all isolates, ranging from approximately 72 in CR5-9 to more than 91 in BP6-25. Fewer BLASTp hits were observed for the other categories. Across all 11 isolates, approximately 400 hits were associated with osmotic, salt, and desiccation responses, with individual isolate counts ranging from 31 in BP6-19 to 45 in LC2-21. Similarly, approximately 300 hits were associated with temperature and oxidative responses across all isolates, with individual counts ranging from 19 in LC2-21 to 31 in BP6-19 for temperature, and from 17 in BP6-7 to 30 in CR2-8 for oxidative responses. The xenobiotic, drug, and metal categories showed individual isolate counts ranging from 30 in BP6-7 to 46 in LI1-17. The distributions of lipid-related metabolic functions were more variable across the yeast isolates (Figure 7D). The highest numbers of BLASTp hits were associated with sterol and ergosterol biosynthesis and metabolism, with the highest counts observed in LC2-23 and NP14-27. Lipid transport, binding, and homeostasis represented the second most frequent category, with hit counts ranging from 13 in LC2-21 to 21 in both LC2-23 and NP14-27. In contrast, relatively few hits were associated with sterol esterification and storage, with counts ranging from 1 in BP6-25, BP6-7, LC2-21, and LI1-17 to 7 in both LC2-23 and NP14-27.

The potential of the studied yeasts as sources of enzymes with commercial and biotechnological applications was assessed by comparing BLASTp results against a set of 41 industrially relevant enzymes reported in the literature. Putative homologs for all 41 enzymes were identified in at least one yeast isolate (Figure 8), with enzymes associated with metabolic and biotransformation processes being particularly well represented. Esterases were among the most common, followed by glutathione S-transferases, β-glucosidases, alcohol dehydrogenases, and aldehyde dehydrogenases. Notably, putative homologs associated with antifreeze proteins, including the antifreeze glycopeptide AFGP polyprotein and antifreeze glycoprotein, were identified in nine yeast isolates, highlighting the occurrence of candidate proteins potentially involved in adaptation to low-temperature environments and of interest for biotechnological applications. In contrast, other industrially relevant enzyme families, including amylases, cellulases, lipases, and ureases, were detected at low frequencies or were absent from most of the analyzed genomes.

2.5. Mitochondrial Genomes Assembly and Functional Annotation

Mitochondrial genomes were assembled and annotated for all yeast isolates, and the resulting putative mitogenomes are summarized in Table 1. A single contig corresponding to a putative mitochondrial genome was recovered for all isolates except BP6-19, for which two contigs were retained as putative mitogenome sequences. The lengths of the assembled mitogenomes ranged from 19,670 nt in BP6-25 to 33,593 nt in CR5-9, while GC content ranged from 31.3% in BP6-25 to 38.8% in LC4-24. Assessment of genome circularity indicated that three mitochondrial sequences were likely circular: contig 2 of BP6-19, LI1-17, and CR5-9. The mitochondrial origin of these assemblies was further supported by the identification of mitochondrial genes and other characteristic mitochondrial features through functional annotation.

The specific mitochondrial genes, structural elements, and intronic features were predicted using MFannot (Figure 9A). Core protein-coding genes involved in oxidative phosphorylation, including apocytochrome b (cob), ATP synthase subunits F0-8, F0-a, and F0-c, cytochrome c oxidase subunits 1–3 (cox1–3), and NADH dehydrogenase subunits 1–6 and 4L (nad1–6 and nad4L), were detected across nearly all isolates and showed a largely conserved distribution. Ribosomal and transfer RNA genes were also widely represented, whereas genes encoding small- and large-subunit rRNAs and ribonuclease P RNA showed more variable or isolate-specific distributions. Intronic features, including group I-derived, group IA1, group IB, group IB 3′ intronic variants, group ID intron 2, and generic introns 1, 3, and 4, displayed greater variation among isolates. BP6-7, CR5-9, LC2-15, and LC2-21 exhibited richer and more diverse intronic profiles, whereas BP6-25 showed notable reductions or absences of several intron categories. MFannot-predicted exons were consistently detected across all analyzed contigs and isolates, with intermediate quantitative values.

In each putative mitogenome, open reading frames (ORFs) were predicted using the yeast mitochondrial genetic code, translated, and functionally annotated. Comparison of the predicted proteins against the NCBI nr_clust database (Figure 9B) identified ORFs with significant similarity to core oxidative phosphorylation proteins, including cytochrome b, cytochrome c oxidase subunit 1 (COX1), and several ATP synthase subunits. The number of ORFs with nr_clust hits varied among isolates, with cytochrome b being represented across multiple isolates. Other core mitochondrial components, including ribosomal protein S3 and NADH dehydrogenase subunits, were also broadly distributed among isolates. In contrast, endonucleases and other mobile genetic element-associated proteins, including homing and intron-encoded endonucleases, showed restricted and isolate-specific distributions. Overall, BP6-25 exhibited the fewest ORFs with nr_clust protein hits among the analyzed mitogenomes. Functional classification of ORFs into Clusters of Orthologous Groups (COGs) showed that most ORFs with COG assignments belonged to the energy production and conversion category across the analyzed isolates (Figure 9C). Other COG categories detected across multiple isolates, but represented by fewer ORFs, included translation, ribosomal structure and biogenesis, and secondary metabolite biosynthesis, transport, and catabolism. KEGG pathway annotation identified oxidative phosphorylation as the only pathway consistently assigned across all isolates. Synteny analysis of the putative mitogenomes revealed distinct patterns of genomic conservation among the isolates (Figure 9D). BP6-25, LC2-23, LC2-15, and CR2-8 showed higher levels of synteny with three other isolates each. The highest levels of synteny were observed between the isolate pairs CR5-9 and LC4-24, BP6-7 and LC2-21, and CR2-8 and LI1-17.

3. Discussion

3.1. Metagenomic Diversity

The predominance of Ascomycota across most sampling sites is consistent with previous studies of fungal communities in high-elevation and arid environments, including the Chilean and Bolivian Altiplano [21,23]. At a broader scale, Ascomycota is also a dominant component of global soil fungal communities, particularly through widespread and environmentally tolerant lineages [41]. However, differences become more apparent at lower taxonomic levels when comparing our results with other high-plateau regions. Studies from the Qaidam Basin, Loess Plateau, and Bolivian Altiplano have reported different dominant genera, including Inocybe, Penicillium, and Humicola [7,9,42]. Within the Chilean Altiplano itself, our results also differ from those reported for other localities, where Coniochaetaceae, Sporormiaceae, and Teratosphaeriaceae were the most abundant families [23]. In contrast, Herpotrichiellaceae and Hydnodontaceae were frequently represented in the present study, highlighting local variation in fungal community composition even within the same broader geographic region. However, the high proportion of unclassified ASVs at the family and genus levels limits taxonomic resolution and should be considered when interpreting differences in community composition among sites. In environmental fungal metabarcoding, incomplete reference database coverage and unresolved taxonomy can constrain taxonomic assignment below higher taxonomic levels, particularly for environmental lineages that remain poorly represented in reference databases [43,44]. Despite this limitation, the co-occurrence networks revealed a densely connected organization of fungal communities at both the family and genus levels, with numerous inferred associations suggesting potential relationships among taxa. Despite this limitation, the co-occurrence networks revealed a densely connected organization of fungal communities at both the family and genus levels, with numerous inferred associations among taxa. These associations may indicate potential ecological relationships among taxa [45], but may also reflect shared environmental constraints or statistical effects [46,47,48]. In the present study, Hydnodontaceae, Venturiaceae, and Pyrenotrichaceae exhibited the highest inferred network connectivity at the family level, whereas Neophaeococcomyces, Subulicystidium, Penicillium, and Solicoccozyma showed the highest at the genus level. This differs from the previous study of the Altiplano at the family level, but Penicillium appeared in high inferred network connectivity in both studies [23].

Fungal α-diversity varied among sampling sites, with differences in Shannon diversity (H′) reflecting contrasting patterns of taxon richness and abundance distribution. Overall, variation in H′ appeared to be more strongly associated with variation in evenness than with variation in richness, suggesting that the distribution of abundances among taxa was an important component of the observed diversity patterns. Similar variation in fungal α-diversity has been reported along high-elevation gradients, where changes in richness and community structure have been associated with environmental variables such as temperature, soil pH, and nutrient availability [9].

Pairwise comparisons of family and genus abundances revealed generally positive associations among sites. Hierarchical clustering broadly reflected geographic grouping, particularly at the genus level, with several Altiplano sites clustering together and Central Mountain sites forming a separate cluster. Similar geographic structuring of fungal communities has been reported in other extreme environments, although differences in community composition may reflect both geographic distance and local environmental conditions [49,50].

3.2. Yeast Genomic Analysis Properties

The combined use of ITS sequences, protein-level taxonomic affiliations, ANI, and genome-level structural and functional analyses provided complementary evidence to differentiate the 11 yeast isolates and identify groups with greater genomic similarity. Overall, these approaches supported the grouping of three S. salmonicolor-associated isolates and eight Naganishia-associated isolates. The discrepancy observed for LI1-17 between ITS- and genome-based assignments highlights the limitations of relying on a single marker for resolving closely related fungal taxa and the value of multilocus and genome-scale comparisons for improving taxonomic resolution [51,52].

The KEGG pathway and functional module profiles were broadly similar among the yeast isolates, although several functions differed in predicted completeness between the S. salmonicolor and putative Naganishia groups. S. salmonicolor showed greater predicted completeness of GPI-anchor biosynthesis, basal transcription factors, glycan degradation, and branched-chain amino acid degradation, whereas the putative Naganishia isolates showed greater completeness of steroid biosynthesis and yeast mitophagy. The largest difference among functional modules was observed for UDP-Gal biosynthesis (M00554), which was incomplete in all putative Naganishia isolates but 50–100% complete in S. salmonicolor; the Leloir pathway (M00632) also showed higher mean completeness in S. salmonicolor. Among these differences, GPI-anchor biosynthesis and galactose-related metabolism provide notable examples of pathways that may warrant further investigation. GPI-anchored proteins contribute to cell-wall organization, adhesion, and nutrient acquisition [53,54], while UDP-galactose and the Leloir pathway are involved in glycan biosynthesis and galactose metabolism, respectively [55]. These results provide bioinformatic evidence of variation in the predicted functional repertoire and point to specific cellular and metabolic processes for future experimental study.

Predicted protein-level analyses showed broadly similar distributions of stress-related functions across the yeast isolates, with general stress, proteotoxic, and starvation responses among the most frequent categories. These patterns were shared by both the S. salmonicolor-associated and putative Naganishia-associated isolates, suggesting a relatively conserved stress-related protein repertoire. In contrast, lipid-associated functions varied more among isolates, with LC2-23 and NP14-27 showing the highest numbers of homologous protein hits related to sterol and ergosterol biosynthesis and metabolism, lipid transport, homeostasis, and sterol storage. Steroid biosynthesis also showed higher mean pathway completeness in the putative Naganishia isolates. Overall, these patterns suggest that stress-related functions are broadly conserved, whereas lipid-associated functions show greater variation among isolates and between the putative taxonomic groups. These differences provide candidates for future studies examining their potential physiological relevance.

The repertoire of putative homologs of enzymes with established industrial applications identified across the yeast isolates indicates a diversity of candidate genes with potential biotechnological relevance. The detection of antifreeze-related protein homologs is also noteworthy and may reflect genetic potential relevant to persistence under low-temperature conditions. More broadly, this repertoire highlights these yeasts as potential sources of candidate genes for heterologous expression in industrially scalable hosts such as Saccharomyces cerevisiae or Pichia pastoris. Similar approaches have previously been used to express pectinase [56], β-galactosidase [57], lipase [58], and cellulase [59] genes from psychrophilic fungi.

The mitochondrial genomes of the yeast isolates contained a conserved set of genes involved in oxidative phosphorylation and mitochondrial translation, together with substantial variation in intronic features, ORF content, and genome organization. This combination of conserved gene content and variable non-core regions is consistent with patterns reported in basidiomycetous yeasts, where mitochondrial genomes retain common core genes but differ in intron content and gene arrangement [60]. Variation in introns and intron-associated endonucleases is also consistent with the known mobility of fungal mitochondrial introns and their contribution to mitogenome structural diversity [61,62]. In contrast, differences in synteny indicate variation in mitochondrial genome organization, with some isolates sharing more similar gene arrangements than others. Notably, the strongest synteny associations were not restricted to isolates from the same locality, suggesting that geographic proximity alone does not explain the observed mitochondrial similarities. These patterns provide a basis for future studies investigating whether mitochondrial genomic variation is associated with differences in environmental or physiological contexts.

4. Materials and Methods

4.1. Soil Samples Collection

Bulk soil samples were collected at a depth of 5–10 cm using sterile 50-mL tubes, avoiding plant material and roots. Samples were immediately sealed and maintained at 4 °C during transport to the laboratory, where they were stored at −80 °C until further processing. Sampling sites were distributed across the Chilean Altiplano, including the Atacama and Antofagasta regions (3,777-4,092 m a.s.l.), and the Central Mountain in the Valparaiso region (2,853–3,817 m a.s.l.) (Figure 10).

4.2. General Molecular Procedures

DNA purification, quantification, characterization, electrophoresis, PCR, and spectroscopic analyses were performed according to standardized molecular biology methodologies [63]. Procedures involving commercial kits were performed according to the respective manufacturers’ instructions.

4.3. Soil DNA Extraction, ITS Sequencing, and ASV Inference

DNA was extracted in triplicate from 0.5–1.0 g of each soil sample using the PowerSoil DNA Isolation Kit (MO BIO Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer’s instructions. DNA concentration was measured using a Qubit 4 fluorometer (Invitrogen, USA), and the amplifiability of the ITS1-5.8S-ITS2 region was verified by PCR using the ITS1 (5′-TCCGTAGGTGAACCTGCGG-3′) and ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) primers. Samples that passed these quality-control steps were submitted to Omics2View Consulting (Kiel, Germany) for library preparation, sequencing, and initial bioinformatics processing. Briefly, ITS2 libraries were prepared by PCR amplification using indexed ITS3/ITS4 fusion primers and purified with Agencourt AMPure XP beads. PCR products were heat-denatured and circularized using a splint oligonucleotide to generate single-stranded circular DNA (ssCirDNA), followed by rolling-circle amplification to produce DNA nanoballs (DNBs). Libraries were sequenced using cPAS technology [64] on an MGI DNBSEQ-G400 platform in paired-end 2 × 300 bp mode. Sequencing results are summarized in supplementary Table S5. Amplicon sequence variants (ASVs) were inferred from the raw paired-end reads using the DADA2 pipeline implemented in the R package dada2 v1.32.0 [65]. Briefly, primer sequences were removed using BBTools v39.19, and reads were quality-filtered based on expected-error thresholds of ≤2 for forward reads and ≤4 for reverse reads. Error-corrected paired-end reads were merged into contiguous sequences using a minimum overlap of 20 bp and allowing a maximum of one mismatch in the overlap region. Chimeric sequences were identified and removed using the consensus method implemented in DADA2, resulting in a set of non-chimeric ASVs. Representative ASV sequences were exported in FASTA format using the Biostrings package [66]. For each ASV and sampling site, ASVs detected in at least two of the three replicates were retained for downstream analysis.

4.4. Taxonomic Assignment

Representative ASV sequences were taxonomically classified using the IDTAXA algorithm implemented in the DECIPHER package for R [67,68]. A custom IDTAXA training set was generated from the UNITE reference database (UNITE_public_19.02.2025) using reference sequences and their associated taxonomic annotations. Reference sequences with valid taxonomic information were retained and used to train the classifier with the LearnTaxa function implemented in DECIPHER. The resulting training set was used to classify representative ASV sequences with the IdTaxa function using a classification threshold of 50, corresponding to a minimum confidence of 50% for assignment at each taxonomic level. Relative abundance of each taxon at each sampling site was calculated from the number of reads assigned to that taxon relative to the total number of reads obtained for the corresponding site.

4.5. Co-Occurrence Network Analysis

Network inference was performed using the SPIEC-EASI framework with Sparse Graphical Lasso regularization to estimate the inverse covariance (precision) matrix [69]. Edge weights and interaction signs, distinguishing positive associations from negative associations (mutual exclusion), were derived from partial correlation coefficients calculated from the estimated precision matrix. Community structures and network modules were identified using the Louvain community detection algorithm [70]. Network topological characteristics, including node degree, node strength, closeness centrality, and betweenness centrality, were calculated to characterize network structure and identify highly connected and central taxa.

4.6. Fungal Community Composition and Diversity Analysis

To evaluate fungal ecological diversity and variation in community structure among sampling localities, taxon richness (S), the Shannon diversity index (H′) [71], and Pielou’s evenness (J′) [72] were calculated. Relationships among richness, evenness, and overall diversity were evaluated using linear regression and variance decomposition analyses to assess the contribution of these components to variation in fungal diversity. Metric values were additionally standardized as Z-scores and ordered spatially by locality to examine patterns and trajectories along ecological gradients across the sampled sites.

4.7. Isolation and Identification of Yeasts

One gram of soil sample was mixed with 1 mL of sterile water and vigorously vortexed for 10 min. After coarse particulate material was settled, 200 μL of the aqueous suspension was spread onto YM 1.5%-agar plates containing 0.3% yeast extract, 0.3% malt extract, and 0.5% peptone, supplemented with 1% glucose, 100 μg/mL chloramphenicol, and 100 μg/mL ampicillin. Plates were incubated at 4, 15, and 22 °C and inspected daily. When a yeast-like colony was observed, it was transferred to a fresh plate and incubated at the same temperature as the source plate. For initial molecular identification, the ITS1-5.8S-ITS2 region was amplified by colony PCR using the ITS1 and ITS4 primers [73] in thermal cyclers (Applied Biosystems). Amplicons were separated by electrophoresis on 1.5% agarose gels, purified using the Wizard^® SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA), and sent to the DNA sequencing service of the Pontificia Universidad Católica de Chile (https://biologia.uc.cl/investigacion/plataformas/secuenciacion-y-tecnologias-omicas). The resulting sequences were analyzed and compared using BLASTn against the NCBI database. Isolates were putatively classified at the species level when the sequence showed ≥99.5% identity to the corresponding database hit.

4.8. DNA Sequencing, Genome Assembly, and Quality Assessment

Yeast isolates grown on agar medium were submitted to SeqCenter (Pittsburgh, PA, USA) for genomic DNA extraction and whole-genome sequencing. Genomic DNA was extracted from approximately 50–100 mg of biomass using the ZymoBIOMICS DNA Miniprep Kit, and DNA concentration was quantified using the Qubit 1X dsDNA Assay with a Qubit fluorometer (Thermo Fisher Scientific). Illumina sequencing libraries were prepared using the Illumina DNA Prep kit, a tagmentation- and PCR-based library preparation protocol, with custom IDT 10-bp unique dual indices (UDIs), targeting an insert size of approximately 280 bp. Libraries were sequenced on an Illumina NovaSeq X Plus platform in one or more multiplexed, shared flow-cell runs, generating 2 × 151-bp paired-end reads. Demultiplexing, quality assessment, and adapter trimming were performed using bcl-convert v4.2.4. Raw paired-end RNA-seq reads were processed using fastp v0.23.2 [74] to remove adapter sequences and low-quality bases. Minimum read-length thresholds were adjusted according to the sequencing read length (≥50 bp for 150-bp reads, ≥40 bp for 100-bp reads, and ≥30 bp for 75-bp reads). Bases with Phred quality scores below 20 were trimmed, and reads failing to meet the minimum length and quality thresholds were discarded. fastp was also used to remove poly-G and poly-X artifacts while maintaining proper paired-end relationships. For genome assembly and quality assessment, high-quality paired-end sequencing reads were assembled independently for each yeast isolate using SPAdes v4.2.0 [75]. Assemblies were generated using multiple k-mer sizes (21, 33, 55, and 77) with the automatic coverage cutoff option (--cov-cutoff auto). Assembly quality was evaluated using QUAST v5.3.0 [76,77]. Assembly statistics, including total assembly length, GC content, number of contigs, and N50, were obtained from the QUAST reports. Genome completeness was assessed using BUSCO v6.1.0 with the fungi_odb10 lineage dataset from OrthoDB [78,79,80]. The proportions of complete BUSCO genes, including single-copy and duplicated genes, as well as fragmented and missing genes, were used to assess genome completeness.

To evaluate pairwise genomic relatedness, Average Nucleotide Identity (ANI) was calculated between query and reference genomes using FastANI [81]. An ANI threshold of 95.0% was used as a reference criterion for species-level similarity, with genome pairs showing ANI ≥95.0% considered consistent with the same species, and pairs below this threshold considered to have lower genomic similarity.

4.9. Gene Prediction and Functional Annotation

For structural gene annotation, ab initio and evidence-informed gene prediction was performed using the BRAKER2 pipeline [82]. Assembled genomes were annotated using pre-existing clade-specific fungal models or, when required, by training species-specific parameter sets with GeneMark [83] and AUGUSTUS [84] in unsupervised training mode (--esmode). Annotations were generated using soft-masked genome sequences (--softmasking), producing standardized GFF3 annotation files and predicted amino acid sequences. Predicted protein sequences were functionally annotated using DIAMOND BLASTp against the NCBI non-redundant (nr) protein database, and through KEGG Orthology assignment using KAAS [85] and functional annotation using eggNOG-mapper [86].

4.10. Synteny Profiling of Genomes

Predicted protein sequences from all genomes were compared against a DIAMOND database constructed from the combined set of predicted proteins from all genomes using an all-against-all DIAMOND BLASTp approach, with an e-value cutoff of ≤1 × 10⁻¹⁰, to identify homologous anchors. For each match, genomic coordinates, e-values, bit scores, sequence identity, and alignment coverage ratios were retained. Homologous loci were evaluated under three nested stringency levels based on sequence identity, bit score, and coverage ratio: low (≥25% identity, bit score ≥40, coverage ratio 0.4–2.5), medium (≥35% identity, bit score ≥50, coverage ratio 0.6–1.8), and high (≥50% identity, bit score ≥80, coverage ratio 0.8–1.2). Redundant matches were removed by retaining the highest-bit-score hit for each genomic locus. Filtered homologous anchors were ordered by genomic coordinates and clustered into syntenic blocks using a gap-aware approach, allowing a maximum distance of 15 kb between consecutive anchors and up to 3 intervening unmatched loci. Block orientation was classified as forward, inverted, or mixed based on the relative orientation of genomic coordinates. For short blocks, orientation was determined from the predominant coordinate direction, whereas for longer blocks, Pearson’s correlation coefficient between genomic positions was used, with |r| ≥ 0.65 indicating a consistent coordinate relationship.

4.11. Assembly and Functional Annotation of Mitochondrial Genomes

Mitochondrial genomes were de novo assembled from paired-end sequencing reads using the GetOrganelle toolkit. Assemblies were generated with get_organelle_from_reads.py using the fungal mitochondrial database profile (-F fungus_mt). Assembly graphs were constructed using a multi-k-mer strategy with k-mer lengths of 21, 45, 65, 85, and 105 (-k 21, 45, 65, 85, 105), with up to 10 extension rounds per sample. Mitochondrial genome annotation was performed using the MFannot toolkit [87], with RNAWeasel [88] used to identify mitochondrial RNA features. Open reading frames (ORFs) were predicted from assembled mitochondrial contigs using the yeast mitochondrial genetic code. The resulting protein sequences were functionally characterized using DIAMOND BLASTp searches against the NCBI nr_clust protein database, and by assigning KEGG Orthology (KO) terms using KAAS [85] and eggNOG-mapper [86].

4.12. Data Processing and Statistical Analysis

Normality was assessed for each dataset using the Kolmogorov–Smirnov test, with p-values ≤ 0.05 indicating a significant deviation from normality. As the normality assumption was not met, statistical comparisons were performed using the non-parametric Kruskal–Wallis test, followed by Dunn’s post hoc test, with a significance threshold of p ≤ 0.05. Data visualization and statistical analyses were conducted using a custom Python workflow based on pandas, NumPy, scipy.stats, Matplotlib, scikit-learn, and seaborn [89,90,91,92].

5. Conclusions

Although Ascomycota predominated in soils across the ten Chilean high-elevation sites, fungal community composition varied markedly at lower taxonomic levels. Herpotrichiellaceae and Hydnodontaceae were among the most ubiquitous and relatively abundant families, with Hydnodontaceae also occupying central positions in a densely connected co-occurrence network. Fungal diversity patterns differed among sites in terms of taxon richness and evenness, with evenness showing a more consistent relationship with Shannon diversity. For the yeast isolates obtained from these localities, integrated ITS and genome-scale analyses supported the identification of three isolates as S. salmonicolor and eight as Naganishia spp. Comparative analyses of synteny, functional genomic profiles, and protein-level profiles revealed genus-associated patterns that clearly differentiated S. salmonicolor from Naganishia. Variation in putative mitochondrial genomes was primarily associated with intron content and mobile genetic elements. Additionally, the identification of potentially biotechnologically relevant enzymes and stress-associated genes highlighted candidate genomic features for future investigation. Overall, these findings reveal the ecological and genomic diversity of high-elevation fungi in Chile and provide a foundation for future studies investigating their environmental adaptation, evolutionary relationships, and biotechnological potential.