Evolutionary origin and population diversity of a cryptic hybrid pathogen

Article https://doi.org/10.1038/s41467-024-52639-1 Evolutionary origin and population diversity of a cryptic hybrid pathogen Jacob L. Steenwyk 1,2,3 , Sonja Knowles 4 , Rafael W. Bastos 5,6 , Charu Balamurugan 2,3 ,DavidRinker 2,3 , Matthew E. Mead 2,3,18 , Christopher D. Roberts 4 ,HuzefaA.Raja 4 , Yuanning Li 7 , AnaCristinaColabardini 5 , Patrícia Alves de Castro 5 , Thaila Fernanda dos Reis 5 , Adiyantara Gumilang 2,3 ,MaríaAlmagro-Molto 8 ,AlexandreAlanio 9,10 , Dea Garcia-Hermoso 9 ,EndrewsDelbaje 5 , Laís Pontes 5 , Camila Figueiredo Pinzan 5 , Angélica Zaninelli Schreiber 11 , David Canóvas 12,19 , Rafael Sanchez Luperini 5 ,KatrienLagrou 13,14 , Egídio Torrado 15 , Fernando Rodrigues 15 , Nicholas H. Oberlies 4 , Xiaofan Zhou 16 ,GustavoH.Goldman 5,17 & Antonis Rokas 2,3 Cryptic fungal pathogens pose disease management challenges due to their morphological resemblance to known pathogens. Here, we investigated the genomes and phenotypes of 53 globally distributed isolates of Aspergillus section Nidulantes fungi and found 30 clinical isolates—including four isolated from COVID-19 patients—were A. latus, a cryptic pathogen that originated via allodiploid hybridization. Notably, all A. latus isolates were misidentified. A. latus hybrids likely originated via a single hybridization event during the Miocene and harbor substantial genetic diversity. Transcriptome profiling of a clinical isolate revealed that both parental subgenomes are actively expressed and respond to environmental stimuli. Characterizing infection-relevant traits —such as drug resistance and growth under oxidative stress—revealed distinct phenotypic profiles among A. latus hybrids compared to parental and closely related species. Moreover, we identified four features that could aid A. latus taxonomic identification. Together, these findings deepen our understanding of the origin of cryptic pathogens. More than 150 million cases of severe fungal infections result in an estimated 1.7 million deaths per year worldwide1.Beyondthe profound toll of lives lost, fungal pathogens pose a substantial economic burden conservatively estimated to be 11.5 billion US$ but may be as high as 48 billion US$2–4. The impact of fungal pathogens on human welfare will likely increase due to global climate change5. The importance of fungal diseases was recently acknowledged by the establishment of the first-ever priority pathogen list by the World Health Organization6, and it has risen in prominence as a prevalent cause of co-infection of COVID-19 patients7–9. Even though a great deal is known about several aspects of fungal pathogenicity, we know surprisingly little about the repeated evolution of pathogenicity in fungi, such as the traits and genetic elements that contributed to it10. Evolutionary approaches have been highly informative both in understanding how pathogens differ from their non-pathogenic close relatives11–14 but also in discovering cryptic species, defined as species that are genetically distinct from known species but difficult to distinguish using morphology15,16.Cryptic pathogens have been identified in several fungal genera that contain pathogens, including Cryptococcus,Aspergillus,Histoplasma,and Fusarium13,15,17–21. Received: 25 June 2024 Accepted: 16 September 2024 Check for updates A full list of affiliations appears at the end of the paper. e-mail: [email protected]r;[email protected] Nature Communications | (2024) 15:8412 1 1234567890():,; 1234567890():,; Recently, we used short-read genome sequencing data from a sample of six clinical isolates to show that Aspergillus latus is a novel cryptic fungal pathogen21.A. latus arose via allodiploid hybridization between A. spinulosporus and a close, unknown relative of A. quadrilineatus, wherein (nearly) the entire genomes of both parents were combined during hybridization. However, the small number of isolates examined, and their highly fragmented short-read genome assemblies left numerous questions unanswered—such as how the genome of the hybrids differs from the parental genomes, the timing and number of hybridization events, and genetic diversity and phenotypic variation in the species. Increased understanding of A. latus genomic and phenotypic trait variation and differences from close relatives could increase accuracy of taxonomic identification of the species and, in the longer term, aid the development of clinical diagnostic techniques. To further understand the biology and evolutionary origins of the cryptic fungal pathogen A. latus, we generated high-quality genome assemblies using long- and/or short-read sequencing technologies on a substantially expanded set of 53 globally distributed isolates (50 clinical isolates; three type strains) from section Nidulantes of the genus Aspergillus. We found that 30 clinical isolates correspond to A. latus and the remaining 20 clinical isolates to A. spinulosporus (8), A. quadrilineatus (1), or A. nidulans (11); one previously published strain of A. latus was included in this dataset21.A. latus likely arose approximately 13 million years ago and exhibits substantial genetic diversity among gene families and biosynthetic gene clusters. Transcriptomic analysis of a clinical isolate revealed that both parental subgenomes are actively expressed and respond to environmental perturbations. Profiling of each species across infection-relevant chemical and physiological traits suggests hybridization contributes to a unique phenotypic profile. Comprehensive genomic and phenotypic characterization of A. latus identified four traits that distinguish A. latus from closely related species. Characterization of A. latus,anunderstudied cryptic fungal pathogen, expands our understanding of the origins offungalpathogenicity and is a key step toward preventing and combating disease. Results Genome sequencing and species identification of 50 clinical isolates Thegenomesofthreetypestrains(A. latus NRRL 200 T;A. spinulosporus NRRL 2395 T;A. quadrilineatus NRRL 201 T) and 40 clinical isolates (Supplementary Data 1) were sequenced and assembled using both long- (Oxford Nanopore) and short-read (Illumina) technologies. The resulting genome assemblies are highly contiguous with mean N50 and L50 values of 3.09 ± 1.27 Mbp and 13.24 ± 21.19 scaffolds, respectively (Supplementary Fig. 1), an improvement of ~10-to-14-fold from previously available genome assemblies21. Nine additional isolates were sequenced using only short-read technology (mean N50: 262.73 ± 40.94 Kbp; mean L50: 74.00 ± 10.98 scaffolds). Species determination via molecular phylogenetic analyzes of the taxonomically informative loci of β-tubulin and calmodulin (Supplementary Fig. 2 and 3) placed the clinical isolates into the clades of four distinct species: A. latus (29), A. spinulosporus (8), A. quadrilineatus (1) and A. nidulans (11). Notably, species determination using matrixassisted laser desorption/ionization time-of-flight (MALDI-TOF) mass spectrometry, a standard method for microbe identification in clinical laboratories22, revealed variable accuracy for each species (p< 0.001, Fisher’s exact test; Supplementary Data 2); all examined A. latus isolates and six out of nine examined A. nidulans isolates were misidentified, whereas all A. spinulosporus isolates were accurately identified. Misidentification of A. latus is largely due to the lack of inclusion of standards for the species fromMALDI-TOF databases. This finding suggests that A. latus infections are likely underreported, and that A. latus data should be included in future updates of MALDI-TOF databases. A previously sequenced and assembled clinical isolate of A. latus,MO46149 21, was added to the dataset resulting in a total of 53 isolates: 50 clinical isolates and three type strains (Supplementary Fig. 4). Clinical isolates were isolated from patients with diverse, but primarily lung-related, pathologies, such as cystic fibrosis and chronic granulomatous disease (Supplementary Data 1). Six misidentified A. latus clinical isolates came from patients co-infected with COVID-19. While some analyzes have been done on A. fumigatus COVID-19- associated pulmonary isolates9,23,identifyingA. latus co-infections adds to the evidence that other Aspergillus species can also establish coinfections with SARS-CoV-224,25. Aspergillus latus is an allodiploid hybrid that has retained the genome organization and content of both parental genomes Genome features, phylogenomics, and macrosynteny provide unequivocal support for the hybrid origin of A. latus. The genomes of A. latus isolates are approximately twice the size, encode twice the number of genes compared to closely related species in section Nidulantes, and contain two copies of the vast majority of nearuniversally single-copy orthologs (or BUSCO genes; specifically, 96.68 ± 2.46% of BUSCO genes are duplicated), all features suggestive of a diploid genome (Fig. 1A; Supplementary Figs. 1 and 5). Macrosynteny analysis between A. spinulosporus NRRL 2395 T,A. quadrilineatus NRRL 201 T,andA. latus NRRL 200 Trevealed large genomic segments that are syntenic with A. spinulosporus and A. quadrilineatus, further supporting an allopolyploid origin (Fig. 1B). Further evidence of hybridization was identified using phylogenomics; specifically, the topologies of 3746 BUSCO genes distributed across all 16 A. latus chromosomes and present in two copies were ohnologs, i.e., they showed evidence of one BUSCO gene originating from A. spinulosporus and the other from a close, but unknown, relative of A. quadrilineatus (Fig. 1C; Supplementary Fig. 6), consistent with previous findings21. Hereafter, we will refer to the gene duplicates that have originated from the allodiploidization that gave rise to A. latus as ohnolog pairs. Reciprocal best BLAST hit analysis between A. spinulosporus NRRL 2395 Tand A. quadrilineatus NRRL 201 Tidentified 9443 orthologs, representing the upper limit of ohnolog pairs in the A. latus hybrids. Reciprocal best BLAST hit analysis between individual A. latus hybrid protein-coding sequences and the combined coding sequences of A. spinulosporus NRRL 2395 Tand A. quadrilineatus NRRL 201 Trevealed an average of 91.63 ± 2.41% (8652.26 ± 227.96) of ohnologs have been retained (Table 1). Among ohnologs where one copy has been retained, 6.26 ± 1.58% (591.58 ± 149.09) originate from the A. spinulosporus subgenome and 1.79 ± 1.16% (169.26 ± 109.79) originate from the A. quadrilineatus-like subgenome. Lastly, 0.32 ± 0.12% (29.90 ± 11.29) of ohnologs pairs were inferred to have been lost. Together, these results indicate the vast majority of ohnologs have been retained, suggestive of genome stability among the A. latus hybrids. Subgenomes have similar evolutionary rates and recombine rarely To further understand the genome organization and diversification of the two subgenomes, we next compared the evolutionary rates between parental subgenomes using tree- and sequence-based measures among 1712 ohnolog pairs. This analysis revealed that the subgenomes are evolving at statistically indistinguishable rates (Fig. 1D; F=3.72, p= 0.08, Multi-factor ANOVA; p= 1.0 for pairwise comparisons for tree- and sequence-based measures, Tukey’s Honest Significant Differences Test). While one parental genome marginally evolved faster than the other at times (Fig. 1E), most ohnolog pairs (N= 1073; 62.68%) evolved at the same rate. This finding suggests that the parental subgenomes are evolving similarly in the A. latus hybrids. The chromosome-level genome assemblies and the high nucleotide sequence divergence (7%) of the subgenomes enabled us to Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 2 examine whether recombination occurs between the subgenomes of A. latus. Using the A. latus NRRL200Tgenome as a reference, we found that its two subgenomes display little evidence of recombination, as measured by scaffolds containing syntenic blocks of 15 or more consecutive genes assigned to one or the other parent (Fig. 1B). Across a broader sampling of isolates with highly contiguous genome assemblies (N= 22), we found that 88.30% (1057/1197) of scaffolds did not show signatures of recombination between the parental subgenomes (Fig. 1F; Supplementary Data 3; Supplementary Fig. 8; Supplementary Data 4). Examination of the relationship between the frequency of recombination and the genome assembly quality ofthe A. latus isolates revealed that the number of scaffolds with evidence of recombination is significantly correlated to the number of scaffolds in a genome Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 3 assembly (ρ=0.57, p= 0.007, Spearman rank correlation) and the genome assembly N50 (ρ=−0.48, p= 0.027, Spearman rank correlation) (Supplementary Fig. 9). Therefore, we cannot exclude the possibility that at least some of the observed recombination is an artifact stemming from genome assembly errors. We conclude that A. latus hybrid genomes display little, if any, evidence of recombination between the parental genomes. The origin of A. latus in time The denser sampling of isolates allowed us to investigate the number and timing of the hybridization event(s) that gave rise to A. latus.The maximum likelihood phylogeny of our phylogenomic dataset, which includes all the parental subgenomes in the A. latus isolates as well as all A. spinulosporus,A. quadrilineatus,andA. nidulans isolates, suggests that A. latus originated via two distinct hybridization events –one for the environmental isolate A. latus NRRL 200 Tand one for all other A. latus clinical isolates (Fig. 2). However, levels of sequence divergence are very low among most A. latus isolates (evidenced by their very short branch lengths in the phylogeny; Fig. 1C) and there is low support for the divergence of A. latus NRRL 200 Tfrom other isolates (SH-aLRT =71.7;UFBoot=97;Fig.1C; Supplementary Fig. 7). Thus,we conducted phylogenetic topology constraint tests of the alternative hypothesis that there was a single hybridization event. We used an approximately unbiased test to determine whether a constrained tree topology where all the A. spinulosporus subgenomes of A. latus isolates were monophyletic had a significantly lower likelihood score than the unconstrained maximum likelihood tree. The likelihood scores of the two trees were not significantly different in both the 3746-gene data matrix ((ΔL = 139.84, p= 0.72) and a 2918-gene data matrix comprised of genes that had at least 90% taxon occupancy (ΔL=29.564,p=0.44). Different methods of topology testing yielded the same result (Supplementary Table 1). Therefore, the most parsimonious explanation is that A. latus originated via a single hybridization event. Relaxed molecular clock analyzes suggest A. latus arose approximately 13.9–13.1 (lower bound: 16.6; upper bound: 7.4) million years ago during the Miocene (Fig. 2). Variation in A. latus gene content is pronounced among biosynthetic gene clusters To begin examining the genetic diversity of A. latus,weassigned orthologous groups of genes (a proxy for gene families) to four different bins of occupancy: a bin of core genes present in all A. latus isolates and three different bins of accessory genes present in <100% of A. latus Isolates. This analysis primarily leveraged genome assemblies generated using long- and short-read technologies (N= 22) since fragmented genomes contribute to error in assessing gene content variation26. Accessory genes were assigned as softcore (present in nearly all isolates; N= 21), shell (present in 5-95% of isolates; 2 ≤N<21), or cloud (present in less than 5% of isolates; N= 1). Across 10,078 gene families, 7485 (74.27%), 1448 (14.37%), 793 (7.87%), and 352 (3.49%) were determined to be core, softcore, shell, and cloud, respectively (Fig. 3A).Wethenusedthesameapproachtoexaminegenecontent variation of biosynthetic gene clusters (BGCs) involved in secondary metabolism (hereafter referred to as the BGCome). Assigning BGCs to BGC families revealed that 46 (48.42%), 9 (9.47%), 25 (26.32%), and 15 (15.79%) BGC families were assigned to the core, softcore, shell, and cloud, respectively (Fig. 3B). Examination of the relationship between the number of gene families or BGCs in the core and accessory (softcore, shell, and cloud) gene families underscored their differing evolutionary dynamics. Core genes are substantially more common than accessory genes among A. latus isolates (Fig. 3C). Incontrast, the accessory BGCome islarger than the core (Fig. 3D). Further sampling of chromosome-level assemblies will enable more preciseestimates of the size of the core and accessory parts of the A. latus genome. Proteins encoded by the core and softcore gene families significantly differ in protein sequence lengths (χ2=165.57,p< 0.001, df = 3; Kruskal–Wallis Rank Sum test) (Fig. 3E). Specifically, the average protein sequence length among core and softcore genes is 515.25 ± 369.06 and 517.45 ± 389.26 amino acids, respectively, and are significantly longer than shell and cloud genes that have an average length of 441.94 ± 412.28 and 401.68±369.06 amino acids, respectively (p< 0.001 for both comparisons; Dunn’stest, Benjamini–Hochberg (BH) multi-test correction). Shell genes are also longer than cloud genes (p= 0.004; Dunn’s test, BH multi-test correction). Core and accessory genes also differed in the proportion of genes that could be functionally annotated, wherein a greater proportionof genes wereannotated in the core and softcore gene families and fewer in the shell and cloud (Fig. 3F). Core and accessory BGCs significantly differed in the number of genes encoded in each (χ2= 166.19, p< 0.001, df = 3; Kruskal–Wallis Rank Sum test) (Fig. 3G). Specifically, softcore BGCs, which have an average of 10.69 ± 3.99 genes, encode significantly fewer genes than core, shell, and cloud BGCs, which have an average of 14.45 ± 4.46, 15.04 ± 6.71, and 15.10 ± 7.04 genes, respectively (p<0.001 for all comparisons; Dunn’s test with BH multiple test correction). The products, if any, of most BGCs are not known, suggesting A. latus isolates may be capable of producing novel secondary metabolites; alternatively, numerous BGCs encoded in the A. latus genomes may be nonfunctional (Fig. 3H). Phenotypic profiles of hybrid isolates are distinct Hybridization can change the phenotypic profile of organisms27. Extensive phenotyping of infection-relevant traits (i.e., growth under cell wall perturbation and integrity stresses, growth under oxidative stress, growthat 44°C, 37 °C, and 30 °C, susceptibility to voriconazole, amphotericin B, and caspofungin, growth in minimal media, asexual spore size, and production of cleistothecia) across all isolates including the synthetic diploid strain R153XR21 of A. nidulans,revealedthat hybrids were phenotypically distinct from the other three species, but were most similar to A. spinulosporus (Fig. 4A). Drug resistance and a cell wall stressor were the largest contributors to the variance along the first dimension in principal component space (Fig. 4B; Supplementary Fig. 11), whereas oxidative stress and asexual spore size, which Fig. 1 | Aspergillus latus originated from hybridization and parental subgenomes evolve at similar rates and remain distinct. A The genomes of Aspergillus hybrids have larger genomes (x-axis), encode more genes (y-axis), and numerous near-universally single copy orthologs (or BUSCO genes) are present in duplicate copies (data point size). The inset table provides averages ± standard deviation for each metric for the diploid genomes of the A. latus hybrids (Alat; purple) and the haploid genomes of A. spinulosporus (Aspi; blue), A. quadrilineatus (Aqua; red), and A. nidulans (Anid; gray). BSynteny analysis between the A. latus,A. quadrilineatus,andA. spinulosporus type strains. Syntenic blocks contain a minimum of 15 genes robustly assigned to one of the parental species. Scaffolds greater than 100 kilobases are depicted. CPhylogenomic tree of the 49 clinical isolates and type strains for each species. Note that the genomes of the hybrid isolates and the A. latus type strain were split into their two subgenomes (blue: A. spinulosporus subgenome; red: A. quadrilineatus-like subgenome). Links connect subgenomes from the same isolate. DComparative analysis of evolutionary rates using treebased and sequence-based measures revealed parental subgenomes are evolving at indistinguishable rates [p=1.0forbothtests;Tukey’sHonestSignificant Differences Test, which was conducted after a multi-factor ANOVA test also revealed no significant differences between the two parental genomes (p= 0.054)]. EThe majority of ohnolog pairs are evolving at the same rate (gray), although sometimes the A. spinulosporus subgenome evolves slightly faster than the A. quadrilineatuslike subgenome (blue) and vice versa (red). FThere is little evidence of recombination between parental subgenomes; specifically, 1057 (88.30%) scaffolds had no evidence of recombination. Source data are provided as Source Data files. Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 4 his correlated with genome size21, contributed most to the variance along the second dimension (Fig. 4C). Notably, 14 A. latus and five A. spinulosporus isolates had high drug resistance to caspofungin compared to reference isolates of A. nidulans (varying from 2.0 to > 4.0 ug/ mL and 0.03 ug/mL, respectively). Variation among the species was also observed in an invertebrate moth model of disease (p=0.04;Log- Rank test) (Supplementary Fig. 12; Supplementary Data 6). Variation among clinically relevant traits suggests that each species maywarrant unique strategies to prevent and combat disease. Examination of secondary metabolite profiles among six isolates— three A. latus isolates and single representatives from A. spinulosporus, A. quadrilineatus,andA. nidulans—revealed that A. latus hybrids were qualitatively more like A. quadrilineatus (Fig. 4D). Correlation analysis revealed that A. latus CI 1908 was most similar to A. quadrilineatus NRRL 201 T(Fig. 4E). Examination of the absolute concentration of a mycotoxin, sterigmatocystin, in the six isolates revealed that A. latus hybrids produced more toxin than either parental species (A. spinulosporus and A. quadrilineatus) but less than A. nidulans (Fig. 4F). Together with other phenotyping data, these results suggest that A. latus hybrids are phenotypically distinct from other species but are more like one parent for some traits, more like the other parent for others, and differ from both parents in yet others. We next determined if A. latus gene expression patterns were biased toward one or the other parent. RNA-sequencing of A. latus CI 1908 at 30 °C and 37 °C revealed nearly equaltranscript abundances in both parental subgenomes (Fig. 4G,H). For example, at 37 °C the A. spinulosporus and A. quadrilineatus-like subgenomes had average transcript per million (TPM) values of 47.72 ± 355.86 and 49.44 ± 369.27, respectively; a similar observation was made for transcript abundances at 30 °C. Examining the estimated codon usage bias of each parental genome (Supplementary Fig. 13) corroborates this finding, revealing that both subgenomes are nearly equally optimized. Differential expression analysis between the two conditions revealed similar numbers of significantly differentially expressed genesineachparentalgenome(Fig.4I; 218 and 213 upregulated genes and 242 and 232 downregulated genes in the A. spinulosporus and A. quadrilineatus subgenomes, respectively; p<0.01 and |log2(fold change)| ≥2 for differential expression analysis). Furthermore, the Table 1 | Genome features and ohnolog information for the parental and hybrid isolates Genome Size (Megabases) Total Number of Genes RBBH with A. spinulosporus RBBH with A. quadrilineatus Retained ohnologs Aspi singleton from an ohnolog Aqua singleton from an ohnolog Ohnolog pairs lost Unique genes A. spinulosporus NRRL 2395 T 31.34 10,496 —9443 —— — —1053 A. quadrilineatus NRRL 201 T 32.94 11,193 9443 —————1750 A. latus hybridNRRL 200 T 63.92 21,496 9938 9677 8339 532 530 42 1881 All A. latus hybrids 63.66 ± 1.54 21,645.9 ± 655.99 10,155.68 ± 128.22 9663.35 ± 197.50 8652.26 ± 227.96 591.58 ± 149.09 169.26 ± 109.79 29.90 ± 11.29 1826.87 ± 596.86 RBBH reciprocal best BLAST hit, Aspi A. spinulosporus,Aqua A. quadrilineatus. Fig. 2 | The origin of A. latus in geologic time. Relaxed molecular clock analyzes of the origin of each parental subgenome suggest the A. latus hybrids first arose approximately 13.7–13.1 million years ago. Divergence times are depicted next to each internode (depicted as: divergence time [lower bound, upper bound]). Species and colors are indicated as follows: A. spinulosporus (Aspi; blue), A. quadrilineatus (Aqua; red), and A. nidulans (Anid; gray). Source data are provided as a Source Data file. Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 5 numbers of differentially expressed parental genes are not significantly different in the two parental subgenomes (Upregulated genes: χ2= 0.058, p= 0.81, df = 1; Downregulated genes: χ2= 0.211, p= 0.65, df = 1; Chi-squared test). Among ohnologs, numerous differentially expressed genes belonged to ohnolog pairs (186 among upregulated genes; 160 among downregulated genes), revealing conserved responses to environmental stimuli (Fig. 4J). Differentially upregulated ohnologs were enriched in functional categories such as glutathione metabolic process, GO:0006749 (Supplementary Data 7); glutathione is involved in many bioprocesses including helping confer resistance against heat shock28. Ohnologs also had divergent, subgenome-specificresponses.Specifically, among ohnologs, 90 and 114 genes from the A. spinulosporus subgenome and 92 and 122 genes from the A. quadrilineatus subgenomes weredifferentially upregulated and downregulated, respectively. Subgenome-specific genes were not enriched in any functional categories. These findings suggest that both Isolates AbsentPresent Core Accessory AbsentPresent Core Accessory Annotated Unannotated Annotated Unannotated 0 3 6 9 12 Number of gene families (x1000) A C EF 50 10 15 20 Number of strains Gene families Biosynthetic gene cluster families 0 250 0 25 50 75 100 Percentage Core Softcore Shell Cloud 236 116 6,640 845 627 166 1,285 163 500 750 *** *** *** *** Average length Core Softcore Shell Cloud ** Core Softcore 7,845 1,448 793 352 Shell Cloud Isolates B D GH 0 2 4 6 Number of BGC families (x10) 510 15 20 Number of strains 0 25 50 75 100 Percentage Core Softcore Shell Cloud Core Softcore Shell Cloud 0 10 20 30 Number of genes *** ** * *** * 46 925 15 Core Softcore Shell Cloud 11 3 4 35 18 4 14 6 Fig. 3 | Differing patterns of diversity among gene families and biosynthetic gene clusters in hybrid A. latus genomes. A Among a total of 10,078 orthologous groups of genes, 7485 are categorized as core (present in 100% of isolates; N=22) and 2593 as accessory (present in <100% of isolates; N< 22); among these accessory genes, 1448 are softcore (present in ≥95% and <100%; N= 21), 793 are shell (5-95% of isolates; 21 > N≥2), and 352 are cloud (present in less than 5% of isolates; N=1). BAmong 95 biosynthetic gene cluster families (BGCFs), 46 are categorized as core and 49 as accessory (9are software, 25 are shell, and 15 are cloud). CThe number of accessory gene families increases as the number of strains increases, suggesting additional sequencing is needed to fully capture A. latus gene content variation. N values varied between 1 and 26,393 depending on the number of strains being analyzed. DThe number of accessory BGCFs substantially increases with additional isolates suggesting the gene content variation of BGCs has also yet to be captured. Notably the core genome is larger among gene families, whereas the accessory portion is larger among BGCFs. Errors bars indicate standard deviation. EProtein sequence lengths differed among gene categories wherein core and softcore genes are longer than shell and cloud genes (N= 10,079). FAs genes were less frequently observed among isolates, they were also functionally annotated less frequently. GThe number of genes in BGCs differed per category wherein softcore BGCs tend to be smaller than BGCs categorizedas core, shell, and cloud (N=2511).HFew BGCs are predicted to make known secondary metabolites. Source data are provided as Source Data files. For panels Eand G, statistical comparisons were made using a Kruskal–Wallis rank sum test (p< 0.01 for both tests); pairwise comparisons were made using the Dunn’s test. One, two, and three asterisks represents a significance threshold of 0.05, 0.01, and 0.001, respectively. In panels C–E,G, average values are depicted and error bars indicate the standard deviation from the mean. Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 6 parental subgenomes contribute to organismal function and respond to environmental perturbations. Genomic and phenotypic traits that distinguish A. latus from close relatives Since all A. latus isolates in the present study were misidentified by traditional methods (Supplementary Data 2), we used the newly generated genomic and phenotypic data to identify traits that taxonomically differentiate A. latus from its closest relatives (Fig. 5). At the genomic level, A. latus isolates have larger genome sizes and gene repertoires than other Aspergillus species (Fig. 1A); the larger genome size of A. latus means that the species can be distinguished from its close relatives through Fluorescence-Activated Cell Sorting (or FACS) analysis of DNA content (Supplementary Data 8). Furthermore, Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 7 amplification and sequencing of single-locus molecular markers, including taxonomically informative loci, is expected to show evidence of two distinct loci that are phylogenetically distinct in a single-locus phylogeny (Fig. 1C). At the phenotypic level, A. latus spores are larger (Fig. 4C) due to their larger genome sizes21. These genomic and phenotypic differences can help taxonomically distinguish A. latus from other pathogens in section Nidulantes and may help elucidate the epidemiology and clinical burden of this cryptic pathogen. Discussion Our study generated and analyzed extensive amounts of genomic, transcriptomic, chemical, and phenotypic data to shed light on a cryptic fungal pathogen’s evolutionary history, genetic diversity, and the signature of hybridization across infection-relevant traits. For example, evolutionary analyzes indicate that the hybrid A. latus first arose approximately 13 million years ago (Figs. 1and 2). While our longread genome assemblies substantially improve previous short-read assemblies, chromosome-level assemblies will help shed light on more aspects of chromosome biology among A. latus hybrids. Moreover, although we are confident one parent is A. spinulosporus,theother parent remains unknown; discovery, sequencing, and phenotyping of the unknown parent will clarify the genomic and phenotypic contribution of the other parental species. Notwithstanding these caveats, our analyzes uncovered several aspects of A. latus genome biology. For example, when one copy of a gene is lost from an ohnolog pair, A. spinulosporus gene copies were more frequently retained (Table 1). Lastly, several isolates are from patients with COVID-19-associated pulmonary aspergillosis. We anticipate that the isolates obtained from patients also infected with SARS-CoV-2 (the virus that causes COVID- 19) and associated data may serve as useful resources to study COVID- 19-related superinfections. Examination of the gene content of A. latus revealed that it harbors considerable diversity, with 74.27% (7485/10,078) of gene families being core (Fig. 3). Although chromosome-level assemblies of a much larger number of isolates will be required to characterize the A. latus pangenome, comparison of A. latus to other species with well-defined pangenomes provides useful insights. For example, the Candida albicans pangenome is estimated to contain approximately 7325 gene families, 5432 (74.16%) of which are core (present in all isolates)29, whereas the A. fumigatus pangenome is more open and only 69.34% (7563/10,907) of pan-genes are core30,31. In contrast to gene families, the accessory BGCs outnumber core BGCs —51.58% of BGCs are accessory —corroborating previous reports that biosynthetic gene clusters evolve rapidly14,32.TheA. latus BGCome is less diverse than the A. fumigatus BGCome; specifically, 46 (48.42%) of A. latus BGCs are core, whereas 30.55% (11/36) of A. fumigatus BGCs are core33. A substantial portion of gene families, especially in the BGCome, remain uncharacterized and may represent an untapped source of novel biology and chemistry, especially in the context of an allodiploid organism. Additional sequencing of isolates will help shed light on the pangenome of A. latus. A. latus has a distinct phenotypic profile compared to parental species and A. nidulans, which is one of the species that A. latus has been previously mistyped as21. We hypothesize that some phenotypes may be additive—for example, A. latus produces the mycotoxin sterigmatocystin in an amount roughly equal to the sum of the amounts produced by its parental species (Fig. 4F). Similarly, both subgenomes Fig. 4 | Phenotypic profiling suggests that A. latus clinical isolates differ from isolates of closely related species in diverse infection-relevant traits. A Principal component analysis of 14 infection-relevant traits reveals each species is phenotypically distinct. Note, some A. latus isolates are like some A. spinulosporus isolates. BDistributions of phenotypic profiles for the three traits that contribute most to the variance along the first and Csecond dimensions. The y-axis for radial growth in the presence of stressors (calcofluor white and menadione).The red dot for conidia size is a measurement from a synthetic diploid isolate. For panels Band C,N=46. DSecondary metabolite profiling of six representative strains and Ecorrelation analysis reveals variation between species and within A. latus.FQuantity of sterigmatocystin production in six strains (N= 24). Errors bars represent one standard deviation. G,HNear-equal RNA abundance between parental subgenomes in A. latus CI 1908 (N= 20,025). IThe numbers of significantly differentially expressed genes between 37 °C and 30 °C are very similar; specifically, among upregulated genes, 218 are from the A. spinulosporus subgenome and 213 are from the A. quadrilineatus-like subgenome (χ2=0.058,p= 0.81, df = 1; Chi-squared test) and, among downregulated genes, 242 are from the A. spinulosporus subgenome and 232 are from the A. quadrilineatus-like subgenome (χ2=0.211,p= 0.65, df = 1; Chisquared test). JAmong significantly upregulated ohnologs genes, 186 belong to ohnolog pairs (93 total ohnolog pairs) were both upregulated. In contrast, 90 A. spinulosporus and 92 A. quadrilineatus-like genes were significantlyupregulated but the corresponding ohnolog was not differentially expressed. Among significantly downregulated genes, 160 belong to ohnolog pairs (80 total ohnolog pairs), whereas 114 A. spinulosporus and 122 A. quadrilineatus-like genes had subgenomespecificresponses.Thisfinding suggests that the A. latus expression profile in response to environmental perturbation contains both conserved elements (both members of ohnolog pairs in the two subgenomes showing the same change of gene expression) and unique ones. Source data are provided as Source Data files. For panels B,C,G,H, box plots depict the 25th, median, and 75th percentiles; whiskers represent 1.5 times the interquartile range from the hinge; and violin plots depict the underlying distribution. Fig. 5 | Genomic and phenotypic traits that taxonomically distinguish A. latus from its close relatives. A. latus (purple), A. spinulosporus (blue), and A. quadrilineatus (red) are indistinguishable in culture. At the genomic level, A. latus isolates have larger genome sizes and gene repertoires than other Aspergillus species and can be distinguished from its close relatives through Fluorescence-Activated Cell Sorting (or FACS) analysis of DNA content. Furthermore, amplification and sequencing of single-locus molecular markers, including taxonomically informative loci, is expected to show evidence of two distinct loci that are phylogenetically distinct in a single-locus phylogeny. At the phenotypic level, A. latus spores are larger than those of other species due to their larger genome size. Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 8 respond to temperature stress with similar numbers of genes, including numerous ohnolog pairs (Fig. 4G). In other cases, the hybrids are more like one parent over the other. For example, across all infectionrelevant traits examined, A. latus is closer to A. spinulosporus than to A. quadrilineatus; however, this observation may partly be explained by the fact that A. spinulosporus is one of the parental species whereas A. quadrilineatus is a close relative of the other parental species (and not the actual parental species). In contrast, A. latus CI 1908 has a secondary metabolite profile more like A. quadrilineatus than A. spinulosporus. Similarity to one parent over the other may be partly due to dominant genetic effects. These observations corroborate previous findings that complex additive and non-additive effects contribute to hybrid phenotypic profiles in fungi and plants34,35. The expression patterns of A. latus subgenomes differ from those previously observed in other allodiploid hybrids. Whereas both A. latus subgenomes are actively transcribed and nearly equally respond to environmental perturbations (Fig. 4G), the plant pathogen Verticillium longisporum displays subgenome-specific gene expression patterns36. Moreover, A. latus genomes tend to maintain allodiploidy whereas V. longisporum undergoes haploidization, a common genomic event in the aftermath of hybridization. Maintaining genetic loci after hybridization may result from selection, as previously proposed for Coccidioides fungi37. Whether selection is acting to maintain considerable portions of both parental genomes in A. latus hybrids remains an open question. Mutation rate analysis revealed that both parental subgenomes evolve at indistinguishable rates (F=3.72, p= 0.08; Multi-factor ANOVA; Fig. 1D,E), suggesting similar selective pressures on both subgenomes. Interestingly,thecoffeeplanthybridCoffea arabica,alsostemsfroman allopolyploid event, retains (near) complete parental genomes, and exhibits no discernable expression biases between subgenomes38. However, the coffee plant originated roughly 0.5 million years ago, whereas A. latus is substantially older, suggesting that allopolyploidy may be, in some cases, a stable evolutionary and genomic state. Accurate taxonomic identification of pathogenic microbes is the first step to their management. Our extensive population-level characterization of genomic and phenotypic diversity in A. latus facilitated the identification of several genomic and phenotypic traits that can reliably taxonomically distinguish A. latus from other species in section Nidulantes, setting the stage for the eventual development of clinical diagnostic tests that enable the medical mycology community to study the epidemiology and clinical burden of this cryptic pathogen. Methods Ethics statement The principles that guide our studies are based on the Declaration of Animal Rights ratified by UNESCO on January 27, 1978, in its 8th and 14th articles. All protocols adopted in this study were approved by the local ethics committee for animal experiments from the University of São Paulo, Campus of Ribeirão Preto (Permit Number: 23.1.547.60.8; Characterization of virulenceand immunopathogenicity of Aspergillus spp in the murine model). Groups of five animals were housed in individually ventilated cages and were cared for in strict accordance with the principles outlined by the Brazilian College of Animal Experimentation (COBEA) and Guiding Principles for Research Involving Animals and Human Beings, American Physiological Society. All efforts were made to minimize suffering. Animals were clinically monitored at least twice daily and humanely sacrificed if moribund (defined by lethargy, dyspnea, hypothermia and weight loss). All stressed animals were sacrificed by cervical dislocation. Genomic DNA extraction Isolates from Aspergillus section Nidulantes were acquired from several facilities and laboratories (Supplementary Data 1). Where possible, a preliminary species identification was performed by matrixassisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF MS; Microflex LT Biotyper, Bruker Daltonics, Bremen, Germany) with the MSI-2 database (https://msi.happy-dev.fr/;2019 release). The high molecular weight DNA was extracted as follows. Fungal conidia were inoculated in liquid media and let to grow at 37 °C for 16 h under agitation. The mycelia were then filtered through Miracloth, freeze-dried, and disrupted by grinding in liquid nitrogen. The ground material was transferred to a 50 mL Falcon tube, lysed in 10 mL of lysis buffer (3.75 mL Buffer A [3.19 g Sorbitol, 5 mL 1 M Tris- HCl (pH 9), 0.5 mL 0.5 M EDTA (pH 8) and ddH 2 Oupto50mL], 3.75 mL Buffer B [10 mL 1 M Tris-HCl (pH 9), 5 mL 0.5 M EDTA, 5.84 g NaCl, 1 g CTAB and ddH 2 O up to 50 mL], 1.5 mL 5 % Sarkosyl; 1 mL 1 % PVP; 100 μL Proteinase K) and incubated 30 min at 65 °C. Subsequently, 3.35 mL of 5 M Potassium acetate (pH 7.5) were added, and the samples were incubated on ice for 30 min. After 30 min centrifugation at 5,000 g and 4°C, the supernatants were transferred to a new Falcon tube. The samples were then treated with Phenol:- Chloroform:Isoamylalcohol (25:24:1), tubes were centrifuged, the aqueous phase (~8 mL) was transferred to a new tube and treated with 100 μL RNase A (10 mg/ml) for 60 min at 37°C. The genomic DNA was precipitated by adding 1/10 volume of 3 M sodium acetate and 1 volume ice-cold 96% ethanol and centrifuged 30 min at 10,000 G-force and a temperature of 4 °C. The resulting pellet was washed with 70% ethanol, dried, resuspended in 500 μL TE, and stored at -20°C until further use. The DNA concentration was estimated using a Qubit fluorometer. All other Aspergillus genomic DNA extraction were performed according to Goldman and colleagues39. The resulting DNA was sequenced using Oxford Nanopore and Illumina technologies at Nextomics Biosciences in Wuhan, China. This resulted in 10,565,680 Nanopore reads (152,924,029,920 total bases) with an average length of 15,492.31 ± 4991.34 base pairs and 1,145,484,780 Illumina reads with a length of 150 base pairs (171,822,717,000 bases). Nine isolates—which have “8954-AG”at the beginning of their identifier (Supplementary Data 1)—were sequenced using 150 bp pairedend short-read technology (NovaSeq6000) at the Vanderbilt Technologies for Advanced Genomics (VANTAGE), generating 250,464,221 reads (37,569,633,150 bases). Library preparation was conducted using the Twist Whole Genome / Metagenomics_LowpassWGS-Twist Miniprep (Twist Bioscience). Genome assembly and annotation Each isolate was assembled using two different approaches. First, a long sequencing read-based assembly was generated by NextDenovo, v2.3.040, using Oxford Nanopore sequencing data alone. The assembly was then polished by Nextpolish, v1.3.041, using both long- and shortsequencing read data. Second, a hybrid assembly was generated by MaSuRCA, v3.4.142, using both long- and short-sequencing data. Before assembly, reads were quality-trimmed using fastp, v0.12.643,with default settings. The two assemblies were compared for their assembly size, scaffold number, and scaffold N50 value. For the A. latus isolates, the hybrid assemblies were selected as they were about twice in size as the corresponding Nanopore-based assemblies. For other isolates, the Nanopore-based assemblies were selected as they were better in continuity. Among the nine genomes that were sequenced using only short-read technology, reads were first quality-trimmed using Trimmomatic, v0.3944 using the following parameters: LEADING:20 TRAILING:20 MINLEN:60. The resulting reads were used as input to SPAdes, v.3.14.045, for genome assembly. For genome annotation, all genome assemblies were first softmasked for repetitive sequences using RepeatMasker, v4.1.0 (http://www.repeatmasker.org;“-species” option set to “Aspergillus”). Gene models were generated using BRAKER, v2.1.446, which combines ab-initio gene predictors (AUGUSTUS, v3.3.347, and GeneMark-ES, v4.5948) and homology evidence (all Eurotiomycetes protein sequences in the OrthoDB, v1049, database). To examine the quality of all genome assemblies examined, diverse metrics that describe genome assembly and gene content Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 9 University of Rio Grande do Norte, Natal-RN, Brazil. 7 Institute of Marine Science and Technology, Shandong University, 72 Binhai Road, Qingdao, China. 8 Max von Pettenkofer-Institut für Hygiene und Medizinische Mikrobiologie, Faculty of Medicine, Ludwig Maximilian University, Munich, Germany. 9 Institut Pasteur, Paris Cité University, National Reference Center for Invasives Mycoses and Antifungals, Translational Mycology Research Group, Mycology Department, Paris, France. 10 Laboratoire de parasitologie-mycologie, AP-HP, Hôpital Saint-Louis, Paris, France. 11 School of Medical Sciences, State University of Campinas, Campinas, Brazil. 12 Departamento de Genética, Facultad de Biología, Universidad de Sevilla, Seville, Spain. 13 Department of Microbiology, Immunology and Transplantation, Katholieke Universiteit Leuven, Leuven, Belgium. 14 Department of Laboratory Medicine and National Reference Centre for Mycosis, University Hospitals Leuven, Leuven, Belgium. 15 Life and Health Sciences Research Institute (ICVS), School of Medicine, University of Minho, 4715-495 Braga, Portugal; ICVS/3B’s-PT Government Associate Laboratory, Braga, Portugal. 16 Guangdong Laboratory for Lingnan Modern Agriculture, Guangdong Province Key Laboratory of Microbial Signals and Disease Control, Integrative Microbiology Research Centre, South China Agricultural University, Guangzhou, China. 17 National Institute of Science and Technology in Human Pathogenic, Fungi, Brazil. 18 Present address: Ginkgo Bioworks, 27 Drydock Avenue, 8th Floor, Boston, USA. 19 Present address: Clinical Microbiology Unit. Synlab Laboratory at Viamed Sta. Ángela de la Cruz Hospital, Seville, Spain. e-mail: ggoldm[email protected]r;[email protected] Article https://doi.org/10.1038/s41467-024-52639-1 Nature Communications | (2024) 15:8412 16