The clinical strain of P. capsulatum (LiaoWQ-2011 = CBS 134186) was isolated from a pulmonary fungus ball in a patient with type 2 diabetes in 2013, and the environmental strain of P. capsulatum (ATCC 48735) was isolated from exposed canvas on Gilbert Island in 1945. The representative isolates were included for whole genome sequence analyses in this study, as described in Supplementary Table S1.

The representative strains for P. capsulatum (i.e., the clinical strain CBS 134186, and the environmental strain ATCC 48735) were grown on potato dextrose agar (PDA) at 25°C for 14 days. Genomic DNA for whole-genome sequencing was extracted from the pure cultures using the DNeasy Plant Mini kit (QIAGEN Co., Ltd, Hamburg, Germany) according to the instruction manual. RNA for transcriptome sequencing was extracted using the RNeasy Plant Mini kit (QIAGEN Co., Ltd, Hamburg, Germany) following the protocol for the isolation of the total RNA.

Our strategy for whole-genome sequencing involved a combination of the Illumina MiSeq (Illumina Inc., San Diego, CA, USA) and the Pacific Biosciences RS II (Pacific Biosciences, Menlo Park, CA, USA) sequencing platforms, and the sequence data from the Illumina platform were used to proofread the PacBio assembly sequence. Illumina libraries were prepared using TruSeq DNA sample prep kits (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s instructions. Two paired-end sequencing libraries were constructed. One was a 250 Pair End (PE) sequencing library with an insert sizes of approximately 400 bp, and the other was a 150PE library with insert sizes ranging from 200 bp to 1 kb. Pacific Biosciences RS II sequencing technology was used as the sequencing platform. A 10-kb Single Molecule Real Time (SMRT) bell library was prepared from sheared genomic DNA via a 10-kb template library preparation workflow. SMRT sequencing was conducted on the PacBio RS II sequencing platform using C3 sequencing chemistry and P5 polymerase with 16 SMRT cells.

The Hierarchical Genome Assembly Process (HGAP) was used to assemble these two sequenced genomes. De novo assembly of the PacBio read sequences was carried out using continuous long reads (CLR), followed by the HGAP workflow (PacBioDevNet; Pacific Biosciences) as available in SMRT Analysis v2.1. HGAP consists of preassembly, de novo assembly with Celera Assembler (CA), and assembly polishing with Quiver. CA software version 7.0 was utilized in the pre-assembly step, and the PacBioRs_PreAssembler with one module and a default minimum subread length of 500 bp, a minimum read quality of 0.80, and a minimum subread length of 7500 bp was used to perform error correction for the raw data generated by the PacBio RS II platform. To polish the assembled sequence from HGAP, the MiSeq read sequences were mapped using BWA v0.7.5a (Li and Durbin, 2009, 2010), and the SNPs and INDELs were called and corrected by SAMtools v0.1.18 (Li et al., 2009) and an in-house script. The genomes of P. capsulatum strains CBS 134186 and ATCC 48735 were deposited in Genbank under the accession numbers JPLR00000000 and JPLQ00000000.

InterProScan was used to obtain functional annotations for all predicted genes and to determine motifs and domains. To obtain high confidence gene models for P. capsulatum, coding genes were predicted with AUGUSTUS v3.1 (Stanke et al., 2008). All of the predicted gene models were functionally annotated based on their sequence similarity to genes and proteins in the NCBI nucleotide (Nt), non-redundant and UniProt/Swiss-Prot protein databases. The gene models were also annotated based on their protein domains using InterProScan (Jones et al., 2014). All genes were classified according to Gene Ontology (GO), eukaryotic orthologous groups (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways. Repeat sequences were masked throughout the genome using Repeat Masker v3.2.9 and the RepBase library v16.08 (Jurka et al., 2005). tRNAs were identified using tRNA scan-SE (Lowe and Eddy, 1997; Schattner et al., 2005).

For analysis of orthologous and unique genes, we used the Markov clustering program OrthoMCL v2.09 (Li et al., 2003) to define a gene family as a group of genes that descended from a single gene present in the last common ancestor of the considered species. The peptide sequences were subjected to an all-versus-all BLASTp search with a threshold value of E ≤ 1e–5 and then clustered by MCL with an inflation value of 1.5. A pipeline was used to cluster individual genes into gene families and to analyze orthologous and unique genes. The pipeline included four main steps: (1) data preparation, (2) BlastP searches for all of the protein sequences against a database containing a protein dataset from all species with E-values under 1E-5 and a 90% Match Cutoff, (3) clustering by MCL, and (4) extracting gene families.

For the phylogenetic analysis based on all 40 strains, alignments of single-copy orthologs were produced using MAFFT v7.221 (Katoh and Standley, 2013). ProtTestv3.4 (Darriba et al., 2011) was utilized to determine the best-fitting model for the data. A maximum likelihood phylogeny was reconstructed using RAxML v8.1.23 (Stamatakis, 2014). For the phylogenetic analysis based on virulence-related genes, 43 virulence-related genes shared among 22 strains were selected for a similar phylogenetic analysis.

The gene clusters related to secondary metabolite biosynthesis were predicted using antiSMASH. All protein sequences were also searched against a local polyketide synthase (PKS)/non-ribosomal peptide synthetase (NRPS) database (a subset of the SwissProt database) using BlastP with a threshold E-value ≤ 1E–5 to identify PKS and NRPS. The putative PKS/NRPS protein sequences were further searched against the NCBI Conserved Domain Database (V3.09) to confirm that the three typical domains (PKS, MT, and ACP) were present. The consensus sequences of PKS/NRPS predicted by antiSMASH (Blin et al., 2013) and manually identified proteins were of high confidence.

Whole genome blast searches were conducted against protein sequences in the pathogen–host interaction database¹ (E-value ≤ 1E-10) (Urban et al., 2015). In the cytochrome P450s analysis, the reference CYP sequences were downloaded from http://drnelson.uthsc.edu/P450seqs.dbs.html. All predicted proteins were then used to search the reference CYP data set using the BLASTP program with a cutoff E-value ≤ 1E–10. SignalP (Petersen et al., 2011) was used to predict signal peptide cleavage sites in amino acid sequences, and putative protein sequences were also examined using the TMHMM method (Sonnhammer et al., 1998) for transmembrane domains. The predicted signal peptide sequences without transmembrane domains were identified as exocrine proteins.

To identify significantly different genes between the genomes of strains CBS 134186 and ATCC 48735, all of the annotated genes were subjected to all-versus-all BLASTP searches with a threshold E-value ≤ 1E–10. Any gene with a coverage rate below 90% was identified as a specific candidate gene. Specific genes were obtained by filtering out genes without a specific annotated function or with any alignments near the gene boundary. Furthermore, comparisons of nucleotide sequences between strains ATCC 48735 and CBS 134186 were carried out for SNP and INDEL detection using the Mauve program with the default values (minimum match length, 20 bp). An in-house script was used to call SNPs and INDELs between strains ATCC 48735 and CBS 134186.

Synteny analysis between ATCC 48735 and CBS 134186 was performed using MUMmer v3.0 (Kurtz et al., 2004) and MAUVE v20150226 (Darling et al., 2010). The two strains of P. capsulatum were visualized using the Circos program v0.68 (Krzywinski et al., 2009).

Poly-A mRNA was isolated using oligo-dT-coupled beads from 40 μg total RNA of each sample, followed by shearing. The isolated RNA samples were used for first strand cDNA synthesis with random hexamers and Superscript II reverse transcriptase. After end repair and the addition of a 3’dA overhang, the cDNA was ligated to the Illumina paired-end adapter oligomix, and approximately 200 bp fragments were size selected by gel purification. After 16 cycles of PCR, the libraries were sequenced using Illumina HiSeq2500 and the paired-end sequencing module. The RNA expression analysis was based on the predicted genes of CBS 134186 and ATCC 48735. First, Tophatv2.0.10 (Kim et al., 2013) was used to map the mRNA reads to the genome, and Cufflinksv2.1.1 (Trapnell et al., 2010) was then used to calculate the expected fragments per kilobase of transcript per million mapped reads (FPKM) as the expression values for each transcript. RNA-Seq data of CBS 134186 and ATCC 48735 are available at National Center for Biotechnology Information (NCBI) under SRR accession numbers of SRR 4031065 and SRR 4051963.

The genomes of both P. capsulatum strains were successively sequenced by a combination of the Pacific Biosciences RS II and Illumina MiSeq platforms. A total of 4.9 and 4.4 Gb of raw data were generated for strains CBS 134186 and ATCC 48735, respectively, using the PacBio RSIISMRT platform and a 10-kb library preparation. A total of 3.2 and 3.0 Gb of raw data from strains CBS 134186 and ATCC 48735, respectively, were generated using the Illumina platform. The genome size of P. capsulatum was estimated to be approximately 35 Mb. Hence, we sequenced the genomes of the two P. capsulatum strains with an average coverage of ∼231 × for strain CBS 134186 and ∼211× for strain ATCC 48735. For strain CBS 134186, the genome size was approximately 34.34 Mb, with a contigN50 size of 3.30 Mb; for strain ATCC 48735, the genome size was approximately 34.37 Mb, with a contigN50 size of 3.19 Mb. The read data from the Illumina platform were used to correct the assembled sequence based on the PacBio platform because the latter normally results in relatively more mistakes (Ferrarini et al., 2013). Fifty SNPs and 23 INDELs were corrected in the genome of strain CBS 134186, whereas 48 SNPs and 29 INDELs were corrected in the genome of strain ATCC 48735. For the corrected SNPs and INDELs, polyN-type errors were the main error type, which accounted for approximately 75 and 73% of the total corrected sides for CBS 134186 and ATCC 48735, respectively. Two different gene footprint coverage methods were applied to validate the coding region coverage of the genome assemblies. The results obtained using the Core Eukaryotic Genes Mapping Approach (CEGMA) showed that 95.56% (237/248) of the core eukaryotic genes that mapped to the assembly genome of CBS 134186 were identified and that 97.18% (241/248) of the eukaryotic genes that mapped to the genome of ATCC 48735 were identified (Table 1).

For the genome of strain CBS 134186, a total of 11,080 gene models were identified based on Gene Ontology (GO) annotations, with an average coding sequence length of 1807.8 bp. For the ATCC 48735 genome, 11,076 gene models were obtained, with a total length of 34.37 Mb and an average coding sequence length of 1808.7 bp (Table 1; Supplementary Table S2). The average GC content in both strains of P. capsulatum was 49.06%. Approximately 1.41 Mb of repeated regions were found in the genome of the clinical strain CBS 134186, which accounted for 4.19% of its genome size, whereas 1.42 Mb of repeated regions were found in the genome of the ATCC 48735 strain. For both strains, interspersed repeats were the predominant type of repeat region, accounting for 84.9% and 85.1% of the repeat regions in the CBS 134186 and ATCC 48735 genomes, respectively. In addition, tRNAScan-SE predicted 185 and 174 tRNAs in the genomes of CBS 134186 and ATCC 48735, respectively.

A phylogenetic analysis based on the whole genome sequence was conducted to determine the relationship of P. capsulatum with other important species belonging to Eurotiales. The included species were pre-divided into four categories: Eurotiales species known to be medically important fungal pathogens; opportunistic fungal pathogens; species newly reported to cause invasive infection; and species unreported as human pathogens. Except for the sequencing data of P. capsulatum, other publicly available and complete genomic data of Eurotiales were downloaded from NCBI, including those data for thirteen Penicillium strains, fourteen Aspergillus strains, five Talaromyces strains, two Neosartorya strains, one Rasamsonia strain and one Byssochlamys strain. Two zygomycete pathogens, Rhizopus oryzae and Mucor circinelloides, were selected as taxonomic out-groups (Supplementary Table S1).

The genomes ranged in size from 26.0 to 35.9 Mb within Penicillium, from 24.2 to 39.5 Mb within Aspergillus, from 28.6 to 37.6 Mb within Talaromyces and from 31.7 to 32.2 Mb within Neosartorya. The genome sizes of Rasamsonia emersonii and Byssochlamys spectabilis, which belong phylogenetically to Eurotiales, were 28.2 and 29.8 Mb, respectively. Within Eurotiales. A. rambellii, and T. stipitatus carried significantly higher percentages of repetitive elements (13.85 and 12.43%, respectively). The number of repetitive sequences of other strains in Eurotiales was relatively low; only 1.01-8.18% of each genome was classified as repetitive. In contrast to the low GC content in the out-group, the genomic GC content of Eurotiales ranged from 46 to 53% with a balanced ratio (Figure 1).

Subsequently, all of the single-copy orthologs were selected for the phylogenetic analysis. Eurotiales were divided roughly into three groups, in which Penicillium had a closer phylogenetic relationship to Aspergillu than to Talaromyces. Fifteen Penicillium strains were grouped together in one clade, suggesting their origin from a single progenitor species (Figure 1). The P. capsulatum strains were confidently placed in a basal position to P. chrysogenum and other species belonging to the subgenus Penicillium. The clinical strain of P. capsulatum showed a very close phylogenetic relationship with the environmental strain, forming a subclade within the Penicillum clade. Talaromyces marneffei, which was previously known as P. marneffei, was grouped in the genus Talaromyces based on its genomic data in this study. In addition, Neosartorya was classified with Aspergillu into one clade according to their genomic characteristics. These results offered us a new perspective based on whole genome data for future phylogenetic classifications.

To further explore the pathogenic factors of P. capsulatum and its closely related strains, we performed a phylogenetic analysis based on 106 virulence-related genes identified from medically important fungal pathogens, such as A. fumigatus. T. marneffei et al. as references for studying the novel fungal pathogens (Supplementary Table S3). Except for two P. capsulatum strains, 18 species were selected according to their medical or economic importance and their phylogenetic position, with the intention to obtain a focused representation of Eurotiales. The following Eurotiales species were included: species known to frequently cause invasive infections as medically important fungal pathogens, including A. fumigates and T. marneffei; species reported as opportunistic fungal pathogens, including A. niger. A. nidulans. A. terreus. A. flavus. A. oryzae. P. chrysogenum. P. oxalicum. P. digitatum, and N. fischeri; and species unreported as human pathogens, including T. stipitatus. A. clavatus. P. camemberti. P. paxilli. P. aurantiogriseum. P. expansum, and P. rubens Wisconsin. Rhizopus oryzae and Mucor circinelloides were included as taxonomic out-groups.

The presence or absence of virulence-related genes was searched in the 22 studied genomes. We found that a total of 43 virulence-related genes were shared the studied genomes. These genes could be divided into the following six categories: toxins, thermo-tolerance, signaling and regulation, resistance to immune response, cell wall, allergens and nutrient uptake (Supplementary Table S4). Generally, the medically important fungal pathogens have more virulence-related genes than non-human pathogens. In this study, the clinical and environmental isolates of P. capsulatum were found to be very similar in terms of their virulence-related genes. Surprisingly, we found in this phylogenic analysis that P. paxilli, which is used as a model to study the biochemistry of indol-diterepene biosysnthesis and is seldom reported to be infectious to humans, had a high proportion of virulence-related genes and very close relationship to P. capsulatum (Figure 2). Thus, we selected P. paxilli and two well-known human pathogens, A. fumigates and T. marneffei, for a more detailed comparison with P. capsulatum.

Using a Venn diagram of four strains, we compared the P. capsulatum genome with a non-human pathogen, P. paxilli, which had a close phylogenetic relationship to P. capsulatum, as well as A. fumigates and T. marneffei, two major human pathogens. By analyzing the homologous genes among the four genomes and enriching the virulence genes among these groups, as expected, the majority of homologous gene groups were in the groups of C1&C2&C3&C4. In this shared group among all four strains, most genes were identified as belonging to the cytochrome P450 (CYP) family and the major facilitator superfamily (MFS), which generally play important roles in the biosynthesis and transportation of metabolites. The second largest group of genes shared by these strains was group C2&C3&C4, indicating that the relationships among P. capsulatum. P. paxilli, and A. fumigate were closer than T. marneffei. Compared with the other three strains, the unique paralog gene groups of A. fumigate in group C3 consisted of typical virulence-related genes, such as genes involved in gliotoxin synthesis (GliZ. Glip) and fumitremorgin synthesis (ftmE-H). Like group C3, group C2 comprised unique paralog gene groups found in P. capsulatum CBS 134186. This group had the most unique paralog gene groups, containing 1865 gene groups (1975 predicted genes in total). Some of these genes may increase the pathogenic capacity of these four strains (Figure 3, Supplementary Table S5).

By further investigating the virulence genes, we found that only 7 out of 15 groups contained virulence-related genes (Supplementary Table S6). The majority of the virulence genes (n = 75) were in the group C1&C2&C3&C4, which had virulence genes among all the four strains. There were 13 virulence factors in group C2&C3&C4, which contained the second greatest number of virulence genes. Interestingly, group C3&C4 contained four virulence genes, which was greater than the number of virulence genes in group C1&C3 and group C2&C3. The enrichment of virulence-related genes in P. paxilli indicated the need for paying increased attention to its potential pathogenicity. However, it is possible that having these virulence-related genes is not sufficient to cause disease. Thus, the underlying difference between fungal pathogens and non-pathogens needs further study.

The genome of strain CBS 134186 was 0.033 Mb smaller than the genome of ATCC 48735 (Table 1). We speculated that the gaps between the large scaffolds in the CBS 134186 genome (n = 104) and the ATCC 48735 genome (n = 101) may be responsible for the different genome sizes observed between these two strains. Pair-wise genome alignments showed that the two P. capsulatum strains were highly syntenic and shared, on average, 99.98% sequence similarity (Figures 4A,C). These results suggested that the pathogenic potential of the clinical strain CBS 134186 and the environmental strain ATCC 48735 may be extremely similar or that the specific genes, SNPs or INDELs present in the CBS 134186 genome may play a key role in causing invasive infection. Among these probable mutations, a total of seven unique genes, SNPs or INDELs in the CBS134186 genome may be involved in the cause of make itself to be a novel fungal pathogen (Table 2, Supplementary Tables S7-S9; Figure 4B). Cytochrome P450, secondary metabolism, and exocrine proteins analyses were also performed on these seven genes. One of the unique genes (g512) was predicted as an exocrine protein gene; there were no positive results in the secondary metabolite and cytochrome P450 analyses. The unique genes coding for vegetative incompatibility protein and glycosyl hydrolases (g512 in contig27), the SNP in the gene coding for circumsporozoite protein (g399 in contig17), and the INDEL in the gene containing the F-actin capping protein beta subunit (g8270 in contig44) in the CBS 134186 genome may be related to its adaptation to stress conditions; a similar phenomenon has been observed in the interaction between A. fumigates and humans (Abad et al., 2010).

To more accurately annotate the P. capsulatum genome, we performed deep transcriptome sequencing of both P. capsulatum strains (CBS 134186 and ATCC 48735), which generated 4.0 and 2.9 Gb of RNA sequencing data, respectively. Compared with ATCC 48735, CBS 134186 contained 481 differentially expressed genes and two up-regulated virulence-related genes, which were observed under in vitro culture conditions by three replicates to give a general estimate. These virulence-related genes were beta-(1,3) glucanosyltransferase (GEL-3) and mannose-6 phosphate isomerase (MPI), which were up-regulated fivefold and threefold, respectively. Moreover, these two genes are involved in cell wall composition and maintenance, and the high expression of these genes may increase the pathogenicity of the clinical strain CBS 134186.

We comparatively analyzed the morphological characteristics between strain CBS 134186 and strain ATCC 48735. Both P. capsulatum strains had smooth walled conidia that were ellipsoid to slightly cylindrical in shape (Figure 5). There were no significantly different morphological characteristics between the two strains, which was in accordance with their highly similar genomes.