香京julia种子在线播放

The wild silkworm with climbing ability and mimicry response and the domestic silkworm P50 strain loss of these two characters were collected from a wild field in Zhejiang Province and maintained as indoor populations. The wild female moth with larger body and weaker flying ability was crossed with the domestic moth P50 (male) to generate the F1 generation. Then, the F1 (male) was backcrossed with the domestic moth P50 strain (female) to generate the BC1 mapping population. Climbing and mimicry behaviors were carefully observed in both F1 and BC1 generations. All the silkworms were fed at 27°C and 70% humidity.

To examine the climbing ability, the BC1 individuals (∼200) were fed in boxes (with insufficient mulberry leaves) under mulberry seedlings at the beginning of the 5th instar for ∼7 days. Individuals with strong climbing ability could climb up ∼50 cm to the mulberry seedlings for foraging. Those climbed up to the branches were put back to the box for repeated observation (at least 3 times/day). Individuals that consistently climb or stay in boxes were considered as individuals with extreme distinct phenotypes in the segregation population and selected for bulked sequencing. We also observed individuals with white or dark body in both segregation groups. Therefore, four independent sets of different larval body colors were generated for re-sequencing, i.e., white body with climbing ability, white body loss of climbing ability, dark body with climbing ability, and dark body loss of climbing ability, 20 individuals for each set.

To estimate the mimicry traits, the BC1 individuals (∼200) were fed on mulberry branches. Then, the branches were shaken artificially to observe the mimicry response. Those extended their head and chest to simulate the dendrite shape and maintained for at least 10 s were considered as extreme mimic individuals (Supplementary Video 1). Those showed mimicry response were also verified by hand touch (Supplementary Video 2). On the contrary, individuals that kept eating even when touched by hand were considered as extreme non-mimic individuals (Supplementary Video 3). Eventually, 18 individuals with strong mimicry response and 20 non-mimicry individuals were selected for re-sequencing, respectively.

The total DNA was extracted for each individual using a traditional phenol–chloroform DNA extraction protocol and quantified by NanoDrop 2000 (Thermo). Equal quantity (1.5 μg) DNA of each individual in each set was pooled to obtain bulked DNA for sequencing on the Illumina Hiseq 4000 platform (PE150, with insert size around 350 bp). The sequencing libraries were constructed according to the manufacturer’s instructions. In addition, both the wild silkworm (P_Wild) and the domestic silkworm (P_P50) were sequenced for comparison.

The reference genome of domestic silkworm was downloaded from the Silkbase¹. The re-sequencing short reads of each bulk were aligned to the reference genome with BWA (Burrows–Wheeler Aligner) (-t 4 -k 32 -M -R) (Li and Durbin, 2009). Alignment files were converted to BAM files using the SAMtools software (-bS -t) (Li et al., 2009). Potential PCR duplications were excluded using the rmdup in SAMtools. If multiple read pairs were mapped to identical coordinates, only the pair with the highest mapping quality was retained. Reliable genome-wide single-nucleotide polymorphisms (SNPs) and insertions/deletions (InDels) were identified with the Unified Genotyper function in GATK (–filterExpression “QD < 4.0 ∥ FS > 60.0 ∥ MQ < 40.0,” -G_filter “GQ < 20,” –cluster WindowSize 4) (McKenna et al., 2010). InDels were further filtered with the Variant Filtration parameter (–filter Expression “QD < 4.0 ∥ FS > 200.0 ∥ Read PosRankSum < −20.0 ∥ Inbreeding Coeff < −0.8 ”). ANNOVAR was used to annotate SNP or InDel based on the GFF3 files of the reference genome (Wang et al., 2010). The gene models were derived from Xiang et al. (2018) that integrative assembled with GLEAN based on a combination of coding sequences from SilkDB, homology, and ab initio sets (AUGUSTUS, SNAP, and GENSCAN).

The homozygous SNPs/InDels between the two parents were extracted as markers for index calculation. Only SNPs/InDels covered by at least seven reads in both parents (∼1/3 of the sequencing depth of parents) were considered. The genotype of the P50 strain was used as the reference, and reads supporting the reference genotype or non-reference genotype were counted in each bulk. The SNP index and InDel index were calculated with reads supporting the non-reference genotype divided by the total reads in each sequenced bulk (Takagi et al., 2013). SNPs/InDels with depth < 7 in both segregation groups or missing in one group were excluded. The genome-wide SNP/InDel indexes were visualized in sliding window of 1 Mb with a step size of 10 kb. The average index of all the SNPs/InDels presented in each window was designated as the index of this window in each group. The ΔSNP/InDel index was calculated to represent the difference of index between the segregation groups, and 1,000 permutation tests were performed with 95% confidence level as threshold. The SNP/InDel index in each segregation group and the ΔSNP/InDel index were shown with ggplot package in RStudio.

Genes related to domesticated behaviors are usually under selection. Based on the published genome-wide SNPs of 137 domestic silkworm and 7 wild silkworm strains (Xiang et al., 2018), we calculated the selection signatures (Fst and π) of the chromosomes harboring the candidate regions related to climbing and mimicry to screen candidate domesticated genes as described previously (Zhu et al., 2019a). The allele frequency for each site was calculated with the reads supporting the reference genotype divided by the total covered reads. The effects of candidate SNPs were then estimated according to the annotation of the reference genome.

The homologous sequences of candidate genes related to climbing were obtained using blastp against the NCBI database non-redundant nucleotide library² and the sequences of Antheraea yamamai in the GigaDB database³ (Kim et al., 2018). Protein sequences with at least 80% identity in other insect genomes were downloaded for comparisons, including the five Lepidoptera species A. yamamai, Danaus plexippus, Chilo suppressalis, Helicoverpa armigera, and Papilio machaon. The downloaded protein sequences were aligned with ClustalW using default settings in MEGA 7 (Kumar et al., 2016). The minimum-evolution tree was constructed to show the evolution of candidate genes in all the collected insect genomes. Branch model (model = 0 and NSsites = 0 for one-ratio model, model = 2 and NSsites = 2 for two-ratio model) and branch-site model (model = 2, NSsites = 2) in PAML software (version 4.8) were used for phylogenetic analysis (Yang, 2007). The likelihood ratio test was used to estimate the fitness of different models.

RNA-seq data from different tissues of B. mori (5th instar day 3 larvae) were downloaded from the NCBI SRA database⁴. The accession numbers are SRR4425245 (ovary), SRR4425244 (testis), SRR4425250 (brain), SRR4425254 (anterior silk gland), SRR4425258 (middle silk gland), SRR4425260 (posterior silk gland), DRR095110 (midgut), SRR4425248 (fat body), SRR7812745 (integument), and DRR095113 (malpighian tubule). Reads were mapped to the silkworm reference genome with TopHat2 (Kim et al., 2013), and the expression levels (FPKM) were further determined with Resm (Li and Dewey, 2011). The expression heatmaps were constructed with normalized Z-scores on the OmicShare online platform⁵.

We also sequenced the RNA of the brains of both parents at three key larval stages (middle larval stage, late larval stage, and wandering stage) to estimate the expression levels of candidate genes related to climbing and mimicry, three replicates for each sample. Reads were mapped to the reference genome of B. mori with TopHat2 (Kim et al., 2013). The expression levels (FPKM) of genes were further calculated with Cuffdiff (Trapnell et al., 2012). The expression heatmaps were constructed as described above. We also normalized the FPKM of each candidate gene according to the total numbers of sequencing depth of the corresponding samples, and Student’s t-test was imported for estimation of the expression differences between wild silkworm and domestic silkworm.

The wild female moth with larger body and weaker flying ability was crossed with the domestic silkworm P50 strain (male) to generate the F1 generation (Figures 1A–C). Interestingly, all of the individuals in the F1 generation showed a strong ability of climbing and mimicry response, indicating that these two behaviors were dominant. The F1 was then backcrossed with the recessive parent (the domestic silkworm P50 strain) to generate the BC1 populations. Both climbing for foraging and mimicry response were segregated in the BC1 populations (Supplementary Videos 1–3). Individuals in BC1 with stronger climbing ability could climb up to over 50 cm from the foot of the mulberry seedling. We also found that the climbing ability for foraging was not linked to the body color (Figures 1D,E), i.e., individuals with strong climbing ability were detected in both individuals of white and dark bodies. Individuals in BC1 with strong ability to mimic could response to artificial shake and hand touch promptly and maintain for at least 10 s (Figure 1F and Supplementary Videos 1, 2). Those unable to mimic showed no response to artificial shaking and hand touching (Supplementary Video 3).

Because the body color was not linked to the climbing for foraging behavior, we therefore generated four bulks to identify the candidate region related to climbing, i.e., white body with climbing ability, dark body with climbing ability, white body unable to climb, and dark body unable to climb (Figures 1D,E and Supplementary Table S1). In total, 14.7–18.2 Gb data were obtained for the four bulks (Supplementary Table S1). The clean reads were mapped to the reference genome for SNP and InDel calling as described in “Materials and Methods” section. More than 93.34% clean reads could be aligned properly, covering 28–34 × of the reference genome in depth (Supplementary Table S1). The two parents, the wild silkworm (P_Wild) and the domestic silkworm P50 (P_P50), were also sequenced to 10.9 Gb (20.39 × in depth) and 10.0 Gb (21.08 × in depth) for comparison, respectively (Supplementary Table S1). Only SNPs/InDels with homozygous genotypes in both parents were extracted as markers for further analyses, resulting in 7,142,097 SNPs and 346,910 InDels in the white set and 7,063,453 SNPs and 340,067 InDels in the dark set, respectively.

To identify the candidate loci related to climbing, the genotypes of the parent P_P50 were used as reference for calculating the SNP/InDel index in sequenced bulks. If all short reads covering one position are identical to the genotype of P_P50, the SNP/InDel index will be 0. In contrast, if all of the short reads support a different genotype from the reference, the SNP/InDel index will be 1. The ΔSNP/InDel index was also calculated to represent the index differences between the climbing ability bulks and the non-climbing ability bulks. Ideally, the candidate loci related to climbing behavior are expected to be heterozygous, and the ΔSNP/InDel index should be 0.5. The candidate region will be outstanding in the ΔSNP/InDel index due to the strong linkages between closely linked SNPs/InDels with the causal mutation.

In comparison of the climbing and non-climbing bulks of white body, two candidate peaks located on chromosomes 14 and 19 were identified (Figure 2A, Supplementary Figure S1, and Supplementary Data). The peak on chromosome 19 nearly reached the 95% confidence interval. Interestingly, this peak was also detected in the comparison of climbing and non-climbing bulks of dark body (Figure 2B, Supplementary Figure S2, and Supplementary Data), indicating a strong linkage to the climbing behavior. By overlapping the SNPs/InDels with high index in both sets, the candidate region could be roughly narrowed down to 3–5 Mb interval on chromosome 19 (Figures 3A,B).

Considering the possible sequencing bias and errors that occurred in the candidate region, we extended the candidate region on 1–5 Mb on chromosome 19 for the identification of candidate domesticated genes related to climbing. We screened the selective sweep of domestication signatures in this region to discover candidate domesticated genes as described in “Materials and Methods” section (Figure 3C). As expected, strong selective sweep signatures were found in this region, especially in 3.7–4 Mb, in which four candidate domesticated genes were identified (Figures 3A–C). These four candidate domesticated genes were annotated as fanconi-associated nuclease 1-like (FAN1), GPI ethanolamine phosphate transferase 3 (Pigo), ATPase ASNA1 homolog (ASNA1), and uncharacterized protein, respectively (Supplementary Table S2). These four genes showed remarkable higher Fst and extremely lower π in the domestic silkworm than in the wild silkworm (Figures 3C,D). In the Silkworm Genome Informatics Database (SGID) (Zhu et al., 2019b), FAN1, Pigo, and ASNA1 were possibly under positive selection (Supplementary Table S2). Consistently, the allele frequencies in this region were also quite different between the domestic silkworm and the wild silkworm, with nearly fixed alleles in the domestic silkworm population (Figure 4A and Supplementary Figure S3). Interestingly, non-synonymous variations divergent in allele frequency between the domestic and the wild silkworm were identified in the coding regions of FAN1, Pigo, and ASNA1 (Figure 4A and Supplementary Figure S3). To further explore the potential significance of these non-synonymous variations, the homologous protein sequences in 5 Lepidoptera species and 62 other insect species were collected for comparative analysis. We found that only the non-synonymous variation on the ATPase ASNA1 was specific in the domestic silkworm, causing glycine (G) to aspartic acid (D) conversion at residue 115 (Figures 4B,C and Supplementary Figures S4, S5). This variation also appeared to be positively selected in the domestic silkworm clade, as detected in PAML analyses (p < 0.0001 using the branch-site model and p = 0.06 using the branch model) (Figure 4D). Of the 68 collected species, ASNA1 was a single copy gene in B. mori and B. mandarina and possibly two copies in only four species (Supplementary Figure S4). The amino acids affected by the non-synonymous variations identified on FAN1 (G₁₀₅D and D₁₅₅E) and Pigo (K₁₀₈E and I₄₆₃V) were common between B. mori, B. mandarina, and some of the other Lepidoptera species, suggesting that these variations were not correlated with the loss of climbing ability in domestic silkworm (Supplementary Figure S6).

To identify candidate regions related to the mimicry phenotype, individuals with/without responses to artificial shaking or hand touching were collected and pooled separately for re-sequencing. The clean data covered about 32- and 37-fold of the reference genome in depth for mimicry and non-mimicry bulks, respectively (Supplementary Table S1). Compared with the genotypes of the recessive parent P_P50, there were 1,877,029 SNPs and 333,064 InDels retained for BSA. Candidate peaks related to the mimicry trait were detected in 13-17 Mb on chromosome 13, 10–14 Mb on chromosome 21, and 9–10 Mb on chromosome 27, respectively (Figure 5 and Supplementary Figure S7). However, all the peaks were not beyond the 95% confidence interval, suggesting the difficulty of accurate mapping of behavior traits by BSA. Screening of the selective signatures resulted in 12 genes in total in these regions (Figures 5B–D and Supplementary Table S2). These genes were predicted to encode proteins with diverse functions, such as neuropeptide receptors, glycine receptor subunit alpha-4, transporters, kinase related genes, and TBC family member (Supplementary Table S2). Of these candidate genes, only six genes were predicted to be under positive selection in silkworm in SGID (Supplementary Table S2). The allele frequencies in these candidate genes were much higher in domestic silkworm than in wild silkworm (Supplementary Figure S8).

By investigating the published RNA-seq data of the domestic silkworm, we noticed that most of the candidate genes related to climbing for foraging and mimicry were expressed at higher levels in the brain than in the other tissues (Figures 6A,B), especially the mimicry-related candidate genes (Figure 6B). These patterns implied that the loci related to domesticated behaviors might primarily function in the brain. Therefore, we further compared the expression of these genes in the brain of the wild and domestic silkworms at three key larval developmental stages (middle larval stage, late larval stage, and wandering stage) as described in “Materials and Methods” section (Figures 6C,D). The expressions of FAN1 and Pigo were higher in the wild silkworm than in the domestic silkworm but no significant differences in ASNA1 (p = 0.2820 for middle larval stage brain, p = 0.2730 for late larval stage brain, and p = 0.1917 for wandering stage brain, respectively) (Figure 6C and Supplementary Figure S9). We also verified the expression levels of candidate genes related to mimicry in the brains of wild and domestic silkworms. It is interesting that most candidate genes related to mimicry except the three non-detected genes were generally expressed at higher levels in the brain of the wild silkworm (Figure 6D and Supplementary Figure S9).