Genome-wide characterization of the Zn(II)₂Cys₆ zinc cluster-encoding gene family in Pleurotus ostreatus and expression analyses of this family during developmental stages and under heat stress

Identification of zinc cluster genes

The zinc cluster domain (PF00172) was downloaded from the Pfam database (http://pfam.xfam.org/) (El-Gebali et al., 2018). Then, an HMM search was performed using the hmmsearch tools embedded in HMMER3.21 (Wheeler & Eddy, 2013) with default parameters based on the protein sequences of P. ostreatus strain CCMSSC0389 (China Center for Mushroom Spawn Standards and Control, Beijing, China) (Qu et al., 2016). The zinc cluster domain region of each potential positive sequence was extracted and multialigned using the Clustal Omega online tool (Madeira et al., 2019). Sequences without complete six cysteine residues in the zinc cluster domain were regarded as false positive sequences and removed. The gene structures of the remaining sequences were corrected by checking the mapped reads positions in the RNA-seq data using the Integrative Genomics Viewer tool (Robinson et al., 2011). After manual removal of false positive sequences and correction, the remaining sequences were submitted to the NCBI Conserved Domain Database (https://www.ncbi.nlm.nih.gov/cdd) and the PFAM database to determine the presence and integrity of the zinc cluster domain. Transcription factor annotation was performed by the hmmsearch tool (E-value cutoff = 0.01) using 52 transcription factor Pfam models (Table S1).

Gene structure and protein analyses of the zinc cluster family

Gene structures were obtained by exon and intron information from the P. ostreatus genome. PoZCPs were classified according to the number of amino acid residues between each two cysteine residues in the typical zinc cluster domain. The isoelectric point (pI) and molecular weight (MW) of PoZCPs were estimated using IPC tools (Kozlowski, 2016). The subcellular localization of each PoZCP was predicted by BaCelLo online tools (Pierleoni et al., 2006).

Sequence alignment, motif analysis and phylogenetic tree construction of the zinc cluster family

The domain region of the zinc cluster sequences was characterized by the NCBI–CDD and Pfam results and classified by the amino acid counts of variable regions. Multiple sequence alignments of the protein sequences of PoZCPs were carried out using the Clustal Omega online tool (Madeira et al., 2019), and the poorly aligned regions were removed by the alignment utility trimAl (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009). The aligned sequences were used to construct a phylogenetic tree based on the maximum likelihood (ML) method by IQ-TREE using the best fitting model and branch support computed with 1000 UltraFast Bootstrap (Trifinopoulos et al., 2016). Conserved motifs in the PoZCP family were assessed by the MEME online program (Bailey et al., 2009). The gene structure and motif distribution of the PoZCP family were displayed using Tbtools (Chen et al., 2018).

RNA-sequencing and analysis

RNA-sequencing (RNA-seq) was conducted to gain insight into the expression patterns of the zinc cluster genes at several developmental stages and heat stress. Fruiting bodies of the dikaryotic strain CCMSSC0389 were grown, and samples, including mycelia, primordia and fruiting bodies, were collected at different developmental stages as described previously (Wang et al., 2017). The mycelia were inoculated on PDA medium and cultured at 28 °C for 6 days and then transferred to 35 °C for 0.5, 1, 2 and 3 h, respectively, for heat stress treatment. Three replicates for each sample were frozen in liquid nitrogen immediately after collection.

Sample collections, RNA extraction and RNA-seq were carried out as previously described (Wang et al., 2018b), except that TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) was used for sequencing library construction and the Illumina HiSeq 2500 platform was used for sequencing. Transcript expression of the zinc cluster genes was analyzed according to the RNA-seq profiles using fragments per kilobase per million mapped reads (FPKM) values. A heatmap was constructed using Tbtools with log and row scales (Chen et al., 2018).

Validation of gene expression by RT-qPCR

To validate the RNA-seq results, we performed an RT-qPCR analysis of 13 genes that showed different expression patterns and large fold changes during various developmental stages or under heat stress. The experimental procedures were described in previous studies (Wang et al., 2017, 2018b). Beta-actin and beta-tubulin were used as reference genes. The primer sequences are listed in Table S2. The relative expression level of each gene was calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Identification of ZCP genes in P. ostreatus

After manual correction, a total of 66 unique genes with conserved zinc cluster domains were identified as ZCP genes in P. ostreatus, accounting for 23.82% of all transcription factors. The 66 genes were named PoZCP01 to PoZCP66, and the coding sequences were submitted to GenBank under accession numbers MN652925–MN652990. Fifty-nine (89.39%) family members contain zinc cluster domains in the N-terminus (front 1/3 part of the protein), three (4.55%) family members have zinc cluster domains in the middle (1/3–2/3 of the protein) of the protein, and four (6.06%) family members possess zinc cluster domains in the C-terminus (last 1/3 part of the protein) (Fig. 1). Fungal-specific transcription factor domains (Fungal_trans, PF04082) are present in 43 (65.15%) family members, and C2H2-type zinc transcription factor domains (PF00096) are also present in three family members. One gene, PoZCP55, encodes three septin domains (PF00735), which are reported to be involved in cytokinesis and septum formation (Longtine et al., 1996), with two zinc cluster domains in the C-terminus.

Classification and characterization of ZCP genes in P. ostreatus

Previous studies in yeast have revealed that in the zinc cluster domain, the numbers of amino acids between C1–C2, C2–C3 and C4–C5 are always 2, 6 and 2, respectively, but that numbers between C3–C4 and C5–C6 vary (Todd & Andrianopoulos, 1997). The PoZCP family was classified into 15 subgroups named types A-O based on the different patterns between C3–C4 and C4–C5 (Table 1), 10 of which belong to the typical zinc cluster motif C–X₂–C–X₆–C–X_5-12–C–X₂–C–X_6-8–C in yeast (MacPherson, Larochelle & Turcotte, 2006). However, five other patterns, including C–X₂–C–X₆–C–X₅–C–X₂–C–X₄–C (type A), C–X₂–C–X₆–C–X₆–C–X₂–C–X₉ (type G), C–X₂–C–X₆–C–X₁₃–C–X₂–C–X₆ (type M), C–X₂–C–X₆–C–X₁₄–C–X₂–C–X₆ (type N) and C–X₂–C–X₆–C–X₁₅–C–X₂–C–X₆ (type O), were also found in the PoZCP family (Table 1). The motifs of 15 types of zinc cluster domains were analyzed, and the motif logos of the seven major types (more than three members) are shown in Fig. 2. For types A, D, F, I and J, amino acids are conserved at position 5 (arginine) and position 8 (lysine); for types A, B, C and F, amino acids at position 12 (aspartate) and position 16 (proline) are conserved (Fig. 2).

The characteristics, including exon and intron length, protein length, isoelectric point (pI), molecular weight (MW) and predicted subcellular localization by BaCelLo online tools are shown in detail in Table S3. The protein length of putative PoZCP genes ranges from 195 to 1,386 amino acids, with MWs ranging from 21 to 152 kDa and theoretical pIs ranging from 4.86 to 9.45, revealing significant differences among PoZCPs in terms of their physical and chemical properties. According to the predicted subcellular localization, 62.12% of PoZCPs are located in the nucleus, 31.34% in the cytoplasm and 6.06% in secretory pathways (Table S3).

Phylogenetic, motif and structure analysis of PoZCPs

A phylogenetic tree was constructed using amino acid sequences of the 66 PoZCP family members. The family members were divided into six well-conserved groups (I–VI) with significant bootstrap support according to phylogenetic relationships (Fig. 3A). PoZCP55 forms a separate clade, which imply its self-dependent origin. The PoZCP family members belonging to types A, B, C and F (Table 1) clustered respectively, suggesting genes in each class may have expanded from a single gene. Conversely, PoZCP family members of types D, I and J (Table 1) were found to be distributed across different clades, indicating that these categories may have different ancestors.

To further investigate the diversity of the ZCP in P. ostreatus, we analyzed PoZCP motifs using the MEME online server. Ten conserved motifs named motif 1 to motif 10 were identified. Motif 1, which forms the zinc cluster domain, was the most common motif found in all family members. The combinations of different motifs differed among the groups (Fig. 3A). For example, the combination of motif 3, motif 7 and motif 9 occurs in most proteins of group II, the combination of motifs 7–9 appeared in most proteins of group IV, the combination of motif 3, motif 7 and motif 8 existed in all members of group V and motifs 2–6 were found in all members of group VI (Fig. 3B). These results suggest that the function of the PoZCPs of different groups may differ.

We deduced the gene structures of individual PoZCP genes to examine the relationship between their evolutionary origins and structural diversity. The results showed that the exon number of PoZCP genes changed over a wide range (from 2 to 23); however, each group in Fig. 3A is relatively conserved in gene structure (Fig. 3C). For example, group II consists of genes of approximately 2,500 bp with similar architecture and 8–10 exons that encode approximately 600 amino acids. Group III is composed of short genes of approximately 1,600 bp that encode approximately 450 amino acids. Group IV and group VI are rich in exons, numbering approximately 14, which is significantly greater than the average exon number of 9 for all family members. Most of the genes in group VI contain many short exons of fewer than 50 bp.

Expression patterns of PoZCPs during different developmental stages and under heat stress

To gain insight into the expression profiles of PoZCP genes during different developmental stages and under heat stress at different time points, RNA-seq data obtained from different developmental stages as well as different heat treatment times were used to analyze expression patterns (Tables S4 and S5). PoZCP14 was not expressed in any expression profile, indicating that it may be a pseudogene.

The expression patterns during different developmental stages were categorized into four clusters according to RNA-seq profiles (Fig. 4A). Genes in cluster I and cluster IV were downregulated and upregulated, respectively, from the mycelial to fruiting body development stages (Fig. 4A), which suggests that these genes may act as negative or positive regulators during development. Genes in cluster II and cluster III were only upregulated in the fruiting body stage and primordia stage, respectively (Fig. 4A), therefore, these genes might be involved in the stage-specific process.

The expression patterns under heat stress at different time points were also categorized into four groups (Fig. 4B). Clusters I and IV were downregulated after heat stress treatment at different time points, which indicates their negative regulatory ability under long-term (2–3 h) and short-term (0.5–1 h) heat stress, respectively (Fig. 4B). Most genes in clusters II and III were upregulated after heat stress treatment at different time points, revealing their positive regulatory ability under short-term (0.5–1 h) and long-term (2–3 h) heat stress, respectively (Fig. 4B).

To verify the gene expression results for the identified PoZCPs, thirteen genes showing distinct variation at various developmental stages or heat treatment in RNA-seq analysis were selected for detection of their expression levels via RT-qPCR. All tested genes showed the same tendency as the transcript profiles (Fig. 5). PoZCP10, PoZCP26, PoZCP31, and PoZCP40 were dramatically upregulated (20–90-fold) in the fruiting body stage, suggesting that they are stage-specific genes. PoZCP15, PoZCP26, and PoZCP31 were markedly upregulated (5–140-fold) under 0.5–1 h heat stress (Fig. 5), which indicates their important roles in the response to short-term heat stress. PoZCP26 and PoZCP31 were upregulated during both developmental stages and under heat stress treatment and thus may have multiple functions in development and the heat stress response.