223 

ISJ 15: 223-233, 2018                                             ISSN 1824-307X 

 
 

RESEARCH REPORT 

 
Composition and diversity analysis of intestinal microbiota in the fifth instar silkworm, 
Bombyx mori L. 

 
C Hou1,2*, Y Shi1, H Wang1, R Li1, M A Nartey1, X Guo1,2 

 
1Jiangsu Key Laboratory of Sericultural Biology and Biotechnology, School of Biotechnology, Jiangsu University of 
Science and Technology, Zhenjiang Jiangsu 212018, China; 
2The Key Laboratory of Silkworm and Mulberry Genetic Improvement, Ministry of Agriculture, Sericultural 
Research Institute, Chinese Academy of Agricultural Science, Zhenjiang Jiangsu 212018, China 

 
 

   Accepted June 12, 2018 

 
Abstract 

Despite the low diversity and high variability at intra- and interspecific levels, gut microbiota of 
insect play important roles in digesting and absorbing nutrients, health and development, even in 
resistance and immunity. In this study, we investigate the composition and diversity of gut microbiota in 
the 5th instar of silkworm larvae since the mount of ingested mulberry leaves of this stage is about 85% 
of the total amount. Intestinal contents were collected at 0, 12, 24, 48 and 72h in the 5th instar. 32, 21, 
29, 24, 14 phyla and 386, 163, 292, 196, 144 genera were detected respectively in the samples 
described above. In general, the dominated phyla in the larvae gut were Firmicutes, Cyanobacteria, 
Proteobacteria, Actinobacteria and Bacteroidetes; the dominated genera were Enterococcus, 
Staphylococcus, Arthrobacter, Pseudomonas, Lactobacillus, Bacteroides, Paenibacillus and Serratia, 
respectively. The range of phyla and genera was different and there were two peaks of microbial 
diversity at 0 and 24 h. Compared with the intestinal bacteria of these time points, there were unique 
bacteria and 56 common ones. In these bacteria, the abundance of Gram-positive bacteria including 
Enterococcus and Staphylococcus were the highest. Our results provided foundation and reference on 
the research of intestinal microorganisms and resistance breeding in silkworm. 

 
Key Words: Silkworm, Gut microbiota, Composition, Pyrosequencing 16S rRNA genes 

 

 
Introduction 

 
The silkworm, Bombyx mori, an important 

economical insects and model organisms of 
lepidopteran, has contributed enormously to the 
study of Lepidopteran genetics and immunology 
(Goldsmith et al., 2005; Hou et al., 2011). The 
silkworm corporeity is closely related to sericulture 
production, the related factors of effected corporeity 
include climatic condition, quality of mulberry leaves, 
silkworm strains, digestion and absorption of 
nutrients. Digestion, absorption and disease 
prevention of silkworm were closely related to the 
microbiota lived in the silkworm digestive tract (Yuan 
et al., 2006). Therefore, understanding the change of 
intestinal bacteria can help us improve the utilization 
of mulberry and health of the silkworm.  
___________________________________________________________________________ 

 
Corresponding author: 

Chengxiang Hou 
The Key Laboratory of Silkworm and 

Mulberry Genetic Improvement 
Ministry of Agriculture 
Jiangsu University of Science and Technology 

Zhenjiang 212018, Jiangsu Province, China 
E-mail: cxhou587@163.com 

 
 
The insects gut inhabits many non-pathogenic 

microorganisms, which is called a special growth 
environment of microorganisms. The microbial 
community structure of insect gut is affected by 
insect species, the morphology of digestive tract, 
feed, feeding habits, insect age, pathological 
conditions and the like (Dillon and Dillon, 2004). 
These intestinal bacteria play important roles in 
metabolism, health, reproduction, immunity and 
pesticide degradation (Maslowski and Mackay, 2011; 
Simon et al., 2011; Lee and Brey, 2013). They are 
also involved in host-parasite interaction, thus 
influence evolutionary and ecological processes of 
insect by affecting host biology (Feldhaar et al., 2011; 
Aptedeshpande et al., 2014). Therefore, it is 
necessary to access the composition and diversity of 
the intestinal bacteria, which can provide some new 
insights in development, immunity, new pesticide 
exploration, even evolution of insect.  

Molecular markers had been widely used in 
analysis of silkworm traits (Li et al., 2007; Xiang et 
al., 2007; Li et al., 2013). 16S rRNA genes are 
important references to identify and analyze 
phylogenetic bacteria, which are common markers 



224 

of bacteria with highly conserved and special 
sequences (Caporaso et al, 2011). The related study 
of intestinal microbiota in silkworm was few. Some 
classical microbiological and molecular techniques 
had been used to investigate the bacterial 
composition of the silkworm gut, such as 
culture-dependent isolation methods and DGGE 
fingerprinting (Yuan et al, 2006; Xiang et al., 2007). 
However, these relatively low-throughput techniques 
cannot provide enough information about the 
microbial composition and diversity in hosts. 
Compared with these methods, the recent rapid 
development of high-throughput sequencing 
techniques provide new methods in analysis of 16S 
rRNA genes, genome sequence and transcriptome 
(Boissière et al., 2012; Hou et al., 2011; 2014; Li et 
al., 2016; Chengxiang et al., 2017; Zhang et al., 
2017). We explored the bacterial compositional 
changes and diversity at different developmental 
times by sequencing 16S rRNA genes to understand 
the host-bacterial relationships in silkworm gut. 
 
Materials and Methods 
 
Silkworm strain 

The silkworm strain Dazao was provided by the 
Sericultural Research Institute, Chinese Academy of 
Agricultural Sciences. The silkworm larvae were 
reared by using fresh mulberry leaves at standard 
temperature and humidity. The newly exuviated 5th 
instar larvae were used for the experiments. 
 
Collection of intestinal contents 

The newly exuviated 5th instar larvae were fed 
with mulberry leaves. The mid guts of 15 larvae were 
dissected in a sterile environment at five time points 
(0, 12, 24, 48 and 72 h after feeding with mulberry 
leaves) respectively. Three repeats of each time 
points were collected. Prior to dissection, the 
dissecting apparatus such as scissors and forceps 
were sterilized by autoclaving and UV treatment. 
The intestinal contents were immediately collected 
into the sterile eppendorf tubes, quickly frozen in 
liquid nitrogen and stored at -80 °C.  
 
DNA extraction, PCR amplification, library 
construction and sequencing 

The genomic DNA of the intestinal contents was 
extracted using SDS method, then the extracted 
DNAs were detected by agarose gel electrophoresis 
and fluorimeter, then the eligible DNA samples were 
diluted to 1 ng/ul. 16S rDNA Amplification 
Sequencing was amplified based on Illumina HiSeq 
sequencing platform and the Paired-End sequencing 
method, using high-fidelity enzyme and specific 
primers with Barcode (515F: 
5ʹ-GTGCCAGCMGCCGCGGTAA-3ʹ, 806R: 
5ʹ-GGACTACHVGGGTWTCTAAT-3ʹ) (Caporaso et 
al., 2011; McHardy et al., 2013). The amplified 
system was 50 μL contained 10 ng DNA, 5 μL 10x 

PCR buffer，0.5 μL 10 mmol/L dNTPs and 5 U/μL 
Plantium Taq, 1 μmol/L 515F and 806R primer 5 μL, 
respectively. Thermocycling conditions were: 94 °C 
for 3 min; 5 cycles of 94 °C for 30 s, 45 °C for 20 s, 
65 °C for 30 s; 72 °C for 90 s; then 20 cycles of 
94 °C for 20 s, 55 °C for 20 s, 72 °C for 30 s; finally 
72 °C for 5 min. The PCR amplification action used 

the Phusion® High-Fidelity PCR Master Mix with GC 
Buffer kit (New England Biolabs, USA). The 
amplification product was tested by 2% agarose gel 
electrophoresis. The aimed amplicons were purified 
by QIAquick Gel Extraction Kit (qiagen, German) 
according to the protocols and their concentrations 
were determined by fluorimeter. The library was 
constructed by TruSeq® DNA PCR-Free Sample 
Preparation Kit (Illumina, USA) and quantified by 
Qubit and Quantitative real-time PCR, then the 
qualified library was sequenced in HiSeq2500 
PE250. Among the sequenced sample libraries, 
there were 3 repeats for each sample of 0, 12, 24, 48, 
and 72 h in the 5th star. 
 
Acquisition of effective tags 

Raw tags of samples were obtained from 
machine data by splitting Barcode sequence and 
PCR amplification primers, then splicing by Flash 
(V1.2.7, http://ccb.jhu.edu/software/FLASH/). Raw 
tags must be conformed to the following procedure 
by strictly filtrating and initially getting the 
high-quality data (clean tags): 1) tags interception: 
the first low-quality base site was truncated in raw 
tags from continuous low quality value (the 
acquiescent quality threshold was ≤ 19 base) to the 
setting length (the acquiescent length was 3 base). 2) 
length filtration of tags: after interception, further filter 
tags in the continuous high quality base with length 
was less than 75% of tags length. 3) acquisition of 
effective tags: the remaining tags after interception 
and filtration need compare and detected chimera 
sequences in the gold database 
(http://drive5.com/uchime/uchime_download.html) 
and remove them 
(http://www.drive5.com/usearch/manual/chimera_for
mation.html), finally obtained the effective tags. 
 
Operational taxonomic units (OTUs) clustering, 
species annotation and diversity analysis 

All the effective tags of samples were clustered 
into OTUs by uparse software (Uparse v7.0.1001, 
http://drive5.com/uparse/) and the threshold of 
acquiescent identity was 97%. The representative 
sequence of OTUs were filtered according to their 
emerging frequencies and analyzed their species 
annotation of these sequence by Mothur method and 
the SSUrRNA database of SILVER 
(http://www.arb-silva.de/), then the information of 
species annotation (the acquiescent threshold was 
0.8-1) and taxonomy were obtained, finally the 
community composition of sample was counted in 
each classification level (kingdom, phylum, class, 
order, family, genus, species). The multiple 
sequence alignment were confirmed by muscle 
software (Version 3.8.31, 
http://www.drive5.com/muscle/) and the systematic 
relationships of OTUs representative sequence were 
obtained. Finally, all the samples data were treated 
by homogenization procedure, the treatment was 
carried out using the least amount in the data as 
standard, the followed Alpha and Beta diversity 
analyses, cluster and phylogenetic analysis were 
based on the normalized data. 

Diversity analysis included Alpha and Beta 
diversity analysis, alpha diversity analysis 
(http://scikit-bio.org/docs/latest/generated/skbio.dive

http://ccb.jhu.edu/software/FLASH/
http://drive5.com/uchime/uchime_download.html
http://www.drive5.com/usearch/manual/chimera_formation.html
http://www.drive5.com/usearch/manual/chimera_formation.html
http://drive5.com/uparse/
http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.html


225 

rsity.alpha.html) was performed to calculate index of 
ACE, Chao1, Shannon and Simpson using Qiime 
software (version 1.7.0), draw rarefaction curves, 
rank abundance and species accumulation curves 
by R software (version 2.15.3). ACE and Chao1 
were used to calculate the bacterial community 
richness, while Shannon and Simpson were used to 
calculate the bacterial community diversity and 
estimate the total number of species. The higher 
index value of ACE, Chao1, Shannon and Simpson 
meant the more community number and diversity of 
microbiota. The bacterial diversity was evaluated by 
rarefaction curves, rank abundance and species 
accumulation curves. Rarefaction curve was a 
representative value of alpha diversity, it was used 
to evaluate whether the sequencing was sufficient to 
cover all taxa or not, which could reflect the 
abundance of species in the samples. Venn 
diagrams were drawn with the online Venn tool 
(http://bioinformatics.psb.ugent.Be/webtools/Venn/), 
and used to represent their intersecting relations. 

Beta diversity analysis was performed to get the 
maps of the phylogenetic tree, Nonmetric 
Multidimensional Scaling (NMDS), principal 
components analysis (PCA) and principal coordinate 
analysis (PCoA). The UniFrac distance was 
calculated by Qiime software (version 1.7.0) and 
analyzed the molecular variance to directly compare 
the structure differences of bacterial community 
between two samples (Schloss, 2008), then the 
phylogenetic tree was constructed by the 
neighbor-joining method of MEGA6.0. The maps of 
PCA and PCoA were drawn by R software (version 
2.15.3). PCoA was used to compare bacterial 
structures based on weighted - UniFrac from each 
library. PCA was utilized to explore the correlation 
between dominant genera. In the maps of PCA and 
PCoA, the difference of community structure was 
determined according to the condition of aggregation 
and separation between samples. NMDS is 
commonly used to compare the differences among 
the samples and evolutionary information. The 
difference of community composition was small if the 
distance between two points was close. 
 

Phylogenetic analysis 
Phylogenetic analysis was performed to get the 

evolutionary and similar relationships. The 
evolutionary relationships of genera Top100 were 
acquired by multiple sequence alignment, the 
relative abundance of these genera was combined 
and integrated, finally the map of phylogenetic 
relationships were gained by GraPhlAn software 
(http://huttenhower.sph.harvard.edu/GraPhlAn) 
(Asnicar et al., 2015). The similar relationships of 
different samples, Unweighted Pair-group Method 

with Arithmetic Mean（UPGMA） was used to cluster 
and analyze samples, construct phylogenetic tree of 
different samples, finally solve classification 
problems. 

 
Results 
 
Analysis of the sequencing data 

After a series of data management, 651,441 
effective tags were obtained from the samples by 
Illumina HiSeq Sequencing platform. 99.43% 
(647,707) of effective tags was assigned to Taxon 
and 3,121 OTUs were obtained (Table 1). There 
were 1200 OTUs in the PC0 group, and 394, 823, 
484 and 220 OTUs in PC1, PC2, PC3 and PC4 
groups, respectively. Among them, there were 109 
shared OTUs (Support Fig 1). After quality controlled, 
the length of total tags were 18,220,501 bp, length 
ranging from 191 to 381 bp and the average length 
was 251.1 bp. 

The index values of alpha diversity were 
calculated by Qiime software (version 1.7.0) (Table 
2). In the index values of ACE and Chao1 for 
calculating bacterial abundance, the values of PC0, 
PC1, PC2, PC3 and PC4 group were gradually 
decreased in the index of acachao1 for calculating 
bacterial abundance, which indicated that their 
bacterial abundance was also decreasing. This 
result was consistent with the data in Table 1. In the 
analysis of alpha diversity, the rarefaction curve 
showed that bacterial richness in the gut contents 
was different in 5 samples and they approached to 

 

 
 
 

Table 1 Sequence data and bacteria in the contents of the healthy silkworm gut 
 

Sample Name Raw PE Effective tags Taxon_Tag OTUs Phylum Class Order Family Genus 

PC0 147,772 139,041 137,271 1200 31 58 120 209 368 

PC1 150,265 146,121 145,682 394 21 35 78 116 163 

PC2 47,301 44,601 43,772 823 29 46 105 169 292 

PC3 161,610 156,628 156,070 484 24 43 87 141 196 

PC4 168,714 165,050 164,912 220 14 28 55 81 144 

total 675,662 651,441 647,707 3,121 119 210 445 716 1,163 

 
The samples of PC0, PC1, PC2, PC3 and PC4 groups were collected from the mid gut contents of five time points 
in the 5th instar silkworm. The five times points were 0, 12, 24, 48 and 72 h after feed with mulberry leaves 

http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.html
http://huttenhower.sph.harvard.edu/GraPhlAn


226 

Table 2 The index values of richness and diversity in the samples 
 

Sample name Shannon Simpson Chao1 ACE Goods_coverage 

PC0 3.771 0.729 1028.908 1021.058 0.995 

PC1 1.016 0.241 383.591 417.949 0.997 

PC2 5.760 0.844 926.030 898.144 0.997 

PC3 1.342 0.428 439.077 469.194 0.998 

PC4 0.240 0.039 227.480 250.154 0.999 

 
PC0, PC1, PC2, PC3 and PC4 are samples mentioned in Table 1 

 
 

 
 
the saturation plateau (Fig. 1). This result indicated 
that sequencing depth has been mainly covered all 
species in the samples and the bacterial richness of 
these time points in the silkworm gut could be 
estimated. It also indicated that species abundance 
of samples and species abundance of PC0 group 
was higher than others. The results of rarefaction 
curve were consistent with the data in Table 1. 

Beta diversity analysis was used to measure the 
changes of species composition on temporal and 
spatial scales by PCoA, NMDS and PCA. The 
differences between individuals and/or groups can 
be observed by PCA and PCoA. NMDS is commonly 
used to compare the differences among the samples 
and evolutionary information. The difference of 
community composition was small if the distance 
between two points was close. 
 
Assignation of OTUs and composition of intestinal 
microbes 

A total of 675, 662 Raw PE, 651,441 effective 
tags and 3121 OTUs were obtained after splitting, 
splicing and filtering. Eliminating tags representing 
chloroplast of mulberry leaves, effective tags with 
97% identity were clustered into OTUs and these 
OTUs were assigned to phylum, class, order, family 
and genus by RDP classifier (Cole et al., 2005). 119 
phyla, 210 classes, 445 orders, 716 families and 
1,163 genera were detected in the intestinal 
contents of 5th instar silkworm larvae (Table 1). 
Among of them, the most abundant phyla were 
Firmicutes (51.41%), Cyanobacteria (32.24%), 
Actinobacteria (3.07%) and Proteobacteria (6.85%), 
respectively. The predominant genera were 
Enterococcus (49.11%), Staphylococcus (6.66%), 
Arthrobacter (3.45%), Lactobacillus (1.04%), 
Pseudomonas (1.07%), respectively.  

 
Changes of intestinal microbes 

To understand the changes of the silkworm 
intestinal microbes, the changes of composition and 
abundance of microbes were analyzed at phylum 
and genus level. The abundance of Firmicutes, 
Cyanobacteria, Proteobacteria, actinobacteria and 
Bacteroidetes was significant in mid-gut as shown in 
the histogram (Fig 2A). In the histogram, the 
abundance of Firmicutes occupied 71.20%, 3.21%, 

13.99%, 72.92%, 98.49% at 0, 12, 24, 48 and 72 h 
time points, respectively; it decreased during 0-12 h 
of the 5th instar, and increased during 12-72 h. At 
these five time points, the abundance of 
Cyanobacteria, Proteobacteria and Actinobacteria 
was 21.44%, 87.32%, 40.04%, 23.56%, 0.3%, 
7.68%, 8.49%, 22.87%, 1.88%, 0.82% and 12.92%, 
0.50%, 6.52%, 1.00%, 0.02%, respectively. Among 
of them, the abundance of Cyanobacteria increased 
rapidly before 12h, then gradually decreased during 
24-72 h; there were no obvious changes about the 
Proteobacteria abundance before 12 h while 
Actinobacteria decreased at this time point, then 
their abundances both increased before 24 h and 
decreased during 48-72 h. 

 
 

 
 

 
 

Fig. 1 Rarefaction curves of the different samples. 
PC0, PC1, PC2, PC3 and PC4 are samples 
mentioned in Table 1. Rarefaction curves of OTUs 
were clustered according to the acquiescent 
sequence identity (97%). In this curve, abscissa axis 
and vertical axis represent the randomly extractive 
sequences number and OTUs number that can be 
constructed based on these sequences number, 
respectively. This curve reflects species abundance 
and sequencing depth. In this curve, species 
abundance of PC0 group was higher than others, 
the quantity of sampling is reasonable for curve 
tends to be flat  



227 

 
 
Fig. 2 The dominant phylum and genera of all the samples at each time point; PC0, PC1, PC2, PC3 and PC4 were 
samples mentioned in Table 1. The map of A and B were the changes of relative abundance of different bacterial 
phylum and genera, respectively. Sequences that could not be classified into the known group were assigned as 
“Unidentified”. In this Fig, the abundance of Enterococcus in Firmicutes varied greatly in five samples, furthermore 
its proportion was higher in PC0, PC3 and PC4 than those of others. The dominant phylum and genera number of 
PC0 and PC2 was more than other samples 
 
 
 
 
 

Enterococcus, Staphylococcus, Arthrobacter, 
Pseudomonas, Lactobacillus, Carnobacterium, 
Bacteroides, Paenibacillus, Serratia, Frondihabitans, 
Brachybacterium and Bacillus were significantly 
abundant (> 1%). Among of them, Enterococcus, 
Staphylococcus, Lactobacillus, Paenibacillus, 
Carnobacterium and Bacillus belong to Firmicutes 
phylum; Arthrobacter, Frondihabitans and 
Brachybacterium are part of Actinobacteria phylum; 
Pseudomonas and Serratia pertain to Proteobacteria 
phylum; Bacteroides belongs to Bacteroidetes 
phylum. The abundance of Enterococcus, 
Staphylococcus and Arthrobacter at five time points 
was 48.32%, 1.99%, 0.44%, 71.91%, 98.03%, 
17.26%, 1.67%, 1.74%, 0.49%, 0.33% and 6.52%, 
0.12%, 1.89%, 0.51%, 0.12%, respectively (Fig 2B). 
In Fig 2B, the abundance of Enterococcus quickly 
decreased before 24 h, then rapidly increased in the 
latter time and reached peak at 72 h; the 
Staphylococcus abundance was gradual downward; 
and the abundance of Arthrobacter decreased 
before 12 h, then increased a little and decreased 
during 48-72 h. 
 
Difference of intestinal microbes at different times 

The gut microbes of silkworm were different at 
five different time points of the 5th instar larvae. At 
phylum level, the number of detected phylum were 
31, 21, 29, 24 and 14, respectively. Among them, the 

number of phylum with significant abundance (> 1%) 
were 5, 3, 5, 4 and 1, respectively (Fig. 2A). The 
abundant changes of Firmicutes, Cyanobacteria, 
Proteobacteria and Actinobacteria were significant. 
The abundance of Firmicutes were approximate at 0 
h and 48 h, their abundances were 22.18, 5.09, and 
0.75 times higher than those of 12 h, 24 h and 72 h, 
respectively. The abundance of Cyanobacteria 
increased before 12 h, then gradually decreased. Its 
abundance of 12 h, 24 h and 48 h was 40.8, 18.68 
and 10.9 times higher than that of 0h, while its 
abundance of 0h was 7.22 times higher than that of 
72 h. The abundance of Proteobacteria was 
basically similar at 0 h and 12 h; compared with the 
abundance of 0 h, it increased at 24 h and was 3.97 
times higher than that of 0 h, then decreased during 
48-72 h, and the abundance of 0 h was 2.88 and 
6.23 times higher than that of 48 h and 72 h, 
respectively. The abundance of Actinobacteria was 
less than the former three phyla, while its abundance 
changed more rapidly. It decreased before 12 h, 
increased at 24 h, then continued to decrease until 
72 h. Based on the abundance of Actinobacteria at 
12 h, the abundance of 0 h, 24 h, 48 h and 72h was 
26.37, 13.27, 2.04 and 0.41 times higher than that of 
12 h, respectively. 

The numbers of detected genera were 368, 163, 
292, 196 and 144, and the number of predominant 
bacterial genera (> 1%) were 8-4-8-2-1 at five time 



228 

 
 
Fig. 3 Common and unique genera analysis of samples at different time points. The analyses results were showed 
by Venn diagram; PC0, PC1, PC2, PC3 and PC4 were samples groups mentioned in Table 1. (A) The Venn 
diagram of PC0, PC1 and PC2 samples groups. (B) The Venn diagram of PC1, PC2 and PC3 samples groups. (C) 
The Venn diagram of PC2, PC3 and PC4 samples groups. (D) The Venn diagram of all samples groups, including 
PC0, PC1, PC2, PC3 and PC4 samples groups 
 
 
 
 
 
points, respectively. Unique genera and common 
ones exist at the samples of different time points (Fig. 
3). Venn diagrams of common and unique genera 
numbers were constructed to understand genera 
changes in the 5th instar. Compared with detected 
intestinal microbes among 0 h, 12 h and 24 h, there 
were 126 common bacteria, 43, 11 and 55 unique 
ones at each time, respectively (Fig. 3A). There 
were 99 shared bacteria, 24, 84 and 16 unique ones, 
among 12 h, 24 h and 48 h (Fig. 3B). 67 common 
bacteria, 101, 20 and 13 unique ones, were present 
in the samples of 24 h, 48 h and 72 h, respectively 
(Fig. 3C). Compared with all detected intestinal 
microbes at these five time points, there were 56 
shared genera, the number of unique ones were 105, 
10, 39, 12 and 16, respectively (Fig. 3D). 

Cluster analysis of intestinal bacteria in the 5th instar 
silkworm 

According to the phylogenetic relationships 
between OTUs, the Unifrac distance between 
samples was calculated by Qiime software (Version 
1.7.0), a heatmap displaying the bacterial similarity 
of different samples by HemI software and the 
analytical map of PCoA by R software (Version 
2.15.3) were constructed (Fig. 4). All the samples 
were assigned to three clusters: (I) PC0.1, PC0.2 
(two repeats of 0 h in the 5th instar) and PC2.2 (the 
sample of 24 h) were grouped into a cluster. (II) 
PC1.1, PC1.2 (two repeats of 12 h) and PC3.1, 
PC3.2 (two repeats of 48 h) were grouped into a 
cluster. (III) PC4.1 and PC4.2 (two repeats of 12 h) 
were grouped (Fig. 4A). 



229 

 
 

Fig. 4 Heatmap and PCoA map of all samples. UniFrac distance was calculated by Qiime software and used to 
construct heatmap and PCoA map. (A) Heatmap was constructed by HemI software. Rows and columns 
represented 5 samples 9 repeats, the similarity represented by the values in the heatmap of microbiota in samples. 
(B) PCoA map was drawn and analyzed by R software, the distance between two dots represented the similarity of 
species structure between two samples 
 
 
 
 
 

In the results map of PCoA, a dot represented a 
sample. The closer the sample distance, the more 
similar the species composition structure was. Our 
PCoA results showed that the samples were clearly 
separated (Fig. 4B); contribution of the principal 
components to samples difference, 72.34% and 
17.14%, were represented by PC1 and PC2 axis, 
respectively. PCoA result suggested that the 
different developmental stage could affect the 
distribution of the silkworm intestinal bacteria. It 
showed that the composition of PC0 and PC2 were 
approached, PC1 and PC3 were nearest. It was also 
showed that the analysis result of PCoA and 
heatmap were concordant. 
 
Phylogenetic tree of predominant microbes phylum 

To understand the similarity of different samples, 
Unweighted Pair-group Method with Arithmetic 
Mean (UPGMA) was used to cluster and analyze 
samples, construct phylogenetic tree of different 
samples, finally solve classification problems. 16S 
rDNA gene sequences of predominant phylum for 
PC0, PC1, PC2, PC3 and PC4 were selected to 
Blast and calculated the Unifrac distance matrix and 
constructed UPGMA phylogenetic tree (Fig. 5). The 
components (phylum) of the phylogenetic tree were 
similar but their abundances were different. 
Compared the abundance of predominant microbes 
phylum, PC0, PC3 and PC4 were assigned to one 
cluster, PC1 and PC2 were assigned to another 
cluster. 

Phylogenetic relationships of genera level 
The relationships of genera Top100 were 

acquired by multiple sequence alignment, the 
relative abundance of these genera was combined 
and integrated, finally the map of phylogenetic 
relationships were gained by GraPhlAn software 
(http://huttenhower.sph.harvard.edu/GraPhlAn) (Fig. 
6). In this map, the evolutionary relationships of 
genera from different samples, information about the 
relative abundance and distribution were displayed. 
 
Discussion 

 
Gut is digestive organ of insect larva, it accounts 

for about 2/3 of the volume of the body cavity at its 
maximum in some insect including silkworm. Insect 
intestinal microbes harbor in digestive tracts with 
feature of low diversity but high variability at intra- 
and interspecific levels (Mathieu et al., 2014), they 
also play important roles in metabolism, 
development, reproduction, protection, resistance 
and immunity (Yuan et al., 2006; Xiang et al., 2007; 
Hedges et al., 2008; Oliver et al., 2009; Vorburger, 
2014). In the same insect species, its composition of 
intestinal microbiota was affected by several factors 
such as age, genetics and environment, especially 
diet (Yun et al., 2014; David et al., 2016). 
Furthermore, the innate immunity and physical 
barriers of insects limited the diversity of gut 
microbes by managing bacteria and formation of 
specific organ to separate the microbes from the host 

http://huttenhower.sph.harvard.edu/GraPhlAn


230 

 
 
Fig. 5 Phylogenetic trees of predominant microbes phylum in different samples. The phylogenetic tree was 
constructed using MEGA6.0 software and the bootstrap value was 1000 replications. Samples groups PC0, PC1, 
PC2, PC3 and PC4 were mentioned in Table 1 
 
 
 
 
 
tissues (Mathieu et al., 2014). The related research 
about other insect was abundant, but few in 
silkworm. 

The information of microbial compositional 
variation in the silkworm gut can help us identify the 
key aspects of microbial composition, and provide 
references for silkworm breeding. In previous study, 
two kinds of methods had been used to study 
intestinal microbes, which were culture-dependent 
and culture-independent methods. The former was 
the traditional low-throughput methods, including the 
methods of isolation and culture-dependent, PCR 
and 16S rDNA- RFLP; they often provide biased 
results and the isolated / detected genera were 
limited (Tian et al., 2007; Xiang et al., 2007). The 
latter based on the analysis of high-throughput 
sequencing and open database, they can provide 
adequate information (Sun et al., 2016). 
Furthermore, rarefaction analysis could judge 
whether the sequencing approach was sufficient or 
not by a plateau and whether the detected bacterial 
diversity was underestimated.  

In the identified intestinal microbes of different 
insects, there were some same phyla in the 
composition, such as Firmicutes and Proteobacteria; 
but there were also some differences with the effect 
of environmental habitat, developmental stage, 
phylogeny of host, especially diet (Xiang et al., 2010; 
Cirimotich et al., 2011). There were also some 
reports about the predominant genera intestinal 
microbes of silkworm. After the silkworm larvae of 
Dongting × Bibo strain reared with tricuspid cudrania 
and mulberry leaves, there were only four shared 
and dominant genera (Brevundimonas, 
Stenotrophomonas, Enterobacter and 
Staphylococcus) in their gut (Xiang et al., 2010). In 
the monophagous C108 and the polyphagous SCN2 
strains, the shared predominant phyla and genera 
were Proteobacteria and Lactobacillales, 
Enterococcus and Thermus, respectively (Xiang et 
al., 2007). In our study some predominant phyla 
included Firmicutes, Cyanobacteria, Actinobacteria 
and Proteobacteria, predominant genera included 
Arthrobacter, Staphylococcus, Enterococcus, 
Lactobacillus and Pseudomonas were identified. 

Compared with the former research, there were 
some same phyla and/or genera, such as 
Staphylococcus and Enterococcus. Maybe like some 
polyphagous lepidopteran insect (David et al., 2016), 
there was a relative stable core microbiota in 
silkworm.  We suspected that silkworm and their 
core intestinal microbes establish a symbiotic 
relationship of co-adaptation and co-evolution to 
maintain survival and reproduction in the 
evolutionary process. 

Some gut microbiotas were related with the 
organism’s resistance, such as Enterobacter 
bacterium (Dillon et al., 2010; Cirimotich et al., 2011). 
In the intestinal bacteria of silkworm, Enterococcus 
is a high frequency microorganism and can reduce 
pH of digestive juice and inhibit some pathogens 
growth in intestinal environment through its products 
(Takizawa et al., 1968). The functions of 
Enterococcus complement each other. Its 
antibacterial function is enhanced in the acidic 
condition, at the same time the pH abatement of 
intestinal digestive juice can protect hosts against 
attacks of toxins (Mead et al., 1988; Broderick et al., 
2004; Li et al., 2010). The appropriate 
supplementary of Enterococcus faecalis in diets can 
increase immunity of organisms (Shi et al., 2015). In 
the prevention and control of silkworm disease, it 
can be used as a possible probiotic for biological 
control since its secreta can strongly inhibit the 
germination of Nosema bombycis spores in vitro (Lu 
and Wang, 2002). These research implied that the 
appropriate supplementary of Enterococcus in vivo 
can enhance immunity of host. In our study, the 
number of Enterococcus increased with the 
prolongation of feeding period. We suspected that 
this may be related with resistance and immunity of 
silkworm. 

The 5th instar is the key stage of silkworm 
development, the amount of ingested mulberry 
leaves of this stage is about 85% of the total amount 
in the whole larval stage. Intake of mulberry leaves 
was continuously increased at the 5th instar before 
gluttonous stage, the numbers of detected genera 
were reduced while the percent of some genera 
were increased in this study, e.g. Enterococcus. It 



231 

 
 
Fig. 6 The map of phylogenetic relationships in genera level. Inside circle was a phylogenetic tree of genera level 
which was constructed by the representative sequences of species. Branch color represented the corresponding 
phylum (Details see the top right corner of legend). Outside circle of each genus represented the relative 
abundance in each group, the colors of relative abundance represented different groups. The thermodynamic 
maps according to the abundance of each genera were then drawn, their values were 10000 times of the original 
value then converted to the value of log 2 
 
 
 
 
may be related to the resistant and metabolic 
function of the silkworm, ingestion and digestion 
substantial mulberry leaves to prepare for the 
development of silk gland and metamorphosis, 
meanwhile it also need increase its resistance to 
defend bad environment and pathogenic 
microorganism. Furthermore, the last instar is also 
the critical period for the growth of symbiotic bacteria 
in host body, the food intake during this period 
provide the conditions for the growth of commensal 
bacteria, and this condition is propitious to the 
common development of insect and its symbiotic 
bacteria (Mathieu et al., 2014). We thought the 
resistance and intake of mulberry leaves may be 
corresponding to the main and/or core microbes of 
silkworm gut.  

Compared with the traditional low-throughput 
methods, much more information of coverage and 
diversity of the gut microbes in hosts was detected 
by high-throughput sequencing techniques 
(Caporaso et al., 2011; Boissière et al., 2012). 16S 
rRNAs were used to evaluate the diversity of 
microbial communities in the gut by high-throughput 
sequencing techniques, microbial taxonomic groups 
were determined and compared by OTUs. The 
changes of proportional OTUs can be observed 
based on the differences of gut microbe. 16S rRNA 
gene is not only a highly conserved and special 
common marker of bacteria, but also can be used to 
predict the functional microbes in the gut for some 
diseases which were associated with changes of 
microbial composition (Turnbaugh et al., 2006; Qin 



232 

et al., 2012). Now some corresponding intestinal 
microbes were used to treat the specific diseases in 
medicine (Zhang et al., 2013). Microbial symbionts 
of insects are potential tools which can improve the 
innate immune homeostasis and contribute to the 
general insect wellbeing (Crotti et al., 2013). 

In this research, we explored the richness and 
dominance of bacterial in the 5th instar silkworm gut 
by sequencing 16S rRNA techniques, and reported 
the differential predominance of bacterial strains in 
five developmental time points. Our results provide 
some theoretical data about the structure and 
change of intestinal microbes in 5th instar silkworm, 
which may be used to develop probiotics, increase 
the conversion rate of leaves and silk, and prevent 
disease in silkworm in the future. Our further studies 
about silkworm gut microorganisms will be focused 
on their composition, acquisition, novel pest control 
strategies and resistant breeding.  
 
Acknowledgments 

This work was supported by the National 
Natural Science Foundation of China (Grant No. 
31672358) and the Natural Science Foundation of 
JiangSu Province (BK20161364). 
 
References 
Aptedeshpande AD, Paingankar MS, Gokhale MD, 

Deobagkar DN. Serratia odorifera mediated 
enhancement in susceptibility of Aedes aegypti 
for chikungunya virus. Indian J. Med. Res. 139: 
762-768, 2014. 

Asnicar F, Weingart G, Tickle TL, Huttenhower C, 
Segata N. Compact graphical representation of 
phylogenetic data and metadata with GraPhlAn. 
PeerJ. 3: e1029, 2015. 

Boissière A, Tchioffo MT, Bachar D, Abate L, Marie 
A, Nsango SE, et al. Midgut microbiota of the 
malaria mosquito vector Anopheles gambiae 
and interactions with Plasmodium falciparum 
infection. PLoS Pathog. 8: e1002742, 2012. 

Broderick NA, Raffa KF, Goodman RM, Handelsman 
J. Census of the bacterial community of the 
gypsy moth larval mid gut by using culturing and 
culture-independent methods. Appl. Environ. 
Microbiol. 70: 293-300, 2004. 

Caporaso JG, Lauber CL, Walters William A, 
Berg-Lyons D, Lozupone CA, Turnbaugh PJ, et 
al. Global patterns of 16S rRNA diversity at a 
depth of millions of sequences per sample. Proc. 
Natl. Acad. Sci. USA 108S1: 4516-4522, 2011. 

Chengxiang H, Han W, Dingding L, Ruilin L, Xijie G. 
Molecular cloning and characterization of High 

Mobility Group box （ HMGB ）  gene from 
Beauveria bassiana- infected silkworm, Bombyx 
mori. Inv. Surv. J. 14: 157-164, 2017. 

Cirimotich CM, Dong Y, Clayton AM, Sandiford SL, 
Souza-Neto JA, Mulenga M, et al. Natural 
microbe-mediated refractoriness to Plasmodium 
infection in Anopheles gambiae. Science 332: 
855-858, 2011. 

Cole JR, Chai B, Farris RJ, Wang Q, Kulam SA, 
McGarrell DM, et al. The Ribosomal Database 
Project (RDP-II): sequences and tools for 
high-throughput rRNA analysis. Nucleic Acids 
Res. 33: 294-296, 2005. 

Crotti E, Sansonno L, Prosdocimi EM, Vacchini V, 

Hamdi C, Cherif A, et al. Microbial symbionts of 
honeybees: a promising tool to improve 
honeybee health. New Biotechnol. 30: 716-722, 
2013. 

David MR, Santos LMBD, Vicente ACP, 
Macieldefreitas R. Effects of environment, 
dietary regime and ageing on the dengue vector 
microbiota: evidence of a core microbiota 
throughout Aedes aegypti lifespan. Mem. Inst. 
Oswaldo Cruz. 111: 577-587, 2016. 

Dillon RJ, Dillon VM. The gut bacteria of insects: 
nonpathogenic interactions. Annu. Rev. 
Entomol. 49: 71-92, 2004. 

Dillon RJ, Vennard CT, Buckling A, Charnley AK. 
Diversity of locust gut bacteria protects against 
pathogen invasion. Ecol. Lett. 8: 1291-1298, 
2010. 

Feldhaar H. Bacterial symbionts as mediators of 
ecologically important traits of insect hosts. 
Ecological Entomol. 36: 533-543, 2011. 

Goldsmith MR, Shimada T, Abe H. The genetics and 
genomics of the silkworm, Bombyx mori. Annu. 
Rev. Entomol. 50: 71-100, 2005. 

Hedges LM, Brownlie JC, O’Neill SL, Johnson KN. 
Wolbachia and virus protection in insects. 
Science 322: 702, 2008. 

Hou C, Qin G, Liu T, Geng T, Gao K, Pan Z, et al. 
Transcriptome Analysis of Silkworm, Bombyx 
mori, during Early Response to Beauveria 
bassiana Challenges. PloS One. 9: e91189, 
2014. 

Hou C, Qin G, Liu T, Mei X, Zhang R, Zhao P, et al. 
Differential gene expression in silkworm in 
response to Beauveria bassiana infection. Gene 
484: 35-41, 2011. 

Lee WJ, Brey PT. How microbiomes influence 
metazoan development: insights from history 
and Drosophila modeling of gut-microbe 
interactions. Annu. Rev. Cell Develop. Biol. 29: 
571-592, 2013. 

Li B, Wang XY, Hou CX, Xu AY, Li MW. Genetic 
analysis of quantitative trait loci for cocoon and 
silk production quantity in Bombyx mori 
(Lepidoptera: Bombycidae). Eur. J. Entomol. 
110: 205-213, 2013. 

Li H, Liu Y, Zhang XW. In vitro evaluation of 
antibacterial effect of 2.5% sodium hypochlorite 
adjusted to PH 12, 9, 7.5 and 6 on 
Enterococcus faecalis. J. Oral Sci. Res. 3: 
323-325, 2010. 

Li J, Qin S, Yu H, Zhang J, Liu N, Yu Y, et al. 
Comparative transcriptome analysis reveals 
different silk yields of two silkworm strains. PloS 
One 11: e0155329, 2016. 

Li M, Hou C, Miao X, Xu A, Huang Y. Analyzing 
genetic relationships in Bombyx mori using 
intersimple sequence repeat amplification. J. 
Econ. Entomol. 100: 202-208, 2007. 

Lu XM, Wang FW. Inhibition of cultured supernatant 
of enterococci strains on germination of 
Nosema bombycis spores in vitro. Acta Sericol. 
Sinica 28: 126-128, 2002. 

Maslowski KM, Mackay CR. Diet, gut microbiota and 
immune responses. Nat. Immunol. 12: 5-9, 
2011. 

Mathieu P, Stephen JS, Fleur P. Towards an 
integrated understanding of gut microbiota 



233 

using insects as model systems. J. Insect 
Physiol. 69: 12-18, 2014. 

McHardy IH, Li X, Tong M, Ruegger P, Jacobs J, 
Borneman J et al. HIV Infection is associated 
with compositional and functional shifts in the 
rectal mucosal microbiota. Microbiome, 1:26, 
2013 

Mead LJ, Khachatourians GG, Jones GA. Microbial 
ecology of the gut in laboratory stocks of the 
migratory grasshopper, Melanoplus sanguinipes 
(Fab.) (Orthoptera: Acrididae). Appl. Environ. 
Microbiol. 54: 1174-1181, 1988. 

Oliver KM, Degnan PH, Hunter MS, Moran NA. 
Bacteriophages encode factors required for 
protection in a symbiotic mutualism. Science 
325: 992-994, 2009. 

Qin J, Li Y, Cai Z, Li S, Zhu J, Zhang F, et al. A 
metagenome-wide association study of gut 
microbiota in type 2 diabetes. Nature 490: 55-60, 
2012.  

Schloss PD. Evaluating different approaches that 
test whether microbial communities have the 
same structure. ISME J. 2: 265-275, 2008. 

Shi Z, Yao Y, Jiang S, Xiao R, Liu Z, Yang F, Dong 
G, et al. Effects of Enterococcus faecalis 
substitute for antibiotic on growth performance, 
diarrhea rate, blood biochemical parameters 
and immune organs of weaner piglets. Chinese 
J. Anim. Nutr. 27: 1832-1840, 2015. 

Simon JC, Boutin S, Tsuchida T, Koga R, Gallic JFL, 
Frantz A, et al. Facultative symbiont infections 
affect aphid reproduction. PloS One. 6: e21831, 
2011. 

Sun Z, Lu Y, Zhang H, Kumar D, Liu B, Gong Y, et al. 
Effects of BmCPV infection on silkworm 
Bombyx mori intestinal bacteria. PloS One. 11: 
e0146313, 2016. 

Takizawa Y, Iizuka T. The aerobic bacterial flora in 
the gut of larvae of the silkworm, Bombyx mori L. 
(I).The relation between media and the numbers 
of living cells. J. Seric. Sci. Jpn.. 37: 295-305, 

1968. 
Tian ZH, Hui FL, Ke T, Kan YC, Wen ZZ. Molecular 

analysis of the bacteria community composition 
in silkworm midgut. Sericultural science 
magazine. 33: 592-595, 2007. 

Turnbaugh PJ, Ley RE, Mahowald MA, Magrini V, 
Mardis ER, Gordon JI. An obesity-associated 
gut microbiome with increased capacity for 
energy harvest. Nature 444: 1027-1031, 2006. 

Vorburger C. The evolutionary ecology of 
symbiont-conferred resistance to parasitoids in 
aphids. Insect Sci. 21: 251-264, 2014. 

Xiang H, Li MW, Zhao Y, Zhao LP, Zhang YH, 
Huang YP. Bacterial community in midguts of 
the silkworm larvae estimated by PCR/DGGE 
and 16S rDNA gene library analysis. Acta 
Entomol. Sin. 50: 222-233, 2007. 

Xiang Y, Wang X, Feng W, Zhou W, Xie H, Wang Y. 
Comparative analysis of the composition of 
dominant intestinal microflora in silkworm 
reared with different forages. Acta Ecol. Sin. 30: 
3875-3882, 2010. 

Yuan ZH, Lan XQ, Yang T, Xiao J, Zhou ZY. 
Investigation and analysis of the bacteria 
community in silkworm intestine. Acta Microbiol. 
Sin. 46: 285-291, 2006. 

Yun JH, Roh SW, Whon TW, Jung MJ, Kim MS, 
Park DS, et al. Insect gut bacterial diversity 
determined by environmental habitat, diet, 
developmental stage, and phylogeny of host. 
Appl. Environ. Microbiol. 80: 5254-5264, 
2014. 

Zhang FM, Wang HG, Wang M, Cui BT, Fan ZN, Ji 
GZ. Fecal microbiota transplantation for severe 
enterocolonic fistulizing Crohn’s disease. World 
J. Gastroenterol. 19: 7213-7216, 2013. 

Zhang J, Blessing D, Wu CY, Liu N, Li J, Qin S, et al. 
Comparative transcriptomes analysis of the 
wing disc between two silkworm strains with 
different size of wings. PloS One. 12: e017956, 
2017.