
==== Front
Nucleic Acids Res
Nucleic Acids Res
nar
Nucleic Acids Research
0305-1048
1362-4962
Oxford University Press

39016185
10.1093/nar/gkae625
gkae625
AcademicSubjects/SCI00010
Narese/24
Methods
CLOCI: unveiling cryptic fungal gene clusters with generalized detection
https://orcid.org/0000-0003-3628-8336
Konkel Zachary Department of Plant Pathology, The Ohio State University, Columbus, OH 43210, USA
Center for Applied Plant Sciences, The Ohio State University, Columbus, OH 43210, USA

Kubatko Laura Department of Ecology and Organismal Biology, The Ohio State University, Columbus, OH 43210, USA
Department of Statistics, The Ohio State University, Columbus, OH 43210, USA

https://orcid.org/0000-0001-6731-3405
Slot Jason C Department of Plant Pathology, The Ohio State University, Columbus, OH 43210, USA
Center for Applied Plant Sciences, The Ohio State University, Columbus, OH 43210, USA

To whom correspondence should be addressed. Tel: +1 614 688 2122; Fax: +1 614 292 4455; Email: slot.1@osu.edu
09 9 2024
17 7 2024
17 7 2024
52 16 e75e75
10 7 2024
01 7 2024
13 11 2023
© The Author(s) 2024. Published by Oxford University Press on behalf of Nucleic Acids Research.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site-for further information please contact journals.permissions@oup.com.

Abstract

Gene clusters are genomic loci that contain multiple genes that are functionally and genetically linked. Gene clusters collectively encode diverse functions, including small molecule biosynthesis, nutrient assimilation, metabolite degradation, and production of proteins essential for growth and development. Identifying gene clusters is a powerful tool for small molecule discovery and provides insight into the ecology and evolution of organisms. Current detection algorithms focus on canonical ‘core’ biosynthetic functions many gene clusters encode, while overlooking uncommon or unknown cluster classes. These overlooked clusters are a potential source of novel natural products and comprise an untold portion of overall gene cluster repertoires. Unbiased, function-agnostic detection algorithms therefore provide an opportunity to reveal novel classes of gene clusters and more precisely define genome organization. We present CLOCI (Co-occurrence Locus and Orthologous Cluster Identifier), an algorithm that identifies gene clusters using multiple proxies of selection for coordinated gene evolution. Our approach generalizes gene cluster detection and gene cluster family circumscription, improves detection of multiple known functional classes, and unveils non-canonical gene clusters. CLOCI is suitable for genome-enabled small molecule mining, and presents an easily tunable approach for delineating gene cluster families and homologous loci.

Graphical Abstract

Graphical Abstract

National Science Foundation 10.13039/100000001 DEB-1638999 Ohio State University 10.13039/100006928
==== Body
pmcIntroduction

Gene clusters are chromosomal loci that contain two or more adjacent genes that are co-inherited with cooperative biological functions. Gene clusters are found across the tree of life and encode diverse metabolic and ecological phenotypes (1–6). Perhaps the most well-studied category of gene clusters is biosynthetic gene clusters, which produce small molecule ‘secondary/specialized metabolites’ that mediate biotic and abiotic interactions. Specialized metabolites have diverse activities with applications in industry and agriculture (7) and a prominent source of new drugs in healthcare (8). These compounds include UV-absorbing carotenoids (9), antimicrobial compounds (10,11) and siderophores (12,13). Beyond specialized metabolite biosynthesis, gene clusters can encode pathways that degrade antagonistic compounds (14), catabolize carbon sources and amino acids (15–17), assimilate nutrients (18), or produce essential vitamins (19) and proteins (20–22). Gene clusters can be thought to span a continuum from generalized to specialized based on how broadly applicable the cluster is to organisms’ niches. Generalized metabolic cluster repertoires present a window into essential metabolic processes (23,24) and widely conserved ecological functions (19,25–27), while specialized cluster repertoires provide insight into ecological niche (28–30,14).

Gene cluster detection is a transformative tool for specialized metabolite drug discovery in the genomic era. Fungal drug discovery rapidly introduced novel drugs in the 20th century by implementing laborious untargeted chemical screening, though the field is increasingly hampered by redundant compound identification due to a focus on easily cultivable biodiversity (31). Additionally, some specialized metabolites are not produced under standard laboratory conditions and myriad expression conditions are often tested to induce silent pathways (32). A reservoir of novel specialized metabolites thus remains hidden in organisms and biosynthetic pathways that are elusive under laboratory conditions (33).

Genomic gene cluster detection unveils specialized metabolite biosynthesis pathways from difficult-to-culture organisms and repressed biosynthetic pathways. Gene cluster detection algorithms have unveiled biosynthetic pathways with novel specialized metabolites (34) and identified novel specialized metabolic pathways that are epigenetically repressed (35). Once gene clusters are identified, it is then possible to predict gene functions and systematically characterize them by heterologous expression (36).

The most widely implemented biosynthetic gene cluster detection algorithms are function-centric (1,37). Function-centric algorithms search for profile models of ‘core’ biosynthetic proteins commonly associated with described biosynthetic gene clusters of bacteria and filamentous Ascomycota (Fungi). Greater than 94% of the fungal clusters in the Minimum Information about Biosynthetic Gene Clusters (MIBiG) database (38) originate from the filamentous fungi in the phylum Ascomycota, depicting bias in gene cluster identification. Furthermore, profile models are largely constructed from a set of ‘canonical’ core genes derived from polyketide, nonribosomal peptide, and terpene biosynthetic clusters (39). Function-centric detection is thus expected to identify clusters derived from canonical functions, though it may have poor precision in lineages with ‘non-canonical’ specialized metabolite gene clusters that lack core genes upon which the models rely (33,40–44). To account for non-canonical classes, the function-centric software antiSMASH incrementally expands its scope by accumulating models of novel biosynthetic gene cluster core functions (1). Each update increases the breadth of future searches to include new biosynthetic classes, though the function-centric paradigm cannot infer novel non-canonical gene cluster classes de novo.

Some gene cluster detection algorithms are function-agnostic, relying instead on using general properties of gene clusters to detect them. Two function-agnostic methods are direct analysis of gene co-expression (44,45) and modeling coordinated expression through identification of shared intergenic motifs (46). Co-expression detection predicts gene clusters from transcriptomic data by identifying significantly coexpressed genes within the same locus. In practice, unsupervised coexpression analysis can be costly and limited to metabolic pathways expressed in the recommended minimum of 15–20 expression conditions required for unsupervised analysis (47). Alternatively, the CASSIS algorithm predicts clusters from genomic data by inferring shared intergenic motifs within loci (43,46,48). Function-agnostic approaches can infer non-canonical gene cluster classes de novo, although CASSIS is primarily implemented in conjunction with function-centric detection to infer cluster boundaries (49).

Another function-agnostic approach identifies gene clusters in genomic data by inferring selection on gene order (28,14,50–54). Gene clusters comprise cooperative genes that are tied to discrete phenotypes, and the clustered state is maintained by selection (55–57). Gene clusters are evidence of selection because it is improbable that multiple distantly related genes convergently co-locate in diverse lineages and endure microsynteny decay over time (14,55). One proxy that identifies selection on gene order is unexpectedly shared microsynteny (synteny) (14), which detects selection for gene colocalization by identifying combinations of gene families that co-occur with one another in an unexpectedly large phylogenetic distribution. Unexpectedly shared synteny is the basis of the algorithms CO-OCCUR and EvolClust (28,52) and can identify gene clusters with large or sparse phylogenetic distributions (11,23,58). However, broadly distributed gene clusters are sometimes embedded in larger regions that also meet unexpectedness thresholds because strong linkage disequilibrium within the cluster may extend to the surrounding regions (59). Additionally, if thresholds of unexpectedness are too low, then shared synteny loci that are not gene clusters, such as regions surrounding the centromere, will erroneously be reported as gene clusters (60). To separate gene clusters from loci that are not clusters, unexpectedly shared synteny algorithms have had to sacrifice either their scope of predicted cluster categories, or restrict the sizes of clusters they can predict and/or their capacity to infer recently evolved clusters (58). These limitations are in part why these algorithms have not been as widely implemented as antiSMASH.

We created CLOCI to test two hypotheses: (i) relaxing thresholds for unexpectedly shared synteny will improve cluster recall; (ii) measuring additional proxies of selection for coordinated gene evolution will increase sensitivity and accuracy. CLOCI lowers unexpectedness thresholds by accurately defining shared synteny locus boundaries and phylogenetic distributions. CLOCI enables separating gene cluster families from unexpectedly shared synteny loci that are not gene clusters using three additional proxies for selection for coordinated gene evolution. This approach yields the greatest recovery and categorical diversity of gene clusters when compared to antiSMASH and EvolClust. The CLOCI detection model has the additional advantage of accurate gene cluster boundary inference, and predicts both novel putative clusters and non-canonical biosynthetic gene clusters.

Materials and methods

CLOCI infers gene cluster families (GCFs) by identifying homologous, unexpectedly shared synteny loci and enriches gene clusters using a model of phylogenetic and alignment-based proxies of gene coinheritance and coevolution. The CLOCI algorithm first extracts loci with more shared microsynteny than expected (Figure 1A–I). It then classifies homologous locus groups (HLGs) by adapting graph clustering methodology from orthologous gene group (orthogroup) classification (61,62) (Figure 1J). CLOCI enriches HLGs for GCFs using proxies of coordinated gene evolution, including: (i) the unexpectedness of the phylogenetic distribution of an HLG; (ii) commitment of genes to the HLG locus; (iii) conservative amino acid substitution bias within an HLG and (iv) the sparseness of the HLG distribution across the species phylogeny (Figure 1K). We implemented CLOCI across a dataset of 2247 fungal genomes assembled by Mycotools (63) (December 2022) and benchmarked metabolic gene cluster recovery against the function-centric software, antiSMASH and the unexpectedly shared synteny algorithm, EvolClust, using a dataset of 57 biosynthetic and 11 non-biosynthetic reference clusters that include catabolic and nutrient assimilation clusters. We benchmark the accuracy of CLOCI boundary detection against the same software by referencing a dataset of 33 well-characterized gene cluster boundaries.

Figure 1. CLOCI pipeline depicting cluster detection, including a hypothetical cluster found in genomes I-III, but not genome IV. GCL and MMI in k. are truncated examples using two shared genes from the full homologous locus group (HLG).

Seeding cluster detection by identifying unexpectedly shared microsynteny

Significantly co-occurring pairs of homologous genes (HG) are the seeds for CLOCI gene cluster detection. HGs are inferred from a reference Mycotools database (63) of complete publicly available genomes using mmseqs2 cluster (64) (Figure 1A, B). The contigs of each genome are then represented as a vector of HGs. All pairs of HGs that co-occur within a five-gene sliding window make up the set of co-occurring HG pairs in a genome (Figure 1C). We default to a five-gene sliding window based on determining which window sizes from 3–9 genes recover and most accurately infer the boundaries of the ergot alkaloid, quinic acid utilization, galactose metabolism, biotin biosynthesis, and sterigmatocystin gene clusters (Supplementary Tables S4 and S5).

In order to quantify the distribution of each HG pair, CLOCI constructs a microsynteny tree of the inputted genomes. The microsynteny tree depicts the gene order similarity among genomes, and thus accounts for shared evolution when quantifying HG pair distribution. To construct a microsynteny tree, CLOCI builds a presence-absence matrix of all HG pairs that include an HG that is a near single-copy ortholog (Figure 1D). We represent microsynteny decay using gene neighborhoods with near single-copy orthologs because these genes are nearly universally conserved and exist in low copy numbers, which limits erroneous identification of orthologous loci. The set of near single-copy ortholog HGs can be determined by BUSCO (65), UFCG (66) or another external method. We used nine UFCG orthologs and five orthologs from a previous study (67) (Supplementary Table S2). Alternatively, CLOCI will attempt to determine a representative set of near single-copy HGs present in all genomes in the lowest median and mean count per genome.

The presence-absence matrix is used to reconstruct a maximum likelihood microsynteny tree (68) via db2microsyntree.py from the Mycotools software suite (63). The microsynteny tree is built using IQ-TREE 2 and a free-rate, general time reversible evolutionary model with ascertainment bias correction (69–71) (Figure 1D). To prevent topological discrepancies with accepted species topology (66,72), the microsynteny topology can be constrained with a species tree or in the absence of an external constraint a consensus tree constructed from 1000 ultrafast bootstrap replicates is generated. We constrained the microsynteny topology with a phylogenomic tree constructed from the same near single-copy orthologs used to detect gene neighborhoods for microsynteny tree construction (Supplementary Table S3). The tree is then rooted on a specified outgroup, and the distribution of each HG pair is calculated as the total microsynteny branch distance (TMD) between the genomes where it is found (Supplementary Equation S1).

To identify unexpectedly distributed HG pairs, CLOCI creates a null distribution of background TMD from randomly sampled HG pairs (Figure 1E). Different lineages have different rates of microsynteny decay, so CLOCI builds local null models for specific lineages or taxonomic ranks. By default, CLOCI builds null models for each genus because microsynteny decay often does not significantly vary within genera (73), and it strikes a balance with computational throughput. CLOCI accounts for species with overrepresented genome samples during null model construction by randomly choosing a species prior to randomly sampling from a genome. The upper 20th percentile from each null distribution is used to define the TMD value above which an HG pair is considered unexpectedly shared (Figure 1f). We selected the 20th percentile to decrease downstream pairwise comparisons while presumably identifying most HG pairs that represent gene clusters.

HG pairs are used as seeds for inferring higher-order HG overlap (Figure 1G). Combinations of three or more co-occurring HGs (HG combination) are identified by pairwise comparison of all loci that correspond to an HG pair. Null distributions are generated for each observed HG combination size, up to the sliding window, as described for HG pairs (Figure 1H). Unexpectedly syntenic HG combinations are extracted from the upper 60th percentile of null TMDs (Figure 1I). We chose the 60th percentile as the threshold below what is necessary to detect the recently evolved psilocybin reference cluster (40,42).

Identifying locus boundaries by assembling shared synteny loci

Shared synteny locus boundaries are predicted by first extracting loci from unexpectedly shared HG combinations (Figure 1J). All loci that correspond to unexpectedly shared HG combinations are extracted. Then any overlapping genes between these loci are added to the locus that contains the HG combination with the greatest TMD. Single gene loci generated by the merging process are fused with an adjacent locus with the most similar phylogenetic distribution (Supplementary Equation S2) if it exceeds a 25% minimum phylogenetic similarity. If a single gene is not merged, then it is aggregated with surrounding singletons, or discarded if no surrounding singletons exist. Spurious loci generated from singleton merging are discarded downstream if they lack homology to other loci.

CLOCI finalizes shared synteny locus boundary inference by identifying groups of similar domains of shared microsynteny and merging adjacent loci with homologous domains. To infer groups of similar domains, CLOCI implements Markov clustering (MCL) (74) on an adjacency graph of locus-locus similarity. CLOCI quantifies the pairwise similarity of all loci following their extraction. We define similarity between two loci as the average amino acid identity of their shared HGs scaled by the overlap coefficient of the shared HG combinations (Supplementary Equation S3). Compared loci must have at least two shared HGs, and if a locus contains multiple gene sequences from a particular HG then the highest identity comparison is incorporated into the average. Locus–locus similarity scores greater than 35% are represented in an adjacency graph, and initial domains are inferred via this graph with inflation value set to 1.1. CLOCI merges extracted loci that are within two genes of one another and the loci belong to the same domain.

Grouping locus homologs by iterative graph clustering

To identify groups of similar loci and families of related gene clusters, by extension, CLOCI groups the finalized extracted loci into homologous locus groups (HLGs). HLGs are inferred using the same locus-locus similarity and graph clustering algorithm as used for detecting initial domains, though in this round CLOCI implements the Sørensen-Dice similarity as the scaling coefficient (Supplementary Equation S3) to penalize against missing HGs in the similarity calculation. HLGs are inferred from the resulting locus-locus adjacency graph with a default MCL inflation value of 1.3. Genes can belong to one HLG, and reported loci can thus comprise directly abutting HLGs. In reality, closely related genes can be binned into separate HLGs, which indicates shared homology between HLGs. The relationships between loci and HLGs can be visualized via a locus-locus similarity network generated from the included hlg2hlg_net.py script.

CLOCI next attempts to repair incomplete or fragmented loci within inferred HLGs. First, CLOCI attempts to complete partial cluster loci by reference to related loci. Each locus within each HLG is extended to include any HG within one gene up/downstream of the locus that is shared across the final HLG, unless the flanking gene is already part of a different HLG. The HG is later removed if it has <35% amino acid identity with all homologous sequences in the HLG. To identify loci that are fragmented across the genome assembly, all contigs are searched for shared HG pairs that are missing from a given locus. CLOCI only searches for two gene combinations because combinations greater than three are already retrieved by identifying unexpectedly syntenic HG combinations.

Characterizing homologous locus groups according to proxies of coordinated gene evolution

Following HLG circumscription, CLOCI characterizes HLGs using proxies of coordinated gene evolution in addition to TMD, including gene commitment to the locus (GCL), locus-locus amino acid similarity, conservative amino acid substitution bias (CSB), and the Phylogenetic Distribution Sparsity (PDS) of HLGs (Figure 1K). To visualize the distribution of HLG proxies and relative contribution of each HLG proxy to this distribution, we reduced the proxy dimensionality using principal component ordination implemented in the scikit-learn package (75). We tested for autocorrelation between proxies using ordinary least squares linear regression (Supplementary Figure S1).

Gene commitment to the locus

Gene Commitment to the Locus (GCL) is a measure of the extent to which an HG has remained in an HLG through time. This is a proxy for the coordinated inheritance of genes in the HLG. The GCL for an HLG is the weighted average of individual HG commitment scores for each HG found in at least two homologous loci of the HLG (Supplementary Equation S4). The HG commitment score is equal to the percent of alignment hits that are sequences within the HLG after discarding paralogs outside the HLG locus in genomes that have it. Sequences are aligned against the entire HG using DIAMOND BLASTp. Hits are determined by sorting the alignments from highest to lowest identity, and retaining sequences up until all homologous sequences within the HLG loci are recovered. Only the maximum scoring sequence is considered if a genome has multiple sequences in the same HG within an HLG.

Conservative substitution bias

Conservative Substitution Bias (CSB) is a proxy for the action of purifying selection on the protein structure that is more computationally tractable than dN/dS. CSB measures the ratio of positive amino acid substitutions to identical amino acids. CSB for an HLG is calculated by dividing the mean minimum observed percent amino acid identity (MMI) by the mean minimum observed percent conservative alignment positions (positives, MMP) across all pairwise comparisons of sequences in an HG. Minimum values approximate the maximum distance between sequences in the HLG. To calculate the minimum identities and positives, all sequences in HGs in two or more homologous loci in the HLG are aligned against the entire HG via DIAMOND BLASTp. The minimum observed positive score for all pairwise alignments is divided by the minimum identity score. The minimum identity and positive scores are set to 0 for alignments that are missing any sequences in the HLG. The HLG CSB is then derived from the ratio of HG minimum identities to minimum positives scaled by the proportion of considered genes in each HG (Supplementary Equation S5).

Phylogenetic distribution Sparsity

Phylogenetic Distribution Sparsity (PDS) measures the prevalence of taxa with a particular HLG in its overall phylogenetic range. PDS is a proxy for the extent of horizontal cluster transfer and/or coordinated cluster gene loss, which are elevated in some categories of gene clusters (25,42,55,58,76,77). CLOCI quantifies PDS for each HLG as the TMD of the HLG divided by the TMD of the most recent common ancestor of the genomes that have the HLG (Supplementary Equation S6). Thus, PDS quantifies the percent of branch length descended from the most recent common ancestor TMD that contains the HLG.

Implementing CLOCI using a dataset of publicly available fungal genomes

We implemented CLOCI across a database of 2247 fungal genomes across the kingdom using default CLOCI parameters (Supplementary Figure S2, Supplementary Table S1) using Mycotools. We additionally constrained the microsynteny tree topology to 14 near single-copy orthologs from Spatafora et al. and Universal Fungal Core Genes (66) and rooted on Rozella spp (Supplementary Table S2). We created GCF networks locus-locus similarity networks using the hlg2hlg_net.py tool packaged with CLOCI. All computation was performed using 26 Quad Core Intel Xeon 6148 Skylake processors with access to approximately 1 terabyte RAM.

Benchmarking CLOCI recovery and boundary accuracy against antiSMASH and EvolClust

We benchmarked detection power with reference to a dataset of known clusters (Supplementary Table S4, Supplementary Figure S3) and the precision and accuracy of cluster boundaries with reference to an independent dataset that primarily contains well-characterized Aspergillus spp. gene clusters (Supplementary Table S5). We cannot directly test the false discovery rate and algorithm precision because gene clusters are not exhaustively sampled in any analyzed genomes. We compared CLOCI reference cluster recovery and cluster boundary precision with those of antiSMASH (v6.1.1) and EvolClust (1,78). We implemented antiSMASH on the genomes containing known clusters and queried the EvolClustDB (accessed 4 April 2023). If the EvolClustDB genome was from a different strain than the query cluster, we first verified that the cluster was present by querying an EvolClustDB genome of the same species using BLASTp. We evaluated the recovery of 57 unique biosynthetic clusters and 11 non-biosynthetic clusters. Biosynthetic clusters were manually selected from MiBIG to maximize diversity of known biosynthetic classes. We included the psilocybin and biotin clusters to increase non-canonical biosynthetic cluster representation (19,42,79). The non-biosynthetic reference clusters were compiled from an extensive search of available literature (23). We tested recovery by aligning a core biosynthetic gene or gene near the center of each reference cluster to the same genome as the reference. Cluster loci were deemed recovered if the top BLASTp hit to the reference gene had a minimum percent identity 90%. If the reference genome was not present in our dataset, we confirmed the presence of the reference cluster in a genome of a species known to have the cluster by identifying two top BLASTp/tBLASTn hits within the same locus with greater than 60% identity to two reference cluster genes. Four clusters were removed from the antiSMASH recovery analysis because the genome annotation format was incompatible with antiSMASH. 19 clusters were removed from the EvolClust recovery calculation because the genome containing the query cluster was not in EvolClustDB.

We assessed cluster boundary accuracy referencing a separate dataset of 33 gene clusters, 25 of which are derived from the CASSIS Aspergillus boundary accuracy dataset (46), because most clusters do not have well characterized boundaries (Supplementary Table S5). We additionally included psilocybin (40) and biotin (19) non-canonical biosynthetic gene clusters, three non-biosynthetic gene clusters (16,25,80), the Penicillium chrysogenum penicillin cluster (11), and two ergot alkaloid biosynthetic gene clusters (30) to broaden the scope of gene cluster categories in boundary detection. We assessed accuracy by the percent missing and percent extraneous genes per cluster. Clusters were submitted to boundary evaluation if at least one gene sequence within the boundaries was identified using BLASTp. All reported loci that contain at least one of the genes within the reference cluster boundaries are considered in boundary accuracy.

Evolutionary analysis of nitrate assimilation gene clusters

We examined the evolution of nitrate assimilation gene clusters by constructing gene phylogenies (81) for cluster genes using the Cluster Reconstruction and Phylogenetic Analysis pipeline packaged in Mycotools (63) with default parameters. We extracted nodes with high support (>0.99 fasttree bootstrap) and generated subsequent phylogenies by aligning with mafft –auto v7.487 (82), trimming using ClipKIT v1.3.0 (83), and constructing phylogenies using IQ-TREE v2.2.0.3 with 1000 ultrafast bootstrap replicates (69,71). To evaluate the alternative hypothesis of vertical evolution compared to horizontal transfer, we compared the likelihood of a Mucoromycota monophyletic constraint using an approximately unbiased test (84) with 10 000 boostrap replicates implemented in IQ-TREE v2.2.0.3 (71) (Supplementary Table S7).

Filtering gene cluster families from homologous locus groups according to proxies of coordinated gene evolution

Gene cluster families (GCFs) are filtered from HLGs by setting thresholds for the minimum proxy values (Figure 1l). In order to evaluate the capacity for coordinated gene evolution proxies to enrich for metabolic gene clusters, we analyzed the effect of variable coordinated gene evolution proxy thresholds on the proportion of metabolic process (GO:0008152) and secondary metabolic process (GO:0019748) gene ontology (GO) terms in the resulting GCFs. To obtain GO terms, we first queried the genes against the Pfam database using hmmsearch v3.3.2 with 50% minimum alignment coverage and 0.001 minimum e-value. Each gene was annotated with its lowest e-value Pfam hit and GO terms were assigned referencing the pfam2go database (current.geneontology.org/ontology/external2go/pfam2go). We incrementally adjusted TMD, GCL, CSB and PDS thresholds by 0.2 from zero to one and quantified the proportion of retained genes annotated as metabolic process or secondary metabolic process. The TMD calculation in this filtering model acts on the distribution of an HLG rather than individual HG pairs/combinations as it does in unexpectedly shared microsynteny detection.

Results

CLOCI recovers functionally diverse families of gene clusters. We recovered all reference gene cluster categories, including both canonical biosynthetic (aflatoxin) and non-canonical biosynthetic (psilocybin), catabolic processes (quinic acid, galactose and proline degradation), and nutrient assimilation (nitrate assimilation) clusters. These gene clusters were recovered from an output of 130 931 homologous locus groups (HLGs) that comprise 334 6973 unexpectedly shared synteny locus domains with a mean size of 4.70 genes per domain (62.2% of overall genes), median size 4 genes, maximum size 53 genes, and minimum size of two genes (Supplementary Data). There are 25.56 locus domains per HLG, and most reference boundary clusters are represented by a single locus domain, though we recovered some clusters as two or three domains. The CLOCI gene cluster detection algorithm has higher cluster recovery than the function-centric algorithm, antiSMASH, and the unexpected shared synteny algorithm, EvolClust. CLOCI also better approximates cluster boundaries when considering all shared synteny locus domains in reference clusters. Filtering using the minimum observed proxy values for reference clusters removes 41 337 HLGs (31.6%). This filtered dataset comprises 89 776 gene cluster families (GCFs), with 33.8 cluster domains per GCF, mean size 4.78 genes per cluster domain (57.2% of overall genes), median size 4 genes, maximum size 48 genes, and minimum size 2 genes. Raising threshold proxy values variably enriches for genes with metabolic process and specialized metabolite process GO terms in the final set of GCFs.

CLOCI gene cluster families recapitulate known gene cluster distributions

To assess gene cluster family (GCF) completeness, we compared the reported phylogenetic distribution of gene clusters to their respective CLOCI GCFs. We recovered known gene cluster distribution in the main categories of gene clusters. We recovered the ergot alkaloid biosynthetic GCF in all Onygenales, Eurotiales, Xylariales, Hypocreales and Helotiales genomes that we previously identified through phylogenetic analysis (70). We identified all known homologs of the non-canonical biosynthesis gene cluster for psilocybin production in our dataset in a single GCF (Figure 1A). Additionally, both parts of the Panaeolus cyanescens psilocybin cluster, which is distributed across two contigs due to a fragmented genome assembly, were recovered in this GCF. We also identified unreported homologs of the clusters that produce the immunosuppressant nonribosomal peptide, cyclosporin, in Dactylonectria estremocensis, and the fungicidal strobilurin polyketides in Mycena spp. (Figure 4C).

Multiple nitrate assimilation GCFs were identified from diverged taxa (25). We recovered at least five GCFs that contain nitrate transporter, nitrate reductase, and nitrite reductase genes associated with nitrate assimilation. These GCFs span three phyla and the clusters are often found directly adjacent to unexpectedly shared synteny regions. GCF #0 consists entirely of Ascomycota clusters, including the characterized Aspergillus nidulans cluster (18). GCF #1 comprises 81.3% Basidiomycota clusters and includes sequences implicated in horizontal transfer between ancestral Ustilagomycotina (Basidiomycota) and Hypocreales (Ascomycota) species (65,69). GCF #2 contains the five-gene Saccharomycotina (Ascomycota) yeast nitrate assimilation gene cluster reported from Ogataea polymorpha (85), as well as unreported homologous clusters in Bifiguratus adelaidae (Mucoromycota). GCF #3 primarily contains Pseudogymnoascus spp. (Ascomycota) (86). Additionally, an unreported putative nitrate assimilation GCF #4 was identified in Basidiomycota and Mucoromycota (Figure 4F), consisting of a nitrate transporter, nitrate reductase, and nitrite reductase. All nitrate assimilation GCFs are connected in the locus-locus similarity network, indicating homology (Figure 4F). Phylogenetic reconstruction of the nitrate transporter, nitrate reductase, and nitrite reductase genes in GCF #4 using crap (63) (Supplementary Table S6) reveals that ectomycorrhizal Amanita spp. (Basidiomycota) are nested within a clade of Mucoromycota spp. with 100% IQ-TREE ultrafast bootstrap support for all three genes (71). Constrained trees that force a Mucoromycota monophyly were rejected via approximately unbiased testing (Supplementary Table S7). Amanita muscaria spp. have close homologs of these genes but lack an intact cluster and saprobic A. thiersii and A. inopinata both lack homologs of the genes (87).

Benchmarking CLOCI against antiSMASH and EvolClust

We benchmarked CLOCI cluster recovery against antiSMASH and EvolClust referencing a dataset of 11 non-biosynthetic gene clusters from literature and 57 biosynthetic clusters, including two non-canonical biosynthetic clusters (Supplementary Table S4). CLOCI recovers 86.0% and 100% of the biosynthetic and non-biosynthetic dataset respectively (Figure 2). In comparison, antiSMASH recovers 76.1% and 0% of the biosynthetic and non-biosynthetic datasets respectively, whereas EvolClustDB recovers 48.7% and 60% respectively. Additionally, CLOCI detects the biotin and psilocybin non-canonical biosynthetic clusters, whereas EvolClust and antiSMASH are missing both. Of the missing clusters, CLOCI did not detect the fumiquinazoline, zearalenone, mycophenolic acid, cephalosporin, ochratoxin A, or the ent-kauren-16alpha-ol clusters.

Figure 2. Comparison of cluster recall and cluster boundary accuracy by CLOCI v0.0b, antiSMASH v6.1.1, CASSIS (implemented in antiSMASH v6.1.1), and EvolClustDB (accessed 2023/04/04). (A) Recovery of 57 biosynthetic gene clusters from (37) and literature, and 11 non-biosynthetic gene clusters from (23). (B) Boundary detection accuracy and precision referencing an independent dataset of 33 characterized cluster boundaries from (45) and literature. Whiskers depict the 1st and 3rd quartiles.

We assessed cluster boundary inference by referencing an independent dataset of well-characterized cluster boundaries (Supplementary Table S5). We determined that CLOCI, antiSMASH and antiSMASH with CASSIS have the lowest median (0) missing genes per recovered cluster, whereas CLOCI has the lowest median (0) extraneous genes per recovered cluster (Figure 2B). 82.1% of reference clusters were predicted as a single domain by CLOCI, 14.3% were reported as two, and one cluster was reported as three domains, with a mean of 6.27 genes per domain and 1.21 domains per reference cluster. For example, aflatoxin was reported as two domains, whereas sterigmatocystin was reported as three domains, one of which is circumscribed in the same GCF as aflatoxin. The pseurotin and fumagillin clusters are directly adjacent and CLOCI reports these clusters as two discrete loci that comprise one and two domains respectively (88). All clusters were recovered as a single domain by antiSMASH, though antiSMASH with CASSIS enabled can distinguish between multiple clusters in a reported region, including the pseurotin and fumagillin clusters. AntiSMASH predicted domains with an average size of 19.04 genes whereas antiSMASH with CASSIS enabled predicted an average of 18.11 genes per domain. AntiSMASH recovers the aflatoxin cluster in the boundary dataset, but did not recover the query gene for the recovery dataset because the recovery dataset cluster references a fragmented assembly. The aflatoxin and cichorine clusters were recovered as two domains by EvolClust, which predicts reference clusters with an average of 13.05 genes per domain and 1.10 domains per cluster

Filtering using proxies of coordinate gene evolution reveals novel gene clusters

CLOCI filters GCFs from shared synteny HLGs using multiple proxies of selection for coordinated gene evolution (Table 1, Figure 3A, B). Proxies have variable power to filter the dataset when their observed minimum values in the reference dataset are used independently. The minimum Total Microsynteny Distance (TMD) reduces HLGs by 23.9%, Gene Commitment to the Locus (GCL) by 7.87%, and Conservative Substitution Bias (CSB) by 8.69%. Phylogenetic Distribution Sparsity (PDS) is 0% for luciferin. However, overall, PDS tends to be greater than 0% (median 93.6%), and PDS also explains most of the variance in the distribution of both principal components following two dimensional ordination of HLG proxies (Table 1).

Table 1. Values of coordinated gene evolution proxies for homologous locus groups (HLGs) that contain known gene clusters from the recovery dataset

Proxy	Minimum	Mean	Median	Standard deviation	Principal component 1	Principal component 2	
Total Microsynteny Distance (TMD)	0.434	0.664	0.683	0.133	0.0531	0.0592	
Gene Commitment to the Locus (GCL)	20.5%	58.0%	55.3%	21.4%	–0.241	–0.0801	
Mean Minimum Identity (MMI)	23.2%	54.3%	53.9%	17.6%	–0.1035	–0.1096	
Mean Minimum Positives (MMP)	30.2%	65.2%	65.5%	16.2%	–0.0803	–0.920	
Conservative Substitution Bias (CSB)	1.76%	17.9%	19.4%	6.57%	0.0422	0.0372	
Phylogenetic Distribution Sparsity (PDS)	0.00%	80.3%	93.6%	26.3%	0.365	–0.117	

Figure 3. Reference gene cluster families (GCFs) are characterized by and extracted from shared synteny homologous locus groups (HLGs) using proxies of coordinated gene evolution: (A) bar graphs of proxies of coordinated gene evolution for 63 recovered clusters. Nitrate assimilation is represented by the GCF containing the Ustilago bromivora cluster. TMD = Total Microsynteny Distance, GCL = Gene Commitment to the Locus, MMI/MMP = Mean Minimum Identity/Positives of amino acids, CSB = Conservative Substitution Bias, and PDS = Phylogenetic Distribution Sparsity. Red font indicates non-biosynthetic clusters, green canonical biosynthetic, and blue non-canonical biosynthetic. (B) Principal component analysis of HLGs highlighting HLGs with reference clusters in red. Loading (feature) vectors correspond to the proxy with the same color as bars in a); Enrichment of (C) metabolic processes (GO:00008152) and (D) secondary metabolic processes (GO:00019748) in GCFs filtered from HLGs with greater than a minimum specified proxy threshold relative to initial HLGs. Significantly enriched data points are in red.

Reference specialized metabolic clusters have the lowest TMDs, while clusters that encode vitamin biosynthesis, catabolic, and cellular protective functions have the highest (Figure 3A). The four clusters with the smallest TMD are the secondary/specialized metabolite clusters omphalotins, tenuazonic acid, cyclosporin, and luciferin. The top four TMD clusters encode valine catabolism, tyrosine degradation/pyomelanin biosynthesis, quinic acid catabolism, and biotin biosynthesis. Biotin is an essential cofactor (19), while the quinic acid cluster catabolizes a carbon source (80,89), the tyrosine degradation/pyomelanin biosynthesis cluster degrade phenolics and produce melanins that may protect fungal cells from the environment (17,90).

TMD is weakly to strongly linearly correlated with other proxies using an arbitrary adjusted R2 cutoff of 0.5, while the other proxies are weakly correlated with one another (Supplementary Figure S1). Log normalized TMD is correlated with GCL (adjusted R2 = 0.816), PDS (adjusted R2 = 0.641), and CSB (adjusted R2 = 0.861). The correlation between GCL and PDS is weak (adjusted R2 = 0.370), while GCL and CSB are correlated (adjusted R2 = 0.614). PDS is also correlated with CSB (adjusted R2 = 0.669).

Filtering using minimum threshold values for proxies of selection for coordinate gene evolution differentially affects enrichment of metabolic process and secondary metabolic process gene ontology (GO) terms. We determined that metabolic process genes are significantly and increasingly enriched as the TMD threshold is raised (Figure 3C). Conversely, increasing TMD reduces secondary metabolic process genes (Figure 3d). GCL reduces metabolic process genes across the thresholding range. However, GCL enriches secondary metabolic process genes from 0.2 to 0.8 with a peak fold change at 0.6 GCL. On its own, raising the CSB threshold to 0.2 enriches metabolic process genes, though larger thresholds lead to diminishment. CSB significantly enriches secondary metabolic process genes across its range, though the effect is small. PDS enriches both metabolic process and secondary metabolic process genes across the range of thresholds, with a peak metabolic process enrichment at 0.2. At 0.8, PDS generates a 209% increase in the proportion of secondary metabolic process genes.

CLOCI putatively identifies an unreported widely shared metabolic cluster and sparsely distributed specialized metabolite cluster. We putatively identified a widely shared iron sequestration gene cluster with an 0.83 TMD that is log normalized relative to all HLG TMDs. This gene cluster contains a ‘widely conserved’ synthetase that produces the Basidiomycete siderophore, basidioferrin (12). The CLOCI reported cluster is relatively large, comprising 23 predicted genes in Psilocybe cyanescens (Figure 4E). The TMD of the basidioferrin GCF is the 9th greatest when compared to reference clusters (Figure 3A), which is congruent with the broad distribution of basidioferrin synthetase homologs (12). We also identified a putative non-ribosomal peptide specialized metabolite GCF (#101062 in the upper 90 percentile of PDS scores. This putative GCF was recovered in four genomes within three taxonomic classes.

Figure 4. Locus-locus similarity graphs and representative synteny diagrams of CLOCI GCFs of multiple gene cluster categories. Nodes represent loci and edge weight/distance depicts log normalized locus-locus similarity (Supplementary Equation S3). Nodes depicted in synteny diagrams in a-e are salmon-colored and ‘*’ indicates clusters recovered from contig edges. Gene arrows are color coded by homology, and grey arrows indicate no inferred homology. (A) The non-canonical psilocybin biosynthetic gene cluster in all reported genomes in the dataset, including across two contigs in Panaeolus cyanescens. (B) An unreported cluster containing a nonribosomal peptide synthetase found in three different taxonomic classes. (C) Strobilurin clusters including unreported homologs in Mycena. (D) The phenolic catabolism and pyomelanin biosynthesis cluster family. (E) A putative basidioferrin gene cluster family that is widely shared across Agaricomycotina (Basidiomycota) genomes. Basidioferrin synteny diagram is centered on the characterized Gelatoporia subvermispora basidioferrin type VI nonribosomal peptide synthetase. (F) Similarity network of five nitrate assimilation GCFs distributed among three phyla.

Discussion

CLOCI generalizes gene cluster detection by inferring selection on coordinated gene evolution

Detecting selection for coordinated gene evolution is an effective tool for identifying gene clusters (28,14,51). Unexpectedly shared synteny (USS) detection identifies gene clusters from diverse categories by identifying selection for gene colocalization (14). However, USS detection does not discriminate between gene clusters and loci that have USS for other reasons, so in order to identify true gene clusters, previous USS algorithms limit their capacity to infer recent, small clusters or generally infer cluster categories. CLOCI identifies recently evolved clusters while remaining function-agnostic in part by accurately inferring locus boundaries in a size-agnostic framework. Detecting other signatures of selection for coordinated gene evolution can then enrich USS loci for metabolic gene clusters.

The signature of selection for gene colocalization can identify gene clusters from diverse categories, including phenolic and amino acid catabolism, nutrient assimilation, and both non-canonical and canonical biosynthetic gene clusters (28,14). Non-biosynthetic clusters have largely been overlooked by function-centric software because they are primarily designed to detect biosynthetic clusters. Non-canonical biosynthesis has been overlooked because these clusters lack screened core gene models, and often have discrete distributions. USS, employed in CLOCI and other methods (28,51,52) is a proxy for selection on gene colocalization, a driver of the assembly and maintenance of gene clusters. Because this general evolutionary signature of gene clusters is not expected to vary based on cluster function, all categories should be detectable by this method, in contrast with function-centric methods (58). Previous USS algorithms, such as CO-OCCUR, indeed can infer gene clusters from diverse categories (28,14). However, limitations in USS cluster recovery and precision have precluded the widespread adoption of these algorithms in comparison to the function-centric software, antiSMASH.

The most prominent challenge of USS algorithms is discriminating between gene clusters and loci that have shared synteny for other reasons. For example, regions near the centromere have shared synteny across relatively large phylogenetic distances because of a reduced rate of gene rearrangement compared to other genomic regions (60). Therefore, low USS thresholds can falsely designate regions with low rearrangement rates as clusters, while high thresholds can exclude recently evolved gene clusters with narrow distributions. Existing USS algorithms have attempted to address this problem in different ways. For example, EvolClustDB omits clusters smaller than five genes and imposes thresholds that remove recently evolved clusters using a global heuristic model of gene cluster size and shared synteny region similarity (78). CO-OCCUR detects shared synteny loci that contain genes with predetermined functions in order to predict gene clusters (28), which limits its scope to targeted searches.

Recently evolved, small clusters with low USS signal are detected by relaxing USS thresholds and increasing the accuracy of USS locus boundaries and phylogenetic distribution. Small clusters have relatively low USS signal because they need larger distributions to reach unexpectedness thresholds. Recently evolved clusters have small phylogenetic distributions, which also require low USS thresholds to identify. However, lowering USS thresholds may simultaneously increase false positives, so the distributions of shared synteny loci need to accurately recapitulate the distributions of the clusters they represent to minimize how much thresholds need to be lowered. CLOCI groups USS loci into homologous locus groups (HLGs) that recapitulate reference gene cluster distributions, including the relatively restricted distribution of the psilocybin gene cluster (42). We attribute the quality of HLG distributions to accurately defining shared synteny locus boundaries which allows for accurate locus-locus similarity calculations. CLOCI builds cluster boundaries from the most widely distributed combinations of gene families that underlay the cluster. This approach is size-agnostic and detects small clusters such as the three gene penicillin cluster, with the caveat that two gene clusters are not part of initial unexpected synteny detection.

The signature of coordinated gene evolution can also help separate gene clusters from USS loci that are not clusters by increasing detection orthogonality. In CLOCI, we implement four proxies of coordinated gene evolution to improve detection of true gene clusters from USS loci. These approximate selection for gene colocalization (TMD), the degree of gene coinheritance within an HLG/GCF (GCL), protein structural conservation (CSB), and loss and horizontal transfer (PDS).

Metabolic gene clusters that perform broadly distributed, or generalized, functions are identified within the highest TMD USS loci, whereas more specialized functions tend to have lower TMD values. The USS proxy is based on detecting unexpected TMD values of groups of co-occurring gene families, though a final TMD threshold is also applied to HLGs since multiple groups of shared synteny loci may represent a GCF. Increasing HLG TMD thresholds increasingly enriches general metabolic process genes. On the other hand, thresholds above 0.4 give diminishing returns for enriching secondary metabolic process genes, a subset of the general metabolic process category. This observation aligns with our finding that specialized metabolite biosynthesis clusters have the minimum observed TMDs, whereas the largest TMD clusters perform more essential metabolic functions (46) (Figure 3A). The limited distribution of specialized clusters may be attributed to their selection by impermanent, specific ecological factors, while clusters that perform essential/broadly distributed functions, such as biotin biosynthesis and iron assimilation, are less affected by environmental heterogeneity (91–94).

Gene Commitment to the Locus (GCL) quantifies the similarity between a gene's distribution and the distribution of its HLG, as a proxy for the degree to which genes are co-inherited within gene clusters. Selection is thought to maintain the clustered state through selection against partial loss of co-adapted genes contributing to the same phenotype (95) and through increased fitness of clustered genes following horizontal transfer (55). We therefore expect clustered genes to be more committed to their clusters than non-clustered genes. GCL enriches secondary metabolic process genes across the tested thresholding range, which suggests genes within specialized gene clusters more often remain clustered. This may result from selection acting on the whole specialized function, which preserves the clustered state by increasing the fitness of clustered genes (77,96). GCL may have diminishing returns in filtering for general metabolic gene clusters because essential metabolic functions are likely to be retained by selection even in the absence of clustering in certain lineages. For instance, lineages may have insufficient effective population sizes to drive clustering through selection for metabolic efficiency (23,24,97).

We implemented the Conservative Substitution Bias (CSB) as a proxy for selection that maintains the coordinated function of genes in gene clusters (97). Genes within clusters coordinately carry-out cooperative metabolic function by co-locating, co-expressing, and binding (97–99). While coordinated function is necessary for clusters to produce their phenotypes, selection for coordinated function may additionally be influenced by the accumulation of toxic intermediate metabolites that are produced from partial metabolic pathways that result from discordance (56). Filtering using CSB thresholds may enrich specialized metabolic clusters in part due to CSB caused by selection for protein-protein interactions that minimize toxic metabolic phenotypes (56).

Phylogenetic Distribution Sparsity (PDS) is the most powerful filtration proxy for secondary metabolic gene clusters, even though many have a low PDS value. PDS approximates the prevalence of horizontal transfer and loss that shape the distribution of metabolic gene clusters (24,55,58,94). Gene clusters may be subject to elevated horizontal transfer when they contain complete selectable phenotypes (55,94,96) and can be rapidly lost in lineages when the cost of the cluster phenotype exceeds its fitness contribution (58). Additionally, the prevalence of horizontal transfer and loss may vary by lineage (96), or type of cluster (Figure 3a), which results in the hypervariability of PDS values (minimum 0%, median 93.6%, standard deviation 26.3%). Gene clusters with specialized functions are often tied to environment-interaction phenotypes, and therefore selected by the ecology of the organisms that contain them (29,30,42,77). This specialization may lead to horizontal transfer among organisms occupying similar niches (42)and cluster loss upon niche switching. Our results support PDS as a proxy for specialized metabolic clusters, as increasing minimum PDS thresholds superlinearly extracts clusters with secondary metabolic process genes. PDS also enriches general metabolic process genes, perhaps in part due to spurious annotation of cryptic specialized metabolic functions as general metabolic processes (100), though PDS selectivity for specialized metabolite clusters becomes more prominent as PDS surpasses 0.2. Horizontal transfer can act across large phylogenetic distances relative to the phylogenetic distribution of the cluster and thus drastically increase specialized cluster PDS signal compared to loss (30,76). Some accessory metabolic clusters that perform functions with broad usefulness, such as nitrate assimilation, may also transfer across large phylogenetic distances (25,101). It is important to note that in situations where TMD approaches a complete distribution, any decrease in PDS could potentially be attributed to sampling limitations (Supplementary Equation S4). However, we did not identify any reference clusters or extracted GCFs that had nearly complete distributions.

Proxies for coordinated gene evolution are partially linearly correlated. We found that TMD is correlated with the other proxies (adjusted R2 > 0.5), which suggests that the distribution of a cluster positively associates with the variance in selection for coinheritance, coordinated function, and prevalence of horizontal transfer and loss. The correlation between TMD and CSB is relatively high (adjusted R2 = 0.861), indicating that broadly distributed clusters have increased selection for structural conservation. This may be because clusters that have large TMDs often perform metabolic functions that are more essential to organismal ecology and are under intense purifying selection to preserve the integrated function, whereas clusters that perform more specialized functions can adapt to new ecological roles and experience positive selection (79–81). Alternatively, the correlation between TMD and CSB may result from increased signal from conservative amino acid substitutions at greater phylogenetic distances. Larger TMD values are associated with larger GCL values (adjusted R2 = 0.816), which could be because GCL is more significantly affected by smaller samples, and a correction for sample size or a phylogenetic approach may be warranted.

CLOCI more accurately identifies gene cluster families when compared to existing algorithms

When compared to the function-centric algorithm antiSMASH, CLOCI recovers more biosynthetic gene clusters and greater cluster diversity. The 11.8% of clusters CLOCI fails to recover can largely be attributed to sparse sampling in our dataset. CLOCI also has greater general cluster recovery when compared to an existing USS algorithm, EvolClust. We attribute CLOCI improvements to USS detection by more precisely inferring cluster boundaries through a robust gene cluster grouping algorithm that globally identifies GCFs. The improvements to recovery and boundary accuracy position CLOCI as a paradigm-shift in gene cluster detection from function-centric algorithms to coordinate gene evolution detection.

CLOCI recovers more gene clusters than function-centric algorithms. CLOCI and the function-centric algorithm, antiSMASH, both detected the majority of biosynthetic gene clusters in the dataset (86.0% and 76.1% respectively). CLOCI additionally inferred 100% of the non-biosynthetic gene cluster categories because USS among gene clusters is a property independent of gene function. In contrast, antiSMASH was initially designed to detect biosynthetic gene clusters, and the capacity for antiSMASH to infer diverse metabolic pathways can be attributed to the versatility of the modeled core genes it searches for. Modeling canonical biosynthetic clusters, such as polyketide and nonribosomal peptide clusters, can reveal most of the known fungal biosynthetic cluster diversity because these clusters can synthesize diverse metabolites from modular core genes (91,92). Non-biosynthetic clusters can be incorporated into antiSMASH detection, though function-centric non-biosynthetic screening may be relatively restricted to the modeled cluster family (27) because non-biosynthetic cluster classes have relatively conserved phenotypes (58). Similar to non-biosynthetic functions, non-canonical biosynthetic gene clusters, such as the psilocybin cluster, may have discrete distributions with conserved phenotypes, so modeling these functions may also have limited capacity for inferring new clusters. In other words, non-canonical clusters may not contain versatile, modular core genes and may have limited capacity for diversified biosynthesis. Therefore, even if antiSMASH and function-centric detection continue to add newly discovered core genes they cannot infer cluster classes de novo because they have not been modeled. CLOCI inherently detects these GCFs, including psilocybin biosynthesis, without having to account for their functions a priori.

While CLOCI successfully identifies most reference clusters, eight of the 68 reference gene clusters (11.8%) were not detected primarily due to genome sample limitations. We failed to recover the zearalenone cluster because it is minimally sampled and has a small distribution within Fusarium (102), which did not meet USS thresholds. Both the zearalenone and fumiquinazoline clusters did not meet USS thresholds in part because their taxonomic orders, Hypocreales and Eurotiales respectively, have relatively highly clustered genomes (24,58,96), so USS thresholds may need to be lowered to identify them. The mycophenolic acid gene cluster was overlooked because the locus was incorrectly assembled due to low support that is partially attributed to excluding Penicillium brevicompactum genomes that failed our genome quality control prior to the analysis. AntiSMASH detects these clusters, which highlights its usefulness compared to CLOCI in detecting canonical gene clusters in lineages with highly clustered genomes because canonical core genes are easily recovered with homology searching and are likely to belong to a gene cluster. Neither antiSMASH nor CLOCI detected the reference gene of the cephalosporin gene cluster in Acremonium chrysogenum. antiSMASH missed the referenced region because it lacks a canonical core gene, and CLOCI overlooked the region because it is a two-gene locus that does not meet CLOCI thresholds for shared synteny. It is worth noting that cephalosporin biosynthesis in A. chrysogenum comprises multiple loci, and both CLOCI and antiSMASH detect the locus that contains the cephalosporin NRPS (103). Both algorithms also failed to detect the ent-kaurene cluster, which CLOCI misses because it primarily comprises multiple tandem gene duplications. Recent tandem duplicates are binned into the same homology group and our CLOCI analysis settings constrained detection to at least three unique homology groups within a five gene sliding window. Missing clusters highlight the benefit of further genome sampling, and suggest CLOCI recovery will improve with further tuning of algorithm variables.

CLOCI improves USS-based detection by implementing an approach that more accurately defines shared synteny loci (Figure 2). The discrepancy in recovery of our cluster recovery dataset between CLOCI (88.2%) and the existing USS algorithm, EvolClust (51.0%), may be partially explained by our increased genome sample because both algorithms rely on sufficiently sampling loci to support cluster detection. However, the EvolClustDB also recovers fewer clusters (75.0%) than CLOCI (84.8%) in the independent cluster boundary dataset (Supplementary Table S5). The cluster boundary dataset is primarily composed of Aspergillus spp., a genus which is well-sampled by both algorithms, and thus we attribute CLOCI’s improvement in recovery to the alternative approach in inferring shared synteny loci. EvolClust implements a global heuristic model of shared synteny locus sizes with a minimum size of five genes to mitigate false positive gene cluster designation. In contrast, CLOCI can detect clusters from as small as two genes up to the chromosome length by assembling gene clusters from domains of shared microsynteny that more accurately recapitulate the boundaries of homologous loci. CLOCI separates the three gene Ustilago nitrate assimilation and Saccharomyces galactose metabolism gene clusters from the surrounding unexpectedly conserved regions, whereas EvolClust infers both gene clusters with 16 extra genes each.

Our results corroborate previous findings (52) that shared synteny detection more precisely infers cluster boundaries than both standalone boundary prediction in antiSMASH and antiSMASH with CASSIS (Figure 2b). According to our functional cluster boundary definition, boundary inference using antiSMASH with or without CASSIS on average infers more than double the reference cluster size. We attribute the similarity between antiSMASH with or without CASSIS enabled primarily to CASSIS adjusting boundaries into surrounding intergenic space, which was unaccounted for in our analysis, or not substantially refining some clusters. CASSIS overestimating cluster boundaries can be attributed to its algorithm operationally defining gene clusters boundaries as entire coregulated regions. In contrast, CLOCI predicts a median zero extraneous and missing gene count per reference cluster. This comparison suggests CASSIS may be more useful in detecting coregulated regions associated with secondary metabolite biosynthesis, while CLOCI appears to more accurately report the size of the cluster directly associated with synthesis of particular metabolites. While the majority of CLOCI clusters are reported as a single shared microsynteny domain (82.1% of reference clusters), we identified some reference clusters as two domains (14.3%) and one cluster as three domains. These domains are directly adjacent, though CLOCI does not report them with other evidence of linkage. The granularity of sub-cluster domains is tunable, and poorly sampled clusters will become reportable as a single domain as the quality of the genome sample improves.

The reference boundary dataset we adopted from CASSIS is biased toward Aspergillus gene clusters, due to the limited research into functional validation of gene cluster boundaries in other taxa. We supplemented with several non-Aspergillus clusters, including the psilocybin gene cluster from Basidiomycota (Supplementary Table S6), though future benchmarking against a more taxonomically diverse dataset of reference boundaries will be necessary to support CLOCI boundary prediction in broader taxonomic datasets.

CLOCI enables unbiased gene cluster analysis

CLOCI does not depend on a priori assumptions of function, which enables unbiased comparative ecological studies on genomes’ cluster repertoires that include ecologically-relevant non-biosynthetic clusters. CLOCI further extends the capacity to globally detect gene cluster categories to GCF circumscription, where GCFs are globally identified using Markov Clustering (MCL) and the granularity of GCF identification is tunable through a single parameter. Examining the network topology of these GCFs can provide insight into the evolution of gene clusters, including evidence of horizontal transfer. We also sift these GCFs for unreported primary and accessory metabolic gene clusters by filtering minimum criteria for coordinate gene evolution proxy values.

CLOCI presents an unbiased framework for ascertaining the ecological functions of genomes as products of their gene cluster repertoires. Gene cluster detection algorithms are extensively implemented to study and compare the ecology of genomes. However, these studies primarily implement function-centric algorithms that are biased toward well-studied lineages and only detect canonical specialized metabolite biosynthesis gene clusters (38,39). Despite this bias, function-centric detection has dominated the literature, in part because previous function-agnostic algorithms do not recover as many gene clusters (Figure 2A). CLOCI recovers more gene clusters without a priori assumptions of function, which presents an unbiased framework for comparatively profiling gene cluster repertoires. This approach has the advantage of cataloging ecologically-significant non-biosynthetic gene clusters that are overlooked by function-centric algorithms, such as nitrate assimilation and galactose metabolism, which presents a broader window into the ecologies of gene cluster repertoires.

CLOCI implements a globally tunable framework for circumscribing GCFs. Robust GCF circumscription accounts for microsynteny, amino acid similarity, and protein domain presence-absence to group gene clusters into ecologically-significant GCFs (29,30), and tie biosynthetic GCF granularity to metabolite structural modifications (54,104,105). GCF circumscription algorithms have had to be tuned for each type of gene cluster because these algorithms use weighted coefficients that must be tailored to the different evolutionary rates of different GCF types (28,106). In contrast, CLOCI implements a globally-tunable model that affects granularity for all types of gene clusters using a single parameter, the MCL inflation value. The MCL random-walk approach to graph clustering is attractive for GCF identification because MCL incorporates normalization at each step, which allows for divergent local graph similarities that reflects the different evolutionary rates of different GCFs. The size-flexibility of this approach is demonstrated by CLOCI inferring the nitrate assimilation GCFs (assimilative) in multiple phyla, whereas the psilocybin GCF (non-canonical biosynthetic) was recovered in only five organisms. Increasing the inflation value will globally increase the granularity of GCF inferences, which can be compared with metabolic data to identify functionally-significant modifications in GCFs.

GCF network topology is consistent with horizontal transfer events and deep phylogenetic divergence in the evolution of gene clusters. We identified at least five nitrate assimilation GCFs that contain characterized nitrate assimilation genes (18,25,107). The locus-locus network topology suggests these divergent GCFs are homologous (Figure 4F). The GCFs correspond with discrete subnetworks because deep phylogenetic divergence results in low sequence similarity and variation in gene composition (25,86,101). For example, some Saccharomycotina nitrate assimilation clusters contain unique transcription factors that result in decreased locus-locus similarity with other nitrate assimilation GCFs (85). Alternatively, convergent evolution may also result in multiple GCFs. Horizontal transfer is suggested by the presence of distantly related taxa within a subnetwork. For example, the aggregation of Hypocreaceae/Bionectriaceae (Ascomycota) and Basidiomycota nitrate assimilation clusters in a common GCF is consistent with known horizontal transfer between these lineages (25). We also identified a nitrate assimilation GCF composed of Mucoromycota and Basidiomycota clusters. Phylogenetic analysis of all genes within this GCF supports the cluster was transferred from Mucoromycota species to Amanita (Basidiomycota). The cluster was only recovered in ectomycorrhizal Amanita spp. (87) and is missing from the saprobic Amanita inopinata and Amanita thiersii. This suggests acquisition of the nitrate assimilation GCF coincided with the emergence of the mycorrhizal ecology in Amanita species, possibly selected by the increased nitrogen demand imposed by host trees.

Proxies of coordinated gene evolution characterize the evolutionary properties of GCFs and can select for clusters of interest. The vast majority of CLOCI output are not yet characterized gene clusters (Figure 3B). The quantity of uncharacterized clusters within our results is consistent with the hypothesis that most clusters have not been studied (4,108,109), and presents an exciting opportunity for mining novel biosynthesis and ecologically important clusters from CLOCI output. To demonstrate CLOCI’s capacity to identify generalized cluster functions, we identified an unreported gene cluster with high log normalized TMD that contains the basidioferrin synthetase throughout the mushroom-forming subphylum, Agaricomycotina. If the function of this cluster is conserved, then it presents a possible conserved mechanism for iron assimilation in Agaricomycotina. We additionally searched the upper 90 percentile of PDS scores to identify putative specialized metabolite clusters. We identified a GCF with an undescribed nonribosomal peptide biosynthetic cluster (Figure 4B) in three different taxonomic classes, which could potentially be explained by multiple horizontal transfer events or convergent origins. The approach we employed to identify this putative specialized metabolite gene cluster can be used as the foundation of specialized metabolite biosynthesis detection, including relating hypothesized non-canonical biosynthetic pathways to CLOCI-predicted clusters.

CLOCI is a lineage-independent foundation for detecting shared synteny loci that are not clusters

CLOCI is tailored toward detecting gene clusters comprised of genes that do not have recent shared ancestry, though future algorithm adjustments may account for these clusters. Clusters that exclusively comprise recent homologs, particularly due to tandem gene duplication, will not be detected by CLOCI. Tandem duplicate clusters are common in animals, plants, and some mushroom-forming fungi. Future work will focus on expanding the capacity to detect such clusters by increasing the sliding window, accounting for tandem duplications in homology group co-occurrences, separating tandem duplicate orthologs into independent homology groups, and further tuning detection parameters with machine learning.

CLOCI can be implemented to classify HLGs and identify unexpectedly shared synteny regions without proxies, or proxies can be adapted to different types of genomic loci that are also under selection for coordinated gene evolution. On its own, HLG inference is an algorithmic adaptation of orthologous gene group (orthogroup) circumscription to the locus level, which thus provides a means for comparing the presence and absence of functionally significant gene neighborhoods. CLOCI HLG detection can even predict loci that are fragmented across multiple contigs, which can be used to inform genome assembly scaffolding (110). Following HLG inference, proxies of coordinated gene evolution can be tuned to accommodate different types of shared synteny loci. For example, defense islands contain genes and gene elements necessary for recognition and immunity against viruses and other antagonistic elements (111) and retroviral DNA is integrated into host genomes in discrete regions (112). CLOCI HLGs provide a powerful foundation for identifying these shared synteny loci in conjunction with gene function assignment.

Recommendations and considerations for implementation

CLOCI is constructed in a comparative genomic framework, and thus requires a large enough sample of sufficiently diverse genomes to identify unexpectedly shared synteny. The inputted genome sample should be based on lineages’ rates of microsynteny decay and the distribution of their gene clusters. CLOCI may best be implemented at least on a subphylum-level dataset to account for the majority of horizontal transfers. Future work will incorporate a single or batch genome module where users can predict HLGs/GCFs from a single genome referencing detected HLGs from a previous large-scale CLOCI analysis.

Improvements to CLOCI’s memory usage and time complexity should be focused on homology group combination inference, CSB quantification, and GCF circumscription, which involve pairwise loci comparisons and self-alignment of entire homology groups. We implement methods to decrease the amount of comparisons, such as homology group pair seeding and gene cluster clan classification; however, the time complexity of these modules remains non-linear. A compiled language implementation of these modules would additionally increase their throughput.

Supplementary Material

gkae625_Supplemental_Files

Acknowledgements

We would like to acknowledge the Ohio Supercomputer Center for providing cutting-edge high performance computing resources. We thank Dr Emile Gluck-Thaler for consultation regarding the development of CO-OCCUR and its application as the conceptual foundation of CLOCI.

Data availability

The data underlying this article are available in the article and in its online supplementary material. CLOCI is available at github.com/xonq/cloci and DOI 10.6084/m9.figshare.24424657.

Supplementary data

Supplementary Data are available at NAR Online.

Funding

National Science Foundation [DEB-1638999 to J.C.S.]; Ohio State University Translational Plant Sciences fellowship (to Z.K.). Funding for open access charge: Ohio State University Department of Plant Pathology.

Conflict of interest statement. None declared.
==== Refs
References

1. Blin K. , ShawS., KloostermanA.M., Charlop-PowersZ., van WezelG.P., MedemaM.H., WeberT. antiSMASH 6.0: improving cluster detection and comparison capabilities. Nucleic Acids Res. 2021; 49 :W29–W35.33978755
2. Burger G. , GrayM.W., ForgetL., LangB.F. Strikingly bacteria-like and gene-rich mitochondrial genomes throughout jakobid protists. Genome Biol. Evol. 2013; 5 :418–438.23335123
3. Ettema T.J.G. , BrinkmanA.B., LamersP.P., KornetN.G., de VosW.M., van der OostJ. Molecular characterization of a conserved archaeal copper resistance (cop) gene cluster and its copper-responsive regulator in Sulfolobus solfataricus P2. Microbiology. 2006; 152 :1969–1979.16804172
4. Keller N.P. Translating biosynthetic gene clusters into fungal armor and weaponry. Nat. Chem. Biol. 2015; 11 :671–677.26284674
5. Mihali T.K. , KellmannR., NeilanB.A. Characterisation of the paralytic shellfish toxin biosynthesis gene clusters in Anabaena circinalis AWQC131C and aphanizomenon sp. NH-5. BMC Biochem. 2009; 10 :8.19331657
6. Nützmann H.-W. , OsbournA. Gene clustering in plant specialized metabolism. Curr. Opin. Biotechnol. 2014; 26 :91–99.24679264
7. Nofiani R. , Mattos-ShipleyK.d., LebeK.E., HanL.-C., IqbalZ., BaileyA.M., WillisC.L., SimpsonT.J., CoxR.J Strobilurin biosynthesis in basidiomycete fungi. Nat. Commun. 2018; 9 :3940.30258052
8. Newman D.J. , CraggG.M. Natural products as sources of new drugs over the nearly four decades from 01/1981 to 09/2019. J. Nat. Prod. 2020; 83 :770–803.32162523
9. Linnemannstöns P. , PradoM., Fernández-MartínR., TudzynskiB., AvalosJ. A carotenoid biosynthesis gene cluster in Fusarium fujikuroi: the genes carB and carRA. Mol. Genet. Genomics. 2002; 267 :593–602.12172798
10. Alberti F. , KhairudinK., VenegasE.R., DaviesJ.A., HayesP.M., WillisC.L., BaileyA.M., FosterG.D. Heterologous expression reveals the biosynthesis of the antibiotic pleuromutilin and generates bioactive semi-synthetic derivatives. Nat. Commun. 2017; 8 :1831.29184068
11. Díez B. , GutiérrezS., BarredoJ.L., van SolingenP., van der VoortL.H., MartínJ.F. The cluster of penicillin biosynthetic genes. Identification and characterization of the pcbAB gene encoding the alpha-aminoadipyl-cysteinyl-valine synthetase and linkage to the pcbC and penDE genes. J. Biol. Chem. 1990; 265 :16358–16365.2129535
12. Brandenburger E. , GresslerM., LeonhardtR., LacknerG., HabelA., HertweckC., BrockM., HoffmeisterD. A highly conserved basidiomycete peptide synthetase produces a trimeric hydroxamate siderophore. Appl. Environ. Microbiol. 2017; 83 :e01478-17.28842536
13. Perrin R.M. , FedorovaN.D., BokJ.W., R.A.C.Jr, WortmanJ.R., KimH.S., NiermanW.C., KellerN.P. Transcriptional regulation of chemical diversity in Aspergillus fumigatus by LaeA. PLoS Pathog. 2007; 3 :e50.17432932
14. Gluck-Thaler E. , SlotJ.C. Specialized plant biochemistry drives gene clustering in fungi. ISME J. 2018; 12 :1694–1705.29463891
15. Arst H.N. , MacDonaldD.W. A gene cluster in Aspergillus nidulans with an internally located cis-acting regulatory region. Nature. 1975; 254 :26–31.1089903
16. Douglas H.C. , HawthorneD.C. Regulation of genes controlling synthesis of the galactose pathway enzymes in yeast. Genetics. 1966; 54 :911–916.5970626
17. Greene G.H. , McGaryK.L., RokasA., SlotJ.C. Ecology drives the distribution of specialized tyrosine metabolism modules in fungi. Genome Biol. Evol. 2014; 6 :121–132.24391152
18. Johnstone I.L. , McCabeP.C., GreavesP., GurrS.J., ColeG.E., BrowM.A.D., UnklesS.E., ClutterbuckA.J., KinghornJ.R., InnisM.A. Isolation and characterisation of the crnA-niiA-niaD gene cluster for nitrate assimilation in Aspergillus nidulans. Gene. 1990; 90 :181–192.2205530
19. Magliano P. , FlipphiM., SanglardD., PoirierY. Characterization of the Aspergillus nidulans biotin biosynthetic gene cluster and use of the bioDA gene as a new transformation marker. Fungal Genet. Biol. 2011; 48 :208–215.20713166
20. Fritsch E.F. , LawnR.M., ManiatisT. Molecular cloning and characterization of the human β-like globin gene cluster. Cell. 1980; 19 :959–972.6155216
21. Forrester W.C. , ThompsonC., ElderJ.T., GroudineM. A developmentally stable chromatin structure in the human beta-globin gene cluster. Proc. Natl. Acad. Sci. 1986; 83 :1359–1363.3456593
22. Lewis E.B. A gene complex controlling segmentation in Drosophila. Nature. 1978; 276 :565–570.103000
23. Slot J.C. Townsend J.P. , WangZ. Chapter four - fungal gene cluster diversity and evolution. Advances in Genetics, Fungal Phylogenetics and Phylogenomics. 2017; 100 :Academic Press 141–178.
24. Wisecaver J.H. , SlotJ.C., RokasA. The evolution of fungal metabolic pathways. PLoS Genet. 2014; 10 :e1004816.25474404
25. Slot J.C. , HibbettD.S. Horizontal transfer of a nitrate assimilation gene cluster and ecological transitions in fungi: a phylogenetic study. PLoS One. 2007; 2 :e1097.17971860
26. Gorfer M. , BlumhoffM., KlaubaufS., UrbanA., InselsbacherE., BandianD., MitterB., SessitschA., WanekW., StraussJ. Community profiling and gene expression of fungal assimilatory nitrate reductases in agricultural soil. ISME J. 2011; 5 :1771–1783.21562596
27. Pascal Andreu V. , AugustijnH.E., ChenL., ZhernakovaA., FuJ., FischbachM.A., DoddD., MedemaM.H gutSMASH predicts specialized primary metabolic pathways from the human gut microbiota. Nat. Biotechnol. 2023; 41 :1416–1423.36782070
28. Gluck-Thaler E. , HaridasS., BinderM., GrigorievI.V., CrousP.W., SpataforaJ.W., BushleyK., SlotJ.C. The architecture of metabolism maximizes biosynthetic diversity in the largest class of fungi. Mol. Biol. Evol. 2020; 37 :2838–2856.32421770
29. Franco M.E.E. , WisecaverJ.H., ArnoldA.E., JuY.-M., SlotJ.C., AhrendtS., MooreL.P., EastmanK.E., ScottK., KonkelZ.et al . Ecological generalism drives hyperdiversity of secondary metabolite gene clusters in xylarialean endophytes. New Phytol. 2022; 233 :1317–1330.34797921
30. Scott K. , KonkelZ., Gluck-ThalerE., DavidG.E.V., SimmtC.F., GrootmyersD., ChaverriP., SlotJ. Endophyte genomes support greater metabolic gene cluster diversity compared with non-endophytes in Trichoderma. PLoS ONE. 2023; 18 :e0289280.38127903
31. Monciardini P. , IorioM., MaffioliS., SosioM., DonadioS. Discovering new bioactive molecules from microbial sources. Microb. Biotechnol. 2014; 7 :209–220.24661414
32. Hewage R.T. , AreeT., MahidolC., RuchirawatS., KittakoopP. One strain-many compounds (OSMAC) method for production of polyketides, azaphilones, and an isochromanone using the endophytic fungus dothideomycete sp. Phytochemistry. 2014; 108 :87–94.25310919
33. Gressler M. , LöhrN.A., SchäferT., LawrinowitzS., SeiboldP.S., HoffmeisterD. Mind the mushroom: natural product biosynthetic genes and enzymes of basidiomycota. Nat. Prod. Rep. 2021; 38 :702–722.33404035
34. Gao S.-S. , LiX.-M., WilliamsK., ProkschP., JiN.-Y., WangB.-G. Rhizovarins A–F, Indole-diterpenes from the Mangrove-derived endophytic fungus mucor irregularis QEN-189. J. Nat. Prod. 2016; 79 :2066–2074.27462726
35. Adpressa D.A. , ConnollyL.R., KonkelZ.M., NeuhausG.F., ChangX.L., PierceB.R., SmithK.M., FreitagM., LoesgenS. A metabolomics-guided approach to discover fusarium graminearum metabolites after removal of a repressive histone modification. Fungal Genet. Biol. 2019; 132 :103256.31344458
36. Yaegashi J. , OakleyB.R., WangC.C.C. Recent advances in genome mining of secondary metabolite biosynthetic gene clusters and the development of heterologous expression systems in Aspergillus nidulans. J. Ind. Microbiol. Biotechnol. 2014; 41 :433–442.24342965
37. Khaldi N. , SeifuddinF.T., TurnerG., HaftD., NiermanW.C., WolfeK.H., FedorovaN.D. SMURF: genomic mapping of fungal secondary metabolite clusters. Fungal Genet. Biol. 2010; 47 :736–741.20554054
38. Terlouw B.R. , BlinK., Navarro-MuñozJ.C., AvalonN.E., ChevretteM.G., EgbertS., LeeS., MeijerD., RecchiaM.J.J., ReitzZ.L.et al . MIBiG 3.0 : a community-driven effort to annotate experimentally validated biosynthetic gene clusters. NucleicAcids Res. 2022; 51 :D603–D610.
39. Medema M.H. , BlinK., CimermancicP., de JagerV., ZakrzewskiP., FischbachM.A., WeberT., TakanoE., BreitlingR. antiSMASH: rapid identification, annotation and analysis of secondary metabolite biosynthesis gene clusters in bacterial and fungal genome sequences. Nucleic Acids Res. 2011; 39 :W339–W346.21672958
40. Fricke J. , BleiF., HoffmeisterD. Enzymatic synthesis of psilocybin. Angew. Chem. Int. Ed. 2017; 56 :12352–12355.
41. Obermaier S. , MüllerM. Ibotenic acid biosynthesis in the fly agaric is initiated by glutamate hydroxylation. Angew. Chem. Int. Ed. 2020; 59 :12432–12435.
42. Reynolds H.T. , VijayakumarV., Gluck-ThalerE., KorotkinH.B., MathenyP.B., SlotJ.C. Horizontal gene cluster transfer increased hallucinogenic mushroom diversity. Evol. Lett. 2018; 2 :88–101.30283667
43. Voigt K. , WolfT., OchsenreiterK., NagyG., KaergerK., ShelestE., PappT. Hoffmeister D. 15 Genetic and metabolic aspects of primary and secondary metabolism of the zygomycetes. Biochemistry and Molecular Biology, the Mycota. 2016; Cham Springer International Publishing 361–385.
44. Wisecaver J.H. , BorowskyA., TzinV., JanderG., KliebensteinD., RokasA A global coexpression network approach for connecting genes to specialized metabolic pathways in plants | plant cell. Plant Cell. 2017; 29 :944–959.28408660
45. Venice F. , DesiròA., SilvaG., SalvioliA., BonfanteP. The mosaic architecture of NRPS-PKS in the arbuscular mycorrhizal fungus gigaspora margarita shows a domain with bacterial signature. Front. Microbiol. 2020; 11 :581313.33329443
46. Wolf T. , ShelestV., NathN., ShelestE. CASSIS and SMIPS: promoter-based prediction of secondary metabolite gene clusters in eukaryotic genomes. Bioinformatics. 2016; 32 :1138–1143.26656005
47. Langfelder P. , HorvathS. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008; 9 :559.
48. Brakhage A.A. Regulation of fungal secondary metabolism. Nat. Rev. Microbiol. 2013; 11 :21–32.23178386
49. Blin K. , WolfT., ChevretteM.G., LuX., SchwalenC.J., KautsarS.A., Suarez DuranH.G., de los SantosE.L.C., KimH.U., NaveM.et al . antiSMASH 4.0—Improvements in chemistry prediction and gene cluster boundary identification. Nucleic. Acids. Res. 2017; 45 :W36–W41.28460038
50. Haas D. , BarbaM., VicenteC.M., NezbedováŠ., GarénauxA., Bury-MonéS., LorenziJ.-N., HôtelL., LauretiL., ThibessardA.et al . SYNTERUPTOR: mining genomic islands for non-classical specialised metabolite gene clusters. NAR Genom. Bioinform. 2024; 6 :lqae069.38915823
51. Winter S. , JahnK., WehnerS., KuchenbeckerL., MarzM., StoyeJ., BöckerS. Finding approximate gene clusters with Gecko 3. Nucleic Acids Res. 2016; 44 :9600–9610.27679480
52. Marcet-Houben M. , GabaldónT. EvolClust: automated inference of evolutionary conserved gene clusters in eukaryotes. Bioinformatics. 2020; 36 :1265–1266.31560365
53. Vignolle G.A. , SchafferD., ZehetnerL., MachR.L., Mach-AignerA.R., DerntlC. FunOrder: a robust and semi-automated method for the identification of essential biosynthetic genes through computational molecular co-evolution. PLoS Comput. Biol. 2021; 17 :e1009372.34570757
54. Louwen J.J.R. , KautsarS.A., BurgS.v., MedemaM.H., HooftJ.J.J.v. iPRESTO: automated discovery of biosynthetic sub-clusters linked to specific natural product substructures. PLoS Comput. Biol. 2023; 19 :e1010462.36758069
55. Lawrence J.G. , RothJ.R. Selfish operons: horizontal transfer may drive the evolution of gene clusters. Genetics. 1996; 143 :1843–1860.8844169
56. McGary K.L. , SlotJ.C., RokasA. Physical linkage of metabolic genes in fungi is an adaptation against the accumulation of toxic intermediate compounds. Proc. Natl. Acad. Sci. U.S.A. 2013; 110 :11481–11486.23798424
57. Price M.N. , HuangK.H., ArkinA.P., AlmE.J. Operon formation is driven by co-regulation and not by horizontal gene transfer. Genome Res. 2005; 15 :809–819.15930492
58. Rokas A. , WisecaverJ.H., LindA.L. The birth, evolution and death of metabolic gene clusters in fungi. Nat. Rev. Microbiol. 2018; 16 :731–744.30194403
59. Rose C.J. , ChapmanJ.R., MarshallS.D.G., LeeS.F., BatterhamP., RossH.A., NewcombR.D. Selective sweeps at the organophosphorus insecticide resistance locus, rop-1, have affected variation across and beyond the α-esterase Gene Cluster in the australian sheep blowfly, Lucilia cuprina. Mol. Biol. Evol. 2011; 28 :1835–1846.21228400
60. Douglass A.P. , ByrneK.P., WolfeK.H. The methylotroph gene order browser (MGOB) reveals conserved synteny and ancestral centromere locations in the yeast family Pichiaceae. FEMS Yeast Res. 2019; 19 :foz058.31397853
61. Emms D.M. , KellyS. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019; 20 :238.31727128
62. Li L. , StoeckertC.J., RoosD.S. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 2003; 13 :2178–2189.12952885
63. Konkel Z. , SlotJ.C. Mycotools: an automated and scalable platform for comparative genomics. 2023; bioRxiv doi: 12 September 2023, preprint: not peer reviewed 10.1101/2023.09.08.556886.
64. Steinegger M. , SödingJ. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat. Biotechnol. 2017; 35 :1026–1028.29035372
65. Simão F.A. , WaterhouseR.M., IoannidisP., KriventsevaE.V., ZdobnovE.M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015; 31 :3210–3212.26059717
66. Kim D. , GilchristC.L.M., ChunJ., SteineggerM. UFCG: database of universal fungal core genes and pipeline for genome-wide phylogenetic analysis of fungi. NucleicAcids Res. 2022; 51 :D777–D784.
67. Spatafora J.W. , ChangY., BennyG.L., LazarusK., SmithM.E., BerbeeM.L., BonitoG., CorradiN., GrigorievI., GryganskyiA.et al . A phylum-level phylogenetic classification of zygomycete fungi based on genome-scale data. Mycologia. 2016; 108 :1028–1046.27738200
68. Zhao T. , ZwaenepoelA., XueJ.-Y., KaoS.-M., LiZ., SchranzM.E., Van de PeerY. Whole-genome microsynteny-based phylogeny of angiosperms. Nat. Commun. 2021; 12 :3498.34108452
69. Kalyaanamoorthy S. , MinhB.Q., WongT.K.F., von HaeselerA., JermiinL.S. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods. 2017; 14 :587–589.28481363
70. Lewis P.O. A likelihood approach to estimating phylogeny from discrete morphological character data. Syst. Biol. 2001; 50 :913–925.12116640
71. Minh B.Q. , SchmidtH.A., ChernomorO., SchrempfD., WoodhamsM.D., HaeselerA., LanfearR. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 2020; 37 :1530–1534.32011700
72. Li Y. , SteenwykJ.L., ChangY., WangY., JamesT.Y., StajichJ.E., SpataforaJ.W., GroenewaldM., DunnC.W., HittingerC.T.et al . A genome-scale phylogeny of the kingdom Fungi. Curr. Biol. 2021; 31 :1653–1665.33607033
73. Li Y. , LiuH., SteenwykJ.L., LaBellaA.L., HarrisonM.-C., GroenewaldM., ZhouX., ShenX.-X., ZhaoT., HittingerC.T.et al . Contrasting modes of macro and microsynteny evolution in a eukaryotic subphylum. Curr. Biol. 2022; 32 :5335–5343.36334587
74. Enright A.J. , Van DongenS., OuzounisC.A. An efficient algorithm for large-scale detection of protein families. Nucleic. Acids. Res. 2002; 30 :1575–1584.11917018
75. Pedregosa F. , VaroquauxG., GramfortA., MichelV., ThirionB., GriselO., BlondelM., PrettenhoferP., WeissR., DubourgV.et al . Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 2011; 12 :2825–2830.
76. Navarro-Muñoz J.C. , CollemareJ. Evolutionary histories of type III polyketide synthases in Fungi. Front. Microbiol. 2020; 10 :3018.32038517
77. Slot J.C. , RokasA. Horizontal transfer of a large and highly toxic secondary metabolic gene cluster between fungi. Curr. Biol. 2011; 21 :134–139.21194949
78. Marcet-Houben M. , Collado-CalaI., Fuentes-PalaciosD., GómezA.D., MolinaM., Garisoain-ZafraA., ChorosteckiU., GabaldónT. EvolClustDB: exploring eukaryotic gene clusters with evolutionarily conserved genomic neighbourhoods. J. Mol. Biol. 2023; 435 :168013.36806474
79. Lim F.Y. , WonT.H., RaffaN., BaccileJ.A., WisecaverJ., RokasA., SchroederF.C., KellerN.P. Fungal isocyanide synthases and xanthocillin biosynthesis in Aspergillus fumigatus. mBio. 2018; 9 :e00785-18.29844112
80. Asch D.K. , ZieglerJ., MinX. Molecular evolution of genes involved in quinic acid utilization in fungi. Comput. Mol. Biol. 2021; 10.5376/cmb.2021.11.0005.
81. Price M.N. , DehalP.S., ArkinA.P. FastTree 2 – Approximately maximum-likelihood trees for large alignments. PLoS One. 2010; 5 :e9490.20224823
82. Katoh K. , MisawaK., KumaK., MiyataT. MAFFT: a novel method for rapid multiple sequence alignment based on fast fourier transform. Nucleic Acids Res. 2002; 30 :3059–3066.12136088
83. Steenwyk J.L. , IiiT.J.B., LiY., ShenX.-X., RokasA. ClipKIT: a multiple sequence alignment trimming software for accurate phylogenomic inference. PLoS Biol. 2020; 18 :e3001007.33264284
84. Shimodaira H. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 2002; 51 :492–508.12079646
85. Siverio J.M. Assimilation of nitrate by yeasts. FEMS Microbiol. Rev. 2002; 26 :277–284.12165428
86. Reynolds H.T. , BartonH.A., SlotJ.C. Phylogenomic analysis supports a recent change in nitrate assimilation in the White-nose Syndrome pathogen, pseudogymnoascus destructans. Fungal. Ecol. 2016; 23 :20–29.
87. Chaib De Mares M. , HessJ., FloudasD., LipzenA., ChoiC., KennedyM., GrigorievI.V., PringleA. Horizontal transfer of carbohydrate metabolism genes into ectomycorrhizal Amanita. New Phytol. 2015; 205 :1552–1564.25407899
88. Yu Y. , BlachowiczA., WillC., SzewczykE., GlennS., Gensberger-ReiglS., NowrousianM., WangC.C.C., KrappmannS. Mating-type factor-specific regulation of the fumagillin/pseurotin secondary metabolite supercluster in Aspergillus fumigatus. Mol. Microbiol. 2018; 110 :1045–1065.30240513
89. Hawkins A.R. , LambH.K., SmithM., KeyteJ.W., RobertsC.F. Molecular organisation of the quinic acid utilization (QUT) gene cluster in Aspergillus nidulans. Mol. Gen. Genet. 1988; 214 :224–231.2976880
90. Schmaler-Ripcke J. , SugarevaV., GebhardtP., WinklerR., KniemeyerO., HeinekampT., BrakhageA.A. Production of Pyomelanin, a second type of melanin, via the Tyrosine Degradation Pathway in Aspergillus fumigatus. Appl. Environ. Microbiol. 2009; 75 :493–503.19028908
91. Gokhale R.S. , SankaranarayananR., MohantyD. Versatility of polyketide synthases in generating metabolic diversity. Curr. Opin. Struct. Biol. 2007; 17 :736–743.17935970
92. Walsh C.T. Polyketide and nonribosomal peptide antibiotics: modularity and versatility. Science. 2004; 303 :1805–1810.15031493
93. Lind A.L. , WisecaverJ.H., LameirasC., WiemannP., PalmerJ.M., KellerN.P., RodriguesF., GoldmanG.H., RokasA. Drivers of genetic diversity in secondary metabolic gene clusters within a fungal species. PLoS Biol. 2017; 15 :e2003583.29149178
94. Wisecaver J.H. , RokasA. Fungal metabolic gene clusters—caravans traveling across genomes and environments. Front. Microbiol. 2015; 6 :161.25784900
95. Slot J.C. , RokasA. Multiple GAL pathway gene clusters evolved independently and by different mechanisms in fungi. Proc. Natl. Acad. Sci. U.S.A. 2010; 107 :10136–10141.20479238
96. Slot J.C. , Gluck-ThalerE. Metabolic gene clusters, fungal diversity, and the generation of accessory functions. Curr. Opin. Genet. Dev. 2019; 58–59 :17–24.
97. Hurst L.D. , PálC., LercherM.J. The evolutionary dynamics of eukaryotic gene order. Nat. Rev. Genet. 2004; 5 :299–310.15131653
98. Chanda A. , RozeL.V., KangS., ArtymovichK.A., HicksG.R., RaikhelN.V., CalvoA.M., LinzJ.E. A key role for vesicles in fungal secondary metabolism. Proc. Natl. Acad. Sci. U.S.A. 2009; 106 :19533–19538.19889978
99. Voland P. , WeeksD.L., MarcusE.A., PrinzC., SachsG., ScottD. Interactions among the seven Helicobacter pylori proteins encoded by the urease gene cluster. Am. J. Physiol.-Gastrointest. Liver Physiol. 2003; 284 :G96–G106.12388207
100. Blei F. , FrickeJ., WickJ., SlotJ.C., HoffmeisterD. Iterative l-tryptophan methylation in psilocybe evolved by subdomain duplication. ChemBioChem. 2018; 19 :2160–2166.30098085
101. Ocaña-Pallarès E. , NajleS.R., ScazzocchioC., Ruiz-TrilloI. Reticulate evolution in eukaryotes: origin and evolution of the nitrate assimilation pathway. PLoS Genet. 2019; 15 :e1007986.30789903
102. Lysøe E. , BoneK.R., KlemsdalS.S. Real-time quantitative expression studies of the zearalenone biosynthetic gene cluster in Fusarium graminearum. Phytopathology®. 2009; 99 :176–184.19159310
103. Coque J.J.R. , MartinJ.F., CalzadaJ.G., LirasP. The cephamycin biosynthetic genes pcbAB, encoding a large multidomain peptide synthetase, and pcbC of Nocardia lactamdurans are clustered together in an organization different from the same genes in Acremonium chrysogenum and Penicillium chrysogenum. Mol. Microbiol. 1991; 5 :1125–1133.1956290
104. Caesar L.K. , MontaserR., KellerN.P., KelleherN.L. Metabolomics and genomics in natural products research: complementary tools for targeting new chemical entities. Nat. Prod. Rep. 2021; 38 :2041–2065.34787623
105. Caesar L.K. , ButunF.A., RobeyM.T., AyonN.J., GuptaR., DainkoD., BokJ.W., NicklesG., StankeyR.J., JohnsonD.et al . Correlative metabologenomics of 110 fungi reveals metabolite–gene cluster pairs. Nat. Chem. Biol. 2023; 19 :846–854.36879060
106. Navarro-Muñoz J.C. , Selem-MojicaN., MullowneyM.W., KautsarS.A., TryonJ.H., ParkinsonE.I., De Los SantosE.L.C., YeongM., Cruz-MoralesP., AbubuckerS.et al . A computational framework to explore large-scale biosynthetic diversity. Nat. Chem. Biol. 2020; 16 :60–68.31768033
107. Willmann A. , ThomfohrdeS., HaenschR., NehlsU. The poplar NRT2 gene family of high affinity nitrate importers: impact of nitrogen nutrition and ectomycorrhiza formation. Environ. Exp. Bot. 2014; 108 :79–88.
108. Keller N.P. Fungal secondary metabolism: regulation, function and drug discovery. Nat. Rev. Microbiol. 2019; 17 :167–180.30531948
109. Keller N.P. , TurnerG., BennettJ.W. Fungal secondary metabolism — From biochemistry to genomics. Nat. Rev. Microbiol. 2005; 3 :937–947.16322742
110. Meleshko D. , MohimaniH., TraccanaV., HajirasoulihaI., MedemaM.H., KorobeynikovA., PevznerP.A. BiosyntheticSPAdes: reconstructing biosynthetic gene clusters from assembly graphs. Genome Res. 2019; 29 :1352–1362.31160374
111. Makarova K.S. , WolfY.I., SnirS., KooninE.V. Defense islands in bacterial and archaeal genomes and prediction of novel Defense systems. J. Bacteriol. 2011; 193 :6039–6056.21908672
112. Fujiwara T. , MizuuchiK. Retroviral DNA integration: structure of an integration intermediate. Cell. 1988; 54 :497–504.3401925
