To compare more genomes of species in Aspergillus section Fumigati, A. viridinutans IFM 47045 and A. pseudoviridinutans IFM 55266 were newly sequenced using the Illumina short read system. The sequenced reads were assembled into 47 and 24 scaffolds, with total length 34.88 and 33.33 Mb, respectively (Table 1). The numbers of predicted proteins were 10,039 and 11,281, and the quality of the genome assembly was confirmed (presence of BUSCO genes: 99.7 and 99.6%), respectively. Our group has reported genome sequences of A. lentulus IFM 54703 and A. udagawae IFM 46973 (Kusuya et al., 2015, 2016). To improve the genome datasets, we here reassembled the sequence of A. udagawae IFM 46973 using ALLPATHS-LG (ver. R52488) and reannotated the genes of A. lentulus IFM 54703 and A. udagawae IFM 46973 (the updated gene IDs are designated Alt_000001-T1 and Aud_000001-T1, respectively). As a result, the sequence assembly was much improved; the genomes of both species were composed of 17 scaffolds. The four abovementioned genome datasets were deposited and updated in the NCBI (https://www.ncbi.nlm.nih.gov/).

In addition to these genomes, genome information for A. fumigatus Af293 (Nierman et al., 2005), N. fischeri NRRL 181 (Fedorova et al., 2008), A. fumigatiaffinis CNM-CM6805 (Dos Santos et al., 2020), A. novofumigatus IBT 16806 (Kjærbølling et al., 2018), and A. thermomutatus HMR AF 39 (Parent-Michaud et al., 2019) was available at the NCBI and was retrieved for this study. Consequently, a phylogenetic tree was constructed using conserved genes in these nine closely-related species in Aspergillus section Fumigati. N. fischeri was the closest relative to A. fumigatus (Figure 1A).

The numbers of predicted proteins in the fungi ranged from 9,702 to 11,534 and are summarized in Table 1. By reciprocal best hit (RBH) search with strict criteria (>80% identity and >80% of length covered), 4,493 genes were determined to be orthologs conserved in the nine Aspergillus section Fumigati species (referred to hereafter as “ASF-conserved genes”), which represented the core-genome (Figure 1B). Among the nine species, there were 948–2,618 species-unique genes (i.e., having no orthologs in other species) (Table 1; Figure 1B).

Synteny analysis of the genomes revealed that most of the A. fumigatus genome is covered by the other species, but that there are several regions unique to A. fumigatus (Figure 1C). The number of A. fumigatus genes that reside in the syntenic region to each other species was 6,874 (A. thermomutatus) to 8,266 (N. fischeri) (69.9–84.0% of A. fumigatus genes). This further suggested that N. fischeri is the closest relative to A. fumigatus among the Aspergillus section Fumigati species.

SM backbone genes encoding PKSs, NRPSs, and PKS-NRPS hybrids were identified from the genome data using antiSMASH software with manual refinement. In total, 34–84 such genes were identified in the nine species (Table 2; Supplementary Table 1). Orthologous SM-related proteins that showed identities of >80% in >80% of the protein were identified. Twenty-seven of the 34 SM backbone genes in A. fumigatus have orthologs in other species, and seven genes are unique to A. fumigatus (Figure 2A). The other species share 17–22 of the SM backbone genes with A. fumigatus. N. fischeri had the most SM genes in common with A. fumigatus (Figure 2A). Notably, 10 genes were shared among all nine studied species of Aspergillus section Fumigati (referred to hereafter as “ASF-conserved SM genes”) (Figure 2B). Meanwhile, there were 18 and 57 genes that were conserved in 5–8 species and 2–4 species, respectively. The species-unique SM backbone genes were identified in each species. A. lentulus has 4 (the fewest), and A. pseudoviridinutans has 41 species-unique SM backbone genes (the most). In total, 510 SM backbone genes were identified and grouped into 247 orthologous types (Figure 2B). A cladogram was generated based on a binary matrix (presence/absence of the SM backbone genes), which revealed that N. fischeri is the most closely related species to A. fumigatus based on the SM backbone protein distribution across species (Figure 2C).

The 10 ASF-conserved SM backbone genes include five NRPSs and five PKSs, among which six were previously characterized in A. fumigatus as being involved in the biosynthesis of fumigaclavine C (O'Hanlon et al., 2012), ferricrocin (Schrettl et al., 2007), fumarylalanine (Steinchen et al., 2013), gliotoxin (Cramer et al., 2006), 1,8-dihydroxynaphthalene (DHN)-melanin (Langfelder et al., 1998), and neosartoricin/fumicyclines (Chooi et al., 2013; König et al., 2013) (Figure 2A). The BGCs containing the ASF-conserved SM backbone genes (hereafter designated BGC1 to BGC10) were compared among the nine Aspergillus section Fumigati species (Supplementary Table 2). BGC1 and BGC3 were almost perfectly conserved across the species in terms of gene composition, gene order, and the similarity of the encoded proteins (Supplementary Figure 1). Although there were a few examples of gene loss or low protein similarity, BGC2, BGC4, BGC5, BGC8, and BGC9 were also well-conserved among the species. Meanwhile, BGC6, BGC7, and BGC10 were diverse due to extensive gene loss and low protein similarity in several components. In particular, BGC7 of A. fumigatus was unique because four of the 10 component genes have no orthologs in the other species (Supplementary Figure 1). BGC2, BGC7, BGC8, and BGC9 contain one or two genes encoding a transcription factor (TF), which are conserved in all nine species (Supplementary Table 2; Supplementary Figure 1). Interestingly, most of the TFs are located at the ends of gene clusters.

To gain insights into SM gene expression, transcriptomic analysis was conducted using five representative species of Aspergillus section Fumigati: A. fumigatus, N. fischeri, A. lentulus, A. udagawae, and A. viridinutans. The fungal strains were cultivated in four different media: potato dextrose broth (PDB), Czapek-Dox medium (CD), Sabouraud broth (SB), and potato dextrose agar (PDA). In A. fumigatus, the median numbers of Transcripts Per Kilobase Millions (TPMs) of ASF-conserved genes (n = 4,493) were 44.6, 28.3, 17.5, and 52.9 for culture on PDB, CD, SB, and PDA, respectively, much higher than the values for A. fumigatus-unique genes (3.6, 1.5, 2.2, and 5.0; n = 1,507) (Figure 3A). This was also the case for the other species (Figures 3B–E). These data indicated that a set of genes that was well-conserved among aspergilli was transcriptionally more active than species-unique genes in all the tested species.

When a gene with an expression level >¹/ 20th of the mean TPM was considered as being expressed, 12.8% of A. fumigatus genes were not expressed in any of the conditions tested (Figure 4A), and were thus considered to be silent genes. Only 3.5% of the ASF-conserved genes were silent in A. fumigatus. In contrast, 37.7% of species-unique genes were silent. This tendency was also observed in the other species (Figure 4A). With regard to SM backbone genes, 32.3% of the genes were silent in A. fumigatus, and the proportion of silent SM backbone genes was 71.7, 76.0, 83.8, and 55.1% in N. fischeri, A. lentulus, A. udagawae, and A. viridinutans, respectively (Figure 4A). These data highlight that the proportion of SM backbone genes that were expressed was quite low in species of Aspergillus section Fumigati.

For ASF-conserved SM genes, the proportion of silent genes ranged from 20 to 60% in the species of Aspergillus section Fumigati (44% in total) (Figure 4B). In contrast, 28.5–100% of the species-unique SM genes were not expressed in any of the culture conditions tested (81.1% in total) (Figure 4B). These results suggest that the ASF-conserved SM genes are transcriptionally more active than the less conserved SM genes. Lists of SM backbone genes with expression values are shown in Supplementary Tables 3–7.

The SM gene clusters that were activated in (a) specific condition(s) were further analyzed. We sought gene clusters whose component genes were coordinately regulated using the MIDDAS-M program (Umemura et al., 2013). Consequently, 14, 13, 12, 12, and 19 sets of genes were found to be expressed as clustered genes in A. fumigatus, N. fischeri, A. lentulus, A. udagawae, and A. viridinutans, respectively (Supplementary Table 8). Of these, 4, 6, 3, 5, and 9 expressed clusters contained SM backbone genes, respectively (Table 3). All these clusters belonged to the ASF-conserved or partly conserved SM genes, and no cluster with species-unique SM genes was found by the MIDDAS-M analysis. Interestingly, most of the expressed clusters have been characterized and are predicted to produce known metabolites. The gliotoxin (gliP) cluster was expressed in A. fumigatus, A. lentulus, and A. viridinutans, and the DHN-melanin (pksP) cluster in A. fumigatus and A. viridinutans. Furthermore, the fusarinine C (sidD), ferricrocin (sidC), hexadehydroastechrome (hasD), and fumitremorgin (ftmA) clusters was expressed in multiple species studied here. In addition to the metabolites whose BGCs are not present in A. fumigatus, the clusters for terrein and viriditoxin were coordinately expressed in A. lentulus and A. viridinutans, respectively.

To gain more insight into variations in the expression patterns of the clusters across the species, expression levels of the component genes in BGCs were depicted using a heat map (Figure 5). This confirmed the result from MIDDAS-M analysis that BGC4 (for pksP) was exclusively expressed in PDA in A. fumigatus and A. viridinutans, and that BGC8 (for gliP) was expressed in A. fumigatus in CD and SB, in A. lentulus in PDA, and in A. udagawae in CD. BGC1 was transcriptionally active in a cluster-wide manner in all conditions in all species.

The transcriptome data revealed coordinate expression of some BGCs, suggesting productions of some metabolites, which have not been reported in the species. First, we examined whether A. lentulus and A. udagawae are capable of producing gliotoxin in CD. No production of gliotoxin was detectable in these species, although A. fumigatus produced the metabolite (Supplementary Figure 2).

The MIDDAS-M analysis also highlighted that the cluster for terrein, containing two PKSs, Alt_007763-T1 and Alt_007763-T2, was markedly activated in A. lentulus. This cluster is very similar to the ter cluster of A. terreus in terms of gene content, order, and directions (Figure 6A; Supplementary Table 9) (Zaehle et al., 2014). The cluster in A. lentulus was coordinately expressed, with the highest expression level observed in PDB (Figure 6B). Culture extracts were examined by comparison with the standard compound terrein, revealing that A. lentulus produced terrein in PDB and on PDA (Figure 6C).

The vdt cluster for viriditoxin, which contains eight genes, has been reported in A. viridinutans (Urquhart et al., 2019) (Figure 7A). In A. viridinutans, the cluster was highly expressed in a coordinated manner in PDB, SB, and PDA (Figure 7B). On the basis of genome comparison, this cluster is also present in A. udagawae and A. pseudoviridinutans, and it was found that the genetic content of the clusters is well-conserved between these species (Figure 7A). It was of interest that no genes of the vdt cluster were expressed in A. udagawae in any of the studied conditions (Figure 7B). This suggested that the vdt cluster was active in A. viridinutans, but silent in A. udagawae.

To gain deeper insight into interspecies variations in the transcriptional regulation of the clusters for gliotoxin, terrein, and viriditoxin biosynthesis, we investigated whether there are cis-elements in the promoter regions of each component gene in each cluster. Each cluster contains a gene encoding a C6-type transcription factor (GliZ, TerR, and VdtR) that may function as a cluster-specific transcriptional regulator (Supplementary Figure 1H, Figures 6A, 7A). Because C6-type transcription factors reportedly bind to inverted CGG triplets spaced by several nucleotides, such as CGG(N_x)CCG (MacPherson et al., 2006), we sought to identify such bipartite motifs conserved in the promoter regions of component genes in the gli, ter, and vdt clusters using BioProspector (Liu et al., 2001). As binding site 5′-TCGG(N₃)CCGA-3′ has been reported for the cluster-specific transcription factor GliZ (Schoberle et al., 2014), a palindromic sequence was sought in the gli clusters of A. fumigatus, N. fischeri, A. lentulus, A. udagawae, and A. viridinutans, setting the gap as three bases (Figure 8A). The consensus sequence was found to be highly conserved among the five closely-related Aspergillus section Fumigati species in terms of position in the promoter region of each gene. No consensus sequence was detected for the gliZ gene in any of the five species; this has been reported previously for A. fumigatus (Schoberle et al., 2014).

For the ter cluster of A. terreus, the consensus sequence for binding of TerR, a cluster specific transcription factor, was proposed to be 5′-TCGGHHWYHCGG-3′ (Gressler et al., 2015). Thus, similar sequences (5′-CGGXXXXXCGG-3′ and 5′-CCGXXXXXCCG-3′) were sought in the ter cluster of A. lentulus. Four ter genes contained either of the sequences in the promoter region, and four of the five sequences corresponded to those in A. terreus (Figure 8B). No consensus sequence has been reported for the vdt cluster in any fungal species. Thus, we sought such a sequence using several combinations of parameters. When the gap was set as three bases, the well-conserved sequence 5′-TCGG(N₃)CCGA-3′ was found in all component genes of this cluster in both A. viridinutans and A. udagawae (Figure 8C). The positions of the consensus sequence were highly conserved between the species. Collectively, cis-elements for the SM gene clusters with a C6-type transcription factor were well-conserved in an interspecies manner.

The strains A. fumigatus Af293, A. lentulus IFM 54703, A. udagawae IFM 46973, A. viridinutans IFM 47045, A. pseudoviridinutans IFM 55266, and N. fischeri NRRL 181 were provided through the National Bio-Resource Project, Japan (http://www.nbrp.jp/) and are preserved at the Medical Mycology Research Center, Chiba University. The genomes of A. lentulus IFM 54703 (Kusuya et al., 2016) and A. udagawae IFM 46973 (Kusuya et al., 2015) were previously sequenced, and the data were retrieved from the NCBI database. A. viridinutans strain IFM 47045 (NRRL 4365) was isolated from rabbit dung in Australia. A. pseudoviridinutans strain IFM 55266 was isolated from a patient in Japan and identified using tubulin and calmodulin partial gene sequences (Lyskova et al., 2018).

Strains were grown in liquid PDB (BD Difco, Franklin Lakes, NJ, USA), SB (BD Difco), or CD (BD Difco) at 37°C for 5 d by inoculating each culture with three 0.5-cm² agar plugs. For asexual stage culture, mycelia cultured in PDB at 37°C for 3 d were harvested using a Miracloth, washed with distilled water, and then placed onto PDA plates (BD Difco) for another 2 d of culture at 37°C.

The genomic DNAs of A. viridinutans and A. pseudoviridinutans were extracted from a 2-d-old culture using phenol-chloroform and NucleoBond buffer set III (TaKaRa, Shiga, Japan). The DNA was fragmented in an S2 sonicator (Covaris, MA, USA), and then purified using a QIAquick gel extraction kit (Qiagen, CA, USA). A paired-end library was constructed using an NEBNext Ultra DNA Library Prep Kit (New England BioLabs, MA, USA) and NEBNext Multiplex Oligos (New England BioLabs) in accordance with the manufacturer's instructions. Mate-paired libraries with insert sizes of 3.5–4.5, 5–7, and 8–11 kb were generated using the gel selection-based protocol of the Nextera Mate Pair Kit (Illumina, San Diego, CA, USA) and a 0.6% agarose gel, in accordance with the manufacturer's instructions. The quality of the libraries was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Paired-end sequencing (100 bp) was performed by HiSeq 1500 (Illumina) using HiSeq reagent kit v1, in accordance with the manufacturer's instructions; 150-bp paired-end sequencing on a HiSeq X system (Illumina) was carried out by GENEWIZ (South Plainfield, NJ, USA).

Adapters and low-quality bases from Illumina reads were trimmed by Trim Galore (ver. 0.6.4) (Krueger, 2015) with default settings (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Mitochondrial genomes were assembled using GetOrganelle (ver. 1.6.4) (Jin et al., 2020) from the trimmed reads. To filter the mitochondrial reads, the trimmed reads were aligned against mitochondrial genomes by BWA (ver. 0.7.17-r1188) (Li and Durbin, 2009), and the mapped reads were filtered by SAMtools (ver. 1.9) (Li et al., 2009) and SeqKit (Shen et al., 2016). The assembly of nuclear genomes of A. udagawae, A. viridinutans, and A. pseudoviridinutans was carried out by ALLPATHS-LG (ver. R52488) (Gnerre et al., 2011).

The annotation of assembled genomes of A. lentulus, A. udagawae, A. viridinutans, and A. pseudoviridinutans was performed by Funannotate pipeline (ver. 1.7.4) (https://funannotate.readthedocs.io/en/latest/). Following identification of repeat sequences by RepeatModeler (ver. 1.0.11) (http://www.repeatmasker.org/RepeatModeler.html) and RepeatMasker (ver. 4.0.7) (https://www.repeatmasker.org), Funannotate ab initio prediction was performed with the option “–busco_seed_species=aspergillus_fumigatus” by Augustus (ver. 3.3.3) (Stanke et al., 2006), GeneMark-ES (ver. 4.38) (Ter-Hovhannisyan et al., 2008), GlimmerHMM (ver. 3.0.4) (Majoros et al., 2004), and SNAP (ver. 2006-07-28) (Ian, 2004) using exon hints from the proteins of A. fumigatus Af293 and N. fischeri NRRL 181 downloaded from the Aspergillus Genome Database (Cerqueira et al., 2014). Functional annotation of predicted genes was performed by using the SwissPROT (Bairoch and Apweiler, 2000), InterPro (ver. 5.42-78.0) (Jones et al., 2014; Mitchell et al., 2019), eggNOG (ver. 4.5.1) (Buchfink et al., 2015; Huerta-Cepas et al., 2016), MEROPS (ver. 12.0) (Rawlings et al., 2018), and dbCAN (ver. 8.0) (Yin et al., 2012) databases. SM gene clusters were predicted by antiSMASH (ver. 4.2.0) (Blin et al., 2019). The completeness of draft genomes and predicted proteins was evaluated by BUSCO (ver. 4.0.6) (Seppey et al., 2019) with the database eurotiales_odb10. Genome synteny analysis between A. fumigatus and other species was conducted per Kjærbølling et al. (2020). Most tools were obtained through Bioconda (Grüning et al., 2018).

A phylogenetic tree of A. fumigatus and related species was constructed using orthologs common in 11 species including A. clavatus NRRL 1 (Fedorova et al., 2008) and A. nidulans FGSC A4 (Galagan et al., 2005) as the outgroup. Orthologous genes among the 11 species were identified by OrthoFinder (ver. 2.3.12) (Emms and Kelly, 2019). Then, 4,465 single copy ortholog genes were aligned by using MAFFT (ver. 7.471) (Katoh and Standley, 2013), and well-aligned regions were extracted using Gblocks (ver. 0.91b) (Talavera and Castresana, 2007) with the options “-t=p, -b4=5, -b5=h.” A total of 4,465 aligned sequences of each strain were concentrated into one long protein sequence. A phylogenetic tree was constructed using multithreaded RAxML (ver. 8.2.12) (Stamatakis, 2014), the PROTGAMMAWAG model, and 100 bootstrap replicates, and visualized by iTOL (Letunic and Bork, 2019).

Identification of SM genes was performed with antiSMASH (Blin et al., 2019). NRPS, PKS, and PKS-NRPS hybrid genes were extracted from the results. The orthologous relationships were determined by RBH with the criteria of BLASTp (ver. 2.7.1+) coverage >80% and identity >80% (Camacho et al., 2009). A cladogram was constructed based on a binary matrix (presence/absence of the PKSs and NRPSs) using Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm#ctv). A tree and heat map were constructed using Tree View (Saldanha, 2004).

Each strain was cultured, harvested, and cells were ground into a fine powder using a mortar. Total RNA was isolated using Sepazol-RNA Super G (Nacalai, Kyoto, Japan) in accordance with the manufacturer's instructions. The RNA isolation was carried out with two biological replicates, and equivalent quantities of the total RNA were pooled for preparation of the RNA-Seq libraries. The RNA-Seq libraries were constructed using a KAPA mRNA Hyper Prep Kit (Nippon Genetics, Tokyo, Japan), in which mRNA was purified by poly-A selection, second-strand cDNA was synthesized from the mRNA, the cDNA ends were blunted and polyAs added at the 3′- ends, and appropriate indexes were ligated to the ends. The libraries were PCR amplified, and the quantity and quality were assessed using a Bioanalyzer (Agilent Technologies). Each pooled library was sequenced using Illumina HiSeq 1500 apparatus. The gene expression levels were estimated using a previously described method (Takahashi et al., 2017). Briefly, the sequencing reads were mapped to the reference genomes using STAR (ver. 2.4.2a) (Dobin et al., 2013). A raw read count was conducted using HTSeq (ver. 0.5.3p3) (Anders et al., 2015), and transcript abundances were estimated as TPMs (Conesa et al., 2016). The expressed PKS and NRPS genes were identified using the 5% mean TPM criterion.

Binary logarithms of TPM values generated in the four culture conditions were analyzed using the MIDDAS-M algorithm to detect gene clusters in which certain condition(s) coordinately expressed or repressed genes (Umemura et al., 2013). Briefly, the induction ratio of each gene was evaluated for every pairwise combination of the four culture conditions. After Z-score normalization, the gene cluster expression scores were calculated for each gene using the algorithm. The maximum cluster size was set as 30. The threshold to detect clusters was set to the value corresponding to a false positive rate of 0, which was evaluated from data in which the original gene order was randomly shuffled. The threshold values were 2.1E5, 2.3E5, 2.4E5, 1.3E5, and 3.5E5 in A. fumigatus, A. lentulus, A. udagawae, A. viridinutans, and N. fischeri, respectively.

For gliotoxin detection, 5 mL of culture supernatant was extracted with an equivalent volume of ethyl acetate. The organic layer was concentrated in vacuo, and the precipitate was dissolved in 2 mL of dimethylsulfoxide (DMSO). The DMSO solution (5 μL) was subjected to HPLC analysis, which was performed using an Infinity1260 modular system (Agilent Technologies, Santa Clara, CA) consisting of an autosampler, high-pressure pumps, a column oven, and a photodiode-array detector with InfinityLab Poroshell 120 EC-C18 column (particle size: 2.7 μM; length: 100 mm; internal diameter: 3.0 mm; Agilent Technologies). Running conditions were gradient elution, 5–40% acetonitrile in water over 18 min; flow rate, 0.8 mL min⁻¹; and detection at 254 nm. Gliotoxin production was identified by comparing the retention times and the UV spectra of natural products with those of a gliotoxin standard that was purchased from Sigma-Aldrich (St. Louis, MO).

For terrein detection, 5 μL of culture supernatant was subjected to HPLC analysis, which was performed using the same HPLC system as described above. Running conditions were gradient elution, 5–100% acetonitrile in water over 30 min; flow rate, 0.8 mL min⁻¹; and detection at 254 nm. Terrein production was identified in the same manner as gliotoxin production. The authentic standard was purchased from Cayman Chemical Company (Ann Arbor, MI).

Sequence 500 bp upstream of each gene encoded in gli and vdt clusters was retrieved as the promoter sequence. Two-block motifs were sought by BioProspector (Release 2) (Liu et al., 2001) with the options “-d 1 -W 4 -w 4 -g 3 -G 3.” The sequence logos were depicted by WebLogo (Crooks et al., 2004).