Lignan is one of the major functional compounds that are mainly found in plants^[1]. Pinoresinol (pin) is a simple lignan in nature and has been found to possess multiple functions, such as anti-cancer^[2], anti-tumor^[3], antifungal^[4] and anti-inflammatory^[5] activities, as well as prevention of hormone-dependent diseases, such as cardiovascular diseases, hyperlipidemia and breast cancer^[6-8]. However, the natural production of lignans, including pin is very low in nature. The identification of lignan biosynthesis pathway and key genes would provide essential information for large-scale production of lignan compounds using microbial fermentation via genetically modified microorganisms. The lignan biosynthesis pathway has been well known as the phenylpropanoid pathway that is widely found in plants^[9], but not currently found in microorganisms^[10-11].

In plants, the phenylpropanoid pathway plays a central role for the biosynthesis of a diverse group of compounds, such as lignin, coumarins, stilbenes, flavonoids, isoflavonoids, and lignans^[12]. In this pathway (Figure 1), D-Erythrose-4P (E4P) from glycolysis and phosphoenolpyruvate (PEP) from pentose phosphate pathway are conversed to chorismate via the shikimate pathway, and then produce phenylalanine, which is transformed into an activated cinnamic acid derivative via the actions of phenylalanine ammonia lyase (PAL), trans-cinnamate 4-monooxygenase (CYP73A), 4-coumarate-CoA ligase (4CL), and specific branch pathways^[13-14].

Many endophytic fungi from plants have been found to possess the capability to produce plant derived compounds that have anti-tumor, anti-inflammatory, antibacterial, and antifungal activities^[15-16]. In recent years, endophytic fungi have attracted great attention of scientists due to their potential application in medicine, agriculture, and other industries. Lignans and their glycosides are the major effective components of Eucommia ulmoide Oliv. (E. lmoides)^[17]. To date, 32 lignans, including pin, pinoresinol-4-O-β-D- glucopyranoside (PMG), and pinoresinol diglucoside (PDG), have been identified and isolated from E. ulmoides^[18]. Several approaches have also been developed to produce pin and its derivatives by chemical or enzymatic methods^[19-21]. Phomopsis sp. XP-8, an endophytic fungus isolated from the bark of Tu-Chung (E. ulmoides), was found to produce pin and its monoglucoside (PMG) and diglucoside (PDG)^[22], similar as the products found in E. ulmoides^[23], providing a potential method to produce pin and its glucoside derivatives by microbial fermentation. Zhang^[24] identified that Phomopsis sp. XP-8 could accumulate the intermediate compounds (cinnamic acid and p-coumaric acid) and had key enzymatic activities (PAL, CYP73A, 4CL) in the phenylpropanoid pathway. However, an understanding of pin and its glucoside derivatives biosynthesis at the genetic level should facilitate greater pin and its glucoside derivatives production^[25].

In order to identify the phenylpropanoid pathway and pin and its glucosides biosynthesis pathway in Phomopsis sp. XP-8 at the gene level, and to help identify more significant genes as well. De novo genome sequencing, feature analysis, were explored in this study.

1 Materials and methods 1.1 Fungal strains and genomic DNA preparation

Phomopsis sp. XP-8, an endophytic fungus previously isolated from the bark of Tu-Chung (Eucommia ulmoides Oliv.) and maintained at the China Center for Type Culture Collection, (Wuhan, China) (code: CCTCC M209291), was selected for genomic analyses in the study. After Phomopsis sp. XP-8 cultivation for 120 h under 28 ℃ and 180 r/min in PDB (containing 200 g potato, 20 g dextrose, and 1000 mL tap water), highly purified fungal DNA extraction was resulted in the construction of libraries using the cetyltrimethyl ammoniumbromide (CTAB) method^[26]. DNA quality and quantity were evaluated using standard 0.8% agarose gel electrophoresis, and a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, USA).

1.2 Genome Sequencing and Assembly

After whole genome DNA extraction, DNA sequencing libraries with 600 bp and 3 kb inserts were constructed using Nextera Mate Pair Library Preparation Kit (Illumina, San Diego, CA, USA), following the standard Illumina protocol. The generated 4418637781 reads were quality control-filtered; that is, low-quality reads were removed. High-quality reads were obtained and then assembled by the de novo method using Velvet v1.2.10^[27]. Final assembly was evaluated for completeness using CEGMA v2.5^[28] and BUSCO v2.0^[29], and the resulting assembled sequences were used for gene model prediction, functional annotation, and downstream information analyses. The whole genome sequencing of Phomopsis sp. XP-8 was performed using the Illumina Hiseq 2500 with 96.8-fold coverage. The sequenced reads were assembled into 330 scaffolds (≥1000 bp).

1.3 Gene Prediction and Functional Annotation

The genes in the Phomopsis sp. XP-8 genome were predicted using various software. Repetitive sequences were de novo predicted using LTR FINDER^[30], MITEHunter^[31], RepeatScout^[32] and PILERDF^[33], and then classified according to PASTEClassifier^[34]. RepeatMasker^[35]was subsequently used against the Repbase^[36] TE library and the de novo repeat library to annotate repeats in the Phomopsis sp. XP-8 genome. miRNAs were predicted by BLASTN of genomic sequence slices against miRBase sequences. tRNAs were annotated using tRNAscan-SEv1.3.1^[37]. rRNA was annotated by aligning the genomic sequence against the Rfam database by using BLASTn with an e-value ≤1e^-5. Protein-coding genes were predicted using Augustus (v3.1 http://bioinf.uni-greifswald.de/augustus/). The predicted CDSs were searched by Gapped BLAST and PSI-BLAST^[38] against the NCBI non-redundant protein (Nr) database, Pfam and Swiss-Prot database, Gene Ontology (GO)^[39], Cluster of Orthologous Groups^[40](COGs), Kyoto Encyclopedia of Genes and Genomes (KEGG)^[41]metabolic pathways. The e-value thresholds were all set at ≤1e^-5. The antibiotics and secondary metabolite analysis shell (antiSMASH) automatic pipeline (http://antismash.secondarymetabolites.org/) was used to predict gene clusters and the core structure of secondary metabolites in Phomopsis sp. XP-8.

1.4 Gene Family Analyses

Orthologous gene clusters of Phomopsis sp. XP-8 and other reference fungal strains were determined using OrthoMCL v1.4^[42] with a default MCL inflation index of 1.5. The result of the OrthoMCL analysis was used to plot Venn diagrams. Protein alignments, of each single-copy ortholog family, were performed using MUSCLE^[43]. All the poorly aligned regions were removed. The maximum likelihood phylogenetic tree of the six species included in this analysis was created based using phyML v3.0^[44]with bootstrapping for 1000 replicates.

1.5 Analysis of Pinoresinol and its Glycosides Biosynthetic Genes

Despite the fact that pin and its glycosides biosynthesis pathway have remained obscure in E. ulmoides and Phomopsis sp. XP-8, our previous analysis indicates that detectable intermediates and corresponding enzyme activities exist, which are similar to those of the plant phenylpropanoid pathway in Phomopsis sp. XP-8^[24]. To reveal the genetic basis of the pin and its glycosides biosynthesis pathway in Phomopsis sp. XP-8, functional annotation of genes using reference plant phenylpropanoid pathway genes was performed according to keyword searches for standard gene names and synonyms or abbreviations. Each search result was further confirmed with BLAST searches.

1.6 Construction of Phylogenetic Trees of Proteins Associated with Pinoresinol Formation and UDP Glucosyltransferase

Dirigent proteins (DIR) or chloroperoxidase (CPO) and UDP Glucosyltransferase (UGT) are the two key enzymes involved in the pin and its glycoside derivatives biosynthetic pathways. However, comparison of CPO, DIR and UGT gene sequences in the genome examined with the sequences in the literature revealed no significant similarity. Therefore, protein searches were performed using the DIR, CPO, and UGT genes in the Phomopsis sp. XP-8 genome by BLASTp against GenBank databases from NCBI. Two phylogenetic trees were constructed using Molecular Evolutionary Genetics Analysis (MEGA) software version 7.0.18. Protein sequence alignment (multiple alignments) was initially conducted by applying ClustalW (Gap opening penalty = 10, Gap extension penalty = 0.1) with default gap penalties. Phylogenetic trees were finally calculated based on the resulting alignments from ClustalW by using 1000 bootstrap steps.

1.7 Nucleotide Sequence Accession Numbers

The sequencing project was registered in the National Center for Biotechnology Information (NCBI) BioProject database with accession number PRJNA376863. Raw sequencing data, assembly, and annotations of Phomopsis sp. XP-8 genome have been deposited in Sequence Read Archive (SRA) database with accession number SRP100812.

2 Results 2.1 Genome Sequencing and General Features

The assembled genome was generated into a total of 1617 contigs with a total length of 54384533 bp and N50 values of 59580 bp. The contigs were assembled into 330 scaffolds (≥1000 bp) with a total length of 55225560 bp and N50 values of 318484 bp. The genome features are shown in Table 1. The overall G+C% content is 53.50%, which is higher than the average G+C% content for fungi, which is 48.96%^[45]. High G+C% content indicates greater thermal stability^[46]. CEGMA analysis showed that 244 of 248 genes were fully annotated (98.39% completeness), and 246 of 248 genes met the criteria for partial annotation (99.19% completeness). BUSCO analysis indicated that 85% and 10% of 290 expected Ascomycota genes were identified as complete and fragmented; meanwhile, only 3% were missed in the assembly. Both assessment methods indicated that the Phomopsis sp. XP-8 genome assembly was fully complete. The genome contains 310 non-coding RNA (ncRNAs) (144 tRNA, 49 rRNA and 117 miRNA). Overall, 0.46 Mb of non-redundant and repetitive sequences were identified, comprising approximately 0.83% of the genome assembly.

Gene prediction was performed using Augustus and yielded 17094 protein-coding genes. For all predicted genes, the functions of 16025 genes (93.75%) were successfully annotated in the different databases. Including 5876 genes > 300 nt and 9766 genes > 1000 nt. The functional annotation of each database is detailed below (Table 2). 15921 genes (93.14%) were aligned using the Nr database, 11893 genes (69.57%) had hits in the Swiss-Prot database, 2907 genes (17.01%) were mapped in the KEGG database, 6239 genes (36.50%) were classified in the GO database, 6788 genes (39.71%) were annotated against the COG database, the e-value thresholds were all set at 1e^-5.

To further demonstrate the functional distribution of all genes, GO, COG and KEGG analyses were used for function prediction and classification. By GO analysis, the gene set was categorized into 45 functional groups (Figure 2). A total of 6319 genes were mapped to GO terms. The assignments were given to cellular components (17.14%), molecular functions (44.67%), and biological processes (38.19%). Among all of the GO terms, the vast majority were related to catalytic activity, metabolic processes, cell process, single-organism and binding. Within the cellular component category, cell part, cell and membrane occupied the majority of the genes. In the molecular function category, catalytic activity and binding proteins were dominant. Under the biological process category, metabolic process, single-organism process and cellular process represented the majorities, indicating that these genes were involved in some important metabolic activities in Phomopsis sp. XP-8.

Further, all genes were searched against the COG database. Overall, 6788 sequences were assigned COG classifications. The COG-annotated putative proteins were classified into 25 functional categories (Figure 3). Among the 25 COG categories, the cluster for general function prediction only represented the largest one, followed by amino acid transport and metabolism, inorganic ion transport and metabolism. Notably, 905 (13.33%) genes were assigned to secondary metabolite biosynthesis, transport and catabolism, and 1074 (15.82%) were assigned to carbohydrate transport and metabolism. To identify the biological pathways in Phomopsis sp. XP-8, the annotated genes were mapped to the KEGG database (Figure 4). A total of 2907 genes were assigned to 116 KEGG pathways and classified into 5 main KEGG categories and 19 subcategories. Among them, the most represented were in the "metabolism pathways", followed by "biosynthesis of secondary metabolites". We focused considerably our attention on the biosynthesis of secondary metabolites (Table 3). 37 genes were identified in the phenylpropanoid pathway to verify that the genes of the pathway are related to lignan biosynthesis. Furthermore, 16 genes were mapped onto terpenoid backbone biosynthesis; 18 for stilbenoid, diaryl heptanoid and gingerol biosynthesis; 12 for plant- pathogen interaction; 16 for flavonoid biosynthesis; 4 for flavone and flavonol biosynthesis, and so on. These annotations would be anextremely valuable resource for further research on gene identification, mechanism for regulating gene expression, structures, functions, and metabolism pathways in Phomopsis sp. XP-8.

2.2 Orthology and Phylogenetic Analysis

Using OrthoMCL, a total of 11818 gene families were identified in the Phomopsis sp. XP-8 genome. Among these gene families, 64 were unique. Compared with other five sequenced Ascomycota genomes, 12635 orthologous genes in Phomopsis sp. XP-8 are shared, whereas 1140 genes are unique for this fungus. 163 genes belong to unique gene families. Furthermore, the majority of 1140 species-unique genes are hypothetical proteins that do not contain recognizable Pfam domains, indicating unique gene functions. The annotated functions of species-unique gene proteins are mainly used in the biosynthesis of amino acids, aromatic amino acids metabolism, secondary metabolism and carbon metabolism processes. A Venn diagram shows that the largest cluster group has 5626 gene families each containing at least one gene from each of the six species. A phylogenomic tree that was generated based on 5042 single-copy orthologs revealed that Phomopsis sp. XP-8 and Diaporth ampelina (D. ampelina) are closely related to the Valsa mali (V. mali), which is consistent with a previous report^[47] (Figure 5).

2.3 Considerable Biosynthetic Capabilities of Secondary Metabolites in Phomopsis sp. XP-8

Endophytic fungi may synthesize various bioactive secondary metabolites^[48]. Phomopsis sp. XP-8 has been found to produce secondary metabolites^[49], including pin and its glucosides^[22].

This finding prompted us to determine the molecular basis of this production by genome sequencing. We detected 80 secondary metabolite gene clusters in Phomopsis sp. XP-8 by submitting the whole genome data to the antiSMASH automatic pipeline. The resulting enzymes included 30 Type Ⅰ polyketide synthases (T1PKS), 2 Type Ⅲ polyketide synthases (T3PKS), 16 non-ribosomal peptide synthases (NRPS), 7 terpene synthase (TS) group, 4 indolel groups, 2 lantipeptide group, 1 siderophore- indole group, 1 NRPS-T1PKS hybrid, 2 T1PKS- NRPS hybrids, 1 NRPS-T1PKS-indole hybrids, and 13 other groups. Among the 80 gene clusters, 53 were identified as PKS and NRPS (> 66% total). In addition, the core structure of secondary metabolites of the Phomopsis sp. XP-8 strain was also analyzed. 14 synthetic compound structures were predicted. These data, along with the obtained results in the numerous novel secondary metabolites, indicated a huge potential for the production of secondary metabolites of Phomopsis sp. XP-8.

2.4 Annotation of Pinoresinol and its Glycosides Biosynthetic Genes

In plants, lignan is normally produced by phenylpropanoid biosynthesis (ko00940). The genes involved in lignan synthesis have been isolated and characterized, such as PAL (DQ115905), C4H (GU014562), 4CL (GU937875), CCR (GQ872418), CAD (GU937874), C3H (JF826963), and CCoAOMT (DQ115905). The transcriptome database, by de novo techniques, provides a pool of candidate genes involved in the biosynthesis of phenylpropanoid in Isatis indigotica^[50-52]. In the present study, we used the genes involved in the phenylpropanoid biosynthetic pathway in plants as a reference, for which 86 genes that encoded up to 11 enzymes were found be involved in this biosynthetic pathway for pin and its glycosides, including PAL, CYP73A, 4CL, shikimate O-hydroxycinnamoyltransferase (HCT), coumaroylquinate (coumaroylshikimate) 3′-monooxygenase (CYP98A), caffeoyl-CoA-O- methyltransferase (CCoAOMT), caffeic acid 3-O-methyltransferase (COMT), cinnamoyl-CoA reductase (CCR), cinnamyl alcohol dehydrogenase (CAD), and CPO, DIR and UGT. All of these enzymes were successfully annotated and had more than one copy but PAL (Figure 6).

2.5 Phylogenetic Analysis of Pinoresinol-Forming Associated genes

In the current study, the phylogenetic tree of the pin-forming associated proteins was generated from alignments of the full-length amino acid sequences from previously published data in GenBank and the proteins of the genome selected from Phomopsis sp. XP-8. As shown in Figure 7, Sesamum indicum (AAT11124.1), Forsythia x intermedia (AAF25357.1), Dysosma tsayuensis (ADD70247.1), Podophyllum peltatum (AAK38666.1), and Dysosma pleiantha (AIF79763.1) were the five pin-forming proteins in GenBank. These proteins had the same large clade, along with Gglean003546.1, which was annotated as the CPO gene in Phomopsis sp. XP-8. Gglean006481.1, Gglean011334.1, and Gglean016114.1 were three CPO genes with the same clade near the aforementioned clade. Gglean016552.1 and Gglean001822.1 comprised the two annotated dirigent protein genes and formed a separate large clade. Analyses indicated that in Phomopsis sp. XP-8, the pin was derived from the coupling of two coniferyl alcohol molecules that were more likely catalyzed by CPO than DIR. Although no reports exist that show that pin is catalyzed by the DIR or CPO in E. ulmoides and Phomopsis sp. XP-8, pin is formed by DIR in other plants.

2.6 Phylogenetic Analysis of Annotated UDP Glucosyltransferase

The phylogenetic analysis of UGTs clearly showed plants, fungi and Saccharomyces formed three separate and larger clades (Figure 8). Gglean016718.1 and several fungi (AQX36236.1, GAN01566.1, EPB92999.1, AAD29571.1), Saccharomyces (AAB67475.1 and AAD29570.1) were in the same large clade; Gglean004896.1, Gglean008180.1, Gglean002755.1, Gglean000988.1, Gglean014280.1, several plants, Eucommia ulmoides (AHX74090.1), Forsythia x intermedia (BAI65911.1 and BAI65912.1), and Sesamum alatum (BAF96581.1) were in another large clade, whereas Gglean007115.1, Gglean008006.1, Gglean013194.1, and Gglean011606.1 formed a separate clade. This analysis showed that Gglean016718.1 was closer to the previously mentioned fungi and Saccharomyces, which could be explored as a candidate gene for a study on glycosylation of pin glucoside derivatives.

The conserved domains database in NCBI was used to align the amino acid sequence of Gglean016718.1. The results showed that Gglean016718.1 was highly homologous to GT1-glycosyltransferase (GT1_Gtf_like, PSSM: cd03784), and the conserved domains were screened. There were 14 core amino acids (active site) involved, six of which were located at the TDP-binding site, one at the acceptor substrate-binding pocket. So, we infer the amino acid backbone of Gglean016718.1 is consistent with the characteristics of NDP- glycosyltransferase, and can be folded to form the corresponding active catalytic center and substrate binding sites.

Gglean016718.1 had multiple sequences aligned with the corresponding sequences of glycosyltransferases in fungal and Saccharomyces proteins (Figure 9). Less than 50% identities were determined between Gglean016718.1 and Mucor hiemalis (AQX36236.1, 45%), Mucor ambiguous (GAN01566.1, 42%), Mucor circinelloides f. circinelloides 1006PhL Mucor ambiguous (EPB92999.1, 42%), Komagataella pastoris (AAD29570.1, 48%), Candida albicans (AAD29571.1, 46%), and Saccharomyces cerevisiae Ylr189cp (AAB67475.1, 46%). The amino acids at the positions of 1125-1174, 1190-1239 and 1351-1400 of Gglean016718.1 were three protein motifs found in the other six glycosyltransferases. These results suggested that Gglean016718.1 has the same biological function as that of fungal glycosyltransferases in the biosynthesis of phenols and sterol compounds.

3 Discussion

Phomopsis was proposed to be an asexual state of Diaporthe^[47]. As the phytopathogenic and endophytic species, Phomopsis have an extremely broad range of plant hosts and are widely distributed, including angiosperms, gymnosperms, pteridophytes, and bryophytes^[53]. Phomopsis sp. are non-systematic fungi lacking host specificity and have been used as model endophytes in numerous studies^[54]. As one class of the most widely distributed endophytic fungi, Phomopsis sp. has recently drawn increasing attention from the research community because of their ability to produce secondary metabolites with various structures and diverse biological activities^{[49, 55-58]}.

In plant, proteins associated with pin formation have been reported, which are classified as DIR. Related genes have been isolated and characterized in Sesamum indicum (AAT11124.1), Forsythia x intermedia (AAF25357.1), Dysosma tsayuensis (ADD70247.1), Podophyllum peltatum (AAK38666.1), and Dysosma pleiantha (AIF79763.1). However, Charles^[19] in 1976 investigated a method of obtaining pin via the utilization of the precursor coniferyl alcohol and the catalyzation activity of CPO produced by Caldariomyces fumago. In plants, glycosylation represents the last step in the biosynthesis of numerous natural compounds. Glycosyltransferases family 1 (GT1), often referred to as a UGT, that is considerably important for chemical stability, water-soluble enhancement, inactivation, and detoxification of natural products. These enzymes generally catalyze the transfer of the glycosyl group from nucleoside diphosphate-activated sugars (UDP-sugars) to a diverse array of substrates, such as secondary metabolites or xenobiotics. The composition of lignans is enriched by glycoslation catalyzed by multiple UGTs. UGT71A18 serves as a (+)-pinoresinol glucosyltransferase in the Forsythia suspension culture that can produce PMG on HPLC analysis^[59].

Pin, PDG and PMG have been detected as bioactive secondary metabolites of Phomopsis sp. XP-8^[22]. DIR or CPO and UGTs are the two key enzymes involved in the pin and its glycoside biosynthesis. 2 DIR, 7 CPO and 14 UGT genes were annotated in the Phomopsis sp. XP-8 genome. To further screen these genes, annotation and identification of genes that code for carbohydrate-active enzymes (CAZymes) in the examined genome were performed using the CAZymes Analysis Toolkit (http://mothra.ornl.gov/cgi-bin/cat/cat.cgi). Combining genomic annotation and phylogenetic analysis, Gglean003546.1 and Gglean016718.1 will be further studied in the biosynthesis of pin and its glycosides as main candidate genes. This study provides genetic basis for the further gene cloning, gene expression, and functional verification, and may also identify more significant genes.

In addition, according to what has been found in plants, pin can be converted to many functional compounds, such as lariciresinol and secoisolariciresinol by pinoresinol/lariciresinol reductase and subsequently to matairesinol via secoisolariciresinol dehydrogenase (SIRD). In this study, 8 pinoresinol reductase genes, 5 pinoresinol-lariciresinol reductase genes, and 5 SIRD genes were found in Phomopsis sp. XP-8, indicating this strain might have the capability to synthesize other lignan compounds using pin as the precursor. Cytochalasins are other secondary metabolites produced by Phomopsis, and are derived from phenylalanine, tyrosine, leucine, tryptophan and alanine^[60]. Schümann and Hertwerck^[61] used RNA silencing to determine that the fungal PKS-NRPS hybrid synthase CheA plays an essential role in cytochalasan biosynthesis in Penicillium expansum. The genomic sequence analysis reveals that Phomopsis sp. XP-8 also has the potential to produce cytochalasan. Overall, the sequencing results indicate a huge potential for the production of a diverse array of industrially important natural products from Phomopsis sp. XP-8.

In conclusion, we report first genome sequencing, assembly, annotation and analysis of candidate genes involved in lignans and other secondary metabolites biosynthesis in the plant endophytic fungus Phomopsis sp. XP-8. This study provides basis for the further metabolic engineering of pin and its glucoside derivatives production, further studies into mining novel bioactive secondary metabolites of plant endophyte, and plant-endophyte interactions. More important, this research will provide important informations on the genes of Phomopsis sp. XP-8 and important reference to similar study of other endophytic fungi.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 31201408). The authors declare that they have no conflict of interests. We thank Biomarker Biotechnology Corporation (Beijing, China) for technical assistance with data analysis and Zhiyuan Yin, Ph.D. for suggestions about the data analysis.