An inter-species hybrid in the emerging fungal pathogen Rhodosporidiobolus fluvialisDownload PDF Download PDF ArticleOpen accessPublished: 05 November 2026Álvaro Redondo-Río1,2 &Toni Gabaldón1,2,3,4 npj Fungal Science volume 1, Article number: 8 (2026) Cite this articleView saved researchAbstractHybridisation between divergent lineages can generate novel trait combinations and has repeatedly been linked to the emergence of opportunistic human pathogens. We re-analysed publicly available genome and read data from two reported clinical isolates of Rhodosporidiobolus fluvialis, a recently described fungal pathogen with high potential to gain virulence and pan-fungicidal resistance. Our whole-genome comparisons showed that one diploid isolate comprised two highly divergent subgenomes with little evidence for recombination and limited sequence identity to the second, haploid isolate. The high divergence between subgenomes and phylogenetic analyses of ITS and protein-coding genes support the first isolate as an allodiploid interspecies hybrid. One parental subgenome is closely related to R. fluvialis, whereas the second seems to derive from an unsampled Rhodosporidiobolus species. These findings highlight the relevance of hybridisation in emerging fungal pathogens and support haplotype-aware genomic analyses for accurate pathogen surveillance and classification.IntroductionThe recurrent emergence of new microbial pathogens poses a global health threat1. This brings about the need to perform constant monitoring and precise identification, as well as to investigate the underlying mechanisms of the emergence of opportunistic pathogenesis. Hybridisation between or within species boundaries has been proposed as a possible route for the emergence of novel human fungal pathogens2,3, where specific combinations of traits in the hybrid lineage may allow adaptation to novel niches. In some cases, such novel niches may include the human host or human-related environments, facilitating the opportunistic infection of immunocompromised patients. Consistent with this idea, hybrid strains have been reported for clinical isolates of a growing number of emerging opportunistic fungal pathogens, including the worldwide distributed Candida orthopsilosis4, Candida metapsilosis5, Candida inconspicua6, and Cryptococcus gattii7,8, but also in cases where the pathogen is rare and more geographically restricted, including single cases of clinical isolations such as Aspergillus latus9, Candida theae10, and the recently discovered Pichia alaskaensis11. These examples suggest that hybridisation may be a common route for the emergence of new fungal pathogens, particularly among rarely isolated species. Notably, genomic evidence indicates that the common yeast pathogen Candida albicans may itself derive from an ancestral hybrid12, raising the possibility that emerging hybrids may later evolve into successful and widespread pathogens. Other emerging pathogens, such as those in the Trichophyton mentanophytes species complex13, may also deserve evaluation from this perspective.Huang and collaborators recently reported the first clinical isolation of a previously undescribed invasive fungal pathogen, Rhodosporidiobolus fluvialis14. This genus belongs to the basidiomycete family Sporidiobolaceae, which also includes Rhodotorula and Sporobolomyces. Pathogenicity within this family is rare, with Rhodotorula mucilaginosa being the most frequent opportunistic pathogen15, although other sporadic human infections have been reported for Rhodotorula glutinis, Rhodotorula minuta15 or Sporobolomyces salmonicolor16. Besides the above-cited report, there are no previous publications of Rhodosporidiobolus species in a clinical context, underscoring the relevance of the new identifications.Huang and colleagues characterised two clinical isolates of R. fluvialis from the China Hospital Invasive Fungal Surveillance Net, designated TZ579 and NJ10314. Both isolates were recovered from blood samples from unrelated elderly patients with systemic infection and multiple comorbidities. Based on ITS and D1/D2 sequences, both isolates were classified as R. fluvialis and used to generate high-quality genome assemblies based on long and short sequencing reads. Notably, the genome assembly of TZ579 was twice the size of NJ103, and Huang et al. divided it into two highly syntenic subgenomes (A and B). This observation drew our attention because, in the absence of explicit phasing, the separate assembly of the two haplotypes usually implies a high divergence between them. Huang and collaborators also reported a peak of approximately 90% sequence identity between homologous proteins in the two subgenomes, a value that appears unexpectedly low for two haplotypes of the same species. Although the TZ579 strain was described as a diploid form of R. fluvialis, we hypothesised that this genomic heterozygosity results from a hybrid nature of TZ579. The original paper ruled out the possibility of a mixed culture with two species, as they confirmed by Fluorescence-Assisted Cell Sorting (FACS) the ploidy of the TZ579 isolate. Here, we re-analysed publicly available genomic data from the two reported clinical isolates of R. fluvialis to test whether TZ579 represents an interspecies hybrid and to infer the evolutionary origin of its two subgenomes.ResultsIsolate TZ579 is an allodiploid hybridTo assess the potential hybrid nature of R. fluvialis TZ579, we first evaluated the reference assemblies generated by Huang and collaborators14. NJ103 and TZ579 assemblies comprised 19 and 36 scaffolds, respectively, and both showed high contiguity (N50 > 1.8 Mb) and identical GC% content (60.18%, Supplementary Data 1). The total assembly size of TZ579 (51.9 Mb) and its predicted gene count (18,537) were approximately twice those of NJ103 (25.8 Mb, 9305 genes), consistent with the presence of two subgenomes in TZ579. Dot-plot analysis confirmed the existence of large syntenic blocks between the two assemblies, together with several inversions (Supplementary Fig. 1A), in agreement with the original report.We next examined ploidy based on mismatch frequencies in read alignments. Consistent with a haploid nature of NJ103 and each of the subgenomes in TZ579, we found no evidence for ploidies higher than 1. This result can be reconciled with the original report’s FACS data indicating TZ579 is diploid14, because alignment-based ploidy inference is performed at each subgenome level due to the low level of cross-mapping. Indeed, we observed very low percentages of mapping reads between TZ579 short reads and NJ103 (58%), as compared to >92% mapped reads when using their respective self-references. This low mapping percentage points to high genetic distance between the two isolates.Depth of coverage analyses further showed homogeneous coverage for NJ103 and a partial aneuploidy (duplication) of scaffold 3 in TZ579 (Supplementary Fig. 1B, C). This scaffold had very low heterozygosity, indicating that the observed duplication was a recent event that had accumulated few mutations.After confirming the overall quality of the assemblies, we compared TZ579 filtered reads against the NJ103 assembly using a k-mer analysis. In brief, this analysis subdivides reads or assembly sequences into fragments of fixed length k (k-mers) and counts the frequency of each k-mer in the sequence dataset of interest. We observed a surprisingly low fraction of k-mers shared between the two isolates (18%, blue area in Fig. 1). This is consistent with the low mapping rate (58%) of TZ579 reads against the NJ103 assembly noted above, and indicates a substantial sequence divergence between the two isolates. The TZ579 k-mer spectrum showed two multiplicity (i.e. k-mer coverage) peaks, with most k-mers concentrated in the lower-coverage peak (Fig. 1). K-mers around the higher-coverage peak, with approximately twice the coverage, were not the result of the duplication in scaffold 3 (Supplementary Fig. 1D), but instead most likely represent shared k-mers between the two subgenomes. The concentration of k-mers in the lower multiplicity peak in TZ579 strongly indicates the existence of two highly divergent subgenomes. This pattern is characteristic of allodiploid hybrids, which are diploid organisms with two chromosome sets from different species3.Fig. 1: K-mer spectrum of TZ579 filtered reads against the NJ103 assembly.Full size imageThe x-axis represents the multiplicity (i.e. coverage) of each k-mer in the TZ579 reads. The primary peak around ~50× represents k-mers unique to one of the two divergent subgenomes, while the secondary peak at ~100× represents conserved k-mers shared across both subgenomes (appearing at double the frequency). The colour of the area indicates k-mer presence in the NJ103 reference assembly: the blue area denotes k-mers in TZ579 reads that appear once in NJ103, while the black area indicates k-mers that do not appear in NJ103. The predominance of black highlights the substantial genetic distance between the two isolates.To further test the hypothesis that TZ579 was an allodiploid hybrid, we sought to separate the two subgenomes. For this, we used minimap2 to align the contigs of the published assemblies of TZ579 and NJ103 and separated the scaffolds based on the average match percentage to NJ103 (Fig. 2). One set of scaffolds showed a mean match percentage of 32% and was assigned to subgenome A, which is genetically more similar to NJ103. Another set of scaffolds had a mean match percentage of 7% and was assigned to subgenome B, genetically more distant from NJ103. We could not find homologous sequences for scaffold 36 in NJ103. Blast results against the NCBI RefSeq database returned some hits to genera Bipolaris and Exserohilum, from Ascomycota. We thus considered this 5 kb scaffold as a contamination and excluded it from our subgenome definition.Fig. 2: Synteny plot and subgenome classification of TZ579 scaffolds relative to NJ103.Full size imageSyntenic blocks were defined with minimap2 alignments and do not represent inversions. Link colour represents the identity (match percentage) in each alignment block. TZ579 scaffolds are colour-coded to indicate subgenome assignment. TZ579 scaffolds have been arranged in the top and bottom tracks accordingly, with scaffold 3 separated into two regions based on identity to NJ103, head (03h) and tail (03t). Scaffolds were arranged so that homologous regions approximately align vertically. Correspondence of scaffold numbers and full names can be found in Supplementary Data 2.Despite their genetic dissimilarity, the two subgenomes remained highly collinear (Supplementary Fig. 1E). Notably, our assignment to subgenomes A and B differs from that proposed by Huang et al. using a k-mer-based method (Supplementary Data 2).The minimap mapping provided little evidence for recombination between the two subgenomes. The only exception was scaffold 3, which contained adjacent regions with markedly different levels of sequence identity to NJ103. The transition point coincided with the ITS locus and presented a very high depth of coverage (>10,000×, Supplementary Fig. 2), consistent with a collapsed ribosomal DNA array. Because repetitive regions frequently lead to assembly artefacts, and considering the abnormal similarity results above, we interpret scaffold 3 as a probable misassembly joining sequences from the two subgenomes, whereas the homologous scaffolds 18 and 29 remained correctly separated. We therefore separated scaffold 3 at the last high-similarity alignment block to NJ103, immediately adjacent to the amplified ITS-containing region (Supplementary Fig. 2), generating fragments 3-head and 3-tail, which were assigned to subgenomes A and B, respectively, for subsequent analyses. Recombination leading to loss of heterozygosity or gene conversion between the two subgenomes would be detected as regions missing in one of the subgenomes and having double depth of coverage in the alternative subgenome. Our mapping and depth of coverage analysis did not detect any such regions. Another possibility is that the recombination event has been correctly assembled, generating identical regions between subgenomes. We did detect two regions where genes were over 98% identical between subgenomes A and B, in sharp contrast with the average 85% identity (Supplementary Fig. 3A). These corresponded to terminal regions of scaffolds 32–35 (Supplementary Fig. 3B). All four scaffolds contained an identical ~35 kb region, even though they are not homologous scaffolds, weakening the possibility that this represents a recombination event. Instead, the small size of these four scaffolds points to a possible misassembly. We also checked repetitive sequences in this region, with no clear pattern pointing to mobile genetic elements or low complexity sequences.Subgenome A in TZ579 is related to NJ103Minimap alignments were useful to observe overall sequence similarity, but they are not designed to precisely quantify sequence divergence between the haplotypes. To obtain a more direct estimate of relatedness between TZ579 subgenomes and NJ103, we predicted genes in their respective assemblies and identified homologous gene pairs based on all-versus-all Blast comparisons. Homologous gene pairs between subgenome TZ579-A and NJ103 showed a mean nucleotide identity of 94% (n = 2533), whereas homologous gene pairs between subgenome TZ579-B and NJ103 showed a mean identity of 85% (n = 2194). These values indicate that subgenome A is substantially closer to NJ103 than subgenome B, although both differ markedly from it.The 15% divergence between TZ579-B and NJ103 strongly suggests they belong to different species. The 6% divergence between TZ579-A and NJ103 is more difficult to interpret, but it is also high enough to question whether both belong to the same species. In addition, the distribution of gene identity between TZ579-A and NJ103 was bimodal (two main peaks at 94.7% and 99.9% identity, Fig. 3A). Mapping these values along TZ579-A revealed several clusters of contiguous genes that were almost identical to NJ103 (Fig. 3B). To further characterise this bimodal identity between NJ103 and TZ579-A, we added coverage and SNP density of NJ103 reads against the TZ579-A scaffolds (Fig. 3B and Supplementary Fig. 4). These results further exposed clear blocks of high identity. This mosaic pattern may reflect a more complex evolutionary history, with a past introgression between the lineage represented by the majority of subgenome A and a lineage more closely related to NJ103. This introgression probably happened before the hybridisation event that gave rise to TZ579.Fig. 3: Bimodal sequence identity between genes in TZ579 subgenome A and NJ103.Full size imageA Distribution of gene identity between TZ579 subgenomes and NJ103. A solid vertical line represents the average gene identity of each subgenome. Dotted lines represent the two local maxima of TZ579-A gene identity (94.7% and 99.9%). B Mapping of NJ103 reads against TZ579-A scaffold 20. The top grey histogram represents relative coverage in 5 kb windows, relative to the median coverage of the scaffold. The next track represents SNP density per kilobase, calculated in 5 kb windows. The last track represents genes predicted in TZ579-A, coloured by their blast identity to their NJ103 homologue. A high-identity region is delimited with dashed lines.As the observed low levels of gene sequence identity suggested that TZ579-A might be considered a different species from NJ103 (R. fluvialis), we next sought to clarify its phylogenetic origin. As no additional R. fluvialis genome assemblies are available, we resorted to phylogenetic reconstruction using publicly available ITS sequences from Rhodosporidiobolus species (Fig. 4, Supplementary Data 3). There are known limitations of this marker for classifying fungal species, and results must be interpreted with caution, but it remains the only marker available for a broader set of R. fluvialis isolates. This analysis identified several sequences closely related to TZ579-A, including three marine R. fluvialis isolates from Thailand17 (Supplementary Data 4). The broader sampling, as compared to the original study by Huang et al.14, provided additional context for the group, suggesting possible taxonomic ambiguity and misclassifications within R. fluvialis and Rhodosporidiobolus azoricus. Note that this analysis did not include TZ579-B, as the assembly does not contain its ITS.Fig. 4: Phylogeny based on the nucleotide sequence of the ITS region from Rhodosporidiobolus species.Full size imageITS sequences from TZ579-A and NJ103 are marked in bold. A scale bar for evolutionary distance in substitutions per site and bootstrap indicators are provided. The left column indicates the species assigned to each sequence in the GenBank database. The right column indicates the broad categories of the source of isolation of R. fluvialis and closely related strains when available, obtained either from GenBank records or scientific literature. More detailed information on the source of isolation is available in Supplementary Data 4. The branch for sequence KJ706794.1 was shortened to ease visualisation.The ITS identity between NJ103 and TZ579-A was 98.94%, which is just above the established thresholds to discriminate between different species in Basidiomycota (98.61%)18. This marker therefore supports the original assignment of both isolates to R. fluvialis made by Huang et al. However, because ITS-based species delimitation has well-recognised limitations19, and because genome-scale comparisons indicate substantially greater divergence (6%), the available data do not rule out that R. fluvialis may represent a cryptic species complex. Public ITS records also revealed additional unreported clinical isolates assigned to R. fluvialis, including two from Brazil (Fig. 4, Supplementary Data 4), suggesting that this pathogen may be more widely distributed and more frequently recovered from clinical settings than currently reflected in the literature.Subgenome B in TZ579 originates from an unsampled Rhodosporidiobolus speciesAfter identifying subgenome A as the TZ579 subgenome most closely related to NJ103, we next investigated the origin of subgenome B. For this, we calculated whole-genome average nucleotide identity (ANI) between all publicly available Rhodosporidiobolus assemblies in NCBI (Supplementary Data 5). All reference assemblies were broadly similar to NJ103 in size (~25 Mb) and GC content (55-65%); all of them had high BUSCO completeness (>87%), and none showed evidence of bacterial contamination (Supplementary Data 1). The contiguity of the assemblies was variable, with some highly fragmented references.ANI between the two TZ579 subgenomes was 82%, which was similar to the ANI observed between TZ579 subgenomes and other species, such as Rhodosporidiobolus microsporus (81%) or R. azoricus (80%). These values are consistent with TZ579 being an interspecies hybrid. TZ579-A shared 92.37% ANI with NJ103, below the 95% threshold previously proposed for fungal species delimitation20. This value is also lower than the 94% mean identity observed across homologous genes, as ANI also captures divergence in intergenic regions. Together, these observations further support the possibility that subgenome A represents a species closely related to, but distinct from, R. fluvialis. In contrast, subgenome B did not show elevated ANI with any available Rhodosporidiobolus assembly (all