mRNA and P-element-induced wimpy testis-interacting RNA profile in chemical-induced oral squamous cell carcinoma mice model.

P-element-induced wimpy testis (PIWI)-interacting RNAs (piRNAs), a novel class of noncoding RNAs, are involved in the carcinogenesis. However, the functional significance of piRNAs in oral squamous cell carcinoma (OSCC) remains unknown. In the present study, we used chemical carcinogen 4-nitroquinoline-1-oxide (4NQO) induced OSCC mouse model. piRNAs and mRNAs were profiled using next-generation sequencing in the tongue tumor tissues from 4NQO induction and healthy tongue tissues from control mice. Furthermore, we analyzed the differential gene expression of human OSCC in Gene Expression Omnibus (GEO) database. According to the common differentially expressed genes in the 4NQO model and human OSCC tissues, piRNAs and mRNAs network were established based on informatics method. A total of 14 known piRNAs and 435 novel predicted piRNAs were differently expressed in tumor tissue compared to healthy tissue. Among differently expressed piRNAs 260 were downregulated, and 189 were upregulated. The mRNA targets for the differentially expressed piRNAs were identified using RNAhybrid software. Primary immunodeficiency and herpes simplex infection were the most enriched pathways. A total of 22 mRNAs overlapped in human and mice OSCC. Moreover, we established the regulatory network of 11 mRNAs, including Tmc5, Galnt6, Spedf, Mybl2, Muc5b, Six31, Pigr, Lamc2, Mmp13, Mal, and Mamdc2, and 11 novel piRNAs. Our data showed the interaction between piRNAs and mRNAs in OSCC, which might provide new insights in the development of diagnostic biomarkers and therapeutic targets of OSCC.

Introduction

Oral squamous cell carcinoma (OSCC) is the most prevalent oral malignant tumor [7, 36]. Despite the rapid development of therapeutics, the 5-year survival rates for OSCC (50% for intermediate cancer patients and <20% for terminal cancer patients) have not altered significantly over the past five decades [12]. Typically, the treatment and outcome prediction depend on the traditional tumor node metastasis (TNM) classification system [14], which is neither precise nor sufficiently useful. Thus, the underlying regulatory mechanisms of OSCC should be elucidated for accurate diagnosis and treatment.

P-element-induced wimpy testis (PIWI)-interacting RNAs (piRNAs) are a class of small noncoding RNAs, firstly identified in MILI protein in mouse testis in 2006 [1]. microRNAs (miRNAs) are usually 21–23 nucleotides (nt) in length, while mature piRNAs are composed of 26–31 nt. Furthermore, piRNAs are distinguished by 2’-O-methylation at the 3’ end [20]. These molecules repress the transposon expression by methylation or chromatin modification in the nucleus [29]. In addition, piRNAs degrade the retrotransposon-associated mRNAs and promote the maturation of these mRNAs by cleavage in the cytoplasm [29]. Reportedly, piRNAs play an essential physical role in germline maintenance and epigenetic regulation [13, 27]. It has been demonstrated that piRNAs function as oncogenic or suppressive regulators in several cancers. piRNAs may also serve as noninvasive markers for the early diagnosis of colon cancers. piR-5937 and piR-28876 are shown to be downregulated in the serum of patients with colon cancer [28]. Intriguingly, two mitochondria-related piRNAs, piR-34536 and piR-51810, are shown to serve as prognostic biomarkers in the tissues of clear cell renal cell carcinoma but not as biomarkers in the serum [35]. piR-1245 is a potential prognostic marker that promotes the progression of colorectal cancer [30]. The overexpression of piRNA-8041 causes cell cycle arrest, apoptosis, and inhibits cell survival pathways as well as cell survival of human glioblastoma in vitro and in vivo [16]. Additionally, a novel piR-823 suppresses the activity of gastric cancer in vitro and in vivo [5]. piR-39980 exhibits anti-tumor effects in fibrosarcoma by repressing ribonucleotide reductase M2 subunit (RRM2) through direct binding at its 3’-untranslated region (UTR). However, this interaction is extensively complementary, which is different from the seed binding of microRNAs [9]. However, the expression profile of piRNAs in OSCC and their role in OSCC development and progression have not been investigated yet. Owing to the importance of piRNAs on cancer prediction and therapy, the functional and clinical significance of piRNAs in OSCC should be evaluated further.

4-nitroquinoline-1-oxide (4NQO) is one of the most intensively studied chemical carcinogens for OSCC induction. It has been demonstrated that the animal model induced by 4NQO can imitate human tumorigenesis at the histological and molecular levels [18]. Several factors are involved in the progression of 4NQO-induced OSCC. MCT4 has the potential as an early diagnostic biomarker, and hence, is essential to drive the progression of OSCC [2]. Tang et al. assessed the gene expression profile to compare the early-stage and late-stage 4NQO-induced OSCC for early diagnosis [25].

Herein, we analyzed the piRNAs and mRNAs expression profile in 4NQO-induced OSCC mouse model. We established the piRNAs-mRNAs regulatory network involved in OSCC. We also used the human OSCC mRNA profile from GEO database to identify the overlapping genes with the mouse OSCC mRNA.

Materials and Methods

Animals and ethics approval

The experimental protocol was approved by the Institutional Animal Care and Use Committee (IACUC) of the First Affiliated Hospital of Guangzhou Medical University of China (Guangdong, China). The present study was conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals.

A total of 36 male FVB mice (6-week-old) were obtained from Southern Medical University of China. All mice were acclimatized to the environment for five days before the experiment. Each cage housed 3–4 animals that were allowed free access to water and food. Every mouse was exposed to 12 h light/dark cycle. The mice were anesthetized (xylazine, intraperitoneal injection, 0.13 mg/kg body weight) and euthanized by exsanguination prior to tissue collection.

Animal model, experimental design, and histological analysis

The experimental flow chart is shown in Fig. 1. Mice were randomly divided into two groups (n=18 each). In the experimental group, mice were fed 50 µg/ml 4NQO (Sigma–Aldrich, Saint Louis, MO, USA) daily in the drinking water, while the control group mice drank water that contained an equivalent volume of propylene glycol. The drinking water was replaced every week. All the animals were weighed weekly. Tongue samples were collected from 8 mice per group after challenging with 4NQO for 25 weeks. The remaining mice (10 mice/group) were subjected to survival analysis. Tongue samples were collected in each group, sectioned in paraffin, and stained with hematoxylin and eosin. Microscopic images were observed, captured, and analyzed using Eclipse Ti-UNIS-Elements software (Nikon Corp., Tokyo, Japan) and Image-Pro Plus version 6.0 software (Media Cybernetics Inc., Rockville, MD, USA).

Fig. 1.

Schematic illustration showing the experimental set up. b: the time for oral and tongue examination. Black arrow: the beginning to 4-nitroquinoline-1-oxide (4NQO) treatment; red arrow: mice sacrificed for tissue collection; green arrow: the endpoint examination of survival rate in mice.

RNA-seq

Total RNA isolation, library preparation, and sequencing were accomplished by BGI (Shenzhen, China). Briefly, total RNAs were extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. Firstly, 5-adenylated, 3-blocked single-stranded DNA adapters were ligated to the small RNA fragments. Reversed transcription (RT) primer hybridization was performed to link with the 3’ adapter. Secondly, 5’ adapter ligated to the 5’ terminus of the small RNAs, and the first-strand cDNA was generated using RT primer for reverse transcription. Subsequently, PCR amplified the cDNA fragments that contain both 3’ and 5’ adapters. The resulting amplicons were purified by PAGE electrophoresis to isolate 100–120 bp fragments from the gel. The double-stranded PCR products were heat-denatured and cyclized to form single-strand circle DNA (ssCir DNA), which formatted as the final library. The library was quantified using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and sequenced using the BGISEQ-500 platform. EBSeq was used for RNA quantification and analysis according to the changed reads counts. Fisher’s exact test was applied to select the significant gene ontology (GO)-term and pathway. The threshold of significance was determined by P<0.05 or FDR<0.05.

piRNA-seq

The small RNA molecules were separated from total RNA using polypropylene acyl amine gel electrophoresis (PAGE). 18–30 nt RNAs were selected and analyzed as small RNA molecules. Firstly, 5-adenylated, 3-blocked single-stranded DNA adapters were ligated to the small RNA fragments. The RT primer hybridization was performed to combine with 3’ adapter. Secondly, the 5’ adapter was ligated to 5’ terminus of small RNA, and the first-strand cDNA was generated using the RT primer for reverse transcription. Subsequently, PCR amplified the cDNA fragments which contain both 3’ and 5’ adapters. The PCR products were purified by PAGE to isolate 100–120 bp fragments from the gel. Quantification and circularization of the PCR products were performed as described above. The library was quantified using the Agilent 2100 Bioanalyzer (Agilent Technologies) and sequenced on the BGISEQ-500 platform. The targeting genes of piRNAs were predicted based on RNAhybrid. The GO and Kyoto encyclopedia of genes and genomes (KEGG) pathway were analyzed according to predicted targets of piRNAs. Also, Fisher’s exact test was applied to select the significant GO-term and pathway. The threshold of significance was determined by P<0.05 or FDR<0.05.

RNA data from Gene Expression Omnibus (GEO) database

The published microarray gene expression profile dataset of OSCC patient cohort was selected from the NCBI GEO (https://www.ncbi.nlm.nih.gov/geo/), accession number GSE31056. It includes 23 tumors and 73 margins from 24 patients with OSCC (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31056). Limma was used to calculate Reads Per Kilobase of exon model per Million mapped reads (FPKM) to detect the altered genes in the GEO data. The log|fold-change|≥1, P<0.05 was considered to be significant. The genes with the same trend between human and mice OSCC were analyzed further.

Integrated analysis of mRNA-seq and piRNA-seq

The common differentially expressed genes between human and mice OSCC were further established the interaction network with piRNAs by Cytoscape. The interaction was based on the previous piRNAs and predicted mRNA data.

qRT-PCR analysis

Total RNA from control or OSCC tongue tissues (n=4) was extracted using TRIzol reagent (ThermoFisher Scientific, Waltham, CA, USA) following the manufacturer’s protocol. The integrity of total RNA was analyzed using denaturing agarose gel electrophoresis. The purity and concentration of RNA was verified using the NanoDrop spectrophotometer. The cDNAs for piRNA was reverse transcribed using Mir-XTM miRNA First-Strand Synthesis Kit (Clontech, Dalian, China). Expression level was evaluated using TB Green® Premix Ex TaqTM II (TaKaRa, Tokyo, Japan) following the manufacturer’s Protocol. The amplification condition was 95°C for 30 s, followed by 40 amplification cycles of 95°C for 10 s, 60°C for 30 s, and 95°C for 30 s, 60°C for 30 s, 95°C for 30 s. U6 small nuclear (U6) was used to calculate relative expression level of piRNAs as an internal control.Relative expression levels were determined by applying the ΔΔCt method using U6 as the endogenous control. Sequences of Primers used for qPCR are showed in Supplementary Table 1.

Statistical analysis

The statistical analyses were conducted using SPSS 22.0 software. Data of each group were presented as mean ± SD. The Kaplan-Meier method was employed for the survival analysis in the OSCC mouse model. P-values <0.05 were considered statistically significant.

Results

Establishment of oral squamous cell carcinoma mice model successfully

To perform bioinformatic analysis of piRNAs-seq, we established the 4NQO induced OSCC model. The mean bodyweight was reduced in the 4NQO group, while significantly increased in the PG group. Also, the mean body weight of the 4NQO group was significantly lower than that of the PG group after week 20 (Fig. 2A). Moreover, the 4NQO group showed a shorter survival as compared to the PG group. Mice in the 4NQO group succumbed to mortality at the 40th week (Fig. 2B). As shown in Fig. 2C, 4NQO exposure successfully induced proliferative lesions on the tongue, ventrum, and dorsum, which were indicated by yellow arrows. Histopathological assessment exhibited epithelial dysplasia in the 4NQO group as compared to the PG group (Fig. 2C). The lesions showed plentiful nuclear pleomorphism and abnormal mitosis as well as lack of keratosis (Fig. 2C), which was consistent with the characteristics of low-differentiated carcinoma. Thus, the 4NQO-induced OSCC model was established.

Fig. 2.

Pathophysiological examination of oral squamous cell carcinoma (OSCC) mice model. (A) Mean body weight curve after treatment, (B) Survival analysis., (C) Gross examination of OSCC developed in tounge (left panel), and histopathological examination of tongue tissue (hematoxylin and eosin staining, right panel). Yellow arrows: OSCC (white spots) on tongue dorsum and papillary, and cauliflower-like lesions on tongue ventrum.

piRNAs-seq was used to detect the piRNAs involved in the carcinogen-induced OSCC. As shown in Fig. 3A, piRNAs clustered in the same group. Here, 246 known piRNAs and 1007 novel piRNAs were identified in the OSCC group, while the control group consisted of 210 known piRNAs and 910 novel piRNAs (data not shown). Also, further analysis of different piRNAs in the induction and control groups revealed one known piRNA and 16 novel piRNAs having significant expression difference between OSCC group and control group (Fig. 3B). Relevant target genes of piRNAs were predicated by RNAhybrid (Energy<−25). According to GO classification, the targeted genes were classified into 123 subtypes. Among these, immune system process represented the biological processes adequately. The external side of the plasma membrane was the mainly cell component. Notably, that hydrolase activity on the acid anhydrides and peptide antigen binding was a vital molecular function (Fig. 3C). KEGG pathway analysis demonstrated that primary immunodeficiency comprised the largest group, followed by herpes simplex infection (Fig. 3D).

Fig. 3.

Bioinformatic analysis of 4-nitroquinoline-1-oxide (4NQO) induced P-element-induced wimpy testis-interacting RNAs (piRNAs). (A) Heatmap of piRNA-seq. (B) Altered piRNA expression in healthy and oral squamous cell carcinoma (OSCC) mice tongue tissue. (C) gene ontology (GO) analysis. (D) Kyoto encyclopedia of genes and genomes (KEGG) analysis.

Integration of RNA-seq and GEO database

RNA-seq has been widely used to detect the changed genes which affect the function of OSCC. Herein, according to the sequencing dataset, 21491046 reads and 21964531 reads were identified in the control and induction groups, respectively. Furthermore, 1727451 and 17675560 unique reads were mapped to the genome in the control and 4NQO groups, respectively. Clustering analysis showed that the gene expression levels displayed a higher clustering among the same group (Fig. 4A). In addition, we identified 424 upregulated genes and 38 downregulated genes in the 4NQO group as compared to the control group (log|fold-change|≥1, P<0.05) (Fig. 5Aa). The GO cluster and KEGG pathway further revealed the functional genes in 4NQO-induced OSCC. Peptide-antigen binding, hydrolase activity, and T cell receptor binding were the main molecular function (Fig. 4B). Antigen processing and presentation and herpes simplex infection were the main enriched pathway (Fig. 4C).

Fig. 4.

Bioinformatic analysis of 4-nitroquinoline-1-oxide (4NQO) induced genes in mice. (A) Heatmap of RNA-seq. (B) gene ontology (GO) analysis. (C) Kyoto encyclopedia of genes and genomes (KEGG) analysis.

To further reveal the common differentially expressed genes of OSCC between human, and mouse, published microarray dataset (GSE31056) from 24 patients was analyzed. A total of 502 genes were altered significantly: 260 were downregulated and 189 were upregulated (log|fold-change|≥1, P<0.05) (Fig. 5Ab). Finally, 26 differentially expressed genes were identified to be shared between human and mouse tongue OSCC tissue. We eliminated 4 genes in subsequent steps due to contrasting changed trend. Among these 22 genes, 3 were downregulated and were 19 upregulated genes in both human and mouse (Fig. 5B and Table 1).

Fig. 5.

Joint analysis of differentially expressed genes in mice and human oral squamous cell carcinoma (OSCC) tissues. (A) Volcanic plot of the degree of differences in the mRNA expression profile of OSCC in (a) mice model and (b) human. X-axis: log2 (fold change); Y-axis: -log10(P value); red plots: significantly upregulated genes; blue plots: significantly downregulated genes; gray plots: no significant difference in gene expression. Each dot represents one gene. (B) Venn diagram of mRNAs in mice and human OSCC tissues.

Table 1.

Differentially expressed genes shared between human and mouse.

Gene symbol	Mouse		Human
	log2FC1	FDR	log2FC2	Pvalue
Agr2	6.82981	6.73E-07	1.93714	0.04887
Muc5b	6.69051	6.64E-06	2.23298	0.049096
Mmp13	5.36666	0.000446	1.99152	0.010736
Cxcl3	5.06897	0.019397	1.56553	0.004034
Azgp1	4.48639	0.00626	2.79622	0.006182
Spdef	4.46587	0.002871	1.26227	0.027927
Smim22	4.17139	0.009688	2.05019	0.002819
Msln	4.11138	0.036795	1.21168	0.003502
Pigr	4.09333	0.024075	1.75301	0.010392
Tmc5	3.91967	0.022035	1.12481	0.020356
Cxcl5	3.84423	0.001197	1.23063	0.034295
Tspan1	3.66396	0.027066	2.43418	0.001052
Galnt6	3.5792	0.017623	2.48925	7.13E-05
Irx2	3.46095	0.000182	1.84839	0.002336
Mybl2	1.59186	0.001786	1.08741	0.000762
Cdc6	1.45573	2.06E-05	1.12566	0.010822
Pthlh	1.2333	0.00058	1.67722	0.022154
Dsg2	1.09247	0.011109	1.87739	0.001406
Lamc2	1.04804	2.85E-08	1.5097	0.008976
Snx31	–3.6658	4.02E-10	–1.7355	0.002078
Mal	–1.5406	0.045645	–4.0425	3.40E-05
Mamdc2	–1.0031	2.84E-06	–1.3796	0.009422

Correlation analysis between piRNAs and mRNAs

To reveal the regulatory network, the correlation analysis of piRNAs and mRNAs was conducted. According to the piRNA-seq data, the targets of piRNAs were predicted by RNAhybrid. The piRNA/ predicted mRNA pairs were established. Then, we acquired 22 common differentially expressed genes from human and mice OSCC. piRNAs could silence the gene expression by “ping-pong” cascade [15], thus, the expression level of piRNAs and their targeting mRNAs were negative correlation. However, 11 genes of total 22 common differentially expressed genes could be paired with piRNAs. The negative regulatory interaction between 11 genes and 11 piRNAs was established according to the piRNA/ predicted mRNA pairs (Fig. 6). The 11 piRNAs were as follows: novel pir354, novel pir415, novel pir494, novel pir443, novel pir444, novel pir446, novel pir984, novel pir162, novel pir1185, novel pir1584 and novel pir832 (Table 2). These 11 genes were as follows: Galnt6, Spdef, Muc5b, Tmc5, Pigr, Lamc2, Mmp13, Mybl2, Mamdc2, Mal, and Snx31. These data presented specific functional network in the development of cancer.The qRT-PCR results of pir1584, pir354, pir415, pir446, pir1185, pir984 and pir494 were in line with the piRNA-seq results.

Fig. 6.

Correlation analysis of common differentially expressed 11 novel P-element-induced wimpy testis-interacting RNAs (piRNAs) and 11 mRNAs.

Table 2.

Differentially expressed P-element-induced wimpy testis-interacting RNAs (piRNAs) in mouse

piRNA Name	piRNA-seq		qRT-PCR
	log2FC	FDR	log2FC	P value
novel_pir494	3.14824	0.02252	4.19470	0.00839
novel_pir1185	2.67348	0.02909	2.36051	0.07796
novel_pir354	–3.57708	0.00602	–4.31661	0.00013
novel_pir1584	–3.21002	3.99E-07	–1.86931	0.08440
novel_pir446	–2.70392	0.01790	–4.46402	0.00625
novel_pir984	–2.13506	8.84E-11	–1.66983	0.00063
novel_pir415	–3.07616	0.00030	–2.76652	0.01861
novel_pir443	–0.90079	0.05417	–	–
novel_pir832	–1.24083	5.57E-07	–	–
novel_pir444	–1.55797	8.77E-05	–	–
novel_pir162	–1.77282	3.41E-07	–	–

-,not detected.

Discussion

OSCC is the sixth most common cancer worldwide. Although the underlying tumorigenesis mechanism has been under intensive focus, the pathogenesis needs to be elucidated further to increase the therapeutic effects. Noncoding RNAs are essential factors for the regulation of cancer progression. 4NQO model is widely used for screening the drug and exploring the mechanism. Interestingly, miR-211, targeting to TCF12 and TGFβRII, promotes the progression of OSCC induced by 4NQO [6, 26]. Furthermore, miR-31 increases the susceptibility of 4NQO-induced OSCC by repressing the expression of Ku80 [4]. A previous study reported that RNA-seq detected the mRNAs changes at early and late stages of post-4NQO treatment [25]. piRNAs are a type of novel noncoding RNAs, which is different from microRNAs and implicated in the cancer progression. piRNAs could inhibit mRNA expression by direct binding. piR-55490 plays an inhibitory role in lung cancer by targeting mTOR signaling [24]. piR-39980 plays an inhibitory role in fibrosarcoma by repressing RRM2 [9]. However, the piRNAs expression profile in OSCC has not been investigated yet. Herein, we elucidated the regulatory network of RNA-seq and piRNA-seq profiles in 4NQO-induced OSCC.

The targets of changed piRNAs were predicted by RNAhybrid. Consequently, 424 upregulated genes and 38 downregulated genes were identified in the 4NQO group as compared to the control group (Fig. 5Aa). In addition, 260 downregulated genes and 189 upregulated genes were identified in the GEO data on human OSCC and adjacent tissues (Fig. 5Ab). We further identified 22 common differentially expressed genes in human OSCC and 4NQO-induced mouse OSCC (Table 1). However, we could find only 11 genes paired with piRNAs. Finally, the regulatory network consisting of 11 genes and piRNAs was established (Fig. 6). Among the 11 genes, GALNT6, SPEDF, and MYBL2 have been reported to involve in the suppression or progression of several tumors [11, 17, 31]. Also, Matse et al. showed that the expression of MUC5B is switched on during the carcinogenesis of mucoepidermoid carcinoma in the parotid gland, and the incomplete glycosylation of MUC5B is responsible for the aberrant glycosylation machinery [21]. TMC5 is a potential biomarker for distinguishing lung adenocarcinoma from squamous cell carcinoma [32], while only a few studies have investigated the role of SNX31 in cancer. As our knowledge, the function of Galnt6, Spedf, Mybl2, Muc5b, Tmc5 and Snx31 on the progression of OSCC has not been reported yet. We found that Galnt6, Spedf, Mybl2, Muc5b, Tmc5 were upregulated while Six31 was downregulated in OSCC (Table 1). The specific function of these genes in OSCC development and progression should be explored further.

Several studies demonstrated that LAMC2 expression promoted the malignant progression of OSCC [10, 22]. MMP13 was found to be involved in the migration and invasion of squamous cell carcinoma cells. Yu et al. suggested that CXCR4 regulated MMP13 via the activation of extracellular signal-regulated kinase (ERK) signaling pathway in OSCC [33]. Furthermore, Smad signaling plays a major role in upregulating MMP13, which in turn, enhances the invasion of head and neck squamous carcinoma cells [19]. Also, MMP13 is upregulated in the human OSCC tissues and 4NQO induced mouse model [25, 34]. miR-196a exerts a pro-tumorigenic effect in OSCC by directly targeting MAMDC2 [8]. Pal et al. confirmed the low expression level of MAL in OSCC, which is a potential biomarker for targeting OSCC [23]. Furthermore, MAL epigenetic silencing significantly increases the tumorigenesis of OSCC in vitro and in vivo [3]. In the present study, the expression of LAMC2 and MMP13 was upregulated, while that of MAMDC2 and MAL was downregulated (Table 1). Furthermore, only seven piRNAs were confirmed by RT-qPCR in this paper, other piRNAs, genes and proteins should be further confirmed by RT-qPCR or westerm blotting. Also, the functional mechanism underlying the interaction of piRNAs with these genes needs to be explored.

In conclusion, we analyzed the expression profile and the interaction of piRNAs and mRNAs in the 4NQO-induced OSCC mice model. The novel piRNAs in mice OSCC and common differentially expressed genes in human and mice were identified. We also established the interaction network of 11 common differentially expressed genes and 11 piRNAs. Primary immunodeficiency and herpes simplex infection were the most enriched pathways. The results of this study provide new insights, which could be helpful in the development of diagnostic biomarkers and therapeutic targets of OSCC.

Author Contributions

LW, YJ, ZZ, HL, MC, ZL and LH organized and conducted this study. JG, KO, and HZ conceived and designed the experiments and analyzed the data. LW, YJ, ZZ, JLP and JG wrote the manuscript and MZ, JLP polish the English writing.