Identification and Functional Analysis of Genes Related to Cocoon Shell Ratio in Bombyx mori
-
摘要:目的 挖掘家蚕茧层率相关基因，以期为家蚕茧层率性状分子遗传改良提供参考依据。方法 以多丝量家蚕品种菁松和中丝量家蚕品种芙蓉为亲本构建BC1代分离群体。在BC1代雄性群体中挑选极端高/低茧层率个体构建子代DNA混池，运用BSA-seq方法对茧层率关联区域进行定位，并运用BLAST软件对关联区域的编码基因进行GO和KEGG等数据库注释及功能预测。结果 重测序数据与家蚕参考基因组平均比对率为98.86%，平均基因组覆盖度为95.79%（1×）和88.63%（5×）；变异检测共获得26 557 646个SNPs；∆(SNP-index)定位到3个与茧层率显著相关的区域，分别为Chr.2:4430～4930 kb、Chr.4:12350～12920 kb和Chr.13:3230～3730 kb，共包含70个编码基因。通过GO与KEGG注释，有58个基因注释到GO数据库，涉及生物过程、分子功能和细胞组分三大类；有19个基因注释到KEGG通路，分布于34个代谢通路中。通过KEGG代谢通路分析，筛选出10个可能对家蚕茧层率有重要调控作用的基因，推测其参与了家蚕丝腺细胞运动、能量代谢和蛋白质合成加工。结论 运用BSA-seq方法在家蚕第2、4和13号染色体上定位到与茧层率关联的区域，筛选到10个可能与茧层率密切相关的候选基因，为茧层率关键调控基因精细定位及克隆奠定基础。Abstract:Objective Genes related to the cocoon shell ratio of Bombyx mori were investigated for improvement on silkworm productivity.Method Segregated BC1 populations of a highly productive silkworm, Jingsong, and a moderately productive Furong were established. DNA pools were constructed by mixing 30 high and low cocoon shell ratio individuals from the BC1 male populations. BSA-seq was applied to identify the regions relevant to the target trait. The coding genes at the candidate regions were annotated with BLAST software in GO and KEGG databases.Result The percentage of clean reads that matched the reference genome of P50 was 98.86%, and those of the average genome coverage 95.79% for 1× and 88.63% for 5×. Variant calling resulted in 26 557 646 SNPs. Three QTLs for the cocoon shell ratio detected by ∆(SNP-index) were in the intervals of 4 430–4930 kb on Chr.2, 12350–12920 kb on Chr.4, and 3230–3730 kb on Chr.13. There were 70 coding genes found in the associated regions. According to the GO database annotation, 58 genes were categorized in the groups of biological process, molecular function, or cellular component. The 19 genes annotated in the KEGG database distributed in 34 metabolic pathways. Of which, 10 might play important roles in regulating the cocoon shell ratio and the functions involving the silk gland cell movement, energy metabolism, and/or protein synthesis.Conclusion The genes related to the silkworm cocoon shell ratio were identified by BSA-seq method. Ten genes were speculated to possibly associate with the silk formation. The present study only managed to locate the regions on chromosomes 2, 4, and 13 of the genes in B. mori that might related to the cocoon shell ratio. Further study will be needed to finely map and clone the key regulatory genes in the highly complex regulatory mechanism.
-
Keywords:
- Bombyx mori /
- silkworm cocoon shell ratio /
- BSA-seq /
- gene annotation
-
0. 引言
【研究意义】蚕丝是家蚕利用最广泛的产物，不仅可用于制作丝绸制品，而且随着丝素蛋白和丝胶蛋白新用途的开发，蚕丝在生物医学、工程材料、食品和化妆品领域也具有较好的应用前景[1]。茧层率是反映家蚕对丝素蛋白、丝胶蛋白合成能力的重要性状，因此挖掘茧层率关键基因，对家蚕分子遗传育种改良和提高家蚕饲养经济效益具有重要意义[2]。【前人研究进展】家蚕饲养至今，茧层率性状经过两次较大的改良，茧层率提高至20%左右，甚至部分高茧层率品种达到30%[3]。但由于家蚕茧层率性状与其生命力、生殖力或抵抗力具有一定的相关关系，部分高茧层率家蚕品种未能在生产上推广应用。采用分子遗传育种手段定向改良茧质性状，是近年来家蚕育种工作的重点[4]。鲁成等[5]、李斌等[6]、司马杨虎等[7]利用AFLP（Amplified fragment length polymorphism，AFLP）标记构建了家蚕分子遗传图谱；侯成香等[8]、Li等[9]构建了SSR（Simple sequence repeat，SSR）和SNP（Single nucleotide polymorphism，SNP）分子遗传图谱，使家蚕茧层率QTL（Quantitative trait loci，QTL）定位研究取得了一定进展。混合群体分离分析（Bulked segregant analysis，BSA）与二代测序相结合的方法（BSA-Seq）是利用目标性状存在差异的2个亲本构建子代分离群体，在分离群体中选取极端性状个体构建DNA混池，通过高通量测序、分子标记挖掘等技术，对目标性状进行基因定位的有效手段[10-13]。BSA-Seq相比构建分子遗传图谱具有测序群体小、成本和测序要求低等优势，已经在油料、水果等作物的农艺性状QTL定位上广泛应用[14-18]。柳海东等[16]利用甘蓝型油菜早花亲本No.4512和晚花亲本No.5246构建成DH系，从中选取早花和晚花单株各20个，利用BSA-seq方法对早花位点cqDTFC8进行定位分析，将早花位点定位在C8染色体上的1.3～6.0 Mb，筛选出2个与光周期调控有关的基因，命名为BnaC08g04860D和BnaC08g04960D。尹明智等[17]利用BSA-Seq法对油菜胞质雄性不育恢复基因进行了定位分析，将恢复基因定位在C09染色体的0～880000 bp区域，共筛选出16个基因与育性恢复有关。欧点点等[18]利用甜瓜黄皮材料和绿皮材料构建了F2分离群体，利用BSA-Seq法将皮色相关基因定位在第4染色体的0.02～5.7 Mb区间和第10染色体的0.08～9.54 Mb区间。BSA-Seq在家蚕茧质性状主效基因定位上也有应用的报道，Li等[2]基于BSA-Seq方法将与丝产量连锁的标记定位在第11、22和23号染色体上，结合候选基因表达模式、丝腺表达分析，筛选出1个调控家蚕茧层量的关键基因，命名为BmRPL18。【本研究切入点】目前有关茧层率相关QTLs的功能验证还有待进一步研究。【拟解决的关键问题】本研究以多丝量家蚕和中丝量家蚕品种为亲本，以构建的BC1群体为研究对象，采用BSA-seq方法对茧层率相关基因进行定位分析，进一步挖掘茧层率关键基因，为高茧层率关键基因克隆及家蚕经济性状分子遗传育种实践提供理论依据。
1. 材料与方法
1.1 试验材料
重测序材料选取多丝量家蚕品种菁松（高茧层率）和中丝量家蚕品种芙蓉（中茧层率）构建的BC1群体，亲本材料由贵州省蚕业研究所提供，均经多代自交，背景纯合。2019年春，在贵州遵义以菁松、芙蓉杂交获得F1代。2019年夏，在贵州贵阳饲养F1代，并以F1代雄蛾回交芙蓉，获得BC1分离群体 芙蓉×（芙蓉×菁松） 。2020年春，在遵义基地同室、同条件饲养亲本和BC1群体。
1.2 茧质性状调查
在上蔟结茧的第7天采集亲本和BC1群体的蚕茧，削茧以鉴别雌、雄蛹。因家蚕茧质性状在雌、雄个体间差异显著，为排除试验误差，本试验在亲本中随机选择雌、雄茧各20粒，在BC1群体中选择所有雄性个体蚕茧（BC1M），调查记录个体的全茧量、茧层量，计算个体茧层率。
1.3 基因组DNA提取及混池构建
根据BC1M群体的茧层率调查结果，从中选择高茧层率和低茧层率极端分离材料各30个，同时取芙蓉、菁松2个亲本，运用DNA提取试剂盒提取样本基因组DNA，用1%的琼脂糖凝胶电泳检测DNA的完整性，用NanoDrop2000检测DNA样品的纯度，运用PicoGreen检测DNA的浓度。质检合格后将30个高茧层率样本DNA和30个低茧层率样本DNA分别等量混合，构建高茧层率子代DNA混池（H池）和低茧层率子代DNA混池（L池），并分别提取2个亲本的基因组DNA构建2个亲本池。
1.4 混池基因组重测序及变异位点检测
将基因组DNA委托给上海美吉生物医药科技有限公司进行建库测序，测序平台为Illumina Novaseq6000。亲本池测序深度为10×，后代混池测序深度为20×。测序数据（Raw data）下机后进行质量控制，过滤得到高质量的Clean data。利用BWA软件将Clean data比对家蚕P50参考基因组序列（http://silkbase.ab.a.u-tokyo.ac.jp/cgi-bin/download.cgi），获得序列的位置归属BAM文件。利用GATK的Best Practices流程对BAM文件进行校正，并进行SNP标记的变异检测，利用SnpEff软件和参考基因组的基因预测信息进行SNP的功能注释，筛选可能影响茧层率的变异位点。
1.5 连锁分析
过滤掉两子代混池间基因型一致的位点、有多个基因型的位点、read支持度小于4的SNP位点。以高茧层率亲本菁松为参考，基于过滤和筛选获得的SNP位点，对H池、L池中每个位点的SNP-index值进行计算；以0.5 Mb为窗口，10 kb为步长，采用滑窗策略对子代SNP-index在染色体上的分布进行作图，并取H池与L池SNP-index的差值计算∆（SNP-index）；进行1000次置换检验，选择99.99%置信水平为筛选阈值，将∆（SNP-index）在置信水平之上的窗口作为候选区域。
1.6 候选区域基因的功能注释
应用BLAST软件对候选区域的编码基因进行NR（Non redundant）、GO（Gene Ontology）、KEGG（Kyoto encyclopedia of genes and genomes）、UniProt（Universal protein）、EggNOG（Evolutionary genealogy of genes: Non-supervised orthologous groups）数据库的深度注释，对候选基因进行功能预测。
2. 结果与分析
2.1 基因组重测序及变异位点检测
对30个高茧层率个体、30个低茧层率个体及2个亲本材料进行全基因组重测序，共获得原始数据量（Raw data）10.32 Gb，过滤后得到高质量Clean reads数在38269934～107624014，Q20≥97.91%，Q30≥93.56%，GC含量为38.66%～38.96%（表1）。将样品高质量测序数据与参考基因组进行比对，平均比对率为98.86%，成功比对到基因组上的reads比例为81.29%～82.72%，平均测序深度为18.35×，平均基因组覆盖度为95.79%（1×）、88.63%（5×）（表2）。说明样本数据量足够，测序质量良好，GC分布正常，建库测序成功；质控数据与家蚕参考基因组同源性较高，比对结果正常，可用于后续变异检测与定位分析。
表 1 测序数据质量Table 1. Statistics on quality of sequencing data样品编号
Sample ID原始数据量
Raw data/bp过滤后数据量
Clean base/bp原测序reads 数
Raw reads过滤后reads 数
Clean readsQ20/% Q30/% GC含量
GC content/%菁松（JS） 380002170 6295123605 42251670 41738170 98.10 94.04 38.83 芙蓉（FR） 5853501108 5771738817 38764908 38269934 97.91 93.56 38.96 H-pool 16080222540 15911340150 106491540 105500096 98.13 94.11 38.89 L-pool 16401318000 16232070628 108618000 107624014 98.13 94.11 38.66 Q20：高质量测序数据中质量值≥20的碱基所占百分比；Q30：高质量测序数据中质量值≥30的碱基所占百分比。
Q20：The percentage of the bases whose Phred value are more than 20; Q30：The percentage of the bases whose Phred value are more than 30.表 2 质控数据与参考基因组比对情况Table 2. Matching between quality control data and reference genome样品编号
Sample ID比对率
Mapped ratio/%比对到基因组上的reads比例
Properly ratio/%平均测序深度
Average depth基因组覆盖度（1×）
Genome coverage(1×) /%基因组覆盖度（5×）
Genome coverage(5×) /%菁松（JS） 98.84 82.72 10.34 94.70 84.71 芙蓉（FR） 98.88 82.28 9.44 94.41 82.61 H-pool 98.90 81.29 26.02 97.04 93.59 L-pool 98.82 82.38 27.59 97.02 93.59 1×覆盖度：1 个碱基覆盖的位点占基因组的百分比；5×覆盖度：5个碱基覆盖的位点占基因组的百分比。
Coverage 1×: the percentage of at least 1 base-covered site in reference genome; Coverage 5×: the percentage of at least 5 base-covered sites in reference genome.根据样品数据与家蚕参考基因组比对结果，利用GATK软件进行SNP变异位点检测，共获得26 557 646个SNP标记。利用SnpEff软件进行变异位点功能预测，结果表明SNP位点主要位于基因间区、内含子区域、基因上游和基因下游，非编码区域的多态性位点显著多于编码区。
2.2 连锁分析
根据SNP过滤原则，对变异检测获得的位点进行过滤，最终得到374 463个高质量SNP位点，使用Circos软件对变异位点在染色体上的分布密度进行可视化（图1）。根据两子代混池的SNP-index进行连锁分析，结果显示∆（SNP-index）在置信水平之上的区域有3个，分别位于第2、4、13染色体上（图2），总长度为1.57 Mb，共包含67个SNP位点、11个InDel位点、70个编码基因（表3）。
![]()
图 1 分子标记及关联信号在染色体上的分布从外到内依次为参考基因组染色体坐标、染色体上基因分布（颜色越深表示基因密度越大）、SNP密度分布（圆点越密集表示SNP密度越大）、InDel密度分布（三角形越密集表示InDel密度越大）、Index值在染色体上的分布。Figure 1. Distribution of SNPs, InDels, and associated signals on chromosomeShown from outside inward: chromosome coordinates of reference genome, genes distribution on chromosome (darker color indicates greater gene density), SNP density distribution (density of dots corresponds to that of SNP), InDel density distribution (density of triangles reflects that of InDel), and distribution of indices on chromosome.![]()
图 2 H池和L池SNP-index、∆（SNP-index）分布情况图中不同颜色表示不同的染色体，横坐标为1～28号染色体上每个window的具体物理位置，纵坐标为位置所对应的Index值。Figure 2. SNP-index and ∆(SNP-index) distribution of H-pool and L-pool.Different colors represent different chromosomes; x-axis is for physical location of each window on chromosomes 1 to 28; y-axis is for Index corresponding to respective locations.表 3 关联区域信息统计Table 3. Statistics of the related genes染色体编号
Chromosome ID关联区域起点
Start of associated
regions/bp关联区域终点
End of associated
regions/bp关联区域长度
Associated region
size/MbSNP数量
SNP number关联区域内基因个数
Gene number in the
associated regions第13染色体 Chr 13 3230000 3730000 0.50 11 13 第4染色体 Chr 4 12350000 12920000 0.50 39 48 第2染色体 Chr 2 4430000 4930000 0.57 17 9 合计 Total 1.57 67 70 2.3 候选区域基因注释及候选基因筛选
利用BLAST软件对关联区域内的编码基因进行多个数据库（NR、GO、KEGG、UniProt、EggNOG）注释，结果显示：70个编码基因中有69个注释到NR数据库中，69个注释到UniProt数据库中，58个注释到GO数据库中，19个注释到KEGG通路，57个注释到EggNOG数据库。
通过对GO数据库中候选区域基因参与的生物过程（Biological process）、分子功能（Molecular function）和细胞组分（Cellular component）进行分类分析，结果表明，在生物过程中候选基因主要富集在上皮细胞发育、脑室系统发育和气管发育等；在细胞组分中，候选基因多富集于细胞核、细胞质、细胞膜；在分子功能中，候选基因富集最多的是组蛋白乙酰转移酶（图3）。
为进一步了解候选区域内基因的生物学功能，对富集到KEGG数据库中的基因进行分析，结果显示候选区域中19个基因分布于34个信号通路中，涉及新陈代谢（Metabolism）、环境信息加工（Environment information processing）、细胞进程（Cellular processes）、遗传信息加工（Genetic information processing）、有机体系统（Organismal systems）（表4）。根据基因参与的代谢通路分析及文献检索，共筛选出10个与家蚕茧层率密切相关的候选基因：KWMTBOMO00853、KWMTBOMO02140、KWMTBOMO02143、KWMTBOMO02145、KWMTBOMO02146、KWMTBOMO02147、KWMTBOMO02152、KWMTBOMO07652、KWMTBOMO07653、KWMTBOMO07659，可能主要参与家蚕丝腺细胞运动、能量代谢和蛋白质合成加工。
表 4 候选基因的KEGG通路分析Table 4. KEGG pathway of genes in candidate regions一级代谢
Primary metabolism二级代谢
Secondary metabolism三级代谢
Tertiary metabolism通路编号
Ko ID基因编号
Gene ID新陈代谢
Metabolism聚糖生物合成与代谢
Glycan biosynthesis and metabolism糖胺聚糖降解
Glycosaminoglycan degradationko00531 KWMTBOMO07657 O-聚糖生物合成
Other types of O-glycan biosynthesisko00514 KWMTBOMO02149 脂质代谢
Lipid metabolism初级胆汁酸生物合成
Primary bile acid biosynthesisko00120 KWMTBOMO02138 有机体系统
Organismal systems免疫系统
Immune systemRIG-I样受体信号通路
RIG-I-like receptor signaling pathwayKo04622 KWMTBOMO02145；
KWMTBOMO07659；
KWMTBOMO02143Toll样受体信号通路
Toll-like receptor signaling pathwayko04620 KWMTBOMO07659 NOD样受体信号通路
NOD-like receptor signaling pathwayko04621 KWMTBOMO07659 内分泌系统
Endocrine system胰高血糖素信号通路
Glucagon signaling pathwayKo04922 KWMTBOMO07653；
KWMTBOMO07652；
KWMTBOMO02146甲状腺激素信号通路
Thyroid hormone signaling pathwayko04919 KWMTBOMO07653；
KWMTBOMO07652胰岛素信号通路
Insulin signaling pathwayko04910 KWMTBOMO02146 神经系统
Nervous system长时程增强效应
Long-term potentiationko04720 KWMTBOMO07653；
KWMTBOMO07652神经营养因子信号通路
Neurotrophin signaling pathwayko04722 KWMTBOMO07659 环境适应
Environmental adaptation生理节律 Circadian rhythm ko04710 KWMTBOMO00853 环境信息加工
Environment information processing信号转导
Signal transductionRas信号通路 Ras signaling pathway ko04014 KWMTBOMO02152 丝裂原活化蛋白激酶信号通路
MAPK signaling pathwayKo04013 KWMTBOMO02145；
KWMTBOMO07659；
KWMTBOMO02143钙离子信号通路
Calcium signaling pathwayko04020 KWMTBOMO02146 低氧诱导因子-1信号通路
HIF-1 signaling pathwayko04066 KWMTBOMO07653；
KWMTBOMO07652Wnt信号通路
Wnt signaling pathwayko04310 KWMTBOMO07653；
KWMTBOMO07652；
KWMTBOMO00853环磷酸腺苷信号通路
cAMP signaling pathwayko04024 KWMTBOMO07653；
KWMTBOMO07652Notch信号通路 Notch signaling pathway ko04330 KWMTBOMO07653；
KWMTBOMO07652Jak-STAT信号通路
Jak-STAT signaling pathwayko04630 KWMTBOMO07653；
KWMTBOMO07652TGF-β信号通路
TGF-beta signaling pathwayko04350 KWMTBOMO07653；
KWMTBOMO00853；
KWMTBOMO07652刺猬信号通路
Hedgehog signaling pathwayko04341 KWMTBOMO00853 细胞进程
Cellular processes细胞生长和死亡
Cell growth and death细胞周期 Cell cycle Ko04110 KWMTBOMO07653；
KWMTBOMO07652；
KWMTBOMO00853细胞凋亡 Apoptosis ko04214 KWMTBOMO07659 卵母细胞减数分裂 Oocyte meiosis ko04114 KWMTBOMO00853 细胞通讯 Cell communication 黏着连接 Adherens junction ko04520 KWMTBOMO07653；
KWMTBOMO07652运输与分解代谢
Transport and catabolism内吞作用 Endocytosis ko04144 KWMTBOMO07659 溶酶体 Lysosome ko04142 KWMTBOMO07657 遗传信息加工 Genetic information processing 折叠、组装和降解
Folding, sorting and degradation泛素介导的蛋白质水解
Ubiquitin mediated proteolysisko04120 KWMTBOMO07659；
KWMTBOMO02114；
KWMTBOMO00853蛋白质在内质网上的加工
Protein processing in endoplasmic reticulumko04141 KWMTBOMO02147；
KWMTBOMO00853蛋白酶体 Proteasome ko03050 KWMTBOMO02148 转录 Transcription 转录因子 Basal transcription factors ko03022 KWMTBOMO02120 剪接体 Spliceosome ko03040 KWMTBOMO02133 翻译 Translation 核糖体 Ribosome ko03010 KWMTBOMO02140 3. 讨论
家蚕的部分经济性状具有性别差异，主要表现在雄蚕比雌蚕体质强健、茧层率和出丝率高，雌蚕比雄蚕全茧量和茧层量高、蛹体重大等，因此学者们推测家蚕性染色体上可能存在控制茧层率、全茧量、茧层量、蛹体重等茧质性状的基因[19-20]。Zhan等[21]以菁松和兰10为亲本构建回交一代群体（BC1M），利用SSR分子标记对家蚕部分茧质性状进行了QTL初步定位，结果在性染色体上检测到了与全茧量、茧层量和蛹体重相关的QTLs，但未检测到与茧层率相关的基因位点。侯成香等[8]在此基础上，同样以菁松和兰10为亲本配置BC1M群体，利用SSR和SNP标记构建家蚕性染色体分子遗传图谱，对茧、丝有关的QTLs进行精细扫描，但同样未扫描到与茧层率相关的QTLs。本研究以菁松和芙蓉为亲本构建BC1M群体，利用BSA-Seq分析方法对控制茧层率的基因进行定位，结果将茧层率相关的QTLs定位在第2染色体（4430～4930 kb）、第4染色体（12350～12920 kb）和第13染色体（3230～3730 kb），也未在性染色体上检测到与茧层率相关的基因位点。鉴于茧层率是茧层量对全茧量的比率，主要由茧层量和蛹体重决定，而前人研究表明性染色体上存在控制茧层量与蛹体重的基因，因此分析认为，茧层率与性别关联可能是由于其决定因素（茧层量、蛹体重）与性别关联。
有关茧层率的数量性状位点定位研究早有报道，鲁成等[5]利用大造和C100杂交的BC1群体，构建了AFLP分子标记连锁图谱，通过复合区间作图将茧层率相关的QTLs定位于第2、11、15连锁群上；李斌等[6]利用大造和C100的BC1群体构建AFLP分子标记连锁图谱，将茧层率相关的QTLs定位在第2、4、14、15、19、25连锁群上；司马杨虎等[7]对湘晖和872构建的F2群体构建了AFLP分子标记连锁图谱，将茧层率性状的QTLs定位在第2、11、13、18连锁群上；Zhan等[21]利用SSR分子标记对菁松和兰10的BC1M群体进行QTL位点检测，将茧层率相关QTLs定位在第18和19连锁群；Li等[9]利用菁松和兰10的BC1M群体，使用STS、SSR和SNP分子标记构建遗传图谱，结果未扫描到与茧层率相关的QTLs。Fang等[22]利用夏芳和野蚕的BC1M群体构建了SNP连锁图谱，对家蚕茧丝产量相关的QTLs进行了定位，将与茧层率有关的QTL定位在第13号染色体上。结合本研究结果可以看出，无论采用相同或不同的亲本材料、作图群体和分子标记，家蚕茧层率相关基因的定位结果均有差异，定位位置重复性较差。分析认为茧层率相关基因定位重复性差与其遗传基础复杂有关，即可能存在多个家蚕茧层率性状关键基因位点，且其效应有差异。本研究筛选出了10个与茧层率相关联的候选基因，其KEGG通路分析显示，候选基因可能参与了丝腺细胞运动、能量代谢和蛋白质合成加工过程，在一定程度上证实了茧层率相关基因调控机制的复杂性。
家蚕茧丝中的丝素蛋白在后部丝腺中合成，后部丝腺细胞为核内复制方式，即细胞周期只有DNA复制时期（S期）和细胞间期（G期），而不具有分裂期（M期）和胞质分裂过程[23]。细胞周期蛋白E（Cycline-E，CycE）、细胞周期蛋白依赖激酶2（Cyclin-dependent kinases 2，CDK2）是S期的主要调节因子，与其他核内周期调节者构成复杂的调控网络，精密协作保证核内周期的顺利完成[24]。前人研究发现，果蝇Ras信号（Ras85D）可能通过调控CycE、CDK2等核内周期调节因子，参与调控细胞的存活、生长等生理过程。Caldwell等[25]通过激活果蝇前胸腺细胞（由核内复制细胞组成）中的Ras信号，使前胸腺细胞体积明显增大。Ma等[26]利用GAL4/UAS技术在家蚕后部丝腺中特异性过表达Ras1CA激活了Ras信号，促进了后部丝腺细胞和细胞核体积增大，全茧量有所增加。马倩等[23]基于RNA-Seq技术比较了野生型与Ras1CA过表达转基因家蚕后部丝腺组织中的差异基因，结果发现Ras信号可使细胞周期依赖蛋白激酶（CDK）、细胞周期蛋白（Cycline）和DP-1,2转录因子等编码基因上调，进而调控细胞周期通路，促进后部丝腺组织细胞核内复制，使丝腺组织增大，其认为在此过程中丝素蛋白编码基因的拷贝数也大量增加，可能有利于丝蛋白的合成。本研究通过BSA-Seq方法定位到10个与茧层率性状密切相关的候选基因，其中KWMTBOMO02152直接参与了Ras信号通路；KWMTBOMO00853、KWMTBOMO07652和KWMTBOMO07653参与了细胞周期通路；KWMTBOMO02143、KWMTBOMO02145和KWMTBOMO07659参与了丝裂原活化蛋白激酶（Mitogen-activated protein kinase，MAPK）信号通路，该通路通过级联方式将细胞外信号逐级扩大并传导到细胞或细胞核内，进而调控细胞周期的运行和基因表达，也是Ras信号通路和细胞周期通路的重要组成部分；KWMTBOMO02146参与了胰岛素信号通路和钙离子信号通路，其中胰岛素信号通路参与调控能量代谢、细胞生长及蛋白质合成过程，钙离子信号通路参与Ras信号通路；KWMTBOMO02140参与了核糖体信号通路，与丝素蛋白的合成密切相关；KWMTBOMO00853、KWMTBOMO02147参与了蛋白质在内质网上的加工。这些基因是否是调控家蚕茧层率的关键基因及其调控机制将是本课题组下一步的研究重点。
4. 结论
运用BSA-seq方法成功地在家蚕第2染色体、第4染色体和第13 染色体上定位到与茧层率关联的区域；通过KEGG功能注释及相关文献检索，筛选到10个可能与茧层率密切相关的基因，参与了Ras信号通路、细胞周期通路、MAPK信号通路、胰岛素信号通路、钙离子信号通路、核糖体信号通路和蛋白质在内质网上的加工过程。研究结果为高茧层率相关调控基因的精细定位奠定了基础。
-
![]()
图 1 分子标记及关联信号在染色体上的分布
从外到内依次为参考基因组染色体坐标、染色体上基因分布（颜色越深表示基因密度越大）、SNP密度分布（圆点越密集表示SNP密度越大）、InDel密度分布（三角形越密集表示InDel密度越大）、Index值在染色体上的分布。
Figure 1. Distribution of SNPs, InDels, and associated signals on chromosome
Shown from outside inward: chromosome coordinates of reference genome, genes distribution on chromosome (darker color indicates greater gene density), SNP density distribution (density of dots corresponds to that of SNP), InDel density distribution (density of triangles reflects that of InDel), and distribution of indices on chromosome.
![]()
图 2 H池和L池SNP-index、∆（SNP-index）分布情况
图中不同颜色表示不同的染色体，横坐标为1～28号染色体上每个window的具体物理位置，纵坐标为位置所对应的Index值。
Figure 2. SNP-index and ∆(SNP-index) distribution of H-pool and L-pool.
Different colors represent different chromosomes; x-axis is for physical location of each window on chromosomes 1 to 28; y-axis is for Index corresponding to respective locations.
表 1 测序数据质量
Table 1 Statistics on quality of sequencing data
样品编号
Sample ID原始数据量
Raw data/bp过滤后数据量
Clean base/bp原测序reads 数
Raw reads过滤后reads 数
Clean readsQ20/% Q30/% GC含量
GC content/%菁松（JS） 380002170 6295123605 42251670 41738170 98.10 94.04 38.83 芙蓉（FR） 5853501108 5771738817 38764908 38269934 97.91 93.56 38.96 H-pool 16080222540 15911340150 106491540 105500096 98.13 94.11 38.89 L-pool 16401318000 16232070628 108618000 107624014 98.13 94.11 38.66 Q20：高质量测序数据中质量值≥20的碱基所占百分比；Q30：高质量测序数据中质量值≥30的碱基所占百分比。
Q20：The percentage of the bases whose Phred value are more than 20; Q30：The percentage of the bases whose Phred value are more than 30.
下载: 导出CSV
表 2 质控数据与参考基因组比对情况
Table 2 Matching between quality control data and reference genome
样品编号
Sample ID比对率
Mapped ratio/%比对到基因组上的reads比例
Properly ratio/%平均测序深度
Average depth基因组覆盖度（1×）
Genome coverage(1×) /%基因组覆盖度（5×）
Genome coverage(5×) /%菁松（JS） 98.84 82.72 10.34 94.70 84.71 芙蓉（FR） 98.88 82.28 9.44 94.41 82.61 H-pool 98.90 81.29 26.02 97.04 93.59 L-pool 98.82 82.38 27.59 97.02 93.59 1×覆盖度：1 个碱基覆盖的位点占基因组的百分比；5×覆盖度：5个碱基覆盖的位点占基因组的百分比。
Coverage 1×: the percentage of at least 1 base-covered site in reference genome; Coverage 5×: the percentage of at least 5 base-covered sites in reference genome.
下载: 导出CSV
表 3 关联区域信息统计
Table 3 Statistics of the related genes
染色体编号
Chromosome ID关联区域起点
Start of associated
regions/bp关联区域终点
End of associated
regions/bp关联区域长度
Associated region
size/MbSNP数量
SNP number关联区域内基因个数
Gene number in the
associated regions第13染色体 Chr 13 3230000 3730000 0.50 11 13 第4染色体 Chr 4 12350000 12920000 0.50 39 48 第2染色体 Chr 2 4430000 4930000 0.57 17 9 合计 Total 1.57 67 70
下载: 导出CSV
表 4 候选基因的KEGG通路分析
Table 4 KEGG pathway of genes in candidate regions
一级代谢
Primary metabolism二级代谢
Secondary metabolism三级代谢
Tertiary metabolism通路编号
Ko ID基因编号
Gene ID新陈代谢
Metabolism聚糖生物合成与代谢
Glycan biosynthesis and metabolism糖胺聚糖降解
Glycosaminoglycan degradationko00531 KWMTBOMO07657 O-聚糖生物合成
Other types of O-glycan biosynthesisko00514 KWMTBOMO02149 脂质代谢
Lipid metabolism初级胆汁酸生物合成
Primary bile acid biosynthesisko00120 KWMTBOMO02138 有机体系统
Organismal systems免疫系统
Immune systemRIG-I样受体信号通路
RIG-I-like receptor signaling pathwayKo04622 KWMTBOMO02145；
KWMTBOMO07659；
KWMTBOMO02143Toll样受体信号通路
Toll-like receptor signaling pathwayko04620 KWMTBOMO07659 NOD样受体信号通路
NOD-like receptor signaling pathwayko04621 KWMTBOMO07659 内分泌系统
Endocrine system胰高血糖素信号通路
Glucagon signaling pathwayKo04922 KWMTBOMO07653；
KWMTBOMO07652；
KWMTBOMO02146甲状腺激素信号通路
Thyroid hormone signaling pathwayko04919 KWMTBOMO07653；
KWMTBOMO07652胰岛素信号通路
Insulin signaling pathwayko04910 KWMTBOMO02146 神经系统
Nervous system长时程增强效应
Long-term potentiationko04720 KWMTBOMO07653；
KWMTBOMO07652神经营养因子信号通路
Neurotrophin signaling pathwayko04722 KWMTBOMO07659 环境适应
Environmental adaptation生理节律 Circadian rhythm ko04710 KWMTBOMO00853 环境信息加工
Environment information processing信号转导
Signal transductionRas信号通路 Ras signaling pathway ko04014 KWMTBOMO02152 丝裂原活化蛋白激酶信号通路
MAPK signaling pathwayKo04013 KWMTBOMO02145；
KWMTBOMO07659；
KWMTBOMO02143钙离子信号通路
Calcium signaling pathwayko04020 KWMTBOMO02146 低氧诱导因子-1信号通路
HIF-1 signaling pathwayko04066 KWMTBOMO07653；
KWMTBOMO07652Wnt信号通路
Wnt signaling pathwayko04310 KWMTBOMO07653；
KWMTBOMO07652；
KWMTBOMO00853环磷酸腺苷信号通路
cAMP signaling pathwayko04024 KWMTBOMO07653；
KWMTBOMO07652Notch信号通路 Notch signaling pathway ko04330 KWMTBOMO07653；
KWMTBOMO07652Jak-STAT信号通路
Jak-STAT signaling pathwayko04630 KWMTBOMO07653；
KWMTBOMO07652TGF-β信号通路
TGF-beta signaling pathwayko04350 KWMTBOMO07653；
KWMTBOMO00853；
KWMTBOMO07652刺猬信号通路
Hedgehog signaling pathwayko04341 KWMTBOMO00853 细胞进程
Cellular processes细胞生长和死亡
Cell growth and death细胞周期 Cell cycle Ko04110 KWMTBOMO07653；
KWMTBOMO07652；
KWMTBOMO00853细胞凋亡 Apoptosis ko04214 KWMTBOMO07659 卵母细胞减数分裂 Oocyte meiosis ko04114 KWMTBOMO00853 细胞通讯 Cell communication 黏着连接 Adherens junction ko04520 KWMTBOMO07653；
KWMTBOMO07652运输与分解代谢
Transport and catabolism内吞作用 Endocytosis ko04144 KWMTBOMO07659 溶酶体 Lysosome ko04142 KWMTBOMO07657 遗传信息加工 Genetic information processing 折叠、组装和降解
Folding, sorting and degradation泛素介导的蛋白质水解
Ubiquitin mediated proteolysisko04120 KWMTBOMO07659；
KWMTBOMO02114；
KWMTBOMO00853蛋白质在内质网上的加工
Protein processing in endoplasmic reticulumko04141 KWMTBOMO02147；
KWMTBOMO00853蛋白酶体 Proteasome ko03050 KWMTBOMO02148 转录 Transcription 转录因子 Basal transcription factors ko03022 KWMTBOMO02120 剪接体 Spliceosome ko03040 KWMTBOMO02133 翻译 Translation 核糖体 Ribosome ko03010 KWMTBOMO02140
下载: 导出CSV
-
[1] 任晓晓, 罗朝斌, 孙运朋, 等. 高原蚕区家蚕茧层率遗传分析 [J]. 农学学报, 2020, 10(6):75−80. REN X X, LUO C B, SUN Y P, et al. The cocoon shell ratio of Bombyx mori from sericultural area of plateau: Genetic analysis [J]. Journal of Agriculture, 2020, 10(6): 75−80.(in Chinese)
[2] LI C L, TONG X L, ZUO W D, et al. QTL analysis of cocoon shell weight identifies BmRPL18 associated with silk protein synthesis in silkworm by pooling sequencing [J]. Scientific Reports, 2017, 7: 17985. DOI: 10.1038/s41598-017-18277-y
[3] 刘娜, 李娟, 秦笙, 等. 家蚕茧丝相关性状的研究进展 [J]. 中国蚕业, 2016, 37(4):6−9. DOI: 10.16839/j.cnki.zgcy.2016.04.002 LIU N, LI J, QIN S, et al. Research progress on cocoon silk related traits of silkworm [J]. China Sericulture, 2016, 37(4): 6−9.(in Chinese) DOI: 10.16839/j.cnki.zgcy.2016.04.002
[4] 栾悦, 李春林, 代方银. 家蚕茧丝性状的遗传基础研究 [J]. 蚕学通讯, 2017, 37(1):21−28. DOI: 10.3969/j.issn.1006-0561.2017.01.005 LUAN Y, LI C L, DAI F Y. Basic researches of the genetics of cocoon traits in silkworm [J]. Newsletter of Sericultural Science, 2017, 37(1): 21−28.(in Chinese) DOI: 10.3969/j.issn.1006-0561.2017.01.005
[5] 鲁成, 李斌, 赵爱春, 等. 家蚕重要经济性状的QTL定位研究 [J]. 中国科学(C辑:生命科学), 2004, 34(3):236−242. LU C, LI B, ZHAO A C, et al. QTL mapping of important economic characters of Bombyx mori [J]. Science in China(SerC), 2004, 34(3): 236−242.(in Chinese)
[6] 李斌, 鲁成, 赵爱春, 等. 家蚕全茧量及重要相关经济性状的多重区间作图分析 [J]. 中国农业科学, 2005, 38(7):1474−1479. DOI: 10.3321/j.issn:0578-1752.2005.07.030 LI B, LU C, ZHAO A C, et al. Multiple interval mapping for whole cocoon weight and related economically important traits QTL in silkworm (Bombyx mori) [J]. Scientia Agricultura Sinica, 2005, 38(7): 1474−1479.(in Chinese) DOI: 10.3321/j.issn:0578-1752.2005.07.030
[7] 司马杨虎, 李斌, 徐海明, 等. 家蚕茧质性状的QTL定位研究 [J]. 遗传学报, 2005, 32(6):625−632. SIMA Y H, LI B, XU H M, et al. Study on location of QTLs controlling cocoon traits in silkworm [J]. Acta Genetica Sinica, 2005, 32(6): 625−632.(in Chinese)
[8] 侯成香, 王修业, 李冰, 等. 家蚕茧丝相关性状的性连锁QTLs定位与分析 [J]. 蚕业科学, 2013, 39(1):35−39. DOI: 10.13441/j.cnki.cykx.2013.01.019 HOU C X, WANG X Y, LI B, et al. Mapping and analysis of Bombyx mori sex-linked QTLs related to cocoon and silk traits [J]. Science of Sericulture, 2013, 39(1): 35−39.(in Chinese) DOI: 10.13441/j.cnki.cykx.2013.01.019
[9] LI B, WANG X Y, HOU C X, et al. Genetic analysis of quantitative trait loci for cocoon and silk production quantity in Bombyx mori (Lepidoptera: Bombycidae) [J]. European Journal of Entomology, 2013, 110(2): 205−213. DOI: 10.14411/eje.2013.031
[10] 张之昊, 王俊, 刘章雄, 等. 基于BSA-Seq技术挖掘大豆中黄622的多小叶基因 [J]. 作物学报, 2020, 46(12):1839−1849. ZHANG Z H, WANG J, LIU Z X, et al. Mapping of an incomplete dominant gene controlling multifoliolate leaf by BSA-Seq in soybean(Glycine max L.) [J]. Acta Agronomica Sinica, 2020, 46(12): 1839−1849.(in Chinese)
[11] 贾秀苹, 卯旭辉, 岳云, 等. 利用BSA-Seq方法鉴定向日葵耐盐候选基因 [J]. 中国油料作物学报, 2018, 40(6):777−784. DOI: 10.7505/j.issn.1007-9084.2018.06.006 JIA X P, MAO X H, YUE Y, et al. Identification of major salt-tolerant genes via BSA-Seq method in sunflower [J]. Chinese Journal of Oil Crop Sciences, 2018, 40(6): 777−784.(in Chinese) DOI: 10.7505/j.issn.1007-9084.2018.06.006
[12] 徐剑文, 刘剑光, 赵君, 等. 利用BSA-seq发掘棉花适宜机采的果枝长度相关QTL [J]. 棉花学报, 2019, 31(4):319−326. DOI: 10.11963/1002-7807.xjwxsh.20190611 XU J W, LIU J G, ZHAO J, et al. The identification of QTL associated with cotton fruit branch length suitable for mechanized harvest utilizing BSA-seq [J]. Cotton Science, 2019, 31(4): 319−326.(in Chinese) DOI: 10.11963/1002-7807.xjwxsh.20190611
[13] YANG T T, AMANULLAH S, PAN J H, et al. Identification of putative genetic regions for watermelon rind hardness and related traits by BSA-seq and QTL mapping [J]. Euphytica, 2021, 217(2): 19. DOI: 10.1007/s10681-020-02758-9
[14] 刘梦雨, 刘小丰, 江东, 等. 利用重测序-BSA分析鉴定金柑油胞发育相关基因 [J]. 园艺学报, 2019, 46(5):841−854. LIU M Y, LIU X F, JIANG D, et al. Identification of genes related to oil gland development in kumquat by using BSA-seq [J]. Acta Horticulturae Sinica, 2019, 46(5): 841−854.(in Chinese)
[15] 许芸梅, 李玉梅, 贾玉鑫, 等. 马铃薯红色薯肉调控基因的精细定位与候选基因分析 [J]. 中国农业科学, 2019, 52(15):2678−2685. DOI: 10.3864/j.issn.0578-1752.2019.15.011 XU Y M, LI Y M, JIA Y X, et al. Fine mapping and candidate genes analysis for regulatory gene of anthocyanin synthesis in red-colored Tuber flesh [J]. Scientia Agricultura Sinica, 2019, 52(15): 2678−2685.(in Chinese) DOI: 10.3864/j.issn.0578-1752.2019.15.011
[16] 柳海东, 赵绪涛, 杜德志. 利用QTL-seq技术定位甘蓝型春油菜早花位点cqDTFC8及其近等基因系构建 [J]. 植物生理学报, 2020, 56(2):219−234. DOI: 10.13592/j.cnki.ppj.2019.0398 LIU H D, ZHAO X T, DU D Z. Mapping of early flowering site cqDTFC8 using QTL-seq technique and construction of its near-isogenic lines in Brassica napus [J]. Plant Physiology Journal, 2020, 56(2): 219−234.(in Chinese) DOI: 10.13592/j.cnki.ppj.2019.0398
[17] 尹明智, 胡燕. 基于BSA-seq法的油菜野芥胞质雄性不育恢复基因的分析 [J]. 西北植物学报, 2020, 40(7):1148−1156. YIN M Z, HU Y. Location analysis of restorer gene of Sinapis arvensis cytoplasmic male sterility in Brassica napus based on BSA-seq method [J]. Acta Botanica Boreali-Occidentalia Sinica, 2020, 40(7): 1148−1156.(in Chinese)
[18] 欧点点, 赵光伟, 贺玉花, 等. 甜瓜果皮颜色遗传分析及基因定位 [J]. 中国农学通报, 2019, 35(13):64−69. DOI: 10.11924/j.issn.1000-6850.casb18110127 OU D D, ZHAO G W, HE Y H, et al. Genetic analysis and gene mapping for melon rind color [J]. Chinese Agricultural Science Bulletin, 2019, 35(13): 64−69.(in Chinese) DOI: 10.11924/j.issn.1000-6850.casb18110127
[19] 祝新荣, 何克荣, 柳新菊, 等. 多丝量雄蚕新品种华菁×平72的育成 [J]. 蚕业科学, 2014, 40(2):248−253. DOI: 10.13441/j.cnki.cykx.2014.02.012 ZHU X R, HE K R, LIU X J, et al. Breeding of new male silkworm variety Huajing × Ping 72 with high silk yield [J]. Science of Sericulture, 2014, 40(2): 248−253.(in Chinese) DOI: 10.13441/j.cnki.cykx.2014.02.012
[20] 司马杨虎, 徐海明, 赵爱春, 等. 性别效应对家蚕茧质性状QTL定位的影响 [J]. 蚕业科学, 2009, 35(4):783−789. DOI: 10.3969/j.issn.0257-4799.2009.04.012 SIMA Y H, XU H M, ZHAO A C, et al. Influence of sex-effects on QTL mapping of silkworm cocoon quality traits [J]. Science of Sericulture, 2009, 35(4): 783−789.(in Chinese) DOI: 10.3969/j.issn.0257-4799.2009.04.012
[21] ZHAN S, HUANG J H, GUO Q H, et al. An integrated genetic linkage map for silkworms with three parental combinations and its application to the mapping of single genes and QTL [J]. BMC Genomics, 2009, 10: 389. DOI: 10.1186/1471-2164-10-389
[22] FANG S M, ZHOU Q Z, YU Q Y, et al. Genetic and genomic analysis for cocoon yield traits in silkworm [J]. Scientific Reports, 2020, 10: 5682. DOI: 10.1038/s41598-020-62507-9
[23] 马倩, 马俐, 李胜, 等. 基于RNA-Seq分析Ras1 CA在家蚕后部丝腺过表达对细胞周期通路基因的影响 [J]. 应用昆虫学报, 2015, 52(2):390−399. DOI: 10.7679/j.issn.2095-1353.2015.043 MA Q, MA L, LI S, et al. RNA-Seq technology based transcriptomic analysis of differentially expressed genes in the cell cycle pathway of Ras1 CA-overexpressed and wild type posterior silk glands of Bombyx mori [J]. Chinese Journal of Applied Entomology, 2015, 52(2): 390−399.(in Chinese) DOI: 10.7679/j.issn.2095-1353.2015.043
[24] 张祥乐, 马俐, 马倩, 等. Ras信号通路通过激活转录因子Myc促进核内复制细胞生长 [J]. 昆虫学报, 2018, 61(8):885−894. ZHANG X L, MA L, MA Q, et al. Ras signaling pathway promotes the growth of endoreplication cells through activating the expression of transcription factor Myc [J]. Acta Entomologica Sinica, 2018, 61(8): 885−894.(in Chinese)
[25] CALDWELL P E, WALKIEWICZ M, STERN M. Ras activity in the Drosophila prothoracic gland regulates body size and developmental rate via ecdysone release [J]. Current Biology, 2005, 15(20): 1785−1795. DOI: 10.1016/j.cub.2005.09.011
[26] MA L, XU H F, ZHU J Q, et al. Ras1 CA overexpression in the posterior silk gland improves silk yield [J]. Cell Research, 2011, 21(6): 934−943. DOI: 10.1038/cr.2011.36
计量
- 文章访问数: 686
- HTML全文浏览量: 229
- PDF下载量: 20
