Genetic characterization of Bombyx mori (Lepidoptera: Bombycidae) breeding and hybrid lines with different geographic origins.

The domesticated silkworm Bombyx mori L. comprises a large number of geographical breeds and hybrid lines. Knowing the genetic structure of those may provide information to improve the conservation of commercial lines by estimating inbreeding over generations and the consequences of excessive use of those lineages. Here, we analyzed the genetic diversity of seven breeds and eight hybrid lines from Eastern Europe and Asia using highly polymorphic microsatellites markers to determine its genetical impact on their use in global breeding programs. No consistent pattern of deviation from Hardy-Weinberg equilibrium was found for most breed and hybrids; and the absence of a linkage disequilibrium also suggests that the strains are in equilibrium. A principal coordinate analysis revealed a clear separation of two silkworm breeds from the rest: one (IBV) originated from India and the other one (RG90) from Romania/Japan. The tendency of the other breeds from different geographic origins to cluster together in a general mix might be due to similar selection pressures (climate and anthropogenic factors) in different geographic locations. Phylogenetic analyses grouped the different silkworm breeds but not the hybrids according to their geographic origin and confirmed the pattern found in the principal coordinate analysis.

The genetic characterization and classification of breeding lines of any domesticated organism is very important for its conservation of the stock material (including local breeds) and for breeding programs creating new, improved hybrid lineages possibly better suited for economical exploitation. In the case of the domesticated silkworm Bombyx mori L., an enormous geographical and genetical diversity exists all over the world, including a large number of traditional genotypes, hybrids, breeding, and commercial lines ( Vijayan et al. 2010 ), whose biodiversity is important for the selection of useful traits ( Staykova 2008 ). After domestication for >5,000 yr ( Ganesh et al. 2012 ), >3,000 silkworm strains can be distinguished ( Nagaraju et al. 2001 ). The major silkworm breeding industries are in China, India, Japan, Thailand, Korea, Italy, Russia, Bulgaria, and Romania ( Nagaraju and Singh 1997 , Goldsmith et al. 2005 , Jingade et al. 2011 ) and four geographical breeds can be identified: Chinese, Japanese, European, and Tropical ( Liu et al. 2010 ). These are characterized by several morphological, physiological, qualitative, and quantitative traits such as origin, voltinism (mono-, bi-, and polyvoltine), larval duration, cocoon shape and color, silk fiber length, number of molts, or resistance to diseases ( Dalirsefat and Mirhoseini 2007 , Kim et al. 2008 ). Voltinism (number of broods per year) has been shown to clearly separate silkworm strains regardless of their country of origin ( Reddy et al. 1999 ). Interestingly, in bivoltine silkworms, temperature fluctuation during development can modify the voltinism of future generations ( Lim et al. 1990 ).

Apart from being economically important for the silk industry and bioreactor production of recombinant proteins ( Tamura et al. 2000 , Tomita et al. 2003 ), the silkworm is also one of the best-studied lepidopteran model systems for insect physiology, immunity, and genetics with properly characterized mutations ( Willis et al. 1995 , Nagaraju et al. 2000 , Ramesha et al. 2010 , Jingade et al. 2011 ). More than 400 mutations have been described in silkworms, and >200 mutants have been placed on linkage maps, which cover 900.2 cM ( Doira et al. 1992 ). Although morphology is widely used to characterize silkworm strains, these traits are often highly variable and environmentally dependent. Molecular markers circumvent these problems, because they are not modified by the environment and can be easily distinguished in any developmental stage ( Gaviria et al. 2006 ). Genetic analyses give fundamental information to distinguish the origin and evolution of breeds and hybrids of silkworms from Eastern Europe and Asia and may be used to improve maintenance programs and future breeds. Several techniques have been developed for the genetic characterization of strains and high-density mapping of the silkworm genome, such as RAPD ( Nagaraja and Nagaraju 1995 , Promboon et al. 1995 , Yasukochi 1998 ), AFLP ( Tan et al. 2001 ), RFLP ( Goldsmith 1991 ), and simple sequence repeat (SSR; Reddy et al. 1999 , Miao et al. 2005 ). Particularly, microsatellites represent a convenient and certain instrument to produce highly polymorphic molecular markers useful in linkage maps ( Miao et al. 2007 ) and have become one of the most popular genetic markers in population genetics ( Sharma et al. 2007 ). Microsatellites (also known as SSRs) are short tandem repeats of 1–6 bp, which can be found in prokaryotes and eukaryotes ( Hamada et al. 1982 , Weber 1990 , Dietrich et al. 1992 ,  Field and Wills 1996 ), in protein-coding and noncoding regions of the genome ( Toth et al. 2000 , Prasad et al. 2005 ).

This study focuses on the characterization and comparison of 15 Eastern European and Asian silkworm breeds and hybrids using highly polymorphic microsatellites markers. We estimate within and among population variation and show the degree of differentiation of the breeds and hybrids to each other in the light of their origin. The results obtained from this study may be useful to provide efficient biodiversity measurements for the management of germplasm conservation.

Materials and Methods

Silkworm Strains and Breeding Conditions

All monovoltine, diapausing eggs of the seven silkworm breeds and eight hybrids (fixed breeding stocks that were initially produced by crossing two stable parental breeds and then subsequently inbred for several generations with selection for desirable characters) used in this study were obtained from the Romanian native genetic sericulture fund, maintained at SC. Sericarom SA. Bucharest, Romania ( Table 1 ).

Table 1.

Genotypes used to evaluate genetic differences among Bombyx mori L. breeds and hybrids

Name	Geographic origin	Voltinism
RG 90	Romania/Japan	Monovoltine
AC/T	China/Japan	Monovoltine
AB (White Băneasa)	Japan	Monovoltine
IBV	India	Polyvoltine a
AC 29 /T (Chinese white)	China	Monovoltine
B 1 (Băneasa 1)	Japan	Monovoltine
S 8 (Saniş 8)	Japan	Monovoltine
S 8 x AC 29 /T (Ana 1)	Japan/China	Monovoltine
AC 29 /T x S 8 (Ana 2)	China/Japan	Monovoltine
AC x B 1 (Cislău 1)	China/Japan	Monovoltine
B 1 x AC (Cislău 2)	Japan/China	Monovoltine
B 1 x Svila 2	Bulgaria	Monovoltine
B 1 x Hesa 2	Bulgaria	Monovoltine
Hesa 1 x Svila 2	Bulgaria	Monovoltine
Vratza 35 x Svila 2	Bulgaria	Monovoltine

All Bombyx strains were received from the Romanian breeding stock center (SC. Sericarom SA., Bucharest, Romania) and were handled as “Romanian” lines for >10 years. Additionally, the original geographic origin and voltinism type are given.

Became monovoltine by acclimatization and selection in Romania.

B. mori larvae were fed ad libitum with mulberry leaves of Morus alba (Ukraina 107 variety) at University of Agricultural Sciences and Veterinary Medicine Cluj-Napoca, Romania, following standard methodologies for larval rearing ( Mărghitaş et al. 2003 ), which ensures constant microclimate (temperatures between 26 and 27°C in the first two instars, 25–26°C in the third, 23–24°C in the fourth, and 22–23°C in fifth instar; humidity of 80–85% in first instars and in the last ones 65–70%) and nutrition conditions for the whole period of larval stage development for all breeds and hybrids according to the larval instars development.

Genomic DNA Extraction, Microsatellite Amplification, and Analysis

DNA was extracted from 10 individuals of each breed and hybrid line from the posterior silk glands of fifth instar larvae, following a modified Chelex extraction protocol after Walsh et al. (1991) , using 150 µl of 5% Chelex (BioRad Laboratories, Hercules, CA) and 3 µl of Proteinase K (20 mg/ml) (Qiagen, Hilden, Germany) following incubation in a PCR cycler (Bio-Rad C 1000, Munich, Germany): 56°C for 60 min, 95°C for 15 min, 37°C for 1 min, 95°C for 15 min, and 4°C.

Five microsatellite pairs designed from the silkworm genome ( Reddy et al. 1999 , Zhang et al. 2005 ) were used in this study ( Table 2 ). Each polymerase chain reaction (10 µl) contained 2 µl 5× MyTaq reaction buffer (including dNTPs and MgCl 2 ), 0.30 µM of each forward and reverse primer, 0.75 U My Taq HS DNA Polymerase (Bioline, Luckenwalde, Germany), and 2 µl of extracted DNA (∼275 ng/µl). Fragment amplification was done in a Bio-Rad C 1000 thermal cycler using fluorescence-labeled primers following the protocol: initial denaturation at 95°C for 1 min; 35 cycles of 95°C for 15s, Tann (°C): 15s, 72°C for 10 s, and final extension for 10 min at 72°C.

Table 2.

Microsatellite loci used in the study by Reddy et al. (1999) and Zhang et al. (2005)

Locus	No. of alleles	Primer label	Range (bp)	T ann (°C)		Sequence (5'–3')
SAT 951	5	D 2	107–122	48	F	ATTGTAACCGATTTGAGAGA
					R	ATTCGCACAATAAGTTCACT
SAT 346	6	Cy 5	139–189	53	F	GAAGACAGAGCGAAGTGGA
					R	ATGGATTCCTGCTGGTAGAT
CA16G03R	10	Cy 5	220–410	55–60	F	ACAGCATCCAGGTCCGTTCC
					R	GCCGAGTAAAGTATTTGCGTCAT
SAT 1423	9	Cy 5.5	130–176	55	F	CTTTCGATCACCGCGTTCTC
					R	CGCTACGAAATACCATTATCTGACA
T01CTA07R	10	Cy 5.5	240–300	55–60	F	GTCAGACCAAATAGCGGAGGAA
					R	TCGCACGCCTTTTGTTTTG

T ann , annealing temperature.

Amplicon lengths were detected using an automated capillary DNA sequencer (Beckman Coulter 8800 Analysis System and CEQ DNA Size Standard Kit 600, Beckman Coulter, Indianapolis, IN), and allele sizes were assigned using the CEQ Fragment Analysis software.

Data Analysis and Statistical Analysis

Effective numbers of alleles ( Ne ), allele frequencies, and RST values, to analyze the degree of genetic differentiation among populations, were calculated using the population genetics software GenAlEx 6.41 ( Peakall and Smouse 2006 ). GenAlEx was also used to detect private alleles and to calculate the fixation index ( F ) to display a potential excess of homozygotes or heterozygotes. Observed ( Ho ), expected heterozygosity ( He ), and the polymorphic information content (PIC) were estimated using CERVUS v. 3.0. ( Kalinowski et al. 2007 ). Tests for deviations from Hardy–Weinberg equilibrium and linkage-disequilibrium for each locus and strain were done using GENEPOP v. 4.0.10 ( Raymond and Rousset 1995 , Rousset 2008 ).

To compare if some of the analyzed populations are more related to each other than others, the genetic distances were estimated among breeds and hybrids using Nei’s distance ( Nei 1972 ) computed with GenAlEx 6.41. Furthermore, we analyzed population substructuring using an analysis of molecular variance (AMOVA) and gene flow ( Nm) implemented in GenAlEx 6.41. AMOVA calculations were based on RST values ( Slatkin 1995 ). A principal coordinate analysis (PCA) of codominant data, implemented in GenAlEx 6.41, was performed to visualize the differentiation of the breeds in two dimensions.

On the basis of genotypic genetic distances for codominant data, as calculated in GenAlEx ( Peakall and Smouse 2006 ), for breeds and hybrids, a highly resolving cluster analysis was performed using Past v. 2.04 ( Hammer et al. 2001 ). Furthermore, Nei’s population distances from GenAlEx 6.41 were used to analyze between population relationships using the program DendroUPGMA ( http://genomes.urv.cat/UPGMA ) and to display similarity among the breeding strains. UPGMA trees are often used to construct phylogenetic trees if the rates of evolution are approximately constant among the different lineages, so that an approximate linear relation exists between evolutionary distance and divergence time ( Nei 1975 ).

A Mann–Whitney U -test was used to compare the within group genetic distances of breeds and hybrids in STATISTICA 8.0 (StatSoft, Tulsa, OK).

Result and Discussion

Genetic Indices

The set of five highly polymorphic microsatellites generated 31 alleles in total among the 15 silkworm breeds and hybrids, with an average number of 6.20 alleles per locus ( Table 3 ). The mean percentage of polymorphic loci over all populations of 82.67% (data not shown), and the average PIC was >0.50 ( Table 3 ). This demonstrates that all microsatellites markers were highly polymorphic, informative, and useful for genetic diversity and conservation studies. The PIC ranged between 0.352 and 0.671 with an average of 0.558 per locus, indicating a high discrimination power for those loci ( Table 3 ). Usually, the PIC value of a locus decreases if the locus has few alleles, which was the case for the locus SAT 346 with only four alleles showing the lowest PIC value (0.352). An exception was seen for locus SAT 1423, even though this locus had the most alleles, the PIC value (0.598) was not the highest. However, both the average and the range of alleles for the breeding and hybrid lines tested in this study were much lower than estimates of others. Several studies showed ranges between 0.299 and 0.919 and an approximate average of 0.717 irrespective of domesticated or wild silkworm lines ( Reddy et al. 1999 , Li et al. 2005 , Zhang et al. 2005 , Hou et al. 2007 ). The absence of any sign of linkage disequilibrium ( Table 4 ) suggests that all five loci segregate independently of each other. Also we found no consistent pattern of deviation from Hardy–Weinberg equilibrium for within-strain analyses (for each breed and hybrid; data not shown and Supp Material 1 [online only]).

Table 3.

Genetic characteristics of the five Bombyx mori microsatellite loci used in this study

Locus	K	N	N e	H O	H E	PIC	F
SAT 951	6	131	2.077	0.155	0.424	0.671	0.557
SAT 1423	8	128	2.162	0.512	0.503	0.598	0.039
TOICTA07R	6	91	2.249	0.390	0.522	0.671	0.198
SAT 346	4	88	1.222	0.101	0.173	0.352	0.267
CA16G03R	7	144	1.792	0.409	0.352	0.502	−0.205
Over all loci	6.2 ± 1.48	116.4 ± 25.3	1.901 ± 0.41	0.313 ± 0.17	0.394 ± 0.14	0.558 ± 0.13	0.171 ± 0.28

K, number of alleles generated by each locus; N, number of individuals; Ne , effective number of alleles; HO , observed heterozygosity; HE , expected heterozygosity; PIC, polymorphic information content; F , fixation index.

Table 4.

Results of the statistical analysis for deviation from Hardy–Weinberg equilibrium (overall loci) and linkage disequilibrium (locus vs. locus), for all silkworm breeds and hybrids

	χ 2	df	P
Hardy–Weinberg equilibrium
SAT 951	1,096,033	260	0.000
SAT 1423	548,057	300	0.003
TOICTA07R	577,969	280	0.000
SAT 346	312,621	80	0.000
CA16G03R	499,021	200	0.000
Linkage disequilibrium
SAT951/SAT1423	1,262,964	22	0.942
SAT951/TOICTA07R	1,800,612	20	0.587
SAT 1423/TOICTA07R	1,176,722	24	0.982
SAT951/SAT346	5,426,226	6	0.490
SAT1423/SAT346	9,630,808	10	0.473
TOICTA07R/SAT346	7,839,478	6	0.250
SAT951/CA16G03R	2,356,224	20	0.262
SAT1423/CA16G03R	3,010,396	18	0.036
TOICTA07R/CA16G03R	1,944,854	14	0.148
SAT346/CA16G03R	1,133,609	8	0.183

Significant P values are in bold.

The observed heterozygosity varied from 0.101 (for SAT 346) to 0.512 (for SAT 1423) and the expected heterozygosity for each locus from 0.173 (SAT 346) to 0.522 (TOICTA07R) ( Table 3 ). Evaluation of the heterozygosity is important to confirm good management and conservation of the genetic resources ( Graner et al. 2004 ).

In the across-population analysis, a highly significant deviation from Hardy–Weinberg equilibrium was observed for each locus ( P  < 0.001, Table 4 ). This result indicates that, as expected, there is a significant subpopulation structure across all silkworm breed and hybrids, instead of one single panmictic population. The heterozygote deficiency associated with this deviation is confirmed by a generally lower observed than expected heterozygosity and a fixation index across all populations of 0.171 ± 0.28 ( Table 3 ). Karlsson and Mork (2005) explained these heterozygosity deficiencies by invoking inbreeding, assortative mating, differential genotypic selection, null alleles, or the Wahlund effect ( Hamilton 2009 ). In our case, the Wahlund effect is most likely causing the observed pattern.

The presence of private alleles can be considered as a measurement of genetic distinctness from others population. Three of five private alleles were found at the locus CA16G03 in the breed IBV (244 bp, 246 bp, and 248 bp), one allele (161 bp) in SAT1423 for the hybrid B 1 x Svila 2 , and one allele (244 bp) in TOICTA07R for the hybrid Vratza 35 x Svila 2 . The existence of private alleles can be applied in conservation genetic cases and has proven useful in the examination of population structure and migration patterns ( Szpiech et al. 2008 ). The majority of private alleles were found in the breeding line IBV from India, this means that this population is genetically different from others.

Differentiation of Breeds and Hybrid Lines

Comparing Nei’s genetic distances between breeds and hybrids ( Supp Material 2 [online only]), respectively, revealed significantly lower distances between  B. mori  hybrids (0.384  ± 0.05) when compared with the within-group genetic distances of the  B. mori  breeds (1.025 ± 0.21) (Mann–Whitney U -test, Z = 3.404214,  P  = 0.0007, Fig. 1 ). The highest values of genetic distance (3.558) were obtained between AB and IBV breeds (Japan and India). The IBV breed (India) also has a greater genetic distance to all other lineages, which is possibly due to its polyvoltine origin. A comparative study between polyvoltine and bivoltine silkworm lines estimated a Euclidian genetic distance of 0.256 between both, using RAPD markers ( Srivastava et al. 2005 ). The lowest genetic distance of 0.018 was recorded between AB and RG 90 breeds (Japan and Romania/Japan) and might be explained as a result of common genetic origin of both breeds from Japan and Romania/Japan ( Supp Material 2 [online only]). In this study, the high genetic variability observed among all tested monovoltine breeding and hybrid lines to the only one breed with polyvoltine origin indicate the potentiality of using this silkworm stock for heterosis breeding.

Fig.

Box plot of Nei’s genetic distance demonstrating much lower within group genetic distances of breeds vs. hybrids ( ntotal  = 150; n  = 10 for each silkworm lineage; MWU test, Z  = 3.404214,  P  = 0.0007; dot: mean; box: standard error; bars: standard deviation).

AMOVA is a powerful tool to analyze genetic variation within and among populations. AMOVA calculations for our microsatellite data were based on RST ( Slatkin 1995 ) and use the stepwise mutation model. If we consider each breed and hybrid line as a population, the AMOVA revealed a large proportion of genetic variation of 20% among all populations; 80% of the variation resided among individuals and 0% within individuals (allelic variation of each individual is unique to the individuals of the same population, given by their phenotype) ( Table 5 ).

Table 5.

AMOVA- RST of the 15 silkworm breeds and hybrids

Source	df	SS	MS	Est. Var.	% Total
Among populations	14	2,196,242.187	156,874.442	5,563.910	20
Among individuals	135	6,155,493.200	45,596.246	22,762.010	80
Within individuals	150	10,834.000	72,227	72,227	0
Total	299	8,362,569.387		28,398.146	100

SS, sum of squares; MS, expected mean squares; Est. Var, estimated variance.

The low genetic variation within individuals revealed by the AMOVA (0%) might be explained as a result that individuals in the populations are homogenous in nature ( Liu et al. 2010 ). Most of the variation is found within the breeds and hybrids (80%) but also among them (20%). This indicates that although some genetic variation within the breeding lines has been maintained, they are moderately distinct. The breeds IBV from India and RG 90 from Romania/Japan are possibly driving this effect as suggested by the PCA.

Gene flow plays a critical role in differentiation between populations, generally Nm  > 1 is sufficient to overcome the effect of the genetic drift. The studies by Slatkin (1981) and Liu et al. (2010) described values of Nm  < 1 as indicator that differentiation is possible to appear between populations due to the effects of genetic drift. An N value of gene flow with 1.026 among the silkworm strains analyzed in this study was above 1, meaning that populations were large enough to prevent intersubpopulation genetic differentiation caused by genetic drift.

PCA analyzes a set of multivariate data in terms of the major axes of variation. Figure 2 shows the PCA results for all five loci. The PCA was used in addition to a cluster analysis because this multivariate approach is more informative on distances between major groups than cluster analysis ( Hauser and Crovello 1982 ). The first two axes explain 57.38% of the total variation (30.66% and 26.72%, respectively). The plot shows a tendency of breeds from different geographic origins to cluster together without any separation except for two silkworm breeds: one IBV originated from India and the other one RG 90 (Romania/Japan). Because of the distinctness of these two strains (IBV and RG 90 ), they are expected to be the most useful sources to avoid inbreeding depression in future crossings. Kumaresan et al. (2007) and Nezhad et al. (2010) found individuals of strains of the same origin in different clusters of the trees indicating the presence of a great genetic diversity among these strains.

Fig.

Principal coordinates analysis of genetic distances of the seven breeding lines ( ntotal  = 70, n  = 10 for each line). The circles mark the polyvoltine individuals IBV from India and monovoltine RG 90 from Romania/Japan.

A phylogenetic clustering analysis of allele frequencies of the breeding lines revealed a clustering of those seven silkworm breeds into two groups in accordance to their geographical origin. The first group contains silkworms of the polyvoltine origin, the breed IBV from India, and few individuals from other breeds ( Supp Material 3 [online only]). This Indian breed became monovoltine by acclimatization and selection in Romania but obviously retains some genetic characteristics of polyvoltine strains. The second group splits in two weak clusters as a mix of three silkworm breeds from China, Japan, and Romania. The monovoltine silkworm breeds RG 90 with origin from Romania and Japan and AB from Japan is weakly separated from China–Japan group due to their origin. This grouping pattern of the seven silkworm breeds confirmed the pattern found in the PCA ( Fig. 2 ), i.e., the separation of the breeds is moderate, with very few exceptions. The genetic relationships among different silkworm strains determined by inter-SSR (ISSR) amplification also showed clustering in monovoltine, bivoltine, and polyvoltine strains ( Reddy et al. 1999 , Li et al. 2007 ).

According to the genetic distances of the breeding lines, a phylogenetic tree was constructed for the hybrids ( Supp Material 4 [online only]). The majority of the silkworm hybrids lines of the same origin (China/Japan and Bulgaria) pooled together and evidenced a close relationship of these strains may stem from geographic isolation. Similar genotypes in mixed populations might be explained by earlier exchanges of genetic material between the countries, ( Chatterjee and Datta 1992 ). To analyze relationship between breeding and hybrid lines on population level, an UPGMA tree was constructed using Nei’s genetic distance. The distribution pattern is shown in Fig. 3 . Again IBV is quite different from the other breeds and hybrids analyzed, and AB and RG 90 are closely related.

Fig.

UPGMA tree using Nei’s genetic distance to analyze relationship between breeding and hybrid lines on the population level.

In conclusion, in this study, five highly polymorphic microsatellites markers were successfully used to investigate genetic diversity and phylogenetic relationship among 15 silkworm breed and hybrid lineages that represent diverse genetic sources with different geographic origins. Microsatellites still present a state-of-the-art molecular tool to study silkworm lineages with respect to diversity and conservation and can be used for future studies to preserve the genetic material for breeding and research. The analyses denote two distinct breeds, one from India, which became monovoltine by acclimatization and selection in Romania, and RG 90 with origin from Romania and Japan, which may be particularly helpful in the future for improvements of silk productivity. The information gained will be useful for their conservation and sustainable use of genetic manipulations and systematic breeding. To maintain the original characters of the silkworm genotypes over years and to know which silkworm stocks are best suited for long-term maintenance, it is important to preserve the genetic variability and the integrity of the breed and hybrid lineages. Practical application of those resources mainly depends on the known information, so breeder and researcher will select desirable genotypes required for selected breeding and other research programs.