Genome-Wide Identification and Analysis of Chitinase GH18 Gene Family in Mycogone perniciosa

Mycogone perniciosa causes wet bubble disease in Agaricus bisporus and various Agaricomycetes species. In a previous work, we identified 41 GH18 chitinase genes and other pathogenicity-related genes in the genome of M. perniciosa Hp10. Chitinases are enzymes that degrade chitin, and they have diverse functions in nutrition, morphogenesis, and pathogenesis. However, these important genes in M. perniciosa have not been fully characterized, and their functions remain unclear. Here, we performed a genome-wide analysis of M. perniciosa GH18 genes and analyzed the transcriptome profiles and GH18 expression patterns in M. perniciosa during the time course of infection in A. bisporus. Phylogenetic analysis of the 41 GH18 genes with those of 15 other species showed that the genes were clustered into three groups and eight subgroups based on their conserved domains. The GH18 genes clustered in the same group shared different gene structures but had the same protein motifs. All GH18 genes were localized in different organelles, were unevenly distributed on 11 contigs, and had orthologs in the other 13 species. Twelve duplication events were identified, and these had undergone both positive and purifying selection. The transcriptome analyses revealed that numerous genes, including transporters, cell wall degrading enzymes (CWDEs), cytochrome P450, pathogenicity-related genes, secondary metabolites, and transcription factors, were significantly upregulated at different stages of M. perniciosa Hp10 infection of A. bisporus. Twenty-three out of the 41 GH18 genes were differentially expressed. The expression patterns of the 23 GH18 genes were different and were significantly expressed from 3 days post-inoculation of M. perniciosa Hp10 in A. bisporus. Five differentially expressed GH18 genes were selected for RT-PCR and gene cloning to verify RNA-seq data accuracy. The results showed that those genes were successively expressed in different infection stages, consistent with the previous sequencing results. Our study provides a comprehensive analysis of pathogenicity-related and GH18 chitinase genes’ influence on M. perniciosa mycoparasitism of A. bisporus. Our findings may serve as a basis for further studies of M. perniciosa mycoparasitism, and the results have potential value for improving resistance in A. bisporus and developing efficient disease-management strategies to mitigate wet bubble disease.

Mycogone perniciosa causes wet bubble disease in Agaricus bisporus and various Agaricomycetes species. In a previous work, we identified 41 GH18 chitinase genes and other pathogenicity-related genes in the genome of M. perniciosa Hp10. Chitinases are enzymes that degrade chitin, and they have diverse functions in nutrition, morphogenesis, and pathogenesis. However, these important genes in M. perniciosa have not been fully characterized, and their functions remain unclear. Here, we performed a genome-wide analysis of M. perniciosa GH18 genes and analyzed the transcriptome profiles and GH18 expression patterns in M. perniciosa during the time course of infection in A. bisporus. Phylogenetic analysis of the 41 GH18 genes with those of 15 other species showed that the genes were clustered into three groups and eight subgroups based on their conserved domains. The GH18 genes clustered in the same group shared different gene structures but had the same protein motifs. All GH18 genes were localized in different organelles, were unevenly distributed on 11 contigs, and had orthologs in the other 13 species. Twelve duplication events were identified, and these had undergone both positive and purifying selection. The transcriptome analyses revealed that numerous genes, including transporters, cell wall degrading enzymes (CWDEs), cytochrome P450, pathogenicity-related genes, secondary metabolites, and transcription factors, were significantly upregulated at different stages of M. perniciosa Hp10 infection of A. bisporus. Twenty-three out of the 41 GH18 genes were differentially expressed. The expression patterns of the 23 GH18 genes were different and were significantly expressed from 3 days post-inoculation of M. perniciosa Hp10 in A. bisporus. Five differentially expressed GH18 genes were selected for RT-PCR and gene cloning to verify RNA-seq data accuracy. The results showed that those genes were successively expressed in different infection stages, consistent with the previous sequencing results. Our study provides a comprehensive analysis of pathogenicity-related and GH18 chitinase genes' influence on M. perniciosa

INTRODUCTION
The button mushroom (Agaricus bisporus) is one of the most widely cultivated and consumed edible mushrooms in the world. Production of button mushrooms in China has rapidly increased in recent years as a result of the expanded area of cultivation and the adoption of technology for improved commercial cultivation (He et al., 2014;Sonnenberg et al., 2017;Li et al., 2018;McGee, 2018). However, diseases caused by fungi, bacteria, and viruses are major constraints to A. bisporus production worldwide, often leading to serious economic losses (Fletcher et al., 1989;Largeteau and Savoie, 2010;Kouser et al., 2013). Wet bubble disease (WBD) is one of the most devastating diseases of A. bisporus, causing yield losses of 15-30% under favorable conditions and up to 75% or total crop loss in the most severe cases (Zhou et al., 2015). WBD is characterized by wet bubbles, malformation, white, fluffy mycelial growth, copious amber droplets (diseased carpophores exuding a brown malodorous liquid), and flocculent mycelia on most substrates (Fletcher et al., 1995;Sharma and Kumar, 2000;Umar et al., 2000). Mycogone perniciosa (teleomorph: Hypomyces perniciosus) is the causal agent of WBD, and it infects a variety of mushrooms (Zhang et al., 2017a,b;Carrasco et al., 2019). M. perniciosa is a fungicolous fungus belonging to the order Hypocreales (Ascomycota) in the family Hypocreaceae.
M. perniciosa is mainly controlled by cultural practices and the application of fungicides (Brankica et al., 2009). Varying levels of WBD resistance have been identified in the A. bisporus germplasm collection in China; however, no major resistance gene has been identified (Fu et al., 2016). The genome of M. perniciosa contains many genes implicated in pathogenicity (Li et al., 2019), but the regulation of these genes in the pathogenesis toward A. bisporus is still unclear. In addition, comparative genomics analysis of M. perniciosa has revealed gene expansion and positive selection of many genes, including GH18 chitinase, peptidase, and secondary metabolite genes (Li et al., 2019). Gene expansion and positive selection contribute to the evolution of virulence genes in microbial pathogens and to the adaptation to different environmental niches through the infection process and via escape from the host defense response (Yoshizaki et al., 2019).
do not contain chitin (Rathore and Gupta, 2015). Chitinases of pathogenic fungi not only play vital roles in spore germination, septum formation, cell division, and morphogenesis, but the enzymes are also important in the host interaction (Elad et al., 1982(Elad et al., , 1983Inbar and Chet, 1995;Chen, 2002;Adams, 2004). In addition to degradation of the host fungal cell wall, chitinases also inhibit hyphae growth and bud tube elongation (Gozia et al., 1993;Stressmann et al., 2004).
Omics and bioinformatics analyses have demonstrated that most fungal chitinases have similar domains that generally contain a signal peptide sequence, a chitinase catalytic domain, a chitinbinding domain, and a short C-terminal domain (van Aalten et al., 2001). Among fungal pathogens, GH18 chitinase gene families are well characterized in Trichoderma species, Fusarium species, and Magnaporthe species (Häkkinen et al., 2012;Han et al., 2019). Currently, 30 chitinase genes have been reported within eight species of Trichoderma, including T. harzianum (30 genes), T. virens (29 genes), and T. atroviride (24 genes) (Kubicek et al., 2011(Kubicek et al., , 2019. The advent of next-generation sequencing technologies has increased the scalability, speed, and resolution of genomic sequencing and reduced genome sequencing cost (Faino et al., 2015). This has rapidly increased fungal genome availability for comparative genomics and genome-wide identification of gene families (Bartholomew et al., 2019). However, few studies have comprehensively examined the structure or the expression of GH18 chitinases of fungal pathogens infecting mushrooms.
In a previous work, we identified 41 GH18 chitinase genes in the genome of M. perniciosa (Li et al., 2019). However, the genome-wide identification, functional characterization, and expression of these GH18 chitinase genes were not considered. In this study, we present the first detailed and comprehensive analysis of the GH18 gene family in the genome of M. perniciosa. The analyses include chromosome location, phylogenetic analysis, protein structure, and motif composition. Furthermore, we performed a transcriptome analysis of the different infection stages of M. perniciosa on A. bisporus to identify pathogenicity-related differential gene expression and expression patterns of the GH18 genes. In addition, we selected five GH18 genes that represented different types and identified their expression levels by RT-PCR and gene cloning to verify the sequencing. The results can provide important information for identifying essential genes as potential antifungal targets in M. perniciosa. Further, the dataset generated in this study may provide a basis for identifying candidate resistant genes in A. bisporus against M. perniciosa and lay a foundation for further research to improve the resistance of A. bisporus to M. perniciosa.

Strains and Culture Conditions
The fungal strains [A. bisporus strain A3 (CCMJ1009) and M. perniciosa Hp10] were obtained from the Engineering Research Center of Edible and Medicinal Fungi, Ministry of Education, Jilin Agricultural University (Changchun, Jilin, China). The M. perniciosa Hp10 used in this study is a highly pathogenic strain (≥90%) able to cause severe disease on all A. bisporus evaluated to date (Li et al., 2019). All the fungal strains were maintained on potato dextrose agar (PDA) at 25 • C. Escherichia coli DH5 α and GV pxt19-t vector were purchased from Beijing TransGen Biotech Co., Ltd., and Beijing Dingguo Changsheng Biotechnology Co., Ltd., respectively.

Mushroom Cultivation, Inoculum Preparation, and Disease Evaluation
Cultivation of A. bisporus and evaluation of the WBD infection process were conducted at the Mushroom Base of Jilin Agricultural University, Changchun, China, using methods described by Li et al. (2019). A. bisporus mycelia were inoculated on autoclaved wheat grains to produce spawn. The spawn was inoculated on a (45 × 33 × 25 cm) basket filled with 7.5 kg compost. After the mycelia overgrew the compost, 4 cm thick casing soil was applied to cover the compost. To induce fruiting, the room temperature, relative humidity, and carbon dioxide (CO 2 ) concentration were set at 15-18 • C, 80-95%, and 1,200-1,500 ppm, respectively.
Spore suspensions of M. perniciosa Hp10 inoculum were from 7-day-old pure PDA cultures grown at 25 • C. The cultures were suspended in 5 ml sterile distilled water (SDW), gently scraped with a glass stick, and filtered through two cheesecloth layers. The spore concentration was determined and adjusted to 1 × 10 5 spores/ml using a hemocytometer.
When the primordial caps reached 0.5 cm diameter after emergence from the casing soil, approximately 50 ml of M. perniciosa Hp10 spore suspension was sprayed on the surface of the caps in each basket. Similarly, 50 ml of SDW was sprayed on the surface of the primordial caps in each basket as a negative control. After inoculation, the mushrooms were observed for changes in disease symptoms every 24 h for 20 days by randomly selecting infected fruiting bodies and observing them under a light microscope (Zhang et al., 2017a). Tissues from the important time-points ( 0-, 3-, 4-, 5-, 10-day old tissues) during M. perniciosa Hp10 infection on A. bisporus were collected, after which the samples were immediately frozen in liquid nitrogen and stored at −80 • C until further use. The test was repeated twice, with four baskets for each test. The disease assessment was recorded for only the first flush. After disease development, the pathogen was reisolated, as previously described.

Identification and Characterization of GH18 Gene Family Members of M. perniciosa Hp10
The annotated protein sequences of M. perniciosa Hp10 were used as queries for a hidden Markov model (HMM) search against the SwissProt 1 (Boutet et al., 2016), InterPro 2 (Mitchell et al., 2019), and carbohydrate-active enzymes databases (CAZy) (Lombard et al., 2013) using HMMER 3.3 3 . The retrieved sequences were searched against the SMART 4 ( Letunic and Bork, 2018) and the National Center for Biotechnology Information (NCBI) Conserved Domain Search Service tool 5 (Lu et al., 2020) to confirm the conserved domains for the GH18 gene family. The number of amino acids, theoretical molecular weight (MW), and isoelectric point (PI) of the GH18 proteins were predicted using ProtParam 6 (He et al., 2019). The subcellular localization was predicted using BUSCA 7 (Savojardo et al., 2018) web server.
distributions of M. perniciosa Hp10 GH18 genes based on the gene starting positions and chromosomal lengths.

Identification of Orthologs and Gene Duplication
OrthoFinder 2 (Emms and Kelly, 2019) was used to determine the orthologous genes and duplicated gene pairs between M. perniciosa Hp10 chitinase GH18 and the other 15 species. GH18 protein sequences of M. perniciosa Hp10 were used in reciprocal BLASTP searches with an E value cutoff of 10e-5 and coverage of ≥80% to give lists of BLAST hits and query/target midpoint positions for each chromosome. Genes on the same chromosome separated by two or more genes in a 100 kb region on a chromosome (Nei and Gojobori, 1986) were considered as tandem array genes. The ratios of non-synonymous to synonymous nucleotide substitution rates (Ka/Ks) of the duplicated genes in M. perniciosa Hp10 were calculated using Ka/Ks Calculator 2.0 (Liberles, 2001;Siltberg and Liberles, 2002).

Expression Pattern Analysis of the GH18 Gene Family of M. perniciosa Hp10
Total RNA was extracted from 100 mg of A. bisporus tissues at each time-point (0-, 3-, 4-, 5-, and 10-day old tissues) using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer's instructions. For each time point, nine tissues were randomly chosen; RNA was extracted, and then equal amounts of RNA from nine tubes were mixed into three new tubes. The purity, concentration, and integrity of RNA samples were determined using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States) and an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). The first-strand cDNA was generated using transcript one-step gDNA removal and cDNA synthesis SuperMix kit (TransGen Biotech Co., Ltd., Beijing, China) following the manufacturer's instructions. cDNA libraries were constructed using a NEBNext R Ultra TM RNA Library Prep Kit for Illumina R (NEB, United States) following the manufacturer's recommendations. The cDNA libraries were sequenced on an Illumina HiSeq X-ten platform with 150 bp paired-end reads at the Novogene Biotech Company (Beijing, China).
Two biological repeats were established for each treatment. After sequencing, low-quality reads and adapter sequences were removed from the raw data using the NGS QC Toolkit 12 (Patel and Jain, 2012). The clean data were then mapped to the M. perniciosa Hp10 and A. bisporus H97 genome sequences using TopHat 13 (Trapnell et al., 2009). The differentially expressed genes (DEGs) analysis was performed with DESeq software (version 1.18.0) (Anders and Huber, 2010). The Fragments Per Kilobase of transcript per Million mapped reads (FPKM) method (Trapnell et al., 2010) was used to obtain the expression levels and to calculate the differential expression multiples among different samples. The false discovery rate (FDR) was used to test the multiple hypotheses of the calculated results, FDR not greater than 0.001, and | log2 ratio | ≥ 1 were defined as the threshold to screen differentially expressed genes (DEGs) and were employed to obtain significantly differentially expressed genes among samples (Padj < 0.05).
All of the DEGs were functionally annotated by mapping to the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) 14 , and InterProScan databases (Ashburner et al., 2000;Jones et al., 2014;Kanehisa et al., 2015; The Gene Ontology Consortium, 2018) using BLASTX program with an E-value cutoff of 10 −5 and identity cutoff of 40%. The Pretty Heatmaps (pheatmap) package (Yao and Liu, 2018) in R software was used to draw the expression pattern of members of the GH18 gene family. These RNA-sequencing data have been submitted to the NCBI SRA database (SRP190007).

Gene Cloning and Analysis of the GH18 Genes
To validate the RNA-seq data, five differentially expressed GH18 genes at different time points (0, 3, 4, 5, and 10 days) of M. perniciosa infection of A. bisporus were analyzed by PCR and RT-PCR. Specific primers ( Table 1) were designed for the five GH18 genes using Primer Premier 5.0 software (Liu et al., 2015). Approximately 200 ng/µL of DNA and cDNA products were used as templates for the PCR and RT-PCR. The reactions were performed in 25 µL containing 2 µL template DNA, 12.5 µL of Premix Taq (TaKaRa, Da Lian, China), 1 µL (10 µM) of each primer, and 10.5 µL of RNase-free water. The reactions were performed in a Bio-Rad T100 thermal cycler (Bio-Rad Lab. Inc., Ltd., California, United States) with the following conditions: 94 • C for 5 min, followed by 30 cycles of 94 • C for 30 s and 60 • C for 30 s, then 72 • C for 30 s, and a final extension at 72 • C for 10 min. The amplified products were detected by 1% agarose gel electrophoresis, purified using an Axyprep DNA Gel Extraction Kit (Axygen Scientific, Inc., California, United States), and cloned into a pXT19-T Vector (Beijing Dingguo Changsheng Biotechnology Co., Ltd., Beijing, China) followed by transformation into E. coli DH5α. The positive clones were screened by LB liquid medium containing ampicillin and confirmed by PCR amplification and agarose gel electrophoresis. Three positive plasmid samples were sequenced using the five GH18 specific primers at Sangon Biotech Co., Ltd. (Shanghai, China), and the sequences were analyzed using DNAMAN 6.0 (Lyu et al., 2017). Two biological replicates and three technical replicates of each sample were used for both the PCR and RT-PCR reactions.

Genome-Wide Identification and Characterization of GH18 Genes in M. perniciosa Hp10
A total of 63 putative GH18 gene sequences were obtained from the genome of M. perniciosa Hp10 after an HMM search of CAZy, InterPro, and SwissProt databases. All candidate GH18 sequences were further analyzed using the CDD and SMART databases to confirm the presence of conserved domains. Forty-one putative GH18 genes were obtained after eliminating short length (100 bp) and low identity sequences ( Table 2). Based on the conserved domains, the 41 GH18 genes were divided into three groups and eight subgroups ( Figure 1A), in which 18 genes contained signal peptides at the N-terminus, and 14 genes had small domain CBM1, ChtBD1, or LysM. Sixteen genes (gene length 1,044-4,577 bp) belonged to group A that contained three subgroups (A-II, A-IV, and A-V); each group contained Glyco_hydro_18. Eight genes (gene length 1,002-1,498 bp) belonged to group B, which contains three subgroups (B-I, B-II, and B-V); each contained GH18_hevamine_XipI_class_III, GH18_CTS3_chitinase, or GH18_chitinase_D-like, and one gene (WH10003001) of subgroup B-II contained a small conserved domain CBM1 at the C-terminus. Seventeen genes (gene length 1,134-5,487 bp) belonged to group C that contained two subgroups (C-I and C-II). Each group contained GH18_zymocin_alpha or GH18_chitolectin_chitotriosidase conserved domain, and some had the ChtBD1 or LysM domain. Five genes of subgroup C-II had two concurrent domains. The gene lengths varied from 1,002 to 5,487 bp, and the coding sequences (CDS) of the genes ranged from 918 to 4,221 bp. The peptide lengths ranged from 305 (WH10006907) to 1406 (WH10002109) amino acids, corresponding to molecular weights (MW) of 33.28 kDa (WH10006907) and 150.78 kDa (WH10002109), respectively. The isoelectric points (pI) ranged from 4.21 (WH10001817) to 8.01 (WH10003001) ( Table 3).
The grand average of hydrophobicity (GRAVY), the fat index, and the instability index of the GH18 genes ranged from −0.564 (10008749) to 0.18 (WH10006436), 63.1 (WH10005062) to 93.15 (WH10006436), and 25.27 (WH10006436) to 48.14 (WH10007469), respectively. The amino acid composition and content analysis showed 13 different kinds of amino acids, and glycine and alanine were the predominant residues. The proteins had variable number disulfide bonds and protein-protein binding regions ranging from 1 to 50 and 5 to 37, respectively. All the genes were predicted to have α-helix (6.72-34.12%), β-fold (8.39-25.68%), and random coil (29.84-71.88%) structures. Most of the GH18 proteins (27, 65.85%) were predicted to be located in the extracellular space (matrix), but some proteins were located in the cytoplasm, plasma membrane, mitochondrion, endomembrane system, and organelle membranes. The amino acid peptide length, molecular weight, isoelectric point, instability index, fat index, GRAVY, and putative in silico subcellular localization predictions of GH18 identified in M. perniciosa Hp10 are listed in Table 3.

Phylogenetic Analysis of GH18 Family Proteins
To predict the functions and better understand the evolutionary relationships of GH18 chitinase proteins among different species, a phylogenetic tree was constructed using the 41 GH18 chitinase protein sequences of M. perniciosa Hp10 and 56 GH18 protein sequences of 15 other fungal species. The phylogenetic tree clustered the 41 GH18 chitinase protein sequences of M. perniciosa Hp10 into three major groups (A, B, and C) and eight different subgroups (A-II, A-IV, A-V; B-I, B-II, and B-V; C-I and C-II) according to sequence similarity of their GH18 catalytic domains (Figure 2). This result was consistent with that using the conserved domains.

Gene Structure, Conserved Motif Analyses, and Chromosomal Location
The exon-intron structure was analyzed to provide further insight into the evolution of the GH18 genes. The numbers of exons/introns in the GH18 gene family ranged from 1 to 15 and 0 to 14, respectively. Although GH18 genes with high similarity were clustered in the same group, the numbers, distribution, and locations of the exons/introns were different ( Figure 1B). Group C contained the highest numbers of exons/introns (1-15 exons and 0-14 introns), and Group B had the lowest numbers of exon and intron (1-5 exons and 0-4 introns). Fifteen distinct motifs with sizes ranging from 15 to 50 amino acids were identified among the 41 GH18 proteins ( Figure 1C). From the Pfam analysis, 11 out of the 15 proteins encoded functional domains. Except for motif 9, which encoded chitin recognition protein, the rest of the motifs (1-5, 7-8, 12, 14, and 15) with functional domains encoded glycosyl hydrolases family 18. Motif 1, followed by motifs 4, 12, and 2, were widely distributed in all of the GH18 genes (31, 30, 30, and 29 genes, respectively). Also, most GH18 proteins in the same group or subgroup shared similar motif types and distributions. For example, except for WH10002282, all the members of groups B, D, and E contained two motifs (motifs 2 and 10).

Identification of Orthologs and Gene Duplication
A total of 35 orthologous clusters were identified among the 16 fungal species. All 41 genes of M. perniciosa Hp10 formed orthologous groups with 13 fungal species (Supplementary S1). In total, 180 gene duplication events were identified within the species tree. Twelve gene duplication events were identified in M. perniciosa Hp10, and these were located on terminal branches of the species tree. Eight genes (WH10000176 and WH10000184, WH10001816 and WH10001817, WH10003001 and WH10002868, and WH10006436 and WH10006445) appeared as tandem repeats on utg195, utg16, utg19, and utg18, but there were not tandem array genes. The estimated Ka/Ks ratios of the genes in the twelve duplication events varied between 0.3842 and 3.5969. The estimated Ka/Ks ratios of four gene duplication events were above 1 (Ka/Ks > 1), whereas the other eight gene duplication events were below 1 (Ka/Ks < 1) (Supplementary S2).  Figure 4A). In total, 78.64 Gb of high-quality clean reads was generated after raw data filtering and trimming, with a Q20 quality score ≥97%. The transcriptome sequencing depths were identified as being close to saturation. The clean reads were mapped to both the pathogen M. perniciosa Hp10 and the host A. bisporus H97 genomes ( Figure 4B). The results revealed that the proportions of total clean reads mapped to the pathogen M. perniciosa Hp10 reference genome were 0.03, 1.45, 10.15, 57.92, and 96.22% at 0, 3, 4, 5, and 10 dpi, respectively. The proportions increased with the post-infection times, reflecting the increase in production of the pathogenic fungal biomass. For the mapping to the host A. bisporus H97 reference genome, the proportions were 82.48, 77.19, 71.18, 32.63, and 0.27% at 0, 3, 4, 5, and 10 dpi, respectively, indicating a decrease with the post-infection times. Therefore, we selected four infection stages (3, 4, 5, and 10) to map to the pathogen M. perniciosa Hp10 genome to perform further identification of the pathogenicity-related genes. The normalized differentially expressed genes (DEGs) were calculated by comparisons between different time points using the log2 fold change (| log2FC [fold change] | > 1, p < 0.005).
The significantly upregulated and downregulated genes were functionally annotated using the GO, KEGG, and PFAM  pathogenicity, peptidase, stress, secondary metabolites, and mating (Table 4).

Expression Pattern of GH18 Genes During M. perniciosa Hp10 Infection of A. bisporus
To understand the possible roles of M. perniciosa GH18 genes in A. bisporus, we characterized the 41 GH18 genes expressed during the disease infection period. The results showed that 23 genes out of the 41 GH18 genes were differentially expressed during the disease infection stage (Figure 5). The results showed no differential expression of GH18 genes at 3 vs. 4 dpi, while 15 GH18 genes were differentially expressed at 4 vs. 5 dpi, and only one gene (WH10009310) was significantly downregulated. The other 14 genes were significantly upregulated, including eight genes in group A (WH10000025, WH10000176, WH10000259, WH10003193, WH10006636, WH10006907, WH10008231, and WH10008749), two genes in group B (WH10002282 and WH10002868), and four genes in group C (WH10002350, WH10002829, WH10005656, and WH10006445). At 5 vs 10 dpi, 13 genes were downregulated, including nine genes in type A (10000025, WH10000176, WH10000184, WH10000259, WH10003193, WH10006636, WH10006907, WH10008231, and WH10009310) and four genes in group C (WH10002350, WH10002829, WH10006445, and WH10009911); six genes had upregulated expression levels, including three genes in type A (WH10003018, WH10003963, and WH10004667), two genes in group B (WH10004309 and WH10006436), and one gene in group C (WH10001817).

Gene Cloning and Analysis of GH18 Genes of M. perniciosa Hp10
Five representative GH18 genes from groups A, B, and C (WH10004667, WH10004309, WH10002868, WH10006436, and WH10001816) (Figure 6) were cloned and sequenced to verify the accuracy of the genome and transcriptome analyses. The results showed that each of the target gene lengths of DNA and cDNA sequences was 1,000-1,500 bp and 1,000-1,300 bp, respectively. The genes showed 100% similarity to the sequences used for the in silico analysis. In addition, the CDS length (969-1,263 bp), amino acids (322-420), and introns (0-4) of conserved domains indicated that these were GH18 genes. From the RT-PCR results, we found that these five genes were involved in the infection process, and some genes were differentially expressed; likely WH10004667, WH10004309, and WH10006436 were upregulated at 5 vs. 10 dpi, while WH10002868 was upregulated at 4 vs. 5 dpi, consistent with the transcriptome analysis.

DISCUSSION
Despite the economic importance of wet bubble disease of A. bisporus, the details of the genes expressed by M. perniciosa during infection of its host remain largely unknown. In our previous work, we sequenced the genome of M. perniciosa and identified 41 GH18 chitinase genes and putative genes related to pathogenicity (Li et al., 2019). The GH18 chitinases are diverse multigene families in a wide range of organisms and are known to play essential roles in biological processes like growth, nutrient acquisition, interspecific interactions, pathogenesis, and defense (Gruber and Seidl-Seiboth, 2012;Nagpure et al., 2013).
In this study, we performed deep genome mining and characterized a total of 41 GH18 genes in M. perniciosa Hp10. Furthermore, we analyzed the transcriptome profile of M. perniciosa Hp10 and the expression patterns of its GH18 genes during infection of A. bisporus. The 41 GH18 proteins' physicochemical characteristics were similar to those of other filamentous fungi, and the genes had detectable transcripts, thereby validating their functionality. However, the number of GH18 genes in M. perniciosa Hp10 was higher than those in F. graminearum, M. oryzae, N. crassa, R. solani, and T. virens (Xue et al., 2012;Zhang et al., 2014;Zapparata et al., 2017), which generally range from 10 to 30 genes. The high number of GH18 genes in M. perniciosa Hp10 is characteristic of mycoparasitic fungi (T. atroviride and T. virens), as it may require several chitinase isozymes acting in synergy to enable degradation of the host chitin cell wall (Lichius et al., 2014).
The pI analysis revealed that most of the GH18 chitinases are acidic enzymes, except two genes that encode alkaline chitinases (Table 3), consistent with the results of previous studies (Seidl et al., 2005). The different subcellular localization patterns of the GH18 proteins suggest that they might be differentially regulated and may have distinct roles in M. perniciosa Hp10. The GH18 proteins secreted into the extracellular space might function as virulence factors and in cell wall remodeling (Zamith-Miranda et al., 2018;Vincent et al., 2019). The finding also suggests that the proteins could be employed as drug targets to control M. perniciosa (Blum et al., 2009). The subcellular localization of proteins is invaluable for understanding their functions and interactions with other proteins (Peng and Gao, 2014).
Phylogenetic analysis revealed that the 41 M. perniciosa Hp10 GH18 genes were clustered into eight subgroups among three groups (A, B, and C), a result that is consistent with the previous classification of GH18 gene families in related genera (Seidl et al., 2005;Karlsson and Stenlid, 2008). The numbers of genes in groups A, B, and C were higher than those reported for Trichoderma species (Kubicek et al., 2011). In addition, M. perniciosa Hp10 GH18 genes in the same group and subgroup had a similar conserved domain and motif distributions to closely related members in the phylogenetic tree, revealing the functional similarity among the same subgroup proteins. Gene structure analysis of the clustered groups showed variation in the intron and exon numbers and lengths, but all had conserved motifs. In general, genes in group C had more introns than those in group A and B. A small number of introns in a gene is the result of genetic evolution and can help the gene be rapidly regulated during a stress response, and the introns also serve as a source of sequence variation (Jacob and Smith, 2017;Naro and Sette, 2017). At the same time, intron retention may also increase the complexity of this protein family (Zhang et al., 2004). Orthology analysis revealed that all the M. perniciosa Hp10 GH18 genes were homologous to the GH18 genes of the 13 species used in the phylogenetic tree. This suggests that M. perniciosa Hp10 GH18 genes and those of the other 13 species evolved from a common ancestor.
The GH18 genes were unevenly distributed on different contigs, and two gene pairs showed tandem repeats on four different contigs. Further analysis revealed the duplication of M. perniciosa Hp10 GH18 genes. We found that the majority of the GH18 genes of 13 species were arranged tandemly. This suggests that tandem duplication has been a major process in the evolution of the GH18 gene family through unequal crossover events (Gruber et al., 2011a). The duplication of genes allows the accumulation of mutations (on independent replication of a single sequence), leading to an increase of divergence and subsequent expansion of the gene family (Bowers et al., 2003;Gu et al., 2003Gu et al., , 2005. The Ka/Ks ratios revealed that M. perniciosa GH18 genes had undergone the process of both positive selection (Ka/Ks > 1) and purifying selection (Ka/Ks < 1). GH18 genes in Trichoderma species have been reported to evolve under positive, neutral, and purifying selections (Ihrmark et al., 2010). The positive selection might change the protein structure and increase the functional divergence of the enzyme, whereas the purifying selection suggests the GH18 proteins preserve the ancestral function of duplicated genes (Ihrmark et al., 2010;Kubicek et al., 2011).
The transcriptome analysis of different stages of M. perniciosainfected A. bisporus revealed that a tiny percentage of clean reads was mapped to the M. perniciosa genome at 1-3 dpi (less than 2%), indicating the low levels of biomass during this period of A. bisporus basidiome colonization (Suberkropp, 2001;Rudd et al., 2015). From the onset of symptoms 4-10 dpi, there was an increased proportion of reads mapped to the pathogen genome (i.e., from 10.12 to 96.22%), which could be attributed to the growth of the fungus completely colonizing the host and the substantial suppression or near destruction of the host RNA by the pathogen (Lai et al., 2014;Rudd et al., 2015). There were some highly upregulated DEGs (96) during 3 vs. 4 dpi; at this stage, the conidia had already attached to the host, and the rapid growth of mycelia was underway. The pathogen uses minimal amounts of host-derived nutrients for the initial phase of colonization. Therefore, there is a need for the rapid growth of the intracellular membrane and transport of components in order for the pathogen to acquire nutrients from the host environment (Xie et al., 2014). At 4 vs. 5 dpi, M. perniciosa Hp10 (Zhang et al., 2017b) overcomes the host defense system and utilizes the host as a source of nutrients, similar to the nutrition style of Hirsutella minnesotensis and Trichoderma species (Lai et al., 2014;Xie et al., 2014).
M. perniciosa Hp10 overcompensates for the host defense by producing more peptidases (such as subtilase, metallopeptidase, lipases, and peptidase family M28), and more glycosyl hydrolases, especially GH 18 and LysM domain and other chitin or carbohydrate-modifying or shielding proteins, such as CFEM (Common in several Fungal Extracellular Membrane proteins) and the WSC (Wall Stress-responsive Component) domain, which protects the pathogen, as well as genes for adaptation to environmental stress (Hsp70 protein) (Kubicek and Druzhinina, 2013;Karlsson et al., 2017). The cerato-platanin and FADbinding proteins of M. perniciosa Hp10 may be associated with hyphal growth and mycoparasitism, as observed in T. harzianum (Gomes et al., 2015;Karlsson et al., 2017). The upregulation of the velvet factor mediates the synthesis of secondary metabolites and initiates sexual reproduction and the production of spores and conidia by regulation of meiotically upregulated genes (Lind et al., 2016).
At 10 dpi, many candidate-effector proteins, protein kinases, secondary metabolites, toxins, sexual reproduction genes [HMG (high mobility group) box and Pheromone A receptor], CAZymes, and transcription factors [CHY zinc finger, Zn (2) Cys (6) and bZIP] (Morán-Diez et al., 2015;Karlsson et al., 2017;Druzhinina et al., 2018) were highly expressed. We suggest these DEGs are associated with the pathogen reproduction and proliferation in the host after the defeat of the host tissues. The high upregulation of sexual reproduction genes is necessary for the pathogen to produce a higher number of conidia in order to persist in compost and for dispersal to continue throughout the infection cycle during subsequent flushes of A. bisporus.
Gene expression analysis showed that 23 of the 41 chitinase GH18 genes of M. perniciosa Hp10 were differentially expressed at different infection stages; 20 genes were significantly upregulated and 16 genes were significantly downregulated during the infection process, indicating that chitinase GH18 genes play different roles at different infection stages. The expression of chitinase genes was significantly different at different infection stages, and the number of chitinase genes with differential expression was less at 3 vs. 4 dpi, but 14 genes were significantly upregulated at 4 vs. 5 dpi, indicating that these 14 genes play important roles in the process of M. perniciosa infecting A. bisporus. Six genes were significantly upregulated at 5 vs. 10 dpi, indicating that these genes are closely related to the pathogen's growth and development and its mass reproduction in the host. The proportions and types of chitinase genes significantly upregulated in different infection or development stages were different, indicating that the chitinase GH18 gene family has undergone functional differentiation during evolution. The genes belonging to different types may have different functions.
In order to verify the accuracy of genome sequencing and to further analyze the function of the chitinase gene family of M. perniciosa, five typical genes were cloned in this study. It was found that all five GH18 genes had α helix structures, β folding structures, and irregular curling structures. The proportion of irregular curling to amino acids of each gene was the largest, while the proportion of α helix and β folding structure to amino acids of each gene was relatively small. In addition, the lengths of these five GH18 genes are all within 1,500 bp, but 322-420 amino acids can be encoded. It is speculated that different RNA splicing methods participate in the expression of chitinase genes, and more chitinase or other proteins can be expressed in the different infection processes and in different host types, allowing the pathogen to expand its host range (Li et al., 2019).

CONCLUSION
In summary, a total of 41 GH18 genes were identified in M. perniciosa Hp 10, with variation in gene structure, protein length, and physicochemical properties. A total of 12 gene duplication events were observed in the GH18 genes and were under positive gene selection. The transcriptome analysis during the infection of M. perniciosa Hp10 infection of A. bisporus revealed that the expression of diverse genes, including those coding for CAZyme, proteases, peptidases, effectors, P450, secondary metabolites, and transcription factors, was involved in pathogenicity and the interaction of the pathogen with its fungal host. The expression patterns of 23 GH18 genes differentially expressed at different infection stages were analyzed, and RT-PCR indicated that the GH18 gene family plays an important role in the infection process. The genes identified by the transcriptome analysis will be valuable targets for functional characterization of the potential pathogenicity factors underlying the molecular mechanism of M. perniciosa, causing wet bubble disease of A. bisporus. Understanding the strategies employed by M. perniciosa to infect A. bisporus and the host's response to the pathogen can serve as the basis for developing efficient disease-management strategies to mitigate wet bubble disease.

DATA AVAILABILITY STATEMENT
The required data type (RNA-seq data) in our manuscript has been submitted to the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra/?term = SRP190007).