Comparative transcriptome analysis of Armillaria gallica 012m in response to ethephon treatment

Fungi growth and culture conditions

A. gallica 012m was derived from G. elata in our previous research. The stock strain was routinely grown on modified Czapek-Dox medium agar (MgSO₄, 0.5 g; FeSO₄, 0.01 g; KCl, 0.5 g; NaNO₃, 3 g; K₂HPO₄, 1 g; sucrose, 30 g; malt extract, 10 g; yeast extract, 10 g; ethanol, 20 g) in the dark at 25 °C for 15 days (Zhan et al., 2020). Exiguous ETH (0.1 ppm, 2.0 ppm, 5.0 ppm) were added to the medium of broth and ager; and set up the BK (blank control) groups. Setting three parallel repeats for each different treatment condition.

The mycelium of A. gallica 012m was grown on liquid media in a 160 r/min shaker at 25 °C for 15 days. Then, the pellets were filtrated with a Buchner funnel, washed with pure water, collected in beaker and weighted with electronic balance. By utilizing independent sample T-test, the data analysis was performed with IBM SPSS v23 (Tanner-Smith & Tipton, 2014). The mycelium was preserved under liquid nitrogen conditions.

RNA extraction and sequencing

The method was as same as our previous work (Cao et al., 2022). RNA was extracted from fungi pellets using the RNeasy mini kit (Qiagen, Hiden, Germany) following the manufacturer’s instructions. The cDNA library was sequenced on the Illumina HiSeq platform (Caporaso et al., 2012) with a double-end sequencing strategy in Novogene Bioinformatics Technology, Beijing, China. The original data were deposited in the National Center for Biotechnology Information database with the accession number PRJNA759758.

Transcriptome sequences data quality control and comparison

To obtain clean transcriptome data, low-quality sequences of RNA-seq were removed with Trimmomatic v0.36 (Bolger, Lohse & Usadel, 2014). Then, the data quality control was evaluated with FastQC v0.11.9 (Brown, Pirrung & McCue, 2017). The reads of RNA-Seq was aligned with the A. gallica 012m genome sequences (https://www.ncbi.nlm.nih.gov/genome/57439?genome_assembly_id=853036) by using Hisat2 v2.1.0 (Kim et al., 2019) to obtain the SAM data.

Differentially expressed gene analysis

Gene expression values were performed according to the reads per kilobase per million mapped reads (RPKM) method. The read count matrix was obtained for expression qualification with StringTie v2.1.0 (Pertea et al., 2015). The read count matrix was imported into R 3.6.3. The read count matrix was imported into R 3.6.3, and the differential gene analysis was carried out with edgeR of R package under an FDR<0.05 and —log2FC—>2. Next, RNASeqPower (DOI: http://dx.doi.org/10.18129/B9.bioc.RNASeqPower), a power analysis calculation software, was used to calculate the statistical power of this experimental design, and the statistical power is 0.81. Then, all transcripts and their corresponding genes were compared by emapper v2.1.3 for functional annotation and classification (Huerta-Cepas et al., 2017). The result of the GO function analysis was performed by using TBtools V0.66836 (Chen et al., 2020). Go terms visualization of DEGs was executed with WEGO (https://wego.genomics.cn/).

Functional analysis of growth-related genes of DEGs

To analyze the function of protein domain, the gene sequences involved biological process functions were extracted and their protein domains were annotated by using the SMART platform (http://smart.embl-heidelberg.de/).

Screened candidate ethylene receptors

Those existing ET receptor protein sequence of fungi were downloaded the from NCBI (https://www.ncbi.nlm.nih.gov/) and compared with the genome protein sequences of A. gallica 012m. Those expressed sequences with E-value ≤ le-5 were selected and detected their conserved domain on SMART (Letunic, Khedkar & Bork, 2021).

Construction of receptor sequence evolutionary tree

An evolutionary tree was constructed with downloaded receptor sequences and screened sequences of A. gallica 012m. ClustalW2 of MEGA7 was used for multiple amino acid sequences alignment (Kumar, Stecher & Tamura, 2016). The phylogenetic tree was constructed by neighbor-joining method and the number of calculations was 1,000.

Morphological characteristics of A. gallica grown under ETH

As shown in Fig. 1, ETH stimulated the mycelial growth of A. gallica 012m in solid media. Moreover, lower concentrations of ETH present better effect on mycelium elongation of A. gallica 012m. The order of influence of mycelium elongation were ETH−0.1 ppm>ETH-2 ppm>ETH-5 ppm.

Effects of plant growth substances on the biomass of A. gallica 012m

The liquid culture was carried out to explore whether or not ETH affected the biomass of A. gallica 012m. As shown in Fig. 2 and Table S1, the biomass of mycelium increased by 88.0 ± 9.1% under ETH−0.1ppm, the biomass of mycelium increased by 66.1 ± 7.9% under ETH-2ppm, the biomass of mycelium decreased by 86.8 ± 5.0% under ETH-5ppm. By utilizing independent sample T-test, it can be considered that the biomass of A. gallica 012m were increased extremely significant (p < 0.01) under ETH−0.1 ppm and ETH-2ppm, while decreased extremely significant (p < 0.01) under ETH-5ppm.

Evaluation of transcriptome sequencing data

Above 10.0 Gb of raw data per sample was obtained by transcriptome sequencing on the Illumina HiSeq platform and could be used to further expression level analysis after quality control. In addition, twelve transcriptome samples of ETH and BK were sequenced by Illumina HiSeq platform to obtain 21,016,459; 25,437,814; 22,314,980; 26,576,772; 24,979,536; 23,865,405; 30,433,600; 21,398,103; 22,659,153; 22,659,153; 27,395,065; 28,784,179 and 26,125,609 pairs of PE reads (Table 1).

After removing the reads with adapter sequences of low quality, an average of 21,916,601 to 30,433,600 pairs of clean reads were retained from ETH and BK, respectively. Read mapped percentage of clean data from all sample were higher than 92.73%, and reads mapped percentages of ETH-2ppm-2 was the highest. What’s more, 92.73% ∼93.82% pure readings were successfully mapped to the reference genome of A. gallica 012m.

Enriched GO terms of up-regulatedand down-regulated at DEGs

In the ETH−0.1 ppm vs. BK comparison, a total of 118 genes were differentially expressed, including 15 up-regulated genes and 103 down-regulated genes (Fig. 3, Table S2). In the ETH-2 ppm vs. BK comparison, a total of 341 genes were differentially expressed, including 224 up-regulation genes and 117 down-regulation genes (Table S3). In the ETH-5 ppm vs. BK comparison, a total of 696 genes were differentially expressed, including 314 up-regulated genes and 382 down-regulated genes (Table S4). Interestingly, the total number of DEGs in the experimental group increased along with ETH concentration.

The gene expression of A. gallica 012m with ETH−0.1 ppm vs. BK, ETH−2.0 ppm and ETH−5.0 ppm vs. BK was analyzed. In the up-regulated genes, the results showed that there were three differential genes shared by ETH−0.1 ppm and ETH−2.0 ppm, and 58 differential genes shared by ETH−2.0 ppm and ETH−5.0 ppm, 11 DEGs in ETH−0.1 ppm group alone, 163 DEGs in ETH−2.0 ppm group alone, and 255 DEGs in ETH−5.0 ppm group alone. In down-regulation genes, the DEGs of the ETH−0.1 ppm group were the same as the ETH−2.0 ppm group more than 15% (Fig. 4).

Enriched GO terms of up-regulated and down-regulated at DEGs

Comparing with the GO database, classification and functional analysis of DEGs were performed for better visualization. The result of the annotation of Level 2 for the GO database was shown in Fig. 5.

At ETH−0.1 ppm vs BK, a total of 90 DEGs were categorized into 22 functional groups under CC (10 functional groups), MF (three functional groups) and BP (nine functional groups). In the CC category, the major terms included ‘cell’ (upregulation: downregulation = 3: 14), ‘cell part’ (upregulation: downregulation = 3: 14), ‘organelle’ (upregulation: downregulation = 3: 14), ‘organelle part’ (upregulation: downregulation = 3:11) and ‘membrane’ (upregulation: downregulation = 2: 9). In the MF category, the GO terms only included three functional groups, which were ‘catalytic activity (upregulation: downregulation = 2: 6), ‘transporter activity (upregulation: downregulation = 0: 4) and ‘binding’ (upregulation: downregulation = 1: 0). In the BP category, the representative GO terms included ‘cellular process’ (upregulation: downregulation = 3: 10), ‘metabolic process’ (upregulation: downregulation = 3: 7), ‘biological process’, ‘response to stimulus’ (upregulation: downregulation = 3: 2) and ‘regulation of biological process’ (upregulation: downregulation = 3: 2).

At ETH-2 ppm vs BK, a total of 269 DEGs were categorized into 34 functional groups under CC (10 functional t groups), MF (nine functional groups) and BP (15 functional groups). In the CC category, the dominant GO terms included ‘cell’ (upregulation: downregulation = 18: 10), ‘cell part’ (upregulation: downregulation = 18: 10), ‘organelle’ (upregulation: downregulation = 14: 9), ‘organelle part’ (upregulation: downregulation = 4: 7) and ‘membrane’ (upregulation: downregulation = 6: 7). In the MF category, the major GO terms included ‘catalytic activity (upregulation: downregulation = 19: 6), ‘transporter activity’ (upregulation: downregulation = 4: 1), ‘binding’ (upregulation: downregulation = 9: 1) and ‘antioxidant activity’ (upregulation: downregulation = 3: 0). In the BP category, the major GO terms included ‘metabolic process’ (upregulation: downregulation = 21: 7), ‘cellular process’ (upregulation: downregulation = 22: 7), ‘response to stimulus’ (upregulation: downregulation = 11: 4), ‘localization’ (upregulation: downregulation = 5: 4) and ‘biological regulation’ (upregulation: downregulation = 7: 3).

At ETH-5 ppm vs. BK, a total of 534 DEGs were categorized into 41 functional groups under CC (11 functional groups), MF (10 functional groups) and BP (20 functional groups). In the CC category, the major GO terms included ‘cell’ (upregulation: downregulation = 48: 31), ‘cell part’ (upregulation: downregulation = 48: 31), ‘organelle’ (upregulation: downregulation = 33: 26), ‘organelle part’ (upregulation: downregulation = 16: 13) and ‘membrane’ (upregulation: downregulation = 24: 18). In the MF category, the representative GO terms included ‘catalytic activity (upregulation: downregulation = 34: 21), ‘transporter activity’ (upregulation: downregulation = 4: 7), ‘binding’ (upregulation: downregulation = 11: 9) and ‘antioxidant activity’ (upregulation: downregulation = 2: 1). In the BP category, the main GO terms included ‘metabolic process’ (upregulation: downregulation = 45: 22), ‘cellular process’ (upregulation: downregulation = 47: 22), ‘biological regulation’ (upregulation: downregulation = 19: 9) and ‘regulation of biological process’ (upregulation: downregulation = 14: 5). These result demonstrated that the main biological functions of the genes expressed in the A. gallica 012m transcriptome.

Identification of genes involved in the growth of A. gallica 012m

To further explore the mechanism of effects on the growth of A. gallica 012m, we respectively compared the biological regulation process of DEGs. The result can be shown as follows: in the ETH−0.1 ppm vs BK, Pyr_redox (pyridine nucleotide-disulphide oxidoreductase) domain was predicted in Armga012mGene24786 which showed up-regulation expressed. This domain is actually a small NADH binding domain within a larger FAD binding domain. This domain exists in NADH oxidases, peroxidases, class I and class II oxidoreductases. In the ETH-2 ppm vs. BK, HMG (High Mobility Group) and GAL4 (C6 zinc finger)/Fungal_trans domain were predicted in Armga012mGene25682 and Armga012mGene16364, respectively and showed up-regulation expressed. HMG-box domains form a large, diverse family involved in the regulation of transcription, replication and strand repair. Gal4 is a positive regulator for the gene expression of the galactose-induced genes of S. cerevisiae. This domain is found in various fungal transcription factors, which regulate cellular and metabolic processes. In the ETH-5 ppm vs. BK, ZnF_C2H2 (C2H2 zinc finger) and ribosomal protein S4 domain were predicted in Armga012mGene07208 and Armga012mGene07627 which showed down-regulation expressed. Znf-containing proteins function in gene transcription, translation, mRNA trafficking, cytoskeleton organisation, epithelial development, cell adhesion, protein folding, chromatin remodelling and zinc sensing. Ribosomal protein S4 is a kinds of small proteins that may be involved in translation regulation.

Prediction and expression of ET receptors in A. gallica 012m

ET receptor is the first element of the ethylene biological effect, and its binding with ethylene can activate downstream ethylene signal transduction. We downloaded the sequence information of ET receptors of fungi from GenBank with a total of 54 sequences, compared it with the genome protein sequence of A. gallica 012m. A total of seven speculated ET receptor domains of A. gallica 012m were predicted by using SMART (Table 2) and expressed. The predicted ET receptor proteins were annotated, and only Armga012mGene13219 has transmembrane region domains. Those predicted ET proteins have similar domain as determined ethylene receptor proteins, and the expression levels of candidate ethylene receptor protein sequences were different (Fig. 6) in A. gallica 012m. Thus, we speculated that exogenous ethylene affected the growth of A. gallica 012 m through ET receptors.

Table 2:

ET receptor protein and domain of A. gallica 012m.

DOI: 10.7717/peerj.14714/table-2

Notes:

ER

endoplasmic reticulum

HisKA

Histidine Kinase A (dimerisation and phosphoacceptor) domain

HATPase_c

Histidine kinase-like ATPase catalytic (Histidine kinase-, DNA gyrase B-, phytochrome-like ATPases) domain

GAF

cGMP-phosphodiesterase, adenylyl cyclase, FhIA domain

P

PAS/PAC (Per-period circadian protein; Arnt-Ah receptor nuclear translocator protein Sim-single minded protein/PAS C terminus)

REC

Receiver domain

GGDEF

diguanylate cyclase domain

EAL

c-di-GMP phosphodiesterase

STYKc

protein kinase–unclassified specificity.

Genes containing ET receptor domain in A. gallica 012m

We downloaded the sequence information of ET receptors in fungi, bacteria, Arabidopsis thaliana from GenBank with a total of 10 sequences (File S3). Finally, by comparing 7 species with the ET receptor domain of A. gallica 012m, an analysis of the phylogenetic relationship is shown in Fig. 7. The result showed that A. gallica 012m (Armgao012mGene04732) had a high homology correlation with Verticillium alfalfae VaMs.102 (GenBank: XP_003000814.1, EEY23199.1). What’s more, A. gallica 012m (Armgao012mGene00417) had homology with Purpureocillum lilacinum (GenBank: XP_018174038.1). The phylogenetic tree analysis inferred that A. gallica 012m might possess ET receptor domains.