Monascus spp. is a highly valuable edible fungus that has been traditionally consumed in China and other Asian countries such as Japan, Republic of Korea, and the Philippines. It has high nutritional value due to the synthesis of a variety of beneficial secondary metabolites, including Monascus azaphilone pigments (MonAzPs), monacolin K (MK), aminobutyric acid, ergosterol, and Hong Qu polysaccharide (Wang et al., 2021), with MonAzPs and MK being of special concern. MonAzPs are a type of chromogenic chemical consisting of a chromophore with a polyketide structure and medium- or long-chain fatty acids. They are highly promising in the food, healthcare, and cosmetic industries due to their excellent coloring properties, good biological efficacy (antioxidant, anti-inflammatory, hypolipidemic, and anti-tumor properties), and non-toxic side effects as natural food additives (Chen et al., 2019). Another important metabolite, MK, is an inhibitor of 3-hydroxy-3-methylglutaryl-coenzyme A (HMG-CoA) reductase, a critical enzyme in endogenous cholesterol production (Zhang et al., 2020). It is a prescription drug with brand names Mevinoline, Lovastatin, or Mevalonate used to treat excessive cholesterol, coronary heart disease, and other disorders (Karthikeyan and Dharumadurai, 2023). However, concerns were raised regarding the safety of Monascus products when citrinin, a mycotoxin hazardous to both humans and animals, was detected in Monascus products as early as 1995 (Shao et al., 2014). Citrinin is nephrotoxic and can cause kidney enlargement, renal tubule expansion, and renal epithelial cell degeneration and necrosis. Fortunately, some citrinin-free Monascus strains and production techniques have been developed recently, including metabolic engineering to eliminate the production of citrinin by disrupting the polyketide synthase gene pksCT (Jia et al., 2010) or dehydrogenase gene citE (Ning et al., 2017), natural screening (Feng et al., 2016), mutagenesis (Kalaivani and Rajasekaran, 2014), genome shuffling (Ghosh and Dam, 2020), low pH (Kang et al., 2014), and genistein addition (Ouyang et al., 2021). With the increased safety of Monascus, people are becoming more interested in the health benefits of this edible fungus.

Despite the economic importance of Monascus, taxonomic research on this genus remains limited (Barbosa et al., 2017). Monascus spp. are categorized under the phylum Eumycota, subphylum Ascomycota, class Plectomycetes, order Eurotiales, and family Monascaceae. Morphologically, this genus generates non-porous perithecia at the top of the stem-like hypha, ascospores distributed throughout the hypha, virtually spherical to wide spherical ascospores, and transparent and oval ascospores detached from the closed capsule. In 1983 Hawksworth and Pitt updated the genus based on physiological and morphological criteria, reducing the number of recognized species to three: M. pilosus, M. ruber, and M. purpureus (Barbosa et al., 2017). Following the discovery of new species, the NCBI taxonomy database now contains more than twenty records on Monascus species.¹ Among these, M. purpureus has the most heterotypic synonyms, such as M. albidus, M. anka, M. erroneous, and M. rubiginosus.

More recently, DNA sequencing and molecular phylogenetics have become essential in microbiological taxonomy. ITS (Dai et al., 2021), LSU (He et al., 2020; Higa et al., 2020), β-tubulin (Tong et al., 2022), calmodulin (He et al., 2020; Ruiz and Radwan, 2021) are frequently used for accurate species identification of close relatives among Monascus spp. However, phylogenetic delineation of this genus was complicated by gene region inconsistency and low support for internal nodes. Phenotype-based identification schemes in Monascus have been difficult to reconcile with the results obtained by ITS, partial LSU and/or β-tubulin gene sequencing (Park and Jong, 2003; Park et al., 2004; Shao et al., 2011, 2014), and the genetic identities of Monascus species are still under debate or are even confusing. Back in Park and Jong (2003) investigated the phylogenetic relationships among the species using sequences from the D1/D2 region of the large subunit (LSU) rRNA genes. They found that M. ruber, M. pilosus, and M. purpureus were closely related and clustered into the same subgroup, but M. ruber and M. pilosus were unable to be distinguished from each other. In Park et al. (2004) amplified and sequenced the ITS and partial β-tubulin genes of 17 ATCC reference strains of Monascus species. Still, they found that M. pilosus and M. ruber could not be differentiated using these sequences (Park et al., 2004), implying that species boundaries in Monascus should be reexamined. In fact, although ITS is the formally recognized fungal barcode, it sometimes does not distinguish among closely related phylogenetic species (Ruiz and Radwan, 2021). Due to insufficient phylogenetic information and gene-specific noise, one or a few loci usually yield incongruent phylogenies, resulting in several weakly supported nodes.

To advance the understanding of phylogeny on these edible fungi, a deeper sampling of larger and identical gene sets across the genome is required (Binder et al., 2013; Ruiz and Radwan, 2021). Whole-genome sequencing (WGS) has provided a significant benefit in establishing phylogenetic relationships, genetic diversity, virulence-related components, and biotechnological features (Carrillo and Blais, 2021). The widespread availability of WGS effectively eliminated data availability as a limiting factor for inferring phylogenetic trees, offering hundreds to thousands of loci for analysis, and eventually replacing molecular phylogenetics with phylogenomics (Nagy and Szollosi, 2017). As the available genomes of Monascus spp., the first WGS project (M. ruber NRRL 1597, SRR1800507, Illumina HiSeq 2000 platform) was published in 2015 as part of the JGI 1,000 Fungal Genomes Project to represent members of the ascomycete family Monascaceae, while the first sequence assembly was deposited in NCBI database in 2017 (M. ruber ASM90018405v1) (Chen et al., 2015), currently up to more than 10 deposits. Only a few have been assembled to the chromosomal level and functionally annotated (Yang et al., 2015).² In the research conducted by Yang et al. (2015) they used 2,053 single-copy orthologs across the genome of M. purpureus YY-1 and other 17 genomes from the class Eurotiomycetes and Sordariomycetes to create a phylogenetic tree by the maximum-likelihood method, revealing the M. purpureus YY-1’s close relationship with the family Aspergillaceae. In another study, Higa et al. (2020) compared the biosynthetic gene clusters (BGCs) of MK, citrinin, and MonAzPs using WGS of M. pilosus NBRC 4520, M. purpureus NBRC 4478, and M. ruber NBRC 4483. They found that M. pilosus and M. purpureus are chemotaxonomically distinct while M. ruber has similar secondary metabolite BGCs to M. pilosus. The genome-scale approach to microbial taxonomy obviously offers improved resolution, stability, and reliability in evolutionary analyses compared to established methods of identifying physiologically and biochemically changeable and monogenic markers. Unfortunately, no phylogenomics studies on Monascus have been conducted within this genus.

On the other hand, WGS enables us to more fully comprehend the intrinsic mechanisms that underlie the phenotypic variations in nutritional profile, pathogenicity, host specificity, secondary metabolite synthesis and so on (San et al., 2019). Genomes within a species or all strains within a clade frequently have a core/conserved component as well as a variable set of genetic material that is referred to as a “pan-genome” among individuals or populations (Barber et al., 2021). The core genome comprises sequences present in all strains and is typically linked to biological functions and key phenotypic traits of the species. The variable/accessory/dispensable genomes contain sequences that are unique to one strain or subset of strains and are related to the adaptability of the species to specific environments or unique biological characteristics, reflecting the characteristics of the individuals. The size of the pan-genome is determined by the effective population size, lifestyle, and niche heterogeneity of the species. The gene pool of a species largely controls its ecological interactions and adaptive capacity, with core and variable genes contributing to the presence-absence variations (PAVs) (Siren et al., 2021). Therefore, a pan-genome of Monascus can allow us to better understand the relationship between individual characteristics and genetic variation within this species.

To expand the understanding of the phylogenetic relationships and metabolic divergence in Monascus species, this study collected 15 genome assemblies from the genus. They were compared at the whole-genome level using average nucleotide identity (ANI), whole-genome alignment (WGA), PAVs of the pan-genome, species tree inferred from all orthologous gene sequences (STAG), and species tree from single-copy orthologous genes (SCOG). Additionally, this study investigated noteworthy characteristics of Monascus, such as secondary metabolite biosynthetic gene clusters (BGCs), carbohydrate active enzymes, secretome, and pathogenicity that occur in this genus’ genomes. This approach moves away from a single reference genome that may not necessarily represent the species as a whole, and allows for better understanding of its metabolic versatility, ultimately leading to better management of these food microorganisms.

Monascus species are widely found in various habitats such as soil, starch, grain, dried fish, surface sediments of rivers, and roots of pine trees (Mohan Kumari, 2009). Due to the production of MonAzPs and MK, these filamentous fungi are frequently used in food and medicine. However, the taxonomy of Monascus has long been a matter of confusion. The phenotype-based identification schemes in Monascus were difficult to match with the results obtained by ITS, partial LSU and/or β-tubulin gene sequencing (Patakova, 2013), especially in the case of single-gene locus-based phylogenetics. For example, Dai et al. (2021) inferred a Neighbor-Joining tree using ITS sequences, which revealed that several isolates from M. pilosus, M. fuliginosus, M. barkeri, M. paxii, M. albidulus, M. ruber, M. purpureus, and M. fumeus were evolutionarily close in the same clade with high bootstrap values. In another clade, several strains from M. purpureus, M. rutilus, M. aurantiacus, and M. kaoliang were clustered. Within these two main clades, the strain-level division was still vague due to the low support level.

To overcome inappropriate taxonomy of fungi caused by a single locus, the genealogical concordance phylogenetic species recognition concept (GCPSR) was proposed as an empirical method for recognizing cryptic speciation (Taylor et al., 2000). GCPSR involves sequencing multiple genes that are then combined in phylogenetic analyses. In a report conducted by He et al. (2020) a phylogenetic tree was constructed using concatenated sequences of five protein genes (BenA, CaM, RPB2, pksKS, and MAT1-1) and two ribosomal RNA genes (ITS and LSU) with a total length of 6,983 bp. Strains from M. ruber and M. purpureus were clustered into different species clades, respectively, with a high Bayesian analysis/bootstrap percentage. Unfortunately, M. pilous was not included in their phylogenetic analysis. In 2017, Patakova divided the Monascus spp. into section Floridani and Section Rubri by concatenated phylogeny based on the sequences of ITS + BenA + CaM + LSU + RPB2 with Bayes/RAXML method (Patakova, 2013). M. ruber, M. pilous and M. purpureus were all in the section Rubri, and M. purpureus was located in a separate subclade from M. ruber and M. pilous with high Bayes/RAXML support. M. ruber and M. pilous were clustered into the same evolutionary branch.

Recently, modern phylogenetic analyses utilize hundreds to thousands of sites from throughout the genome, which are orders of magnitude larger than traditional sequencing datasets. Thus, the size of these datasets significantly reduces the impact of random error and data availability, making them promising for addressing historically recalcitrant nodes in the tree of life. In this study, using the Monascus assembly deposits from the NCBI genome database, a genome-level phylogenetic analysis was conducted. These 15 genomes were predominantly descended from two clades, the M. purpureus clade and the M. ruber-M. pilous clade, from either the SCOG tree according to the concatenated phylogeny based on 4,589 single copy protein orthologs or the STAG analysis with all the 5,565 single-/multiple- copy protein orthologs. Both of the phylogenomic analysis strategies were reliable with 100% support values for the two major clades. Furthermore, evidence from the nucleic acid level was proposed to support this conclusion, including analyses of average nucleotide identity and genomic collinearity. In particular, there was considerable genomic collinearity between M. ruber and M. pilous, indicating a high degree of similarity between these two species, rather than with M. purpureus. However, genome collinearity analysis, which frequently relies on genome assembly quality, makes it difficult to distinguish between strains when using fragmented genome assemblies. Hence, it is essential to employ a combination of several methods from various perspectives for evaluating the similarity between genomes.

Through identification and characterization of the BGCs of Monascus, it was found that all genomes contain the MonAzPs BGCs. However, the citrinin BGCs were only discovered in M. purpureus, while the BGCs of MK were only present in M. pilosus and M. ruber. The production of the mycotoxin citrinin, was originally described in M. purpureus and M. ruber in Blanc et al. (1995). It was subsequently shown that the M. ruber used in that study was, in fact, M. purpureus (Chen et al., 2008) because the pksCT gene for citrinin polyketide synthase was only present in M. purpureus and M. kaoliang (a synonym for M. purpureus), but not in M. pilosus, M. ruber, M. floridanus, M. sanguineus, M. barkeri, or M. lunisporas. Despite this, citrinin production in M. ruber was later demonstrated by other authors (Li et al., 2010, 2021). However, whether or not this was due to incorrect strain classification requires additional validation. Another possibility is that this strain-specific citrinin synthesis might originate from a horizontal gene transfer of the BGC among fungi, because the citrinin pathway belongs to the general pathway shared by many Penicillium, Aspergillus, and Monascus species (Higa et al., 2020).

Additionally, when compared to the reference BGCs, the BGC sequences of MonAzPs and citrinin revealed the highest degree of conservation. Interestingly, seven additional inserted genes in the MonAzPs BGCs could well differentiate M. purpureus apart from M. pilosus or M. ruber. Furthermore, the BGCs of MK were more varied in M. pilosus but were conserved in M. ruber. Although this suggested that the three strain classes might be distinguished by divergence from the MK gene cluster sequence, further sequencing evidence from more strains is necessary to support this.

According to the findings of this study, the investigated Monascus species can be classified into two groups: the M. pilosus-M. ruber clade and the M. purpureus clade. This classification may have significant implications in Monascus-related industries. Typically, commercial Monascus products are divided into two categories, those intended for the production of MonAzPs for food coloring, and those for the production of MK, a secondary metabolite to lower cholesterol and treat hypolipidemia. Among them, M. purpureus is a prominent red-colored mold species (Yang et al., 2015), but a number of strains including M. purpureus, M. pilosus, M. sanguineus and M. ruber were reported to produce MK (Wen et al., 2020). This is confusing because it is unclear whether the synthesis of MK is due to individual differences among Monascus strains or incorrect classification, as the genome of all M. purpureus strains used in this study lacked the complete gene cluster for MK biosynthesis. Therefore, further identification of these strains at the phenotypic and genotypic levels is necessary. For example, in a recent study by Higa et al. (2020) the metabolite analysis based on liquid chromatography-mass spectrometry revealed significant differences in the MonAzPs and related metabolites produced by the three species (M. pilosus, M. ruber, and M. purpureus) in liquid media, despite M. ruber had similar biosynthetic and secondary metabolite BGCs to M. pilosus.

Moreover, genome-level analysis provides insights into metabolic differences among individual strains. Comparative genomic analyses showed that the genome size of M. purpureus was smaller than that of M. pilosus/M. ruber and had undergone significant gene losses, particularly in cellular transport and maintaining homeostasis, as a result of specialized adaptation to the environment. Compared to A. oryzae, Monascus also displayed gene losses in both enzyme species and quantities involved in carbohydrate metabolism, which might be due to strain degradation resulting from prolonged domestication of the starch-rich matrix (Dufossé, 2018). Furthermore, Monascus’ genomes were predicted to contain a number of secreted proteins that could act as allergens or be involved in virulence and pathogenicity, although they exhibit interindividual variability. One such protein is CARP_ASPFU, which is present in most of the genomes and has a high identity ≥ 80%, and may play a crucial role in processing or signaling to allergens, causing rare human infections such as lung aspergillosis and mycotic keratitis. Given the importance of food safety, it is necessary to confirm whether the toxins produced by particular Monascus strains are actually produced by gene expression or just exist in the genome.

In conclusion, Monascus is a widely consumed commodity strain due to its good coloring and health care efficacy. However, its taxonomic profile remains under debate and is difficult to classify using a single locus sequence alignment and comparison, as is the case with other fungi. This study utilized phylogenomics, which draws information from comparing entire or large portions of genomes, to gain more detailed insights into evolutionary relationships, as compared to traditional phylogenetic methods that rely on a smaller number of genetic markers. Our findings clearly demonstrate differences between M. purpureus and M. pilosus/M. ruber, as well as a high degree of similarity between the genomes of M. pilosus and M. ruber. By comparing the genomes of different Monascus strains, we were able to identify differences in metabolism for environmental adaptation, carbohydrate-active enzymes, secretome, fungal pathogens, as well as in secondary metabolite gene clusters. Genome-level research provides insights into the species classification of Monascus, and as the cost of genome sequencing decreases, this high-resolution phylogenetic method will become an important means for evaluating the safety of Monascus and other edible fungi.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

ZZ, HL, and XL conceived and designed the study. ZZ, MC, PC, JL, ZM, YM, ZL, QG, and CW performed the data analysis. All authors wrote the manuscript and approved the final manuscript.