Genome Sequence, Comparative Genome Analysis and Expression Profiling of Chitinase GH18 Gene Family in Cordyceps javanica Bd01

1. Introduction

Entomopathogenic fungi (EPF) have a broad host range and represents one of the largest groups of entomopathogenic microorganisms. They are also among the earliest entomopathogenic microbes applied in pest control. The fungus Cordyceps javanica (Formerly known as Isaria javanica) belongs to the family Cordycipitaceae (Ascomycota: Sordariomycetes: Hypocreales) [1]. According to the existing reports, it has shown a great biocontrol potential; it can control 11 species across 10 genera of insects. Pathogenicity tests have demonstrated that this fungus exhibits strong lethal effects on various insect species across different orders. It has shown great toxicity against several insect pest of the order Hemiptera and induced 100% mortality in Diaphorina citri [2], the Bemisia tabaci [3] and Bemisia argentifolii [4], as well as aphids such as the Hyalopterus pruni and the Aphis pomi [5], and the Empoasca vitis [6]. Additionally, it is effective against the Thrips palmi (Thysanoptera) [7] and red imported fire ant Solenopsis invicta (Hymenoptera) [8]. In the case of the order Isoptera, it is also lethal to the Coptotermes gestroi [9]. Furthermore, the fungus has demonstrated good insecticidal effects on several Lepidopteron species, including the Mamestra brassicae [10], and the Lymantria dispar [11].

Entomopathogenic fungi initiates the fungal infection through insect cuticles adhesion by secreting various enzymes, mainly chitinases, proteases, and lipases. These enzymes are crucial determinants of pathogenicity to the insect hosts [12,13,14,15,16]. The chitinases of pathogenic fungi play important roles not only in spore germination, septum formation, cell division, and morphogenesis but also in host interactions [17,18,19,20]. Based on amino acid sequence similarity, chitinases are classified into two families: glycoside hydrolase 18 and 19 [21]. The GH18 chitinase gene family is widely distributed among bacteria, fungi, viruses, animals, and higher plants [22,23,24,25]. GH18 chitinase play a role in deglycosylation of glycoproteins [26,27]. Until now, 75 chitinases of entomopathogenic fungi have been characterized from M. anisopliae, B. bassiana and I. fumosorosea. The chi2 and chi3, identified from M. anisopliae, were involved in the pathogenicity of the fungus [28,29]. Overproduction of Bbchit1 found in B. bassiana enhanced the virulence of B. bassiana for aphids [30]. The Ifchit1 of I. fumosorosea was involved in fungal growth and virulence functions, while the Ifu-chit2 may play a role in the early stage of pathogenesis [31,32,33]. However, little is known about chitinases and other entomopathogenic-related genes in the fungus C. javanica, and there is no information available in the literature about the structural and functional properties of C. javanica chitinases.

In this study, we aimed to reveal the chitinase GH18 gene family evolution based on sequencing and assembling the genome of C. javanica Bd01 isolated from the larvae of Xylotrechus quadripes. Furthermore, we performed a comparative genome analysis of C. javanica Bd01 genome with other available genomes of entomopathogens. Moreover, the phylogenetic relationships, structural domains, motifs, and protein models of the C. javanica Bd01 GH18 chitinase were analyzed. This would provide additional information on the genomic background of entomopathogenic fungi, facilitating a deeper understanding of their biological functions. These findings would serve as a reference for future protein manipulation and can be used for further studies on their mechanisms of action.

2. Results

2.1. Genome Sequencing, Assembly and Annotation

Complete genome sequencing of C. javanica Bd01 using Novaseq 6000 platforms (Illumina, USA) revealed a total number of 78.1 million clean reads (~229.6 × coverage). The assembled genome represents 98.28% completeness (based on BUSCO analysis) with a contig N50 length of 5.4 Mb. The resulting genome assembly was based on 8 contigs with a 53.18% of GC content, 9590 protein-coding genes, 43 rRNA genes, and 131 tRNA genes. The repetitive sequence 0.34Mb, represents 1.04% of the total genome (Figure 1; Table 1).

Pfam domains were assigned to 9,590 proteins based on InterProscan program. Subsequently, a total of 5,236 proteins were assigned to the eukaryotic orthologous groups (KOG) databases (Supplementary Table S1; Figure 2A). The abundance of metabolism was about 50%, which was the highest in these four categories. 3111 genes (32.7%) were annotated using KEGG map, and KEGG annotations contained four major pathways including metabolism (41.2%), genetic information processing (24.5%), cellular processes (8.2%), environmental information processing (1%) (Supplementary Table S2; Figure 2B). Gene ontology (GO) terms were divided into three major function categories: molecular functions (51.1%), cellular components (30.6%) and biological processes (18.3%). The 6,990 genes (73.7%) were annotated and assigned to the three categories, 35 genes were classified into the secondary metabolite biosynthetic process (Supplementary Table S3; Figure 2C), and 1096 genes were unknown function in KOG, GO, and KEGG.

2.2. Comparative Genomic Analysis

To evaluate the genome evolution across C. javanica Bd01 and its related entomopathogenic fungal strains and to identify factors underlying lineage diversification, we evaluated changes associated with fungal evolution in relation to Expansion and contraction of gene families, as well as gene duplications and deletions that occurred within the gene families. By OrthoFinder analysis, among the ten fungal genomes (C. javanica Bd01, C. javanica IJ1, C. javanica IJ2, C. fumosorosea ARSEF2679, C. militaris CM01, C. militaris ATCC34164, Akanthomyces lecanii RCEF1005, B. bassiana ARSEF2860, M. libera RCEF2490 and M. anisopliae JEF290), a total of 11,781 orthologous groups were constructed, and 42.68% (5028) of which were shared among the ten strains (Figure 3A, 3B), indicating their close evolutionary relationships.

Then, based on the amino acid sequences of 4,096 single-copy orthologous groups, we obtained the maximum likelihood phylogenetic tree of the ten fungal strains. It was shown that C. fumosorosea and C. militaris formed a branch, independent from the C. javanica branch (Figure 3C; Supplementary Table S4), indicating that the three species have a closer relationship with each other than with other species. The number of gene family expansion and contraction was mapped to the phylogenetic tree (Figure 3C). 74 and 413 gene families were identified as expansion and contraction in the fungus C. javanica Bd01, respectively (Figure 3C). To investigate the function of gene family expansion in C. javanica Bd01 during the evolutionary process, the present study was conducted to annotate the KOG function of the gene family expanded at the differentiation node of C. javanica Bd01 and other related species of the genus Cordyceps. Based on KOG annotation (Figure 3D), in addition to the unknown functions, the gene families for expansion in C. javanica Bd01 had the highest proportion of gene families annotated to the ABC transporter. Using CAZy database models, the number of genes for glycoside hydrolases in the C. javanica Bd01 genome was markedly greater than that in other entomopathogenic fungal genomes (Supplementary Table S5), especially in family GH18 and family CBM50.

2.3. Genome-Wide Identification and Characterization of GH18 Genes in Cordyceps javanica Bd01

A total of 27 putative GH18 gene sequences were obtained from the genome of C. javanica Bd01 after an HMM 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. Twenty-seven putative GH18 genes were obtained after eliminating short length (100 bp) and low identity sequences (Supplementary Table S6).

2.4. Conserved Motif and Chromosomal Location of GH18 Genes in Cordyceps javanica Bd01

In addition to the predictions of subcellular localizations and the analysis of trans-membrane domains, the evaluation of conservative motifs is one of the important means for the functional prediction of GH18 proteins. In this study, a sum of 15 conserved motifs (motifs 1–15) were predicted using MEME online software (Figure 5A). Several motifs were found in most GH18 proteins, such as motifs 1, 2, 3, 7, and 10. Additionally, Cordyceps0G049700.1 contained motifs 1–15 and motif 3 existed in all GH18 members, whereas Cordyceps0G070700.1 only included one member (Figure 5A).

Most of GH18 proteins within the same clan exhibited similar motif components while a high discrepancy was observed among different clans, indicating that the GH18 members within the same family may undertake semblable functions, and that some motifs may play a vital role in the family-specific functions. The chromosomal location of each GH18 gene in the C. javanica Bd01 genome is shown in Figure 5B. A total of 27 GH18 are unevenly distributed on 7 contigs. Out of the 27 GH18, seven were mapped onto contig 1, followed by 6 on contig 4, whereas only 1 GH18 were located on contig 3. The GH18 with relatively high densities were detected on the top and middle arms of the contig.

2.5. 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 27 GH18 chitinase protein sequences of C. javanica Bd01 and 63 GH18 protein sequences of 4 other fungal species (C. sinensis, C. militaris, C. cicadae, and B. bassiana). The phylogenetic tree clustered the 27 GH18 chitinase protein sequences of C. javanica Bd01 into three major groups (A, B, and C) (Figure 5). This result was consistent with that using the conserved domains (Figure 4A).

2.6. Characterisation of Physicochemical Properties

The predicted physicochemical properties as computed by the ProtParam tool showed that the chitinases under investigation from C. javanica Bd01 had a varied number of amino acids (315 - 1842) and molecular weights (14.03 kDa – 81.41 kDa) (Table 2). The theoretical pI of all the chitinases ranged from acidic pH (4.40) to slightly alkaline pH (7.92). Analyses of the sequences showed that the number of negatively charged amino acid residues (aspartate and glutamate) was higher than the number of positive residues (arginine and lysine) in all the strains, except Cordyceps0G055830.1 in C. javanica Bd01 (Table 2). The instability index, the measure of the in vivo half-life of a protein, of most of the analysed chitinases were under 40. GRAVY values for all the chitinase were predicted to be negative (Table 2). This variable is used for the prediction of the hydrophobicity or hydrophilicity of a protein over the entire amino acid sequence. The aliphatic index of the chitinases in this study ranged between 59.48 to 89.41. Aliphatic index of 30% chitinases was above 75, while another 40% chitinases were between 65-75 and the rest 30% was around 64. Due to the wide variation in the aliphatic index of chitinases from the same source, C. javanica Bd01, it is predicted that the enzymes had varying temperature optima and thermostability. The amino acid composition of chosen chitinases from C. javanica Bd01 was evaluated and profiled as depicted in Figure 6A. The results substantiated that alanine and glycine were the most widely distributed amino acids in the selected chitinases. These protein sequences were also rich with amino acids such as leucine, serine, threonine, aspartic acid, valine, proline and isoleucine. Secondary structures of proteins are governed by the inclination of amino acids to form either helix or sheet form.

2.7. Secondary and Tertiary Structure Analysis

The secondary structures of all the analysed chitinases are dominated by a random coil region as this conformation has the highest mean value of approximately 58.3% (Figure 6B). This was followed by the alpha helix region (27%) and Extended strand (14.7%).

The tertiary structure of proteins is critical in distinguishing the function of proteins, their interactions with other compounds such as ligands, other proteins, and nucleotides as well as to understand the phenotypical effects of mutations. The Swiss-Model server was used to predict the three-dimensional structure of C. javanica Bd01 GH18 chitinases based on the known crystal structures of homologous proteins. For the 3D structure prediction, Most of the GH18 protein sequences have very similar three-dimensional structures, and their structural models are shown in Figure 7, with a higher degree of similarity among members within subfamilies than among subfamilies.

2.8. Analysis of cis-Regulatory Elements in the GH18 Gene Family of Cordyceps javanica Bd01

To further investigate the cis-acting elements in the promoter region of the GH18 gene, the 2000 bp region upstream of the transcription start site of the GH18 gene family was analyzed using the PlantCARE online tool. In terms of element types, hormone responsiveness was the most diverse and was found in all GH18 gene promoter regions. Regarding the distribution of individual elements, 90% of the family members contained CGTCA motif, G-BOX and ABRE elements in the promoter region. It is suggested that these genes may bind to various transcription factors and be involved in the regulation of pathogenicity, and some of the elements in the promoter region of the GH18 gene are shown in the Figure 8.

3. Discussion

Complete genome sequencing of C. javanica Bd01 revealed a total number of 78.1 million clean reads (~229.6 × coverage) and BUSCO analysis represents 98.28% completeness with a contig N50 length of 5.4 Mb. The resulting genome assembly was based on 8 contigs with a 53.18% of GC content, 9,590 protein-coding genes, 43 rRNA genes, and 131 tRNA genes. The repetitive sequence 0.34Mb, which represents 1.04% of the total genome. These findings are in line with Yu et al. [34]; their findings revealed a genome size of 47,239,278 bp, comprising 27 contigs, with a GC content of 51.16% from sequencing of Fusarium solani KMZW-1. The genome completeness was reported as 97.93% (BUSCO analysis). Moreover, the DFVF sequence identifier was Fusarium 0G092560.1, and AntiSMASH analysis identified 35 gene clusters associated with secondary metabolite biosynthesis. Similarly, the genome of C. cicadae strain CCAD02 (a sexual type of I. cicadae) was also sequenced using asexual fruiting bodies [35]. Both of them have close genome size (33.8~33.9 M).

Delineating gene family emergence and extinction within phylogenetically related organisms can identify molecular determinants that underlie species adaptation and evolution [36]. By the comparative genomic analysis, we found that the Chitinases (GH18) is significantly expanded in genes in the genomes of C. javanica Bd01. 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 defence [37,38]. Interestingly, chitinases from B. bassiana and other entomopathogenic fungi have been identified as important virulence factors in their pathogenicity towards arthropods and nematodes [39]. Many fungi feed on chitin or chitin-containing organisms. Subgroup B chitinases appear mainly to be involved in nutritional functions and more aggressive functions such as invasion and pathogenesis [40]. Most studied subgroup B chitinases in mycoparasitic and entomopathogenic fungi are inducible by nutritional stimuli including chitin or host-specific carbon sources. The expression of subgroup B chitinases, such as Chit33 in Trichoderma harzianum, Ech30 in T. atroviride and Chi2 in M. anisopliae was upregulated during starvation but repressed by glucose or other easily metabolizable carbon sources [40]. Subgroup B chitinases also play a role in insect infection. Knockout of Chi2 from M. anisopliae resulted in decreased virulence toward insects, while Chi2 overexpression constructs showed higher efficiency in host killing [41]. The other two subgroups (subgroup A and C) of fungal chitinases are also related to pathogenesis [42,43,44,45]. The compartmentation and localisation of different proteins, including enzymes have been well linked to their general biological functions and the reactions they catalyse specifically. In C. cicadae genome, Lu and colleagues identified 135 carbohydrate-degrading enzymes, including 16 GH18 proteins [35]. This observation is in accordance with the previous reports that entomopathogenic fungi had more GH18 chitinases than plant pathogenic fungi and mammalian pathogenic fungi [46]. We identified 27 GH18 chitinase genes and putative genes related to the pathogenicity. Furthermore, it is believed that information on the subcellular localisation may provide useful insights into the specific enzymatic pathways of proteins and serve as a guide in subsequent wet-lab experiments.

Thus, in the present study, a variety of readily available computational tools were employed to gain insights into the overall physical parameters; primary, secondary and tertiary structures; functional analysis, domains and motifs, and protein model of C. javanica Bd01 chitinase. However, the function of GH18 in the pathogenesis of C. javanica Bd01 is unknown and needs to be further investigated, and we are only making predictions based on the existing conditions.

The molecular weight of the enzyme is a critical characteristic for the selection of appropriate purification protocol. The molecular weight range observed in this study suggested the use of gel exclusion chromatography with Sephadex G-100 SF having a fractionation range of 2-120 kDa, best suited for the purification of chitinases under investigation. Orthology analysis revealed that all the C. javanica Bd01 GH18 genes were homologous to the GH18 genes of the 4 species used in the phylogenetic tree. This suggests that C. javanica Bd01 GH18 genes and those of the other 4 species evolved from a common ancestor. This indicates the divergence in the relative amount of basic/acidic amino acids in their sequences. Hence isoelectric focusing mediated purification should be tailor-made for the individual chitinases as they are expected to precipitate in buffers with different pH values. However, the results showed that most of the C. javanica Bd01 GH18 chitinases contain more acidic amino acid residues than basic residues. Results of pI analysis showed that most of the chitinases are acidic enzymes, which is in agreement with previous studies [47,48].This characteristic can also be used for the purification of the said chitinases by employing anion exchange chromatography with weak and strong anion exchangers such as DEAE-sepharose and Qsepharose. These results are in correlation with the earlier reports where chitinases from B. bassiana has been purified with DEAE-cellulose and Mono Q column [49,50]. Recent findings have suggested the effects of nascent peptide charge on protein translation efficiency and protein expression [51]. Hence, it can be speculated that most of the chitinases in this study, being more acidic will be polysome-translated proteins. Proteins with instability index below 40 were predicted to possess in vivo half-life of above 16 hours, whereas, those with instability index above 40 have in vivo half-life of below 5 hours [52]. Therefore, it could be inferred that besides ten chitinases, all the other C. javanica Bd01 chitinases analysed, possess remarkable stability in vivo. It is posited that the ten proteins might be products of gene duplications in the organism as strain possess other chitinases which were predicted to be relatively stable [53]. The results from this prediction could specifically indicate the significant thermodynamic and kinetic stabilities of the chitinases which are believed to be necessary for their natural function as entomopathogenic enzymes as well for their potential applications in industrial processes [54]. Alanine, leucine, aspartic acid and proline have a high propensity to form helix confirmation while valine, isoleucine, threonine, glycine, and serine favours the formation of sheet confirmation [53]. The dominance of random coils and alpha-helix regions among the secondary structures indicate the stability of the enzymes as well as a significant amount of conservation within chitinases GH18 from C. javanica Bd01 [55].

In order to verify the accuracy of genome sequencing and to further analyse the function of the GH18 gene family of C. javanica Bd01, 27 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. 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.

In filamentous fungus, only the TATA box, CCAAT motifs, and CT rich sections of core promoter elements have been functionally studied. Other elements, such as Inr, DPE, MTE, BRE, and CpG islands, mostly present in higher eukaryotes, are not characterized in filamentous fungi. Further, multiple copies of a cis-acting region were added to Aspergillus to enhance its promoter activity [56]. The majority of fungal promoters, viz., gpdA, oliC, trpC, citA, and agaA, include a region rich in pyrimidines (CT) upstream of the TSP. The presence of these CT-rich sequences was more evident in highly expressed genes (especially housekeeping 400 Page 6 of 24 World Journal of Microbiology and Biotechnology (2024) 40:400 genes) and in genes missing the TATAAA and CCAAT motifs, indicating a potential role for these sequences as promoter elements in filamentous fungi [57]. G-boxes are a class of cis-acting elements, including the ACGT family, that are capable of responding to a variety of environmental cues, including ABA, light, UV, injury, and pathogen signaling [58,59].

In summary, a total of 27 GH18 genes were identified in C. javanica Bd01, with variation in gene structure, protein length, and physicochemical properties. The lower GRAVY score and higher aliphatic index substantiated the hydrophilicity and thermostability of the chitinases. Collectively from the phylogeny and structure analysis, it was observed that the twenty-seven sequences are highly conserved and might also share a similar mechanism for chitinase activity. approach afforded some useful information about the chitinases from C. javanica, which could be helpful for isolation and characterisation of the enzyme in vitro and to understand the possible structure and functions of unknown proteins and hence could be subsequently exploited for various applications.

4. Materials and Methods

4.1. Culture of C. javanica Bd01

The strain C. javanica Bd01 was isolated from the dead larvae of Xylotrechus quadripes and the fungal isolate belongs to the collection of the Institute of Microbiology, Chinese Academy of Sciences, Conservation No. CGMCC23078. The isolate was identified morphologically and by means of the nucleotide sequence of the ITS1-5.8S rRNA-ITS4 region. The internal transcribed spacer sequence of C. javanica Bd01 was deposited in the GenBank nucleotide sequence database under the accession number MZ831846 [60]. Routinely maintained on sabouraud dextrose agar medium with yeast extract (SDAY) [61]. Bd01 was grown on SDAY at 27°C temperature, 70% ± 5 relative humidity and 12:12 (dark: light) photoperiod for 14 days.

4.2. DNA Extraction and Sequencing

DNA extraction was carried out on pure cultures of C. javanica Bd01 grown on SDAY. The experimental procedure was performed according to the standard protocol provided by Oxford Nanopore Technologies (ONT, Shanghai, China), including sample quality testing, library construction, library quality testing, and library sequencing.

Genomic DNA was extracted using the Omega Fungal DNA Kit D3390-02 according to the manufacturer’s instructions (Illumina, USA). Purified genomic DNA was quantified by TBS-380 fluorometer (Turner BioSystems Inc., Sunnyvale, CA). High quality DNA (OD260/280=1.8~2.0, > 15ug) was used to do further research.

Genome was sequenced using a combination of PacBio Sequel Single Molecule Real Time (SMRT) and Illumina sequencing platforms. The Illumina data was used to evaluate the complexity of the genome. For Illumina sequencing, at least 5μg genomic DNA was used for each strain in the sequencing library construction. DNA samples were sheared into 400-500 bp fragments using a Covaris M220 Focused Acoustic Shearer following the manufacture’s protocol. Illumina sequencing libraries were prepared from the sheared fragments using the NEXTflex^TM Rapid DNA-Seq Kit. Furthermore, in sum 5’ prime ends were first end-repaired and phosphorylated. After that, the 3’ ends were A-tailed and ligated to sequencing adapters. The third step was to enrich the adapters-ligated products using PCR (Polymerase chain reaction). The prepared libraries were then used for paired-end Illumina sequencing (2 × 150 bp) on an Illumina HiSeq X Ten machine. For Pacific Biosciences sequencing, an aliquot of 8 μg DNA was spun in a Covaris g-TUBE (Covaris, MA) at 6,000 RPM for 60 seconds using an Eppendorf 5424 centrifuge (Eppendorf, NY). DNA fragments were then purified, end-repaired and ligated with SMRTbell sequencing adapters following manufacturer’s recommendations (Pacific Biosciences, CA). Resulting sequencing library were purified three times using 0.45 x volumes of Agencourt AMPure XP beads (Beckman Coulter Genomics, MA) following the manufacturer’s recommendations. Next, a ~10kb insert library was prepared and sequenced on one SMRT cell using standard methods.

4.3. Genome Assembly and Annotation

The genome sequence was assembled using both the PacBio reads and Illumina reads. The original image data was then transferred into sequence data via base calling, which is defined as raw data or raw reads and saved as FASTQ files. A statistic of quality information was applied for the quality trimming, by which the low-quality data can be removed to form the clean data. The reads then were assembled into contigs using CANU. Finally error correction of the PacBio assembly results was performed using the Illumina reads.

The predicted protein were blast (e-value: 1e-5) against Nr [62], Swiss-Prot [63], TrEMBL [63], KEGG [64], KOG [65]. Blast2go [66] was used for GO [67] annotation. Hmmer [68] was used for Pfam [69] annotation.

4.4. Orthologous and Phylogenomic Analysis

We used OrthoFinder version 2.5.5 [70] to analyze gene families of ten entomopathogens, including C. javanica Bd01, C. javanica IJ1, C. javanica IJ2, Cordyceps fumosorosea ARSEF2679, Cordyceps militaris CM01, C. militaris ATCC34164, Akanthomyces lecanii RCEF1005, Beauveria bassiana ARSEF2860, Moelleriella libera RCEF2490 and Metarhizium anisopliae JEF290. For each of the 4,096 single-copy gene families, alignment of amino acid sequences from six strains was performed by MUSCLE version 5.0 [71]. Sequences from all single-copy gene families were concatenated by an in-house Perl script. Then, the poor alignment regions were removed using Gblocks version 0.91b [72]. The phylogenetic trees were constructed based on the alignment using the Neighbor-Joining method (1,000 repeats) with the parameters of the Jones-Taylor-Thornton model, uniform rates among sites, and partial deletion of gaps was done in MEGA version 11.0.13 [73].

4.5. Identification and Characterization of GH18 Gene Family in C. javanica Bd01

The annotated protein sequences of C. javanica Bd01 were used as queries for a hidden Markov model (HMM) searched against the SwissProt [63], InterPro [74], and carbohydrate-active enzymes databases (CAZy) [75] using HMMER database. The retrieved sequences were searched against the SMART [76] and the National Center for Biotechnology Information (NCBI) Conserved Domain Search Service tool [77] to confirm the conserved domains for the GH18 gene family. The number of amino acids, theoretical molecular weight (MW), instability index, aliphatic index, grand average of hydropathicity (GRAVY) and isoelectric point (PI) of the GH18 proteins were predicted using ExPasy-ProtParam tool (http://web.expasy.org/protparam/) [78]. The subcellular localization was predicted using WoLF PSORT (https://wolfpsort.hgc.jp/) [79] web server.

4.6. Conserved domain, Conserved Motif Analyses, and Chromosomal Location

The conserved domain structures were collected by aligning protein domain and protein family models constructed as multiple sequence using the NCBI Batch CD-Search [80] Server. The conserved motifs of the genes were predicted using MEME version 5.5.7 [81] with default parameters. The secondary structure and tertiary structure of the chitinase GH18 gene family were predicted using SOPMA version 2.16.0 [82]. server and SWISS-MODEL version 1.0.2 [83] workspace, respectively.