Bioinformatic comparison of Kunitz protease inhibitors in Echinococcus granulosus sensu stricto and E. multilocularis and the genes expressed in different developmental stages of E. granulosus s.s.

Background Cystic and alveolar echinococcosis caused by the tapeworms Echinococcus granulosus sensu stricto (s.s.) and E. multilocularis, respectively, are important zoonotic diseases. Protease inhibitors are crucial for the survival of both Echinococcus spp. Kunitz-type inhibitors play a regulatory role in the control of protease activity. In this study,we identified Kunitz-type domain protease inhibitors(KDPIs) present in the genomes of these two tapeworms and analyzed the gene sequences using computational, structural bioinformatics and phylogenetic approaches to evaluate the evolutionary relationships of these genes. Hi-seq transcriptome analysis showed that E. granulosus s.s. KDPIs were differentially expressed in the different developmental stages. We validated some of the genes expressed in adult worm, protoscolex and cyst germinal membrane of E. granulosus s.s. and E. multilocularis by quantitative PCR. Results A total of 19 genes from E. multilocularis and 23 genes from E. granulosus s.s. were predicted to be KDPIs with the most containing a single Kunitz-domain. A maximum likelihood method phylogenetic tree indicated that the E. granulosus s.s. and E. multilocularis Kunitz domain peptides were divided into three branches containing 9 clusters. The ratio of positively charged residues and neutral residues are different between E. multilocularis and E. granulosus s.s. KDPIs. We also found that E. multilocularis had higher percentage of sequences containing signal peptides (17/19, 89.47%) than that of E. granulosus s.s. (14/23, 60.87%). Transcript analysis showed all the E. granulosus s.s. KDPI genes were expressed differentially in four developmental stages of the worm. Transcription analysis showed that 9 KDPIs (including EG_07244,EGR_08716 and EGR_10096) were highly upregulated in adult worm, and 2 KDPIs (EG_09268 and EG_09490) were highly expressed in the cyst germinal membrane. Quantitative gene expression analysis(qPCR) of four genes confirmed the expression of these genes. EGR_08716 and its homologous gene (EmuJ_001137000) were highly and specifically expressed in adult worms of the two worms. Conclusions A total 19 and 23 KDPIs were identified in the genomes of E. multilocularis and E. granulosus s.s. , respectively. The differential expression of these KDPIs in different stages may indicate their different roles in the different hosts. The difference in characterization of KDPIs may be associated with the different pathology of metacestode stage of these two parasites. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08219-4.


Background
Cyst echinococcosis (CE) and alveolar echinococcosis (AE) are both medically and economically important diseases caused by the metacestode stages of Echinococcus granulosus sensu stricto (s.s.)and E. multilocularis respectively. The diseases impact on hundreds of millions of people in Asia, Europe, American and Africa [1]. The control and treatment of echinococcosis are difficult. High frequency of dosing dogs with the drug praziquantel has played a key role in the control of the disease [2,3], but undertaking the control measure is challenging in remote areas. A vaccine against adult worms in dogs is urgently needed [4].
The life-cycle of these two tapeworms involves four developmental stages including adult worm, oncosphere, cyst and protoscolex present in their definitive and intermediate hosts. The major definitive and intermediate hosts of E. granulosus s.s. are dogs and sheep respectively, whereas, the natural definitive and intermediate hosts of E. multilocularis are fox/wolf and rodent small mammals respectively. Humans are the intermediate hosts of these two tapeworms. The survival of these tapeworms relies on evading host immune responses and avoiding attack by proteases; this is especially important for the adult parasites which reside in the gastrointestinal duct, a location where high concentration of proteases are present which are harmful and toxic for the worms.
Eukaryote proteases including serine (trypsin/chymotrypsin-like), cysteine (thiol) and aspartic (pepsin/ cathepsin/rennin) proteases play a fundamental role in the regulation of protein function. Their functions and activities are controlled largely by protease inhibitors which play crucial roles in the regulation of proteases involved in a range of biological processes system including cell proliferation, inflammation, cell homeostasis and immune mechanism [5][6][7]; protease inhibitors act mainly through the control of potentially disadvantageous, excessive or inopportune proteolytic activity. Protease inhibitors including aspartic, cysteine, metallo, serine, and threonine inhibitors are super-families based on their similarities at the amino acid sequence level and tertiary structure [8]. Similarities in primary structure and tertiary structure support the common ancestry of many protease inhibitor families.
In the present study, we identified all KDPI sequences predicted in the E. granulosus s.s. and E. multilocularis genomes and used computerized programs to characterize these Kunitz domain protease inhibitors. We show that the majority of E. granulosus s.s. KDPIs were differentially expressed in different life cycle stages and some have a range of GO numbers indicating these inhibitors likely function in different ways in the tapeworm's development.

General characterizations of Kunitz domain protease inhibitors
InterproScan and Motif scan identified 19 and 23 genes encoding KDPIs from the E. multilocularis and E. granulosus s.s. genomes, respectively ( Table 1). The KDPI family has a typical Kunitz domain of about 60 amino acids in four genes confirmed the expression of these genes. EGR_08716 and its homologous gene (EmuJ_001137000) were highly and specifically expressed in adult worms of the two worms.

Conclusions:
A total 19 and 23 KDPIs were identified in the genomes of E. multilocularis and E. granulosus s.s. , respectively. The differential expression of these KDPIs in different stages may indicate their different roles in the different hosts. The difference in characterization of KDPIs may be associated with the different pathology of metacestode stage of these two parasites.
Keywords: Kunitz inhibitors, E. multilocularis, E. granulosus sensu stricto, Bioinformatic analysis, Stage expression size ( Fig. 1) with a special secondary structure formed by three disulphide bonds or bridges (Additional file 3: Fig.  S1). The echinococcal Kunitz domains contain an average of 52.85 aa (range 47-55 aa) with the majority comprising 53 aa (Fig. 1 Table S1). The majority of single Kunitz domain proteins comprise less than 100 amino acid (Additional file 1: Table S1).
We used the instability index to estimate the stability of the KDPIs. An instability value >40 is an unstable protein.
The index value representing rigidity/flexibility of each peptide varied (26.51-94.41 Table S1).
Sequence analysis showed that these echinococcal KDPI sequences contain a high percentage of hydrophobic residues including alanine (A), valine (V), leucine (L) and isoleucine (I En-t-In(T/C) 9/2 9/4 (See figure on next page.)

Cluster and phylogenetic analysis of Kunitz protease inhibitors
Multiple sequence alignment (Fig. 1) and phylogenetic analysis (Fig. 2) of the amino acid sequences were used to infer the evolutionary relationships between the E. multilocularis and E. granulosus s.s. KDPIs and to make a comparison with other species. Figure 2 shows the different evolutionary distances of the genes containing single Kunitz domain of the KDPIs using the maximum likelihood method. The maximum likelihood tree indicated that these echinococcal KDPIs were divided into three branches with 9 clusters. One branch contains several Echinococcus KDPIs such as EGR_3480, EGR_3481, EGR_05316 and EGR_07242 are close to the KDPIs from human, sheep, cattle and other species (Fig. 2 and Additional file 6: Fig. S4). Whereas other two branches are relatively Echinococcus unique (Fig. 2).
It showed that two KDPIs of EGR_05482 and EmuJ_000077700.1 from E. granulosus s.s. and E. multilocularis shared a high homology with mamba venom toxin P00982.1 and P00981 from Dendroaspis angusticeps and Dendroaspis polylepis respectively (Fig. 2), indicating these echinococcal KDPIs may have similar pathophysiological functions of mamba venom blocking ion channels and membrane receptors.

Comparison of KDPI genes predicated from the E. granulosus s.s. and E. multilocularis genomes
We compared the KDPI genes predicted from the genomes of E. granulosus s.s. and E. multilocularis and found that some genes are species-specific. E. multilocularis does not have homologues of E. granulosus s.s. sequences EG_07242, EG_07266, EG_07243, EG_09006 and EG_09008; whereas EmuJ_001136700.1 and EmuJ_001137100.1 are specific to E. multilocularis (Additional file 5: Fig. S3).
The specificity of a protease inhibitor against a protease is mainly determined by the nature of the amino acid residue at position P1 of its active site. It has been shown that Lys(K) and Arg(R) mutants of bovine pancreatic trypsin inhibitor (BPTI) bind to bovine trypsin about 10 5 -fold stronger than BPTI with P1 Tyr(T) [16]. In addition, it has been shown that typical trypsin inhibitors have Arg(R) or Lys(K) at P1, and chymotrypsin inhibitors have Leu (L) or Met (M) at the P1 position [17]. Therefore, the sequence analysis shows that the E. multilocularis has 10 KDPIs containing R at P1 and E. granulosus s.s. has 11 KDPIs containing R at P1, which belong to typical trypsin inhibitors. Furthermore, the two tapeworms have 3 or 4 sequences containing L at P1 respectively, which are chymotrypsin inhibitors ( Fig. 1; Table 1).

D and 3 D structure of Kunitz domain protease inhibitors
The majority of E. multilocularis and E. granulosus s.s. single KDPIs are small proteins sized 16-kDa and contain a relatively high percentage of Lys and Arg residues at the C-terminus. Like most Kunitz domain protease inhibitors, the Em-and Eg-KDPIs contain a conserved Kunitz type sequence with 6 cysteine residues forming three disulphide bridges (C1-C6, C2-C4 and C3-C5) (Additional file 3: Fig. S1) and these bridges play a key role in the formation of the 2D and 3D structure of these KDPIs. For the single Kunitz domain sequences, the secondary structure prediction revealed 19.01-52.71% and 18.6-60.35% of α-helix and random coil structures in Eg-KDPIs, followed by extended strands and β-turn structure, accounting for 13.1-26.67% and 1.89-10.84% respectively. Em-KDPIs α-helix and random coil structures account for 19-40.45% and 32.5-55.99% of the protein sequence respectively, followed by extended strands and β-turn, accounting for 8-36.25% and 0-10.71% (Table 3).
It is accepted that there is a close relationship between the structure and function of a protein. Therefore, we used SWISS-MODEL to predict 3D structures based on the homology modeling of KDPI templates from PDB (protein database) including single and multiple Kunitz domain proteins (Fig. 3 and Additional file 4: Fig. S2).
Three-dimensional structure analysis showed that a single Kunitz domain sequence with three disulphide bridges has a similar structure containing an α-helix and random coils with similar structures (Fig. 3). The structure of some single Kunitz domain sequences lose the second cysteine (C2) which may impact the 3D structure of these KDPIs. (Fig. 1; Table 1).

Expression of E. granulosus s.s. KDPIs in different developmental stages
To estimate expression of the KDPIs, Hi-seq techniques were employed to obtain the transcript reads of these genes from total RNA from each of 4 developmental stages of E. granulosus s.s. . The transcript read information was published in a previous paper of ours [18]. The transcriptome analysis showed that these Kunitz peptides were differentially expressed in the different developmental stages of E. granulosus s.s. ( Table 2). All the inhibitors, except EG_09006, were expressed in one or 4 stages of E. granulosus s.s. with some being highly and differentially expressed in one or two stages. Transcription analysis showed that 9 KDPIs including EGR_03480, EGR_03481, EGR_07242, EGR_07243, EG_07244,EGR_08716, EGR_08720, EGR_09269 and EGR_10096 were highly up-regulated in adult worm, and two KDPIs (EG_09268 and EG_09490) were highly expressed in the cyst germinal membrane. Sequence analysis showed that some of these adult worm up-regulated genes are extracellular including EG_03480, EG_07242, EG_08716, EG_08720, EG_09490 and EG_10096, and some are intracellular such as EG_03481, EG_07243 and EG_07244 (Table 2 and Additional file 1: Table S1).
EG_08716 is an extracellular protease inhibitor and has 42 predicted GOs, including cytoplasmic vesicle for neuromuscular process controlling balance, ionotropic glutamate receptor signaling pathway, regulation of the activity of epidermal growth factor receptor and synapse, regulation of mitotic cell cycle and translation and cellular copper and calcium ion homeostasis (Additional file 1: Table S1). The expression analysis indicated that this gene may play an important role in adult worm development and against host protease attack. It is interesting that EG_07244 is predicted having endopeptidase activity, indicating that the protein has two functions, as a peptidase and as a protease inhibitor in adult worms. This needs to be identified. EG_08721 is an extracellular inhibitor and was differentially highly expressed in the oncosphere compared with the other stages, indicating this protease inhibitor plays an important role in oncosphere biology, the only stage for primarily infection and EG_08721 may play an important role in oncosphere against host protease attack which may be a candidate for vaccine development.
Although we activated PSC with pepsin, only three KDPIs (EG_01779, EG_05317 and EG_07944) were slightly elevated in this stage. Importantly, we found that EG_09268 and EG_09490 were highly expressed in the cyst germinal membrane and the proteins expressed by these genes may be potential targets for drug development.
To validate the expression of these genes, a qPCR assay was employed to quantify the transcript levels of 4 genes of E. granulosus s.s. (EGR_08716, EGR_07244, EGR_07944 and EGR_09490), which confirmed the expression in cyst germinal membranes, protoscoleces and adult worm stages of the tapeworm. We also identified the homologous KDPI genes expression from E. multilocularis (EmuJ_001137000, EmuJ_000534800, EmuJ_000302900, and EmuJ_000929500) in three stages of the tapeworm (Fig. 4). Only one gene, EmuJ_001137000 showed a similar expression pattern as E. granulosus s.s. EG_08716, is highly and specifically expressed in the adult worms of E. multilocularis (Fig. 4).

Table 3 The secondary structure prediction of the single Eg-KDPIs and Em-KDPIs
Note: aa/%, number/percentage of amino acids in each secondary structure Other three genes were different from those homologous genes of E. granulosus s.s. , which were highly expressed in protoscoleces of this tapeworm (Fig. 4). The qPCR showed that E. multilocularis adult worm EmuJ_001137000 was 4.43 folds higher than EGR_08716 in the adult worms of E. granulosus s.s. .

Discussion
KDPIs occur in almost all living organisms from bacteria to plants and animals. Kunitz peptides show diverse biological activities including inhibition of proteases and/or blocking or modulating ion channels. Helminth parasites have been reported expression KDPIs [19], such as in Schistosoma mansoni [12], S. japonicum [20] Fasciola hepatica [21,22], Ancylostoma caninum [23], Ancylostoma ceylanicum [24], and Steinernema carpocapsae [25]. Beside their role in inhibition of proteases such as pancreatic elastase, neutrophil elastase, chymotrypsin and trypsin, KDPIs play an important role in helminth's immune evasion and development, modulating the inflammatory responses, especially impairing Th1/Th17-associated inflammation response and reduction of neutrophil recruitment [19]. In this study, based on the genomic information available for E. granulosus s.s. and E. multilocularis we identified 23 and 19 KDPIs, respectively. These genes were differential expressed in different developmental stages, indicating these KDPIs may play different role in both parasite development and regulation the interface responses between parasites and hosts.
A remarkable difference between the larval stages of E. multilocularis and E. granulosus s.s. is the difference in the lesion pathology in the intermediate hosts. The metacestode of E. multilocularis is a tumor-like, infiltrating structure consisting of many small vesicles embedded in the stroma of connective tissue. The continual growth of parasite vesicles in a proliferative style causes damage of liver tissues, which results in a high mortality of AE patients. In contrast, E. granulosus s.s. cysts develop in internal organs (mainly liver and lungs) of humans and other intermediate hosts as unilocular fluid-filled bladders with clear edge between cyst and host tissue. CE causes mortality in very few patients and there is a relatively good prognosis after surgical removal of the cystic lesion. Contrastingly, AE causes severe damage to the liver and patients require extensive treatment with albendazole to prevent relapse. However, little is known about the molecular mechanisms underpinning biological differences between the two parasites and the diseases they cause.
The differential expression of these KDPI genes between E. granulosus s.s. and E. multilocularis may be associated with the differences in pathology caused by the metacestodes of the two species. It would be informative to determine whether these genes play a role in determining the different pathologies resulting from infection by the two cestodes in their intermediate hosts.
Signal peptide analysis showed that 89.47% of E. multilocularis KDPIs contain a signal peptide compared to only 60.87% of E. granulosus s.s. KDPIs containing signal peptide sequences. It may indicate that the proportion of secretion protein in E. multilocularis KDPIs is relatively higher than that of E. granulosus s.s. KDPIs. It is not known whether the high percentage KDPIs of E. multilocularis containing signal peptides is associated with the virulent pathology of AE lesion.
E. granulosus s.s. has 5 genes EG_07242, EG_07266, EG_07243, EG_09006 and EG_09008, that E. multilocularis does not have. Whereas, these two genes, EmuJ_001136700.1 and EmuJ_001137100.1 are only existed in E. multilocularis genome. These differentially presented genes may play a role in the difference of pathology between the two parasites.
Gastrointestinal helminths survive in an environment containing proteases and these parasites must have mechanisms to control protease activation. Therefore, Kunitz domain inhibitors are important for parasite survival, especially intestinal dwelling helminth parasites, to counteract protease attack. Two Echinococcus stages, the oncosphere and adult worm, are found in the gastrointestinal duct. The oncosphere is activated in the stomach and penetrates through the intestinal wall before being passed into the internal organs, whereas the adult worm spends its whole life in the gastrointestinal duct which contains high concentrations of proteases such as pepsin, trypsin and chymotrypsin. We previously showed that two KDPIs, EgKI-1(EG_08721) and EgKI-2 (EG_7242) function as protease inhibitors. EgKI-1 (also has accession number EUB56407.1) is highly expressed in the oncosphere and EgKI-2 (GenBank: EUB57880.1) is highly expressed in the adult worm [14]. These KDPIs are differentially expressed and stage-specifically protect E. granulosus s.s. from protease attack [13]. In this study, we showed that 11 out 23 Eg-KDPIs were highly expressed in adult worms. These Eg-KDPIs likely protect against protease attacks in the gut during adult worm development. EG_05483 and EG_08721 were relatively highly expressed in oncospheres, suggesting their expressed products might be potential vaccine We did not find any KDPIs that were differentially and highly expressed in protoscoleces in this study, although a previous study described a multigene family of eight (EgKU1-EgKU8) secreted Kunitz proteins from E. granulosus s.s. protoscoleces preferentially expressed by pepsin/H (+)-treated worms [15].
The secondary structures of proteins, especially the α-helix and β-strands play key roles in molecular function, cell stability, mechanical signaling, and tissue constitution as random coils are easily folded and exposed to the protein surface [26]. The basic structure of a Kunitz peptide domain contains a typical sequence with six highly conserved cysteine residues connecting three disulphide bridges (C1-C6, C2-C4 and C3-C5) which stabilizes the protein structure. Among the disulphide bridges, the C1-C6 and C3-C5 bridges are required for the maintenance of native confirmation [27],whereas the C2-C4 bond stabilizes the folded structure [28]. We found 10 sequences had lost the #2 cysteine, including 5 from E. granulosus s.s. , indicating no C2-C4 bridge in these proteins. The reduction of disulfide bonds may affect the stability of the protein. It is not known whether these 5 proteins formed different bridges impacting on the function of these KDPIs, indicating that these genes may have a different functional role.
Hydrophilicity analysis showed that the Em-and Eg-KDPIs have high hydrophobicity, which is a typical characteristic of membrane proteins. The transmembrane regions consist of 20 hydrophobic amino acids, which could have an anchoring effect on cell membranes.
We previously showed that EgKI-1 is highly expressed in the oncosphere, indicating this protein helps protect this stage from digestion by trypsin, chymotrypsin and pancreatic elastase before it penetrates the intestinal wall.
In this study, qPCR results showed that 3 out 4 KDPI genes of E. multilocularis homologous to those E. granulosus s.s. KDPIs were expressed in different pattern in three stages of this tapeworm including adult worm, protoscolex and cyst germinal membrane. It will be interesting to identify whether the expression of KDPIs is associated with the pathological difference in both intermediate hosts and definitive hosts of these two tapeworms.
It is predicted that EgKU-1 and EgKU-4 functionally blocks voltage-dependent potassium channels (Kv). The similar sequences of E. multilocularis EmuJ_001136700 and EmuJ_001137100(Additional file 8: Fig. S6), may likely function as a Kv blocker [29], which may be a clue for designing drug targets against alveolar echinococcosis. In addition, Eudiplozoon nipponicum EnKT1 (a KDPI) is an antihemorrhagic snake venom factor-like protein which exhibited a higher activity against plasmin and Factor Xa which can act as a C3 and C5 convertase and impaired both haemostasis and complement activation [30]. The different expression of similar KDPIs in E. granulosus s.s. and E. multilocularis may be associated with the difference in cellular pathology of these two echinococcal diseases. It needs to be uncovered in future studies.
It is shown that parasite KDPIs with immunomodulatory activity can be vaccine candidates. Fasciola hepatica FhKTM decreased dendritic cell activation and may be involved in the immune evasion mechanisms of the parasite [21]. FhKTM induced mice protection against F. hepatica challenge by preventing liver damage and improving survival, likely through eliciting potent IFNgamma and IL-17 A with high levels of antigen-specific IgG1, IgG2a, and IgA serum antibodies [31]. The KDPI EGR_05316 of E. granulosus s.s. has sililar sequences with Smp_147730 of S. mansoni [12,32]. Two KDPIs EGR_05483 and EmuJ_000077800.1 of E. granulosus s.s. and E. multilocularis have similar sequences with FhKT1.1, FhKT1.2 and FhKT2 of F. hepatica [33] (Additional file 7: Fig. S5).
A study showed that Schistosoma mansoni Kunitz peptides were highly protective in vaccinated BALB/c mice in terms of reduction in recovery of adult females (89~91%) and in the numbers of eggs trapped in the livers (77~81%) and guts (57~77%)of mice [32]. In addition, SmKI-1 showed 23~33% of reductions in adult worm [13]. EGR_03481, EGR_07242 and EGR_05316 of E. granulosus s.s. KDPIs have very similar sequences with Smp_147730 ( Fig. 2), this indicates that these KDPIs of Echinococcus likely can be candidates for vaccination against echinococcosis.

Conclusions
In conclusion, based on whole genome analysis, 19 and 23 Kunitz domain protease inhibitors were identified in E. multilocularis and E. granulosus s.s. . The differential expression of these KDPIs in different developmental stages of E. granulosus s.s. suggests that they may have different functions in regulation of host immune responses. The difference in characterization of KDPIs may be associated with the different pathology of metacestode stage of these two parasites. These should be further illuminated to determine their roles in echinococcal development and interface interactions between host and the tapeworms and such information may provide new insights for the prevention and treatment of cystic and alveolar echinococcosis.

Identification of E. granulosus s.s. and E. multilocularis Kunitz domain sequences
The E. granulosus s.s. and E. multilocularis genomes were previously completed by the Chinese National Human Genome Centre in Shanghai (CHGC) and the Wellcome Sanger Institute, United Kingdom in 2013 [18,34]. Genome data is available from http:// www. sanger. ac.uk/ resources/downloads/helminths/ (E. multilocularis, E. granulosus). The complete genome annotation is available at www. genedb. org.
E. granulosus s.s. , E. multilocularis and Fasciola hepatica sequence data reported in this manuscript are accessible at WormBase ParaSite and the corresponding accession numbers indicated in Additional file 9:
The multiple sequence alignment was analyzed and ordered by Clustal omega (http:// www. ebi. ac. uk/ Tools/ msa/ clust alo/) and then visually edited with BioEdit software v7.1.3. The phylogenetic tree was constructed by MEGA version 7 (http:// www. megas oftwa re. net/) and Interactive Tree of Life iTOL v6 (https:// itol. embl. de) with the maximum likelihood method. Bootstrapping analysis was performed and the bootstrap values which display support values (1000 bootstraps) are shown on the nodes. The clades are represented by different color. The complete KDPI protein sequences of single domain KDPIs were used for phylogenetic tree analysis.

Expression of Kunitz domain inhibitors in E. granulosus s.s. and E. multilocularis
Transcript reads were obtained for each of the KDPI genes expressed in the adult worm, oncosphere, protoscolex and cyst (cyst germinal membrane) of E. granulosus s.s. using Hiseq techniques as described [18]. Based on these data, the expression of 8 KDPI genes from the E. multilocularis and E. granulosus s.s. expression was validated for adult worm, protoscolex and cyst stages using quantitative reverse transcription PCR (qRT-PCR). Isolated protoscoleces were washed 10 times with saline and stored at -80 ℃ until used. TRIzol reagent (Invitrogen, USA) was used for extraction of total RNA according to the manufactory's instructions. RNA concentration was determined by a spectrophotometer (NanoDrop 2000, Thermo Scientific, USA) at 260 nm and the purity of RNA was considered satisfactory if the ratio of absorbance at 260 nm and 280 nm(A260/280) is ranged from 1.9 to 2.0.
Total RNA was reverse transcribed into cDNA using Reverse Transcriptase Kit (Takara, Dalian, China) following the instructions of manufacturer. The cDNAs were stored at −80 °C until used.

Quantitative PCR (qPCR)
Gene-specific primers were designed with Primer 3 software and listed in Table 4, eif3 was used as the internal reference [35].The qPCR reactions were carried out in a CFX96 Real-Time PCR Detection System (Bio-Rad, USA). Experiments were performed with Quanti-Nova SYBR Green PCR Kit (QIAGEN). Each individual sample was run in triplicate. The qPCR cycling reactions were initially denatured at 95 °C for 2 min followed by 39 cycles of 95 °C for 5 s, 60 °C for 30 s. The melting curve analysis of 5 s per step from 65 to 95 °C after amplification was conducted to assess primer specificity. The mRNA expression level of the target cytokines relative to the reference gene were analyzed using the 2 −ΔΔCt method [36].

Statistical analysis
Data are presented as means or median. All data are presented as the means± standard deviation (SD) of three individual experiments unless otherwise stated. Group comparisons were assessed by Two-tailed Student's t test, Mann-Whitney U test and one-way analysis of variance (ANOVA) for statistically significant differences using GraphPad Prism software (Version 8.01). Chi square test followed by Fisher's Exact Test was used to compare the sample rate (or constituent ratio) between the two groups. P values of <0.05 was considered significant in statistical analysis. (* P value ≤ 0.05; ** P value ≤ 0.001; *** P value ≤ 0.0001).