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]利用甜瓜黄皮材料和绿皮材料构建了F₂分离群体，利用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的功能验证还有待进一步研究。【拟解决的关键问题】本研究以多丝量家蚕和中丝量家蚕品种为亲本，以构建的BC₁群体为研究对象，采用BSA-seq方法对茧层率相关基因进行定位分析，进一步挖掘茧层率关键基因，为高茧层率关键基因克隆及家蚕经济性状分子遗传育种实践提供理论依据。

1.   材料与方法

1.1   试验材料

重测序材料选取多丝量家蚕品种菁松（高茧层率）和中丝量家蚕品种芙蓉（中茧层率）构建的BC₁群体，亲本材料由贵州省蚕业研究所提供，均经多代自交，背景纯合。2019年春，在贵州遵义以菁松、芙蓉杂交获得F₁代。2019年夏，在贵州贵阳饲养F₁代，并以F₁代雄蛾回交芙蓉，获得BC₁分离群体 芙蓉×（芙蓉×菁松） 。2020年春，在遵义基地同室、同条件饲养亲本和BC₁群体。

1.2   茧质性状调查

在上蔟结茧的第7天采集亲本和BC₁群体的蚕茧，削茧以鉴别雌、雄蛹。因家蚕茧质性状在雌、雄个体间差异显著，为排除试验误差，本试验在亲本中随机选择雌、雄茧各20粒，在BC₁群体中选择所有雄性个体蚕茧（BC₁M），调查记录个体的全茧量、茧层量，计算个体茧层率。

1.3   基因组DNA提取及混池构建

根据BC₁M群体的茧层率调查结果，从中选择高茧层率和低茧层率极端分离材料各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，Q₂₀≥97.91%，Q₃₀≥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 reads	Q20/%	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 

| 显示表格

根据样品数据与家蚕参考基因组比对结果，利用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 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.

下载: 全尺寸图片 幻灯片

表  3  关联区域信息统计

Table  3.  Statistics of the related genes

染色体编号 Chromosome ID	关联区域起点 Start of associated regions/bp	关联区域终点 End of associated regions/bp	关联区域长度 Associated region size/Mb	SNP数量 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 

| 显示表格

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）。

图  3  候选基因的GO注释

Figure  3.  GO annotation of candidate genes

下载: 全尺寸图片 幻灯片

为进一步了解候选区域内基因的生物学功能，对富集到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 degradation	ko00531	KWMTBOMO07657
		O-聚糖生物合成 Other types of O-glycan biosynthesis	ko00514	KWMTBOMO02149
	脂质代谢 Lipid metabolism	初级胆汁酸生物合成 Primary bile acid biosynthesis	ko00120	KWMTBOMO02138
有机体系统 Organismal systems	免疫系统 Immune system	RIG-I样受体信号通路 RIG-I-like receptor signaling pathway	Ko04622	KWMTBOMO02145； KWMTBOMO07659； KWMTBOMO02143
		Toll样受体信号通路 Toll-like receptor signaling pathway	ko04620	KWMTBOMO07659
		NOD样受体信号通路 NOD-like receptor signaling pathway	ko04621	KWMTBOMO07659
	内分泌系统 Endocrine system	胰高血糖素信号通路 Glucagon signaling pathway	Ko04922	KWMTBOMO07653； KWMTBOMO07652； KWMTBOMO02146
		甲状腺激素信号通路 Thyroid hormone signaling pathway	ko04919	KWMTBOMO07653； KWMTBOMO07652
		胰岛素信号通路 Insulin signaling pathway	ko04910	KWMTBOMO02146
	神经系统 Nervous system	长时程增强效应 Long-term potentiation	ko04720	KWMTBOMO07653； KWMTBOMO07652
		神经营养因子信号通路 Neurotrophin signaling pathway	ko04722	KWMTBOMO07659
	环境适应 Environmental adaptation	生理节律 Circadian rhythm	ko04710	KWMTBOMO00853
环境信息加工 Environment information processing	信号转导 Signal transduction	Ras信号通路 Ras signaling pathway	ko04014	KWMTBOMO02152
		丝裂原活化蛋白激酶信号通路 MAPK signaling pathway	Ko04013	KWMTBOMO02145； KWMTBOMO07659； KWMTBOMO02143
		钙离子信号通路 Calcium signaling pathway	ko04020	KWMTBOMO02146
		低氧诱导因子-1信号通路 HIF-1 signaling pathway	ko04066	KWMTBOMO07653； KWMTBOMO07652
		Wnt信号通路 Wnt signaling pathway	ko04310	KWMTBOMO07653； KWMTBOMO07652； KWMTBOMO00853
		环磷酸腺苷信号通路 cAMP signaling pathway	ko04024	KWMTBOMO07653； KWMTBOMO07652
		Notch信号通路 Notch signaling pathway	ko04330	KWMTBOMO07653； KWMTBOMO07652
		Jak-STAT信号通路 Jak-STAT signaling pathway	ko04630	KWMTBOMO07653； KWMTBOMO07652
		TGF-β信号通路 TGF-beta signaling pathway	ko04350	KWMTBOMO07653； KWMTBOMO00853； KWMTBOMO07652
		刺猬信号通路 Hedgehog signaling pathway	ko04341	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 proteolysis	ko04120	KWMTBOMO07659； KWMTBOMO02114； KWMTBOMO00853
		蛋白质在内质网上的加工 Protein processing in endoplasmic reticulum	ko04141	KWMTBOMO02147； KWMTBOMO00853
		蛋白酶体 Proteasome	ko03050	KWMTBOMO02148
	转录 Transcription	转录因子 Basal transcription factors	ko03022	KWMTBOMO02120
		剪接体 Spliceosome	ko03040	KWMTBOMO02133
	翻译 Translation	核糖体 Ribosome	ko03010	KWMTBOMO02140

下载: 导出CSV 

| 显示表格

3.   讨论

家蚕的部分经济性状具有性别差异，主要表现在雄蚕比雌蚕体质强健、茧层率和出丝率高，雌蚕比雄蚕全茧量和茧层量高、蛹体重大等，因此学者们推测家蚕性染色体上可能存在控制茧层率、全茧量、茧层量、蛹体重等茧质性状的基因^[19-20]。Zhan等^[21]以菁松和兰10为亲本构建回交一代群体（BC₁M），利用SSR分子标记对家蚕部分茧质性状进行了QTL初步定位，结果在性染色体上检测到了与全茧量、茧层量和蛹体重相关的QTLs，但未检测到与茧层率相关的基因位点。侯成香等^[8]在此基础上，同样以菁松和兰10为亲本配置BC₁M群体，利用SSR和SNP标记构建家蚕性染色体分子遗传图谱，对茧、丝有关的QTLs进行精细扫描，但同样未扫描到与茧层率相关的QTLs。本研究以菁松和芙蓉为亲本构建BC₁M群体，利用BSA-Seq分析方法对控制茧层率的基因进行定位，结果将茧层率相关的QTLs定位在第2染色体（4430～4930 kb）、第4染色体（12350～12920 kb）和第13染色体（3230～3730 kb），也未在性染色体上检测到与茧层率相关的基因位点。鉴于茧层率是茧层量对全茧量的比率，主要由茧层量和蛹体重决定，而前人研究表明性染色体上存在控制茧层量与蛹体重的基因，因此分析认为，茧层率与性别关联可能是由于其决定因素（茧层量、蛹体重）与性别关联。

有关茧层率的数量性状位点定位研究早有报道，鲁成等^[5]利用大造和C100杂交的BC₁群体，构建了AFLP分子标记连锁图谱，通过复合区间作图将茧层率相关的QTLs定位于第2、11、15连锁群上；李斌等^[6]利用大造和C100的BC₁群体构建AFLP分子标记连锁图谱，将茧层率相关的QTLs定位在第2、4、14、15、19、25连锁群上；司马杨虎等^[7]对湘晖和872构建的F₂群体构建了AFLP分子标记连锁图谱，将茧层率性状的QTLs定位在第2、11、13、18连锁群上；Zhan等^[21]利用SSR分子标记对菁松和兰10的BC₁M群体进行QTL位点检测，将茧层率相关QTLs定位在第18和19连锁群；Li等^[9]利用菁松和兰10的BC₁M群体，使用STS、SSR和SNP分子标记构建遗传图谱，结果未扫描到与茧层率相关的QTLs。Fang等^[22]利用夏芳和野蚕的BC₁M群体构建了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技术在家蚕后部丝腺中特异性过表达Ras1^CA激活了Ras信号，促进了后部丝腺细胞和细胞核体积增大，全茧量有所增加。马倩等^[23]基于RNA-Seq技术比较了野生型与Ras1^CA过表达转基因家蚕后部丝腺组织中的差异基因，结果发现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信号通路、胰岛素信号通路、钙离子信号通路、核糖体信号通路和蛋白质在内质网上的加工过程。研究结果为高茧层率相关调控基因的精细定位奠定了基础。