Genome sequencing of Aspergillus glaucus ‘CCHA’ provides insights into salt-stress adaptation

Microorganisms, culture medium and inoculum development

The sampling site for the isolation of salt-tolerant A. glaucus ‘CCHA’ was the air-dried surface of wild vegetation from a solar salt field in Northeast China. For fungal isolation, potato extract plus 20 g · L⁻¹ glucose, supplemented with 250 g · L⁻¹ NaCl was used. The production of an intense greyish-green colour under and around the fungal colony in the salt-supplemented agar was considered a halophilic positive reaction. The strain was identified based on morphological characterisations and subsequently confirmed by 5.8S rDNA sequencing of fungal genomes from NCBI. The strain was deposited under the accession name: CCHA. A defined medium, composed of 20 g · L⁻¹ glucose and 200 g · L⁻¹ potato extract prepared in 50 g · L⁻¹ NaCl, was used for growth and maintenance. Fully sporulated slants were stored at 4 °C and subcultured once every 2 weeks. After being dislodged using a sterile inoculation loop under precise sterile conditions, the spores in the suspension were characterised by microscopy. The effects of salinity on the growth of the fungal isolate were determined by preparing the medium supplemented with NaCl to a final saline concentration of 1–8% (in 1% increases, e.g. 1 for 1%, 2 for 2% and 8 for 8%) and 9–21% (in 2% increases, e.g. 9 for 9%, 10 for 11% and 15 for 21%) for 6 days under 30 °C.

Genome sequencing, assembly and validation

High-quality genomic DNA was extracted from the fruiting bodies using a modified Benzy method (Min et al., 1995). RNase A and proteinase K were separately used to remove RNA and protein contamination, respectively. The quantity and quality of the isolated DNA were separately checked by electrophoresis on a 1.0% agarose gel and on a Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Three Solexa genomic sequencing libraries (500 bp, 2 KB and 15 KB) were separately constructed using Cre-Lox recombination following the Illumina instructions (Van Nieuwerburgh et al., 2012) and subsequently sequenced using an Illumina HiSeq 2000 platform. Paired-end and mate-paired reads from Solexa sequencing were assembled by SOAPdenovo to construct contigs (Li et al., 2010). The contigs were assembled using the GS DE novo Assembler (Roche, 2011). A total of 2,595.2 M Illumina reads were selected to perform the genome size estimation. The distribution of 17-kmer had a major peak at 82x. Based on the total number kmers and the corresponding kmer depth of 84, the A. glaucus ‘CCHA’ genome size was estimated using the formula: Genome size = kmer_number/Peak_Depth.

Genome annotation

Putative protein-coding genes of A. glaucus were predicted by combining several different ab initio gene predictors and sequence evidence, including protein sequences from closely related species and ESTs assembled in this study. The quality validation of gene models was evaluated by aligning the transcriptome, Ascomycota BUSCOs and homologous peptides to our gene predictions. A popular software tool RepeatMasker was used to identify and classify the repetitive elements (Tarailo-Graovac & Chen, 2009).

Comparative genomic analysis

A comparative genomic analysis was performed using the MAUVE programme (Version 2.3.1). Default scoring and parameters were used to generate the alignment. Unique and shared gene families among A. glaucus ‘CCHA’ and other Aspergillus strains (AN513, AG516 and AA106) were clustered using the OrthoMCL method.

Transcriptome sequencing and analysis

The strains of every sample were frozen in liquid nitrogen, and the total RNA was extracted using TRIzol and then treated with Dnase I according to the manufacturers’ protocols. mRNA purified by Sera-Mag magnetic oligo (dT) beads (Illumina, San Diego, CA, USA) from 20 µg total RNA per sample was reverse transcribed into double-stranded cDNA to generate the DNA libraries. The cDNA ends were repaired, and the cDNA was amplified, denatured and then sequenced on an Illumina Genome Analyzer IIx using special reagents.

The RNA-seq reads were aligned to the genome of A. glaucus ‘CCHA’ using the Hisat2 programmes (Trapnell, Pachter & Salzberg, 2009). Gene expression levels were measured in terms of FPKM values, for which the expression level of each gene is the sum of the values of its isoforms (Toung et al., 2011). The criteria FDR ≤ 0.05 and |log2FC| >= 1 were used to verify the significance of gene expression differences using edgeR.

Isolation of salt-stress-related genes

Construction of a size-fractionated cDNA library was performed using the SMART cDNA Library Construction Kit (Clontech, Mountain View, CA, USA) according to the manufacturer’s instructions, with minor modifications (Fang et al., 2014). E. coli BL21 was transformed separately using cDNA obtained from plasmid pools and selected for salt stress on a 0.75-M NaCl-containing medium that led to the full repression of the empty vector. The resulting transformants were replicated onto media containing different salt concentrations and screened for strains that grew a rate similar to the empty vector-containing transformants on normal media but grew more readily on the salt-stress media.

Phenotypic analysis of transgenic Arabidopsis thaliana

The recombinant pBI121 vector, into which CCHA-2114 cDNA was inserted, was transformed into A. thaliana (Columbia ecotype) plants by Agrobacterium tumefaciens ‘EHA105’. Wild type (WT) and T3-homozygous CCHA-2114-overexpressing A. thaliana plants were sown in soil under long-day conditions (16 h light, 8 h dark). After 30 days, transgenic plants and control plants were soaked in Hoagland’s medium for 1, 3 and 5 days with or without 200 mm NaCl. Total SOD activity and MDA levels, a biomarker for lipid peroxidation, were assayed using previously described methods (Giannopolitis & Ries, 1977; Li et al., 2017). The degree of ion leakage was calculated as follows: conductivity (before boiling)/conductivity (after boiling) × 100%. The total chlorophyll, chlorophyll a and chlorophyll b contents were measured after extraction with 80% acetone (Fang et al., 2014). All the experiments were conducted at least three times independently in four lines (WT and three transgenic plant lines).

Strain isolation and the salt-tolerance assessment of A. glaucus ‘CCHA’

The halophilic fungal strain was isolated from the surface of wild vegetation growing around a saltern, and it was identified as A. glaucus strain CCHA based on morphological properties, an ITS sequence comparison (Liu et al., 2011), and its salt-tolerance level (Figs. 1A–1C). The strain was extremely tolerant to salt levels over 25%, which is much greater than the levels tolerated by halotolerant A. niger and A. flavus (Figs. 1B and 2A–2E). The phylogenetic tree of the CCHA 5.8S rDNA sequence and the corresponding sequences of seven other reported Aspergillus strains provided insights into the evolution of the A. glaucus genome (Fig. 1D).

Genome sequencing, assembly and annotation

Because microbial genomes adapt to environmental conditions (Chan et al., 2014; Gibbons & Rokas, 2013), to better use the RNA-seq data to confirm the candidate gene of saltern-isolated A. glaucus CCHA, we denovo sequenced the A. glaucus ‘CCHA’ genome using Illumina Solexa sequencing technology (Table 1). We generated approximately 17.5 GB of raw sequencing read data, resulting in a 358× genome coverage. The obtained reads were assembled using the Berry Genomics Inc. platform, which yielded a high-quality genome assembly. The N50 length of the contigs was 209 KB and the longest contig was 474 KB. The complete genome of A. glaucus ‘CCHA’ contained 76 chromosomal scaffolds of approximately 26 MB (82% coverage of the 31.6 MB genome estimated by the k-mer method in Fig. 3), with an overall G + C content of 47.9%. The genome contains 10,066 predicted open reading frames. The general characteristics of vthe complete genome are summarised in Table 2. The whole-genome sequence data reported in this paper have been deposited in the Genome Warehouse of the BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under accession number GWHAAFV00000000, which is publicly accessible at http://bigd.big.ac.cn/gwh (BIG Data Center Members, 2018). Genomic comparisons between A. glaucus ‘CCHA’ and other fungal strains revealed that strain CCHA has a genome size comparable to those of Magnaporthe grisea (28.6 MB) (Dean et al., 2005) and Aspergillus fumigatus (29.4 MB) (Nierman et al., 2005).

The clusters of orthologous groups (KOG/COGs) programme is extremely useful for compiling the annotated gene data of a complete genome, for readily describing the data, and for systematically deducing the functions of protein families (Tatusov, Koonin & Lipman, 1997). The predicted 10,066 genes of A. glaucus ‘CCHA’ were sorted according to COG classifications (Fig. 4A) and 1,647 of the genes were associated with metabolism, 10.02% with information storage and processing, and 13.82% with cellular processes and signalling. Some of the genes (2.71%) could not be categorised into COG classes because the functions and features of these genes were poorly characterised (Fig. 4B).

Comparative genome analysis

The Aspergillus strains shared most of the genes assigned to general cellular functions. As shown in the Venn diagram constructed for four representative Aspergillus genomes (Fig. 5), all the strains shared 6,309 coding sequences (CDSs). In addition, strain AG516 shared the most (1,079) additional CDSs from the core collection of genes with ‘CCHA’ (Fig. 5A). The fewest unique CDSs (21) were detected in ‘CCHA’, while ‘AG516’ had the greatest number of unique genes.

An ortholog analysis revealed that the four Aspergillus species had 11,923 clusters, 11,522 ortholog clusters that included at least two species, and 5,986 single-copy gene clusters (Fig. 5B). To better understand the relationships among taxa based on the genomes, a MAUVE progressive alignment was performed for ‘CCHA’ with other genomes selected based on the phylogenetic tree of concatenated genes. Strain CBS_516 had a similar gene cluster distribution to that of ‘CCHA’, supporting the assignment of strain CCHA to A. glaucus (Fig. 6). To compare the genomes of A. glaucus ‘CCHA’ with those of widely recognised A. glaucus ‘CBS_516’, a dot-blot analysis was conducted using the RAST programme. The dot-blot analysis showed that the genome of ‘CCHA’ was very similar to that of ‘CBS_516’, as indicated by the diagonal line (Fig. 7A). When the whole-genome sequences were compared using MAUVE software, the locations and sizes of the genes in ‘CCHA’ were similar to those in ‘CBS_516’ (Fig. 7B). Among the genes annotated by a BLAST algorithm-based analysis, conserved genes between ‘CCHA’ and ‘CBS_516’ had high-query coverages.

Differentially expressed genes

To investigate the gene expression profiles and identify the critical genes involved in A. glaucus responses to salt stress, we compared the DEGs between different treatment groups (Figs. 8A and 8B; Table 3). There were 889 DEGs in ‘36H’ exposed to a 3-M salt stress compared with ‘36H’ exposed to a 1-M salt stress (Table S1). Here, ion transport-related genes, including sodium P-type ATPase (g4579), voltage-gated K⁺ channel beta subunit (g5829) and MFS transporters were up-regulated in response to salinity. In addition, redox-related oxidoreductases and superoxide dismutases were also enriched. There were 582 DEGs in ‘84H’ exposed to a 0-M salt stress for compared with in ‘84H’ exposed to a 1-M salt stress. Some genes involved in basic metabolism, such as sugar transporters and electron carriers were down-regulated, indicating that the normal development of A. glaucus ‘CCHA’ requires a certain amount of salt. Moreover, there were 574 DEGs in ‘84H’ exposed to a 3-M salt stress and 1,731 DEGs in ‘84H’ exposed to a 4-M salt stress compared with those in ‘84H’ exposed to a 1-M salt stress. Long-term stress mainly influenced ion transmembrane transport and redox adjustment-related genes, while under high-salinity stress conditions, secondary metabolism became the driving force behind cellular homeostasis. These results showed that salt stress affects A. glaucus ‘CCHA’.

Gene ontology and pathway classification enrichment analyses of DEGs

On the basis of the gene ontology (GO) analysis, an internationally standardised gene functional classification system (Ashburner et al., 2000), DEGs were classified into the three major functional categories of biological process, cellular component, and molecular function (p-value < 0.05; Table S2). The GO enrichment analysis based on Sample36h_3M vs. Sample36h_1M showed that in the biological process classification, oxidation–reduction progress, melanin metabolic progress and transmembrane transport were the highly represented categories. For molecular function, oxidoreductase activity was the most highly represented category. Among the DEGs in Sample84h_0M vs. Sample84h_1M, for the biological process classification, oxidation–reduction progress was the most highly represented category. For cellular component, intrinsic component of membrane figured prominently. In Sample84h_3M vs. Sample84h_1M, for the biological process classification, transmembrane transport was the most highly represented category. For molecular function, oxidoreductase activity was the most highly represented category. In Sample84h_4M vs. Sample84h_1M, for the biological process classification, transmembrane transport was the most highly represented category. After comparing these groups, we inferred that the salt condition affected the redox state and transmembrane transport of A. glaucus.

To evaluate the functions of DEGs in responses to salt stress, the KEGG database was used to analyse pathways. The results of the KEGG pathway enrichment analysis are shown in Table S3. For Sample36h_3M vs. Sample36h_1M, the DEGs were most highly represented in ‘starch and sucrose metabolism’ and ‘tyrosine metabolism’ pathways (Table S2). For Sample84h_0M vs. Sample84h_1M, the most significant pathways were identified as ‘pentose and glucoronate interconversions’ and ‘starch and sucrose metabolism’. Salt stress had significant effects on multiple pathways in A. glaucus, suggesting that these pathways and associated processes might be critical for developmental and metabolic variations in responses to salt stress.

To verify the effectiveness of annotation, a COG classification of up regulated genes for Sample36h_3M vs. Sample36h_1M was performed, resulting in 20 COG groups. Of these, carbohydrate transport and metabolism (10.8%) was the dominant group, followed by translation, ribosomal structure and biogenesis (9.39%) and Amino acid transport and metabolism (8.92%). These results provided multiple candidate genes implicated in responding to high-salt stress in A. glaucus ‘CCHA’.

Isolation and characterisation of salt-tolerance-related genes

We previously observed that a 0.75-M NaCl treatment caused a severe growth inhibition in E. coli ‘BL21’ Star (DE3) cells. Therefore, we hypothesised that one way to identify genes involved in the activation of salt-tolerance in A. glaucus ‘CCHA’ was to identify genes that caused growth-inhibition recovery when heterologously expressed in E. coli ‘BL21’ (Fang et al., 2014). We hypothesised that mRNAs for tolerance inducers most likely would be expressed highly after exposure to stress and that mRNAs for many regulatory genes would be present at high levels relative to other genes. Consequently, we concluded that the best way to assure randomness in our expression library would be to fuse salt-induced cDNAs to pYES2.0. In this manner, we examined all of the transformants and identified five conditional growth-inhibited recovery mutants (Fang et al., 2014). Then, we verified the salt-stress tolerance of candidate genes in A. thaliana using Agrobacterium-mediated transformations. Ag_2114, which encodes a DJ-1/PfpI family protein, confers salt tolerance to A. thaliana in transgenic plants (Fig. 9). Three transgenic lines and a WT control line were used for physiological studies under salt-stress conditions. In total, greater damage was observed on the leaves of the WT A. thaliana than on the leaves of the transgenic lines when treated with 200 mm NaCl. At all the salt-stress durations, the SOD activity and total chlorophyll, chlorophyll a and chlorophyll b levels, in CCHA-2114-overexpressing plants were significantly greater than those in control plants. Moreover, CCHA-2114’s overexpression led to reductions in the MDA level, lipid peroxidation and ion leakage compared with the untransformed line and thus likely increased the tolerance to salt stress (Figs. 10A–10F).