The insects, living under natural light, often have ingenious color patterns that provide great aid to them, such as in behavior1, immunity2, thermoregulation3, UV protection4, and especially in antagonist to visual predators by either aposematic5, mimicry6, cryptic color patterns or combination among them7. Impressively, the color patterns of insects, such as cryptic color, can deceive visual animals, based on their clever designed which highly similar the surrounding environment. In addition, the color patterns are divergent and convergent among populations8. Due to this gorgeous visual feature and highly active adaptive evolutionary phenotype, lead to the genetic basis and evolution mechanism of color patterns have been highly concern for a long time.

The body color of insects is actually main formed by the absorption, reflection or/and diffraction of natural light by the body surface, or/and some appendage (such as bristles and scales)9. Thus, the difference in epidermal structure and the difference in pigment deposits can form divergent color patterns. Cuticular proteins (CPs) are the main structural proteins that form the epidermis. There are large number genes encoding cuticular proteins in arthropods genome. For example, in Bombyx mori over 220 genes belong to the family of CPs10. These cuticular proteins belong to more than 10 subfamilies, including CPR, CPF, CPFL, CPH, Tweedle, CPT, CPG and et al11. However, the fine localization molecular mechanism of cuticular proteins in epidermis is still unclear. Moreover, the biochemical metabolic pathways of pigments, which used for color patten of insect, such as melanins, pterins, ommochromes, have been reported 12. The deficient of enzymes in pigments metabolism pathway can lead to the change of color pattern13. However, some pleiotropic factors, such as wnt114, Apontic-like15, clawless, abdominal-A, and Abdominal-B16, bric à brac (bab)17, Distal-less (Dll)18, play an important role in the regulation of color pattern. Nevertheless, the molecular mechanism by which these factors participate in the regulation of body color patterns is unclear. In addition, there are many undiscovered factors involved in the gene regulatory network (GRN) involved in the regulation of insect color patterns. Hence, the identification of new color pattern regulatory genes and the study of molecular mechanisms are helpful to further understand the color patterns plan.

Silkworm, Bombyx mori (B. mori), is a completely domesticated Lepidoptera, for more than 5000 years, is famous in silk fiber production 19. Because of long-term domestication, artificial selection and genetic research, its genetic background is clear, and about 200 strains of mutants involved in color pattern have been preserved20, which provides a good resource for color pattern research. The black dilute (bd) mutant, complied with the recessive Mendelian inheritance law, showed a dark gray larvae body color and its female is sterile. The bdf (black dilute fertile) is one of its alleles, the body color is lighter than that of bd and the female is fertile (fig.1). The two mutants were mapped to a single locus on 22.9 centimorgan (cM) of linkage group 9 20, 21. Because there are no pigmentation related genes reported at this locus. For it, research may discover a new color pattern related gene, which have aroused our concern.

Larval phenotype of the wild type Dazao, bd and bdf.

The epidermis of bd is dark black gray and the epidermis of bdf is light black gray at 5th instar day 3 of silkworm. The bar indicates 1 centimeter.


Candidate gene of bd allele of silkworm

To identify the genomic region responsible for the bd alleles, the positional cloning of bd locus was performed. Due to the female infertility of the bd mutant and the female of its allele bdfis fertile. So, we used bdf and wild type Dazao strain as the parents for mapping analysis.

The 1162 back-crossed filial 1 (BC1M) generation individuals, between bdf and Dazao, were used for fine mapping by using molecular markers (fig 2). It found that the genome region of about 390 kb (kilo base) was responsible for bd phenotype (table s1). This region included five predicted genes (BGIBMGA012516, BGIBMGA012517, BGIBMGA012518, BGIBMGA012519 and BGIBMGA014089). And then, qRT-PCR (Fluorescence quantitative PCR) was performed, which found that BGIBMGA012517 and BGIBMGA012518 were significant down-regulation expressed at day 3 of five instar of larvae in the bd phenotype individuals (fig s1). In addition, the two genes are consistent expression trend in those strain, respectively. To determine the gene structure of the BIGBMGA012517 and BIGBMGA012518, we using forward primers of the BIGBMGA012517 gene and reverse primers of the BIGBMGA012518 gene to amplify the cDNA from wild type strain Dazao. By gene cloning, the predicted genes BGIBMGA012517 and BGIBMGA012518 were proved to be one gene. For the convenience of description, we temporarily call it 12517/8 gene. The 12517/8 gene produce two transcripts, the long transcripts contain an open reading frame (ORF) with 2397 bp, and the ORF of short one was 1824 bp, which share the same 5’-terminal in wild type strain Dazao (fig s2). The 12517/8 is significant low expression in bdfmutant, and were found multiple variations in nearby region of this gene by analysis of comparative genome between bdf and Dazao (table s2). In addition, the 12517/8 is silence completely due to the deletion of DNA fragments (about 168 Kb) from first intron to upstream, and insert 3560 bp in bd mutant (fig s3).

Positional cloning of the bd locus

(A) Using 532 BC1 individuals to map the bd locus between the PCR marker C and I. The numbers above DNA markers indicate recombination events. (B) Using 1162 BC1individulas to narrow the bd locus to ∼400 kb region. (C) The partial enlarged view of responsible region of bd. In this region contains 4 predicted genes, BGIBMGA012517, BGIBMGA012518, BGIBMGA012519 and BGIBMGA014589. (D) The analysis of nucleotide differences of bd responsible region. The green block indicated the deletion of the genome in bd mutants. The grey lines indicated the SNP and the red lines indicated the Indel of bdfmutants.

To estimate the function of the gene 12517/8, we performed a BLAST search using its full-length sequence and found a transcription factor, mamo (maternal gene required for meiosis), in D. melanogaster, which had a highly sequence similarity to that of the 12517/8 gene (fig s4). So, we named the 12517/18 gene Bm-mamo.

The gene mamo, belong to Broad-complex, Tramtrack and bric à brac/poxvirus zinc finger protein (BTB-ZF) family. The BTB-ZF family, which use the zinc finger as recognition motif of DNA specific sequences, and the BTB domain promotes oligomerization and the recruitment of other regulator factors23. Most of these factors are typically repressors of transcription, such as N-CoR (nuclear receptor corepressor) and SMRT (silencing mediator for retinoid and thyroid hormone receptor)24, but also some are activators, such as p30025. Therefore, these displaying features of commonly serve as regulator of gene expression. The mamo is enriched in embryonic primordial germ cells (PGCs) in D.melanogaster. Deficiency of mamo are able to undergo oogenesis, but fail to execute meiosis properly and leading to female infertility in D.melanogaster26. However, the lack of the Bm-mamo gene not only leads to sterility, but also causes significant melanism, in caterpillar.

Expression pattern analysis of Bm-mamo

To analyze the expression profiles of Bm-mamo, we examined the quantitative RT–PCR. The different developmental stage expression level of Bm-mamo gene was investigated in the whole body from embryonic stage to adult stage in Dazao strain. The gene is highly expressed at molting stage of caterpillar and the expression was up-regulated in the later pupa and adult stage (fig 3A). This suggests that the Bm-mamo gene may respond to ecdysone titer and participate in the process of molting and metamorphosis of silkworm. In the investigation of tissue expression levels in the 5th instar 3rd day larvae of Dazao strain, the midgut, head, epidermis, have highly expression level, and the trachea, nerve, silk glands, testis, ovary, muscle, wing disc and malpighian tubules have moderate expression level, in blood and fat body were low expression level (fig 3B). This result suggests that the Bm-mamo may be involved in the development regulation of multiple tissues of silkworm. Due to the melanism of the epidermis of the bd mutant and the high expression level of Bm-mamo gene in the epidermis, we detected the expression level of this gene in the epidermis of the 4th to 5th instar of Dazao. In the temporal expression of epidermis, the Bm-mamo gene was up regulated in the molting period, and the highest expression was at the beginning of molting (fig 3C).

Spatiotemporal expression pattern of Bm-mamo.

(A) Temporal expression of Bm-mamo. The molting stage and adult stage this gene is (B) M: moulting, W: wandering stage, P: Pupal stage, A: adult stage. The 1th to 5 th denote the first instar of larvae to fifth instar of larvae, respectively. (B) Tissues expression of the 4th instar molting larvae. The Bm-mamo have relatively high expression level in midgut, head and epidermis. (C) Detailed analysis of Bm-mamo at the 4th larval stage in epidermis of Dazao strain. The Bm-mamo is up-regulation expression during molting stage. HCS indicates head capsule stage.

Functional analyses of Bm-mamo for pigmentation

In order to study the function of the Bm-mamo gene, we carried out RNAi experiment. Five microliter of siRNA were injected into haemolymph form stigma of day 3 of fourth instar, before the expression peak of Bm-mamo in larval stage, and conducted the electroporation experimental27, immediately. We found that there was significant melanin pigmentation in the epidermis of the newly molted 5th instar larvae. This indicate that deficient of Bm-mamo can cause melanin pigmentation (fig 4).

Effects of Bm-mamo siRNA injection on pigmentation formation.

(A) The siRNA was introduced by microinjection followed by electroporation. The “+” and “-” indicate the positive and negative poles of the electrical current. Scale bars, 1 cm (B) Partial magnification of siRNA experimental group and negative control group. Scale bars, 1 cm (C) and (D) Relative expression levels of Bm-mamo in a negative control and RNAi group as determined by quantitative RT–PCR analysis. The means ± s.d. Numbers of samples are shown at the upper right in each graph. **P<0.01, paired Student’s t-test. (NS, not significant).

In addition, the gene knockout was performed. Two sgRNAs were designed to target exon 4 of Bm-mamo by CRISPRdirect ( online program. The ribonucleoproteins (RNPs) formed by gRNA and recombinant Cas9 protein were injected into 450 silkworm embryos. Then, the phenotype identification and DNA cloning of the individuals of third generation were performed. In the G0 generation, individuals with mosaic melanization phenotype were found. These melanistic individuals were raised to moth and then crossed. The homozygous line of gene knockout was obtained through generations of screening. It is revealed that the gene edited individuals with significantly melanism body color and the female moth is sterile (fig 5).

Disruption of Bm-mamo using CRISPR/CAS9.

(A) Larval phenotype in G3 forth instar larva of Dazao targeted for Bm-mamo. Scale bars, 1 cm (B) The after 36 hours of laying eggs, the pigmentation of the eggs indicates that the eggs produced by knockout homozygous females cannot undergo normal pigmentation and development. (C) Genomic structure of Bm-mamo. Open reading frame (blue), untranslated region (black). The gRNA 1 and gRNA 2 sequences are shown. (D) Sequences of the Bm-mamo knock out individuals. Line 1 and line 2 deletion 15 and 71 bp, respectively.

It is indicated that the Bm-mamo gene has negative regulation on melanin pigmentation in caterpillar and participated in the reproductive regulation of female moths.

Downstream target gene analysis

Because Bm-mamo belongs to zinc finger proteins family, which specifically recognize downstream DNA sequences based on their zinc fingers. We have identified homologous genes of mamo in multiple species and conducted phylogenetic analysis (fig s5). In addition, by sequence alignment, found that the amino acid residues of the zinc finger motif of mamo-S protein are much conserved among 57 species, and the mamo-L protein are conserved among 30 species (fig s6). Then, we used a software ( to predict the DNA recognition sequence of Bm-mamo, especially the core sequence, GCGTG, which highly consistent with the binding sequence of mamo in those species (fig s7). Using the predicted position weight matrices (PWMs) of recognized sequences of Bm-mamo and FIMO (Find Individual Motif Occurences) software of MEME29 suit to perform whole genome scanning for possible downstream target gene. The 2 Kb genomic regions, which locate in upstream and downstream of the total 14623 predicted genes, were investigated. It is found that there are 10622 genes contain the recognition sites of Bm-mamo protein in silkworm genome (fig 6). Moreover, we used a comparative transcriptome data between homozygotes and heterozygotes of bd mutant in 4th instar beginning molting, which is high expression level stage of Bm-mamo30. In the integument tissue of this period, there were 182 genes with significant expression differences between homozygotes (bd/bd) and heterozygotes (+/bd) (table s3). Further, we found 147 genes that are differentially expressed in the transcriptome data and have binding sites of Bm-mamo protein near the upstream and/or downstream.

The DNA recognize sequence logo of Bm-mamo protein and the analysis of downstream target genes.

(A) The DNA recognize sequence logo of Bm-mamo-s protein. (B) The DNA recognize sequence logo of Bm-mamo-L protein. (C) According to MEME’s FIMO program package, potential downstream target genes of Bm-mamo protein in the silkworm genome were identified.

The protein function annotation was performed, and found that a group of cuticular protein genes were significantly up-regulation expressed in homozygotes mutant. We detected the expression levels of 20 cuticular protein genes and eight key genes of melanin synthesis were investigated in heterozygous and homozygous gene knockout line of Bm-mamo. Among them, yellow, tan, DDC (DOPA decarboxylase) (fig 7) and 19 cuticular protein genes were significantly up-regulated, and cuticular protein gene CPR11 were significantly down-regulated, in the homozygous Bm-mamo knockout individuals (fig 8).

The melanin metabolism pathway and Fluorescence quantitative PCR of relative genes.

(A) The melanin metabolism pathway. The blue indicates amino acid and catecholamine, the red indicates the name of genes, the black indicates the name of enzymes. (B) Relative expression levels of eight genes in heterozygous type of Bm-mamo knockout group (blue) and homozygous Bm-mamo knockout group (red) as determined by quantitative RT–PCR analysis. The means ± s.d. **P<0.01, paired Student’s t-test.

The expression level of cuticular protein genes were detected by qRT-PCR. The means ± s.d. **P<0.01, paired Student’s t-test.

The up-regulation of gene expression level of yellow, DDC and tan can promote the generation of melanin. Moreover, some cuticular protein gene, for example CPH24, when it is defective can lead to the remarkable decreased of pigmentation 31.

Therefore, the Bm-mamo gene can affect the color pattern of caterpillar by regulating the key genes of melanin metabolism pathway and cuticular protein genes.


Insects are known as one of the most successful populations in evolution, which evolved many important phenotypes in the process of adapting to the environment. Among these, color pattern is one of the most magnetic traits. The biochemical metabolic pathways of pigment, such as melanin, ommochrome, pteridine have been identified in insects12. However, the regulation of pigment metabolism related genes, and the process of transport and deposition of pigment substances is not clear. In this paper, we discover the Bm-mamo gene negative regulate melanin pigmentation in caterpillar. When this gene is deficient, the body color of silkworm is signification melanism, changes in the expression of some melanin synthesis genes and some cuticular protein genes are also significant change.

Deficient of some genes in the melanin synthesis pathway can lead to variations in the color pattern, such as TH, yellow, ebony, tan, aaNAT 3234. However, those often only lead to localized pigmentation changes in later instar caterpillar. For example, the yellow mutant of silkworm, which eye spot, lunar spot, star spot, spiracle plate, head, foot hook and some sclerotization area appears reddish brown, and the other epidermis is consistent with the wild-type35. Why not exhibit a yellow phenotype in other parts of the larval epidermis? The possible reason is that, the expression of yellow is limited to a certain region by transcription regulator, or\and other factors, such as transporter, cuticular protein in the epidermis limit the pigmentation of melanin. Because, the biochemistry synthesized pigment substances need to be transported from epidermis cells and embedded into the cuticle to form coloring. Therefore, the correct location of cuticular protein in epidermis may be one of the important factors of pigmentation. Previous studies have shown that the lack of cuticular protein, CPH24, can lead to defects in pigmentation which supports this argument31.

How does cuticular protein affect pigmentation? Previous research finding that some cuticular proteins can form “pore canals” to transport some macromolecules36. In addition, some cuticular proteins can be crosslinked with catecholamines which synthesized by melanin metabolism pathway37. Because of no live cell in cuticle, melanin precursor substances may be transported by pore canals, formed by some cuticular proteins, and fixed to specific positions through cross-linking with some cuticular proteins. The cuticular protein, CPR4, high expression in black individuals of Bm-mamo knockout line, is required for formation of pore canals in Tribolium castaneum 38. In addition, the vertical pore canal is lacking in less pigmented cuticles of T. castaneum 39. This suggestion that the pore canals constructed by CPR4 may transport the pigments and contribute in epidermis pigmentation. Consequently, to maintain the accuracy of the color pattern, the localization of cuticular proteins in the epidermis may be very important. In one study of microarray analysis, it was found that there are different kinds of cuticular proteins in different colored areas of epidermis, in papilio xuthus 40. These ten cuticular protein genes were conserved in different insects (table s4). In addition, ten cuticular protein genes, which are highly expressed in the black stripes of papilio xuthus larvae, are also significantly up-regulated in Bm-mamo knockout lines of silkworm (fig s8). Therefore, we believe that the location of cuticular proteins and architecture of epidermis, the transport and cross-linking of pigments in cuticle, is the important bases for the formation of color patterns. Hence, the cuticular proteins and pigments metabolism are close cooperated by the uniformly regulated of transcription factors, such as Bm-mamo, thus forming color patterns in caterpillars.

In D.melanogaster, the mamo is required for functional gamete production. In silkworm, the Bm-mamo has conservative function in the female reproduction. Furthermore, the Bm-mamo play an important role in regulation of color pattern of caterpillar. It is generally believed that the main reason for a conserved transcription factor to gain new functions is the evolution of cis-regulatory elements (CREs) of target genes. Because CREs have important function in spatiotemporal expression pattern regulation of genes, and large number of them are located in the non-coding region, which are prone to sequence changes compared with the sequences of coding genes41. In addition, transcription factors (TFs), because of them recognize relatively short sequences, generally between 4 and 20 bp in length, can have a large number of binding sites on the genome42. Therefore, one member of TFs has the potential to regulate almost all genes in one genome. For example, according to the PWM sequence of mamo protein recognize sequence, 10622 genes (∼73% of the total number of genes) contain the binding sites in 2 kb of their up/down stream region, in silkworm. If the core recognition sequence (GCGTG) of mamo protein is used to compare these genome regions, there were 13567 genes (∼93% of the total number of genes) contain the binding sites in 2 kb of their up/down stream region (table s5). And, such binding sites also exist in intron or/and beyond the 2 kb of their up/down stream region of the remaining 1056 genes. The shorter the recognition sequence, the more widespread it can be in the genome. However, in fact, there are not as many genes whose expression is significantly affected by the change of one TF as expected. Therefore, there are some transcription factors binding sites in a gene, but this gene not necessarily regulated by these factors. The possible reason is gene regulation network restricts the target gene set of regulatory factors by such as, whether the target chromatin is open or not43, interplay of TFs44, and the participation of cofactors45.

We have noticed that a lepidopteran caterpillar has multiple visual factors in its skin, and even when its compound eyes are covered, it can rely on the skin to feel the color of the environment, so as to change its body color to adapt to the environment. This visual gene is common in many Lepidoptera46. In the industrial melanism event, in a short period of several decades, the body colors of multiple of insects have changed significantly47. If we only rely on random mutations in the genome, it cannot be so swift.

We speculate that under the stimulation of the natural environment, caterpillars can distinguish the light wave frequency of the background, so that this signal can guide the epigenetic factors and transcription factors, turn on or off the expression of corresponding gene, such as pigment synthesis gene and cuticular protein genes, thus forming a body color pattern similar to the background. In addition, the BTB-ZF family member can participate in chromatin remodeling, and their zinc finger motifs can recognize DNA sequences48. During the opening of chromatin, the probability of mutation of exposed DNA sequence will increase49. The transposon insertion occurs in a timely manner upstream of the cortex gene in the melanism Biston betularia50, may be caused by the similar binding of transcription factors and opening of chromatin.

Therefore, we speculate that evolution of color pattern may be programmed events of the system, rather than relying on random mutations. Based on this assumption, animals can rapidly undergo phenotypic changes under the influence of the environment, and thus the speed of species evolution may be much faster than previous cognition.


Silkworm strains

The bd and bdf mutant strains, wild-type strain Dazao and N4 were from the bank of genetic resources of Southwest University. Silkworms were reared on mulberry leaves or artificial diets at 25°C and 73% relative humidity, dark 12 hours and light 12 hours.

Positional cloning of the bd locus

For mapping of the bd locus, the F1 heterozygous individuals were obtained from a cross between a bdf strain and a Dazao strain. Then the female of F1 was crossed with bdf male (BC1F), and the male of F1 was backcrossed with a bdf female (BC1M). A total of 1,162 BC1M individuals were used for recombination analysis. Genomic DNA was extracted from the parents, Dazao, bdf and F1, and each BC1 individuals using phenol chloroform extracting method. Search for available DNA molecular markers through polymorphism screening and linkage analysis. The primers used for mapping is listed in supplementary Table s6.

Phylogenetic analysis

To determine whether Bm-mamo orthologues existed in other species, we performed the Blastp program on the NCBI nr database ( The sequences of multiple species were executed a multiple sequence alignment of predicted amino-acid sequences by MUSCLE. And then the phylogenetic tree was constructed using the neighbor-joining method with the MEGA7 program (person model). The confidence levels for various phylogenetic lineages were estimated by bootstrap analysis (2,000 replicates).

Quantitative RT–PCR

Total integument RNA was isolated from the silkworm bd mutant, bdf, Bm-mamo knockout line and strain Dazao using TRIzol reagent (life), purified by phenol–chloroform extraction and then reverse transcribed with a random primer (N6) using a PrimeScript™ RT reagent Kit (TAKARA), according to manufacturers’ protocol.

qRT-PCR was performed using a CFX96™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA) with a Real-Time PCR System with iTaq Univer SYBR Green Supermix System (Bio-Rad). The cycling parameters were as follows: 95 °C for 3 min followed by 40 cycles of 95 °C for 10 s and annealing for 30 s. The primers used for target genes are listed in supplementary Table s6. Expression levels for each sample were determined with biological replicates, and each sample was analyzed in triplicate. The gene expression levels were normalized against the expression levels of the ribosomal protein L3 (RpL3). The relative expression levels were analyzed using the classical R=2−ΔΔCt method.

siRNA for gene knockdown

siRNAs for Bm-mamo were designed by siDirect program ( The target siRNAs and negative control were synthetized by Tsingke Biotechnology Co., Ltd. The siRNA (5 μl, 1μl/μg) was injected from the abdominal spiracle into haemolymph at the fourth instar day 3 larval stage. Immediately after injection, phosphate buffered saline (PH 7.3) droplets were placed nearby and the 20-voltage pulse was applied 3 times. The phenotype was observed for fifth instar larvae. The left and the right of epidermis were separately dissected from the injected larvae, and the RNA were extract, respectively. Then the cDNA was synthetized and the expression level of gene were detected by quantitative RT–PCR.

sgRNA synthesis and RNP complex assembly

We used CRISPRdirect (, an online software originally, to screen appropriate guide RNA (sgRNA) target sequence. The gRNAs were synthetized by Beijing Genomics Institute. sgRNA templates were transcribed using a T7 polymerase with the RiboMAX™ Large-Scale RNA Production Systems (Promega, Beijing, China) according to the manufacturer’s instructions. RNA transcripts were purified using 3M Sodium Acetate (pH 5.2) and anhydrous ethanol (1:30) precipitation, were washed using 75% ethanol and were eluted in Nuclease-Free Water. All injection mixes contained 300 ng/µL Cas9 Nuclease (TAKARA Bio, Dalian, China) and 300 ng/µL purified sgRNA. Before injection, the injection mixtures were incubated for 15 min at 37 °C to reconstitute active ribonucleoproteins (RNPs) directly before use.

Microinjection of embryos

For embryo microinjection, microcapillary needles were manufactured using a PC-10 dual stage glass micropipette puller (Narishige, Tokyo, Japan). Microinjection was performed using Eppendorf’s electronic microinjector TransferMan NK 2 (Eppendorf, Hamburg, Germany) with suitable needles. The female moth of wild type silkworm was mated with male moth for 4 hours and laid eggs. The newly laid eggs were adhered onto a clean glass slide. CRISPR/Cas9-messenger RNP mixtures with volumes of approximately 1 nL were injected into the middle of the eggs, and glue the wound. All injected embryos were allowed to develop in the artificial environment of 25 °C and humidity 80%.

Comparative genomics

The reference genome of Dazao was extract from Silkbase ( The bdf genome from silkworm pan-genome project22.

DNA recognition sequence analysis and Downstream target gene screening

The software ( was used to predict DNA binding site for Cys2His2 zinc finger proteins. Chose the confident ZF domains, which scores higher than 17.7. The prediction model chooses RF regression on B1H. And then the sequence logo and position weight matrices (PWM) of DNA binding site was Obtained. The FIMO (Find Individual Motif Occurences) package of MEME suite program was used to search the binding site in silkworm genome.


The research was supported by the National Natural Science Foundation of China (No. 32002230), the Fundamental Research Funds for the Central Universities in China (No. SWU120024).

Competing interests

The authors declare that they have no competing interests.

Author contributions

FY Dai conceived and designed the experiments. SY Wu and CX Peng performed the study. SY Wu, XL Tong, KP Lu, JW Luo, CL Li, X Ding, YR Lu, XH Duan, CH Zhang, D Tan and H Hu analyzed the data. SY Wu wrote the paper. XL Tong and FY Dai edited and revised the manuscript.