Mutation Rate and Spectrum of the Silkworm in Normal and Temperature Stress Conditions

Mutation rate is a crucial parameter in evolutionary genetics. However, the mutation rate of most species as well as the extent to which the environment can alter the genome of multicellular organisms remain poorly understood. Here, we used parents–progeny sequencing to investigate the mutation rate and spectrum of the domestic silkworm (Bombyx mori) among normal and two temperature stress conditions (32 °C and 0 °C). The rate of single-nucleotide mutations in the normal temperature rearing condition was 0.41 × 10−8 (95% confidence interval, 0.33 × 10−8–0.49 × 10−8) per site per generation, which was up to 1.5-fold higher than in four previously studied insects. Moreover, the mutation rates of the silkworm under the stresses are significantly higher than in normal conditions. Furthermore, the mutation rate varies less in gene regions under normal and temperature stresses. Together, these findings expand the known diversity of the mutation rate among eukaryotes but also have implications for evolutionary analysis that assumes a constant mutation rate among species and environments.


Introduction
The mutation rate (single-base substitution rate) is a central parameter of evolutionary biology [1]. Many crucial evolution inferences, such as diversity within species and divergence among species, rely on this parameter. The mutation rate of many model species has been explored, including Escherichia coli [2], yeast [3], Caenorhabditis elegans [4][5][6], Drosophila melanogaster [7], mouse [8], and human [9]. However, the number of species in which the rate of mutation has been studied remains limited.
The silkworm, Bombyx mori, is the only completely domesticated insect involved with historical trade and cultural exchange between countries of the East and West (Ancient Silk Road,~2000 years ago), and occupies a special position among domesticated species. Thus, the detailed domestication history of B. mori has long been a concern to silkworm geneticists. Currently, many studies have suggested that B. mori was domesticated from

Sample Source and DNA Extraction
The inbred domesticated silkworm (Bombyx mori) strain Dazao used in this study was obtained from the Silkworm Gene Bank, Southwest University, China. This strain was previously used to produce the B. mori reference genome [38,39]. One normal condition and two temperature-stress-inducing silkworm families were prepared ( Figure 1). The recently laid eggs (aged one hour) were treated with different conditions including a normal condition (25 • C), high temperature (32 • C, 10 h), and low temperature (0 • C, 10 h). B. mori in larval stage is sensitive to temperature, and temperature stresses can easily lead to silkworm death. Therefore, we performed temperature stresses in the egg stage. After treatment, all individuals were reared at 25 • C under 12:12 h (L:D) photoperiod. Both parents and 30 randomly selected offspring (15 males and 15 females) in the normal family were selected to extract genomic DNA. For each stress-inducing family, DNA from both parents and 10 randomly selected offspring were extracted. The whole genomic DNA of individual moths (whole body, including all tissues) was obtained by phenol-chloroform extraction.

Genome Sequencing
DNA library and genome sequencing were performed by Novogene(Beijing Novogene Technology Co., Ltd., Beijing, China) For each moth, a DNA fragment library with a mean insert size of~350 bp was constructed using the TruSeq Library Construction kit (Beijing Novogene Technology Co., Ltd., Beijing, China). Quantity and insert size of each library were checked by Qubit2.0 and Agilent 2100, respectively. Paired-end (2 × 150) sequencing reads with a mean read depth of~30 times for each individual were obtained by Illumina HiSeq 4000 platform.

Identification of the Mutations
Single nucleotide variations (SNVs) were called by standard Genome Analysis Toolkit (GATK) pipeline (version 4.1.2.0) [40,41]. Reads of each individual were mapped onto Genes 2023, 14, 649 3 of 13 the new B. mori reference genome (http://silkbase.ab.a.u-tokyo.ac.jp/cgi-bin/index.cgi (accessed on 12 July 2021) using BWA-MEM algorithm (Burrow-Wheeler Aligner version 0.7.15-r1140) [42]. Sam file was transformed into bam file by SAMtools (version 1.5) [43]. Picard tool (http://broadinstitute.github.io/picard/ (accessed on 15 July 2021)) was used to sort the bam file. GATK was used to mark duplicates and carry out local realignment around indels. We then called SNPs using GATK HaplotypeCaller(version 4.1.8.1). A similar criterion of the previous study was used to reduce false mutations [16]. We disregarded sites: (1) with read depth less than 10 or mapping quality less than 20, (2) marked as low quality by GATK, (3) the proportion of supported reads of mutation base less than 30%, (4) with more than 2 non-reference sites present at same read and (5) the sites of the parents that have an alternate allele or read depths less than 10. Then, the aligned reads near each mutation site were checked using an Integrated Genomics Viewer(version 2.11.9) [44]. Finally, a mutation site only present in one offspring but absent in both parents was defined as a candidate de novo mutation site (DNMs). SAMN13671397, SAMN13671410-SAMN13671421, and SAMN13671422-SAMN13671433, respectively.

Identification of the de Novo Mutations
We sequenced 56 samples (mean depth = 34.2, SD = 7.4) from three full-sib families treated with different environmental conditions including normal conditions as control, high temperature, and low temperature using parents-progeny Illumina sequencing (Figure 1). We obtained a high proportion of mapped reads ( >99.6%) with an average mapping quality of 21.2 (SD = 0.2, including MQ = 0) (Supplementary Table S1). After mutation calling and filtering, we identified 106 candidate mutations from 30 offspring of the normal family. For 10 offspring from each stress treatment, we obtained 48 mutations in the family treated at high temperature and 47 mutations for the low-temperature family (Supplementary Table S2). We found that all of the mutations were heterozygous (Supplementary Table S2), suggesting that these mutations are newly generated in offspring. Further, we analyzed the distribution of the de novo mutations on the different chromosomes and found that the observed distribution was consistent with the expected distribution (assuming the mutations had a random distribution on each chromosome, pvalue > 0.2, K-S test) (Supplementary Figure S1).

Estimation of False Positives
The de novo mutations (DNMs) were verified using PCR amplification and Sanger sequencing. For each DNM, nongenome-amplified genomic DNA of parents and mutation individuals were used as templates for PCR amplification and Sanger sequencing. We verified 105 out of 106 detected DNMs in the normal family. We could not design a proper primer paired with the excluded DNMs. We also randomly selected 20 DNMs in each stress-inducing family for verification.

Estimation of False Negatives
We used a previous synthetic mutation protocol to estimate false negatives [7,16]. For each family, we generated 1000 synthetic point mutations. Each synthetic point mutation was randomly selected from genomic positions and offspring. If the read depths of the selected site were x, we randomly selected y reads to change the reference base in selected points to a randomly non-reference base. The y value was randomly sampled from a binomial distribution with size = x and prob = 0.5. We then used modified reads to substitute original reads in fastq files. Finally, we realigned the changed fastq files to the reference genome and recalled SNPs using the same pipeline and criteria as the above original data.

Statistical Analysis and Silkworm Domestication Time Estimation
R package 4.0.2 was used to perform all statistical analyses in this study. We used previous method (compared the observed numbers of mutants versus non-mutant sites across comparisons to null expectations of discrete uniform mutation distributions, X 2 test) to evaluate the significance of mutation differences among five insects [4]. To evaluate the 95% confidence interval for the mean value of the mutation rate across different treatments, genome regions, and specific bases, we used the bootstrap method with 1000 bootstrap estimates. For example, we randomly resampled 30 values from the rates of offspring of silkworms reared under normal conditions, and then the mean value of the 30 numbers was recorded as one bootstrap estimate. Finally, we sorted the 1000 bootstrap estimates from small to large and used the values at the 25th and 975th positions as the minimum and maximum values of the 95% confidence interval, respectively. Student's t-test was used to evaluate the significance of mutation rate differences of the whole genome, different genomic regions, and specific bases between each stress-inducing family and normal condition family.
To estimate the domestication time of B. mori, the parameter estimates of population mutation parameter (θ = 0.034) and the time (τ = 0.0035) of the wild silkworm and domesticated silkworm diverged were obtained from the previous study [12]. Then, both parameters can be converted into effective population size (Ne) and time in years using the formulation of Ne = θ/4µ and in 4Ne units (divergence time T = 4τNe). The parameter µ represents the mutation rate per site per year (generally, the silkworm is raised 3 generations a year).

Data Access
All sequencing data were deposited to NCBI's short read archive (project ID: PR-JNA597265). The data of normal condition, low-temperature, and high-temperature treatment are accessible with BioSample accession numbers SAMN13671366-SAMN13671397, SAMN13671410-SAMN13671421, and SAMN13671422-SAMN13671433, respectively.

Identification of the de Novo Mutations
We sequenced 56 samples (mean depth = 34.2, SD = 7.4) from three full-sib families treated with different environmental conditions including normal conditions as control, high temperature, and low temperature using parents-progeny Illumina sequencing ( Figure 1). We obtained a high proportion of mapped reads (>99.6%) with an average mapping quality of 21.2 (SD = 0.2, including MQ = 0) (Supplementary Table S1). After mutation calling and filtering, we identified 106 candidate mutations from 30 offspring of the normal family. For 10 offspring from each stress treatment, we obtained 48 mutations in the family treated at high temperature and 47 mutations for the low-temperature family (Supplementary Table S2). We found that all of the mutations were heterozygous (Supplementary Table S2), suggesting that these mutations are newly generated in offspring. Further, we analyzed the distribution of the de novo mutations on the different chromosomes and found that the observed distribution was consistent with the expected distribution (assuming the mutations had a random distribution on each chromosome, p value > 0.2, K-S test) (Supplementary Figure S1).
To estimate the true positive rate of the identified mutations, PCR amplification and Sanger sequencing were used to check the two parents and the mutation offspring. For the 106 de novo mutations identified in the normal family, we verified 84 loci except for 21 loci that gave no PCR product and for one site for which a specific primer could not be designed. We found that 79 out of 84 candidates (positive rate =~94%) were genuine mutations (Supplementary Table S2). We also randomly selected 20 sites for verification in each stress-inducing family. The results showed that all verified 40 candidates were genuine mutations (Supplementary Figure S2 and Supplementary Table S2). To estimate the false negative rate, we randomly generated 1000 simulated mutation sites for each family. We then used the same protocol with the real data to call simulated mutations. We obtained 914 mutations out of 1000 sites in the normal family (false negative rate =~8.6%), 887 sites in the high-temperature-treated family (false negative rate =~11.3%) and 884 sites in low-temperature-treated family (false negative rate =~11.6%).

The Base Mutation Bias
We found that the number of transition mutations (Ts: C↔T and G↔A) was higher than the number of transversion mutations (Tv: C↔A, C↔G, T↔A, and T↔G) in all three families. Assuming the base mutation has no bias, the ratio of Ts/Tv should be 0.5. Yet, for the three silkworm families, Ts/Tv ranged from 1.35 to 1.52 (Supplementary Figure  S3). In addition, the proportion of GC→AT mutation (including C:G→T:A and C:G→A:T mutations) was greater than AT→GC mutation (containing A:T→C:G and A:T→G:C) in the three silkworm families (Supplementary Figure S4A). We found that the GC→AT mutation bias is a common phenomenon in eukaryotic genomes by summarizing previously investigated species (Supplementary Figure S4B). This bias could be the reason why most eukaryotic genomes are rich in A/T bases. Further, the mutation rate of C/G bases is susceptible to temperature. There was no significant difference (p > 0.79, Student's t-test) in the mutation rate of A/T bases among the three families ( Figure 2A, Table 1, Supplementary  Table S3). However, the mutation rates of C/G bases of the two temperature-stress-inducing families were significantly higher (p ≤ 0.05, Student's t-test) than that of a normal condition family ( Figure 2B, Table 1).
the 106 de novo mutations identified in the normal family, we verified 84 loci except for 21 loci that gave no PCR product and for one site for which a specific primer could not be designed. We found that 79 out of 84 candidates (positive rate = ~94%) were genuine mutations (Supplementary Table S2). We also randomly selected 20 sites for verification in each stress-inducing family. The results showed that all verified 40 candidates were genuine mutations (Supplementary Figure S2 and Supplementary Table S2). To estimate the false negative rate, we randomly generated 1000 simulated mutation sites for each family. We then used the same protocol with the real data to call simulated mutations. We obtained 914 mutations out of 1000 sites in the normal family (false negative rate = ~8.6%), 887 sites in the high-temperature-treated family (false negative rate = ~11.3%) and 884 sites in low-temperature-treated family (false negative rate = ~11.6%).

The Base Mutation Bias
We found that the number of transition mutations (Ts: C↔T and G↔A) was higher than the number of transversion mutations (Tv: C↔A, C↔G, T↔A, and T↔G) in all three families. Assuming the base mutation has no bias, the ratio of Ts/Tv should be 0.5. Yet, for the three silkworm families, Ts/Tv ranged from 1.35 to 1.52 (Supplementary Figure S3). In addition, the proportion of GC→AT mutation (including C:G→T:A and C:G→A:T mutations) was greater than AT→GC mutation (containing A:T→C:G and A:T→G:C) in the three silkworm families (Supplementary Figure S4A). We found that the GC→AT mutation bias is a common phenomenon in eukaryotic genomes by summarizing previously investigated species (Supplementary Figure S4B). This bias could be the reason why most eukaryotic genomes are rich in A/T bases. Further, the mutation rate of C/G bases is susceptible to temperature. There was no significant difference (p > 0.79, Student's t-test) in the mutation rate of A/T bases among the three families ( Figure 2A, Table 1, Supplementary Table S3). However, the mutation rates of C/G bases of the two temperature-stressinducing families were significantly higher (p ≤ 0.05, Student's t-test) than that of a normal condition family ( Figure 2B, Table 1).  Mutation rates of A/T bases and C/G bases of the silkworm under different conditions. The mean value and 95% confidence intervals of the mutation rate of (A) the A/T bases and (B) the C/G bases. The 95% confidence interval of the mutation rate was estimated using bootstrap method with 1000 bootstrap estimates. Student's t-test was used to test differences in the mutation rate (observed mutation rate per genome) between each stress-inducing family and normal condition family. The number of G:C→A:T mutations was higher than the other five mutation types (including A:T→C:G, A:T→G:C, A:T→T:A, G:C→T:A, and G:C→C:G) in all three silkworm families ( Figure 3A, Table 2, Supplementary Table S4). Assuming no mutation bias, the number of G:C→A:T mutations constituted about 16.7% of all de novo mutations. Nonetheless, the proportion of the G:C→A:T mutations in each silkworm family was more than 32% ( Figure 3A). This tendency is a common phenomenon in eukaryotes (Supplementary Figure S5). A previous study showed that the C:G→T:A transition bias could be caused by the methylation of CpG islands, which is very common in vertebrates [45]. After oxidative deamination of methylated cytosines, the cytosines are changed into thymine [46][47][48]. However, we found that the C:G→T:A transition bias in B. mori is not related to the methylation of cytosines at CpG islands. The C:G→T:A mutations in B. mori presented no CpG island distribution bias or CG dinucleotide bias ( Figure 3B,C). In addition, a previous study showed that there are fewer methylated cytosines in the B. mori genome and that most methylcytosines occur in CG dinucleotides [49]. Therefore, the C:G→T:A transition bias in B. mori could be caused by a methylation-independent mechanism.

The Mutation Rates under Different Environments
To estimate the mutation rate of the three silkworm families, callable sites and the mutation rate of each offspring of each family were calculated (Supplementary Table  S5). Then, the mutation rate was estimated to be 0.41 × 10 −8 (95% confidence interval = 0.33 × 10 −8 -0.49 × 10 −8 , bootstrap method, see Materials and Methods) in the normal family (Table 1). This value is significantly higher (p = 5.1 × 10 −5 , X 2 test) than the mutation rate observed in four other insect species (0.28 × 10 −8 in Drosophila melanogaster, 0.29 × 10 −8 in Heliconius melpomene, 0.34 × 10 −8 in Apis mellifera and 0.36 × 10 −8 in Bombus terrestris) [7,16,23,24], suggesting that the spontaneous mutation rate of insect genomes are not as constant as previously claimed. Further, we compared the single-base mutation rate of silkworms under normal and stressed conditions and found that the mutation rate of the high-temperature family (0.58 × 10 −8 ) was significantly (p = 0.04, Student's t-test) higher than the mutation rate of the normal condition family ( Figure 4A, Table 1). The mutation rate of the low-temperature family (0.57 × 10 −8 ) was marginally significant (p = 0.06 using Student's t-test; p = 0.04 by Wilcoxon rank-sum test) higher than that of the normal family ( Figure 4A, Table 1). mine [46][47][48]. However, we found that the C:G→T:A transition bias in B. mori is not related to the methylation of cytosines at CpG islands. The C:G→T:A mutations in B. mori presented no CpG island distribution bias or CG dinucleotide bias (Figures 3B,C). In addition, a previous study showed that there are fewer methylated cytosines in the B. mori genome and that most methylcytosines occur in CG dinucleotides [49]. Therefore, the C:G→T:A transition bias in B. mori could be caused by a methylation-independent mechanism.  To investigate the mutation rate in different genome regions, for each offspring of each silkworm family, we calculated the callable sits and the mutation rate in intergenic and gene regions (combining exons and introns and including 5 -and 3 -UTRs) (Supplementary Table S6). We found that the mutation rate of the gene region is constant (P = 0.76, Student's t-test) in normal conditions and temperature stresses ( Figure 4B). The mutation rates of the gene region in the normal condition and high-temperature-and low-temperaturetreated families were 4.52 × 10 −9 (95% CI = 3.23 × 10 −9 -5.92 × 10 −9 ), 4.95 × 10 −9 (95% Genes 2023, 14, 649 8 of 13 CI = 2.65 × 10 −9 -7.23 × 10 −9 ) and 4.92 × 10 −9 (95% CI = 3.94 × 10 −9 -5.90 × 10 −8 ), respectively (Table 1). In contrast, the rate (6.05 × 10 −9 ) of the intergenic region in the hightemperature family was significantly (p = 0.04, Student's t-test) higher than the mutation rate (3.84 × 10 −9 ) of the intergenic region of the normal condition family ( Figure 4C, Table 1). The mutation rate of the intergenic region in the low-temperature family (5.85 × 10 −9 ) was marginally significant (p = 0.06 by Student's t-test; p = 0.02 by Wilcoxon rank-sum test) higher than that of the normal condition family ( Figure 4C, Table 1).
Genes 2023, 14, x FOR PEER REVIEW 8 of 13 Figure 4. The mutation rates in whole genome, gene region, and intergenic regions of the silkworm under different conditions. The mean value and 95% confidence intervals of the mutation rate of (A) the whole genome, (B) gene region, and (C) intergenic region were shown. The 95% confidence interval of the mutation rate was estimated using bootstrap method with 1000 bootstrap estimates. Student's t-test was used to evaluate the difference in the mutation rate between each stress-inducing family and normal condition family, and p-value is shown at the top of the figure. The data for the test of significance are the observed mutation rates (per genome) rather than the bootstrap estimates.

The Mutation Rate Varies in Different Species
In this study, we obtained the spontaneous mutation rate (0.41 × 10 −8 ) of a model insect, Bombyx mori. This is the fifth direct estimate of the single-nucleotide mutation rate in insects. The rate of B. mori is similar to the estimates in Daphnia pulex [50], Canis lupus [51], and Ficedula albicollis [52] and is lower than the rates of humans [18][19][20][21][22]. To the best of our knowledge, the human genome has the highest mutation rate (per site per generation) of any species for which mutation rates have been detected. That could be because Figure 4. The mutation rates in whole genome, gene region, and intergenic regions of the silkworm under different conditions. The mean value and 95% confidence intervals of the mutation rate of (A) the whole genome, (B) gene region, and (C) intergenic region were shown. The 95% confidence interval of the mutation rate was estimated using bootstrap method with 1000 bootstrap estimates. Student's t-test was used to evaluate the difference in the mutation rate between each stress-inducing family and normal condition family, and p-value is shown at the top of the figure. The data for the test of significance are the observed mutation rates (per genome) rather than the bootstrap estimates.

The Mutation Rate Varies in Different Species
In this study, we obtained the spontaneous mutation rate (0.41 × 10 −8 ) of a model insect, Bombyx mori. This is the fifth direct estimate of the single-nucleotide mutation rate in insects. The rate of B. mori is similar to the estimates in Daphnia pulex [50], Canis lupus [51], and Ficedula albicollis [52] and is lower than the rates of humans [18][19][20][21][22]. To the best of our knowledge, the human genome has the highest mutation rate (per site per generation) of any species for which mutation rates have been detected. That could be because the generation time of humans is much longer than other species. If we use the mutation rate per site per year, the mutation rate of silkworms (~3.74 × 10 −8 ) will be much higher than that of humans (~0.048 × 10 −8 , assuming one generation time equals 25 years). In addition, the rate of B. mori is significantly higher (by 1.5-fold, p-value = 5.1 × 10 −5 , X 2 test) than the rate of four other insects as estimated by a parents-offspring sequencing approach [7,16,23,24] (Supplementary Figure S6A), which could be caused by two reasons. First, the genome size of the silkworm is higher than that of the prior four insects. We found a significant positive correlation between the mutation rate and genome size (r = 0.83, p < 0.01) among the species whose mutation rate has been estimated (Supplementary Figure S6B). Another possible explanation for the mutation rate differences among insects is that the generation time varies among the species. For instance, the generation time of the silkworm is~40 days, but~10 days for Drosophila melanogaster. A previous study in yeast indicated that the mutation rate per generation increases with the generation time [3]. That conclusion is partially consistent with the hypothesis of time-dependent mutagenesis [20,53,54].
Our finding suggests that prior estimates of the silkworm domestication history relying on the mutation rate of Drosophila or Heliconius requires revision. Our results suggest that the mutation rate per site per generation in the silkworm is higher (1.4~1.5-fold) than that in Drosophila and Heliconius. Based on the mutation rate estimated in this study and parameter estimates of divergence time from previously study [12], the start time of the silkworm domestication can be traced back to 9700 years ago (detailed inference see Materials and Methods). Nonetheless, we also found that the mutation rate varies in different environments (discussed below). The living environment varies significantly between wild silkworms and domestic silkworms. While the domestic silkworm lives protected, the wild silkworm lives in the wild. Is the mutation rate constant between wild silkworm and domestic silkworm? Moreover, the domesticated silkworm has a multitude of different local strains based on geographical distribution including torrid, temperate and frigid zones. Is the mutation rate constant among these local strains? The answers to these questions will help us to obtain an accurate estimate of the silkworm domestication history.

The Mutation Rate Varies in Different Environments
We compared the mutation rate of B. mori under temperature-stress-inducing environments. Compared with the normal family, the mutation rate was significantly higher in temperature-stress-inducing environments. This implies that the mutation rate can be affected by the environment, although this study does not clearly distinguish between somatic mutations and germ cell mutations. An increase in the mutation rate under stress conditions has also been observed in unicellular organisms. For instance, the mutation rate increased when bacteria experienced antibiotic and other stress treatments [28][29][30][31]. The mutation rate of yeast varies significantly under different nutrition conditions [3]. In yeast, the yearly mutation rate is less variable than the mutation rate per generation because there is a negative correlation between the yeast growth rate and the mutation rate per generation [3]. However, we did not find this trend in B. mori. The generation time (~40 days) did not change in the different temperature-stress-inducing silkworm families. Generally, a low temperature will lead to a longer generation time of B. mori, while a high temperature will speed up development and shorten the generation time [55]. However, long-time exposure to high temperatures (32 • C) or low temperatures (0 • C), especially in the larval stage, can lead to the death of B. mori. In contrast, the silkworm has a relatively strong tolerance to temperature in the egg stage or embryonic stage, so we observed that the generation time of B. mori was not affected by short high (32 • C, 10 h) and low temperature (0 • C, 10 h) treatment in the egg stage. Currently, the mechanism involving the mutation rate variation in different stresses remains unclear, which is an interesting topic for future research.

The Mutation Rate Is Constant in Gene Regions
Our results showed that the mutation rate in gene regions is more constant than the rate in intergenic regions. We found the mutation rate variation (differ by 1.09-fold) of gene regions among the normal and temperature-stress-inducing silkworm family was smaller than the corresponding variations (>1.5-fold) in intergenic regions. This phenomenon may be attributed to replication time-dependent mutagenesis. Generally, gene regions replicate at the early S-phase, and intergenic regions replicate at the late S-phase [56]. Single-strand DNA will appear with a higher probability at the late S-phase due to the slowing down or stalling at the replication with a decreasing dNTP pool and is more vulnerable to environmental damage [53,[57][58][59]. Currently, many studies have used wholegenome sequence information to carry out evolution inference or phylogenomic analyses with sequencing technology development and reduced sequencing costs. However, we found that the mutation rate of the intergenic region that constitutes most of the genome in multicellular eukaryotes can be affected by external temporary temperatures including low temperature (0 • C) and high temperature (32 • C), which occur frequently in nature. This finding implies that using the sequence of gene regions to carry out evolution analysis assuming a constant mutation may give results that are more reliable. This finding may also help explain why there is a positive correlation between the mutation rate and genome size.  Figure S4). The phylogenetic tree was generated by TimeTree (http://www.timetree.org/). (B) Correlation analysis between the mutation rate and genome size was performed using the R program with Pearson's method. Supplementary Table S1. Sequencing depth and coverage of all samples. Supplementary Table S2. The information of all de novo mutations. Supplementary Table S3. The callable sites, de novo mutations, and mutation rate of A/T bases and C/G bases of each genome. Supplementary Table S4. The callable sites, de novo mutations, and mutation rate of each mutation type of each genome. Supplementary Table S5. The callable sites, de novo mutations, and mutation rate of each genome in the silkworms under different conditions. Supplementary Table S6. The callable sites, de novo mutations, and mutation rate in gene region and intergenic regions of each genome.  Data Availability Statement: All sequencing data were deposited to NCBI's short read archive (project ID: PRJNA597265). The data of normal condition, low-temperature and high-temperature treatment are accessible with BioSample accession numbers SAMN13671366-SAMN13671397, SAMN13671410-SAMN13671421 and SAMN13671422-SAMN13671433, respectively.