Key points

  1. 1.

    RNA-Seq and iTRAQ were used to identify AQ synthesis-related genes and enzymes.

  2. 2.

    Ubiquinone and other terpenoid-quinone biosynthesis pathway was upregulated.

  3. 3.

    Genes such as metE, CoQ2, and ubiA are important for AQ production.

Introduction

Antrodia camphorata, a unique basidiomycete, is indigenous species of Taiwan. It has been used as a traditional medicine for treat diverse discomforts such as abdominal pain, hangover, and diarrhea (Wu et al. 1997). The fruiting body of A. camphorata is rich in bioactive metabolites, but the fruiting body can only grow on the rare Cinnamomum kanehirae Hayata at an extremely slow rate. The mycelium of A. camphorata can easily be cultured on a large scale via submerged fermentation, but it has been found to be deficient in specifically bioactive metabolites, such as antroquinonol (AQ) (Lu et al. 2013). AQ, an ubiquinone derivative, was identified from the solid-state fermented mycelium of A. camphorata in 2007 (Lee et al. 2007). The remarkably therapeutic activities of AQ in many diseases have been documented, including liver and kidney diseases (Angamuthu et al. 2019), Alzheimer’s disease (Chang et al. 2015), and cancer (Chiang et al. 2010).

In a previous study, it is successful to stimulate the biosynthesis of AQ in submerged fermentation by addition of precursor coenzyme Q0 (Xia et al. 2018). Therefore, we hypothesize that the addition of coenzyme Q0 would stimulate the expression of genes involved in AQ synthesis. However, a lack research on the molecular genetics of related metabolites in A. camphorata, such as transcriptomic and proteomic resources, currently hinders such studies on AQ biosynthetic genes and enzymes. Large-scale sequencing techniques have been widely used, which improve the efficiency of understanding the differential gene expression patterns within microorganism (MacLean et al. 2009; Juan et al. 2010).

In the present study, we focused on the gene transcriptional and protein expression features of A. camphorata S-29 in submerged fermentation with and without the precursor coenzyme Q0 using second-generation sequencing on the Illumina HiSeq™ 2000 platform. We further validated functional genes associated with AQ synthesis through quantitative polymerase chain reaction (q-PCR). Two “omics” levels of analyses indicated that the formation of AQ constitutes a highly complicated and genetically programmed process that requires the participation of multiple regulators. This investigation advances our understanding of genes and enzymes involved in AQ synthesis during FM.

Materials and methods

Microorganism and cultivation

The microorganism used in this study was A. camphorata S-29, which was deposited in the China General Microbiological Culture Collection Center(CGMCC No. 9590). Seed medium was perpared according to the method of Xia et al. (2018).

The fermentation medium was prepared according to our previous report (Liu et al. 2020). For the investigation on the regulatory mechanisms, 0.3 g/L coenzyme Q0 was filtered and added into the fermentation broth. The experiments were carried out with three biological repeats. Mycelia in the fermentation broth were collected by centrifugation at 3000g for 10 min at 4 °C. The mycelia sediment was washed twice with 40 mL of potassium phosphate buffer, rapidly frozen in liquid nitrogen, and stored at − 80 °C for further use.

The absence of precursor (coenzyme Q0) during fermentation was regarded as conventional submerged fermentation (KB). Adding coenzyme Q0 during fermentation was named as coenzyme Q0-induced fermentation (FM). KB4 (or FM4), KB5 (or FM5), and KB10 (or FM10) representing the mycelia of A. camphorata S-29 were collected on days 4, 5, and 10, respectively.

Differential transcriptome analysis

The days 4, 5, and 10 samples were sequenced by Shanghai Majorbio Biopharm Technology Co., Ltd (Shanghai, China), with the Illumina HiSeq™ 2000 platform (Bai et al. 2015). The raw data were filtered by removing the low quality reads, sequences of length less than 20 bp, and reads containing adapters to obtain clean data. A rapid comparison of the sequencing data with a reference genome was performed by Tophat2 software, which can also be used to detect events such as variable shear and gene fusion (Trapnell et al. 2012). Analyses of sequence saturation, gene coverage, and duplicate reads in the transcriptome were performed using RSeQC-2.6.3 software (Vera Alvarez et al. 2019). Analysis of gene expression was performed using RSEM software (Li and Dewey 2011). Furthermore, edgeR software was used to calculate the differential expression based on gene read count data (Robinson et al. 2010). The differentially expressed genes (DEGs) were classified by the Gene Ontology (GO, http://geneontology.org/) database according to the biological process, cellular components, and molecular functions. KOBAS software (Peking University, Beijing, China) was used to test the enrichment of the DEGs; in particular, the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) pathways (Xie et al. 2011).

Protein quantification by isobaric tags for relative and absolute quantitation (iTRAQ)

Protein quantification by iTRAQ was carried out according to the method of Lu et al. (2017). The mixed peptides were sequenced by Majorbio Biopharm Technology Co., Ltd (Shanghai, China). All identified proteins were annotated by GO and KEGG. Differential proteins were screened based on differential multiples and P-values. The data were analysed by the Majorbio I-Sanger Cloud Platform (https://cloud.majorbio.com/).

The enrichment ratio (ER) was calculated as follow: ER = CN/BN, where CN is the number of proteins (or genes), which were enriched to this KEGG in the protein (or gene) set, and BN is the number of proteins (or genes), which were enriched to this KEGG among all annotated proteins (or genes).

Validation of AQ synthesis-related genes using q-PCR

Reverse transcription of cDNA was performed according to the manufacturer’s instructions (PrimeScript 1st strand cDNA synthesis kit, Takara, Japan). Q-PCR was performed using LightCycler® 96 system (Roche, Switzerland) detection. The 18S rRNA gene of A. camphorata S-29 was used as the internal standard. Primers used in this research for q-PCR (Additional file 1: Table S1) were designed by Oligo 7.0 software. Relative expression level of gene was quantified based on the 2−∆∆Ct method. It was carried out with three biological repeats. Values are given as the means ± standard deviations (n = 3).

Sequencing data

The FASTQ format raw reads were deposited to the National Center for Biotechnology Information Short Read Archive (NCBI SRA) database (accession: PRJNA543624 and PRJNA622907). The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the iProX partner repository with the dataset identifier PXD018004 (Ma et al. 2019).

Results

Identification of DEGs via transcriptome sequencing

The highest number of the expressed genes of A. camphorata S-29 was from FM5 (7478), followed by KB5 (7470), KB4 (7456), FM4 (7436), KB10 (7409), and FM10 (7356), respectively (Additional file 1: Figure S1). KEGG map shows that the highest number of genes were involved in metabolism. Over 200 genes annotated to the KEGG functional categories of translation, carbohydrate metabolism, folding, sorting and degradation, and amino acid metabolism (Additional file 1: Figure S2). Scatter plots and Venn diagram were used to analyse of DEGs between KB and FM (Fig. 1). As shown in Fig. 1a–c, the DEGs between KB and FM increased gradually from 4 to 10 days. The number of DEGs were 51, 151, and 732 on days 4, 5, and 10, respectively (Fig. 1d). A total of 35 DEGs were only expressed between KB4 and FM4, with 101 DEGs expressed between KB5 and FM5. In comparison, 677 DEGs were only expressed between KB10 and FM10.

Fig. 1
figure 1

DEGs between KB and FM. a–c a scatter plot of the DEGs during FM compared with KB on days 4, 5, and 10, respectively. Significantly upregulated or downregulated genes were labeled with red dots or blue dots, respectively. d Venn diagram of the DEGs on days 4, 5, and 10, respectively

Full size image

Functional classification of genes by GO and KEGG analysis

These genes were classified into 20 GO categories, and all GO categories were then assigned into three main categories: biological process, cellular component, and molecular function. As shown in Fig. 2, in the biological processes, 9 significant enrichment items were found. Particularly, a large number of genes involved in metabolic process, single-organism metabolic process, and cellular processes were significantly enriched in FM10 compared to KB10. In the cellular component, 7 significant enrichment items were found. In the molecular functionality, 4 significant enrichment items, especially catalytic activity and binding, were identified in FM10 compared to KB10.

Fig. 2
figure 2

GO annotation of DEGs. Each bar represent the percentage and number of genes in each GO category. Yellow, blue, and red colouring represent DEGs on days 4, 5, and 10, respectively

Full size image

In the early phase of coenzyme Q0 addition (on day 4), only the SNARE interactions in vesicular transport pathway was significantly enriched in FM4 compared to KB4 (Additional file 1: Figure S3A, P < 0.05). The arginine and proline metabolism, tryptophan metabolism, and C5-branched dibasic acid metabolism were significantly enriched in FM5 compared to KB5 (Additional file 1: Figure S3B, P < 0.05). As shown in Fig. 3, 20 metabolic pathways were significantly enriched on day 10 (P < 0.05). Arginine and proline metabolism, lysine degradation, cyanoamino acid metabolism, citrate cycle (TCA cycle), and valine, leucine and isoleucine degradation pathways were significantly enriched (P < 0.001). Amino sugar and nucleotide sugar metabolism pathway had the largest number of enriched the DEGs among 20 most significant enrichment pathways (Table 1). The genes of atoB, ALDH, and ECHS1 participated in lysine degradation, valine, leucine and isoleucine degradation, tryptophan metabolism, and fatty acid degradation pathways at the same time. It is noteworthy that the DEGs of ubiquinone and other terpenoid-quinone biosynthesis pathway were CoQ2 and CoQ5 (Table 1). Additional file 1: Table S2 shows that the CoQ2 gene, which encodes 4-hydroxybenzoate polyprenyltransferase, was upregulated, while the CoQ5 gene, which encodes 2-methoxy-6-polyprenyl-1,4-benzoquinol methylase, was downregulated on day 10.

Fig. 3
figure 3

Analysis of the DEGs between KB10 and FM10 based on the KEGG enrichment map. *P < 0.05, **P < 0.01, ***P < 0.001

Full size image
Table 1 Annotations of the KEGG pathway for DEGs on day 10
Full size table

Identification of differentially expressed proteins (DEPs) using iTRAQ

The iTRAQ is a technique in quantitative proteome that accurately measures large-fold changes in protein expression within dynamic ranges of protein abundance (Juan et al. 2010). Here we used this technique to perform proteomic analyses of the key enzymes involved in AQ biosynthetic pathways. The iTRAQ data show that 44,751 unique peptides, which were assembled into 3987 proteins (Additional file 1: Figure S4). All proteins were then assigned into three main GO categories: biological process (BP), cellular component (CC), and molecular function (MF). As shown in Additional file 1: Figure S5, a large number of proteins were assigned into the GO categories of biological process and cellular component. Particularly, the proteins of more than 1000 involved in metabolic process, cellular process, single-organism metabolic process, cell, cell part, organelle, catalytic activity, and binding. In total, 3914 proteins did not change significantly, 24 proteins were downregulated, and 49 proteins were upregulated in FM10 compared to KB10 (Additional file 1: Figure S6).

GO analysis was also used to further elucidate the putative functions of upregulated DEPs on day 10 (Fig. 4). According to the GO analysis, the upregulated DEPs in FM were mainly grouped into the following functions (number of proteins): sulfur amino acid metabolic process (4), one-carbon metabolic process (3), methionine metabolic process (3), sulfur compound biosynthetic process (4), serine family amino acid metabolic process (3), homoserine metabolic process (2), sulfur amino acid biosynthetic process (3), 5-methyltetrahydropteroyltri-l-glutamate-dependent methyltransferase activity (2), 5-methyltetrahydropteroyltriglutamate-homocysteine S-methyltransferase activity (2), and S-methyltransferase activity (2) (P < 0.001, Fig. 4 and Additional file 1: Table S3).

Fig. 4
figure 4

GO annotation of upregulated proteins on day 10. *P < 0.05, **P < 0.01

Full size image

KEGG pathway enrichment analysis was also performed to investigate the enriched pathways of DEPs on day 10. Seven KEGG pathways were significantly upregulated in FM: cysteine and methionine metabolism (P < 0.001), selenocompound metabolism, ubiquinone and other terpenoid-quinone biosynthesis (P < 0.01), betalain biosynthesis, glutathione metabolism, glycine, serine and threonine metabolism, and cyanoamino acid metabolism (P < 0.05, Fig. 5a). One KEGG pathway was significantly downregulated: glycolysis/gluconeogenesis (P < 0.001, Fig. 5b).

Fig. 5
figure 5

Enriched pathways of upregulated (a) and downregulated (b) proteins on day 10 based on KEGG annotation. *P < 0.05, **P < 0.01, ***P < 0.001

Full size image

Among these upregulated DEPs in FM, we further identified key enzymes that might be related to the biosynthesis of AQ (Table 2 and Additional file 1: Table S4). Five enzymes (ec:2.5.1.6, ec:2.5.1.18, ec:2.5.1.39, ec:2.5.1.21, and ec:2.5.1.54) are relate to transfer alkyl or aryl groups. Two enzymes (ec:2.1.1.14 and ec:2.1.1.41) are methyltransferases. Seventeen of them are oxidoreductases. The ec:2.1.2.1 and ec:2.1.2.10 are relate to hydroxymethyltransferase and aminomethyltransferase, respectively. Four enzymes (ec:2.3.3.1, ec:2.3.3.8, ec:2.3.3.13, and ec:6.2.1.1) can generate acetyl-CoA during the catalytic reaction.

Table 2 Critical upregulated proteins based on proteomic analysis (10-day samples)
Full size table

Integrated analyses of transcriptomic and proteomic datasets of A. camphorata S-29

The correlation between the proteome and transcriptome of A. camphorata S-29 was demonstrated by scatter plot analysis (Fig. 6). The distribution of the dots with different colors shows that a large number (3589) of mRNAs and proteins did not exhibit changes in expression (black dots, Fig. 6 and Additional file 1: Table S5). On the other hand, a different trend was observed between protein and transcription levels. For instance, 253 genes were upregulated at the transcriptional level but unchanged at the protein level (sky blue). Moreover, 72 genes were downregulated at the transcriptional level but unchanged at the protein level (slate blue). Due to the regulations of transcription and translation, the mRNA expression levels are not always consistent with protein levels (Washburn et al. 2003). It is noteworthy that 18 unigenes were upregulated at the transcriptional level and the protein level (red).

Fig. 6
figure 6

Correlation analysis of proteomic and transcriptomic data (10-day samples). Differentially expressed mRNAs or proteins are represented by solid dots

Full size image

For the comparative analysis of proteomic and transcriptomic data, we focused on the 9 unigenes that were both upregulated in transcriptomic and proteomic sequencing (Table 3). The ubiA and CoQ2 genes simultaneously participate in the formation of 4-hydroxybenzoate polyprenyltransferase. Of note, among the genes for synthesized benzenoids, PKS63787 gene was the DEG that exhibited the greatest upregulation. Ac-8 cytochrome P450 (CYP) monooxygenase synthesis-related gene P450 and methyltransferase-related gene metE were significantly upregulated in the proteome and transcriptome data. E1.14.13.1, E1.14.13.7, and ADK were associated with the electron transport respiratory chain and energy metabolism.

Table 3 Significantly upregulated genes and proteins on day 10
Full size table

Validation of gene expression data using q-PCR

The mRNA expression levels of 5 genes in potential AQ biosynthesis on days 4, 5, 8, and 10 were studied (Fig. 7). The mRNA expression levels of 5 genes (ubiA, CoQ2, PKS63787, P450, and metE) in KB and FM were not significantly increased on day 4 (Fig. 7a, P > 0.05). Compared with the KB, the expression levels of the ubiA and P450 genes increased significantly in FM on day 5 (Fig. 7b, P < 0.05). The expression levels of the ubiA, CoQ2, PKS63787, P450, and metE genes were significantly higher in FM than in KB at 8 and 10 days (Fig. 7c, d, P < 0.05). Therefore, the gene expression levels with the transcriptome sequencing analyses were consistent with those of the q-PCR analyses (Table 3 and Additional file 1: Table S6). The production of AQ increased significantly on 5–10 days compared with 4 and 5 days (Additional file 1: Figure S7). The AQ synthesis-related genes expression levels are in line with its production.

Fig. 7
figure 7

The q-PCR analysis of mRNA expression of genes potentially related to the biosynthesis of AQ on days 4 (a), 5 (b), 8 (c), and 10 (d). Means with different letters differed significantly (P < 0.05). The experiments were carried out with three biological repeats. Values are given as the means ± standard deviations (n = 3)

Full size image

Discussion

AQ is considered one of the most bioactive metabolites in A. camphorata (Lu et al. 2014). The addition of the AQ precursor (coenzyme Q0, orsellinic acid, etc.) via submerged fermentation is considered the most efficient method for the industrial production of AQ. In this study, we used transcriptomic and proteomic data to reveal AQ synthesis-related genes and enzymes during FM. GO enrichment and KEGG pathway analyses of transcriptomic data provided insights into the secondary metabolite biosynthesis pathways. The GO assignment system was used to obtain functional information for the DEGs, further assisting in understanding the distribution of DEG functions at a macro level. Given that fewer differential genes were noted between KB and FM groups in the early stage of adding coenzyme Q0, the number of genes on days 4 and 5 annotated in the GO categories were significantly less than that on day 10 (Figs. 1 and 2). Numerous genes involved in metabolic processes were significantly enriched in FM compared to KB on day 10. This finding indicates that the addition of the precursor coenzyme Q0 significantly increased the synthesis of metabolites.

KEGG pathway analysis provides classifications that are valuable for studying the complex biological functions of genes. KEGG analysis of transcriptomic data showed that the pathway of fatty acid degradation was significantly enriched in FM (Fig. 3). The ECHS1 gene, which encodes enoyl-CoA hydratase, and the atoB gene, which encodes acetyl-CoA C-acetyltransferase, were upregulated. In fungi, the key metabolic pathway for fatty acid degradation is beta-oxidation. Fatty acids are broken down in a repeating process involving four steps in mitochondria, and each round of processing removes two carbons, in the form of acetyl-CoA, from the fatty acid chain (Shen and Burger 2009). In addition, pyruvate is oxidized to acetyl-CoA in mitochondria during pyruvate metabolism. ACACA, E2.3.3.10, OXCT, IVD, and PAAH genes, which are associated with acetyl-CoA, were also upregulated. The formation of AQ may be closely related to the ubiquinone biosynthesis pathway (Hu et al. 2016). The initial biosynthetic pathway of ubiquinone involves the conversion of three acetyl-coA units into farnesyl diphosphate (FPP) (Ericsson et al. 1992). Therefore, acetyl-CoA may enter either the mevalonate pathway or polyketide pathway, which form the farnesyl tail precursor or ring precursor of AQ in A. camphorata S-29, respectively (Chou et al. 2017). The CoQ2 and ubiA genes encode p-hydroxybenzoate polyprenyl transferase, which is required in the biosynthetic pathway of ubiquinone (Forsgren et al. 2004; Liu et al. 2017). P-hydroxybenzoate polyprenyl transferase, which is also referred to as the ‘CoQ2 enzyme’, mediates the second step in the final reaction sequence of coenzyme Q biosynthesis, namely, the condensation of the polyisoprenoid side chain with 4-p-hydroxybenzoate (Forsgren et al. 2004). The coupling of the aromatic substrate and isoprenoid chain is presumed to be the rate-limiting step in ubiquinone synthesis. Chou et al. (2019) proposed that the farnesylation of coenzyme Q0 at C-3 forms 5-demethoxy-coenzyme Q3 or coenzyme Q3, which is further reduced to form AQ. Therefore, we hypothesize that overexpression of CoQ2 and ubiA genes promoted the linkage of coenzyme Q0 and the polyisoprene side chain, thus increasing AQ synthesis during FM.

The iTRAQ has the advantages of high coverage, accuracy, and sensitivity. Proteomic analyses of A. camphorata S-29 using the iTRAQ technique could provide useful information for understanding the AQ biosynthesis pathway. Based on the numerous unique proteins identified in GO functional categories, the largest categories for the functional group were sulphur amino acid and sulphur compound biosynthetic processes (Fig. 4). In microorganisms, sulphate is used in the synthesis of organic sulphur metabolites, mostly cysteine, methionine, and S-adenosylmethionine. The sulphur amino acids (SAA) include methionine and cysteine (Glenda et al. 2012). Methionine is metabolized via three major metabolic pathways: transmethylation, remethylation, and transsulfuration (Vigneaud et al. 1994). Fungi utilize the methyl group of methionine for biological methylation. Methionine serves as the major methyl group donor. The expression of S-adenosylmethionine synthetase is significantly increased during FM of A. camphorata S-29 compared with control (Fig. 4). S-adenosylmethionine synthetase is the transferring enzyme with S-methyltransferase activity. In addition, upregulated metE genes could encode an enzyme (ec:2.1.1.14) with methyltransferase activity. According to the reported results, it is important to have a sufficient supply of methyl groups for AQ structure formation (Wijayasinghe et al. 2014; Hu et al. 2016). These results suggest that S-adenosylmethionine synthetase may play a vital role in the biosynthesis of AQ.

According to results from Yu et al. (2016), the biosynthesis of aromatic compounds in A. camphorata occurs via the polyketide pathway. PKS63787 plays a critical role in the biosynthesis of benzoquinones and benzenoids. Three metabolites, including benzenoids, benzoquinones, and AQ, are related to PKS63787. Yu et al. (2017) used phylogenetic analysis to demonstrate that PKS63787 likely encodes an orsellinic acid synthase. The biosynthetic route of AQ follows a synthetic sequence similar to that of yeast coenzyme Q6. The pathway starts with the farnesylation of orsellinic acid to form 3-farnesyl-orsellinic acid, which undergoes further ring modification to form 5-demethoxy-coenzyme Q3, coenzyme Q3, and finally AQ (Chou et al. 2017). Table 3 and Additional file 1: Table S2 show that PKS63787 gene expression is significantly increased in FM on day 10 (P < 0.001). After adding coenzyme Q0 on day 4, its content gradually decreased from days 4–10. Coenzyme Q0 was completely utilized after day 10 (Additional file 1: Figure S7). Therefore, we hypothesized that excess polyisoprene side chain induced the overexpression of the PKS63787 gene to synthesize the cyclohexanone moiety of AQ. However, further investigations are required.

CYP genes are ubiquitously distributed genes in the genomes of almost all organisms. CYP genes encode a superfamily of haeme-containing monooxygenases that metabolize a wide variety of endogenous and xenobiotic compounds (Bernhardt 2006). Many fungal species have evolved complex and diverse CYP systems that play an important role in primary and secondary metabolism, allowing these fungi to exhibit unique and superior metabolic functions (Brink et al. 1998). CYP enzyme systems of fungi also include a wide range of catalysed reactions, such as reduction, heteroatom oxygenation, epoxidation, dealkylation, aromatic hydroxylation, and carbon hydroxylation (Brink et al. 1998; Ling et al. 2017). Hsu et al. (2011) identified 10 CYP genes (ac-1 to ac-10) from A. camphorata. According to phylogenetic tree of fungal CYP amino acid sequences, ac-8 CYP was assigned to a new CYP family. Since ac-8 CYP did not exhibit close evolutionary relationships with known CYPs, the function of ac-8 CYP is unknown. In the transcriptomic and proteomic analysis, the expression of ac-8 CYP monooxygenase increased significantly on day 10. In addition, mRNA expression levels of P450 were significantly increased in the FM compared to the KB at 5, 8 and 10 days (Fig. 6b–d). Therefore, we hypothesized that ac-8 CYP monooxygenase plays an important role in AQ synthesis.

According to several reports, coenzyme Q0 is associated with the electron transport respiratory chain (Redfeam and Burgos 1966). Several DEGs and DEPs that might be closely related to the electron transport respiratory chain were noted in this study. E1.14.13.1, E1.14.13.7, and ADK are three energy-related genes that exhibit increasing levels of expression upon the addition of coenzyme Q0. One of the main functions of these genes is to provide energy for cells and electron transport. Farnesylation is a chemical reaction that involves covalent bonding, requiring higher energy. Overexpression of these genes may facilitate oxidative phosphorylation and accelerate electron transport through cytochrome c within the inner mitochondrial membrane, and finally activate ATP synthase to synthesise the energy carrier ATP. In addition to farnesylation steps, the conversion step of coenzyme Q3 to AQ involves the elimination of one double bond and the reduction of a keto group. These chemical reactions require several protons that are derived from NADH. When NADH is transformed into NAD+, this reaction releases two electrons, which can be accepted by electrophile substances, such as oxygen (Wei et al. 2020). This reaction may be beneficial for the conversion of coenzyme Q3 to AQ.

In summary, transcriptomic and proteomic analyses showed that specific genes and enzymes associated with ubiquinone and other terpenoid-quinone biosynthesis pathways as well as various transcription factors are involved in AQ formation. Two “omics” levels of analyses identified several important candidate genes and enzymes that regulate or participate in AQ biosynthesis (Additional file 1: Figure S8). This study is beneficial to analyse AQ biosynthesis pathways and provide strategies for targeted increases in AQ production by genetic means in the future.