任曉曉，孫運朋，卿 卓，楊萬軍，羅朝斌

（貴州省農業(yè)科學院蠶業(yè)研究所，貴州 貴陽 550006）

0 引言

【研究意義】蠶絲是家蠶利用最廣泛的產物，不僅可用于制作絲綢制品，而且隨著絲素蛋白和絲膠蛋白新用途的開發(fā)，蠶絲在生物醫(yī)學、工程材料、食品和化妝品領域也具有較好的應用前景[1]。繭層率是反映家蠶對絲素蛋白、絲膠蛋白合成能力的重要性狀，因此挖掘繭層率關鍵基因，對家蠶分子遺傳育種改良和提高家蠶飼養(yǎng)經濟效益具有重要意義[2]?！厩叭搜芯窟M展】家蠶飼養(yǎng)至今，繭層率性狀經過兩次較大的改良，繭層率提高至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相比構建分子遺傳圖譜具有測序群體小、成本和測序要求低等優(yōu)勢，已經在油料、水果等作物的農藝性狀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～880 000 bp區(qū)域，共篩選出16個基因與育性恢復有關。歐點點等[18]利用甜瓜黃皮材料和綠皮材料構建了F2分離群體，利用BSA-Seq法將皮色相關基因定位在第4染色體的0.02～5.7 Mb區(qū)間和第10染色體的0.08～9.54 Mb區(qū)間。BSA-Seq在家蠶繭質性狀主效基因定位上也有應用的報道，Li等[2]基于BSA-Seq方法將與絲產量連鎖的標記定位在第11、22和23號染色體上，結合候選基因表達模式、絲腺表達分析，篩選出1個調控家蠶繭層量的關鍵基因，命名為BmRPL18?！颈狙芯壳腥朦c】目前有關繭層率相關QTLs的功能驗證還有待進一步研究?！緮M解決的關鍵問題】本研究以多絲量家蠶和中絲量家蠶品種為親本，以構建的BC1群體為研究對象，采用BSA-seq方法對繭層率相關基因進行定位分析，進一步挖掘繭層率關鍵基因，為高繭層率關鍵基因克隆及家蠶經濟性狀分子遺傳育種實踐提供理論依據。

1 材料與方法

1.1 試驗材料

重測序材料選取多絲量家蠶品種菁松（高繭層率）和中絲量家蠶品種芙蓉（中繭層率）構建的BC1群體，親本材料由貴州省蠶業(yè)研究所提供，均經多代自交，背景純合。2019年春，在貴州遵義以菁松、芙蓉雜交獲得F1代。2019年夏，在貴州貴陽飼養(yǎng)F1代，并以F1代雄蛾回交芙蓉，獲得BC1分離群體 芙蓉×（芙蓉×菁松） 。2020年春，在遵義基地同室、同條件飼養(yǎng)親本和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委托給上海美吉生物醫(yī)藥科技有限公司進行建庫測序，測序平臺為Illumina Novaseq 6000。親本池測序深度為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）；進行1 000次置換檢驗，選擇99.99%置信水平為篩選閾值，將Δ（SNP-index）在置信水平之上的窗口作為候選區(qū)域。

1.6 候選區(qū)域基因的功能注釋

應用BLAST軟件對候選區(qū)域的編碼基因進行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數 在 38 269 934～ 107 624 014， 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分布正常，建庫測序成功；質控數據與家蠶參考基因組同源性較高，比對結果正常，可用于后續(xù)變異檢測與定位分析。

表1 測序數據質量Table 1 Statistics on quality of sequencing data

表2 質控數據與參考基因組比對情況Table 2 Matching between quality control data and reference genome

根據樣品數據與家蠶參考基因組比對結果，利用GATK軟件進行SNP變異位點檢測，共獲得26 557 646個SNP標記。利用SnpEff軟件進行變異位點功能預測，結果表明SNP位點主要位于基因間區(qū)、內含子區(qū)域、基因上游和基因下游，非編碼區(qū)域的多態(tài)性位點顯著多于編碼區(qū)。

2.2 連鎖分析

根據SNP過濾原則，對變異檢測獲得的位點進行過濾，最終得到374 463個高質量SNP位點，使用Circos軟件對變異位點在染色體上的分布密度進行可視化（圖1）。根據兩子代混池的SNP-index進行連鎖分析，結果顯示Δ（SNP-index）在置信水平之上的區(qū)域有3個，分別位于第2、4、13染色體上（圖2），總長度為 1.57 Mb，共包含 67個 SNP位點、11個InDel位點、70個編碼基因（表3）。

表3 關聯區(qū)域信息統(tǒng)計Table 3 Statistics of the related genes

2.3 候選區(qū)域基因注釋及候選基因篩選

利用BLAST軟件對關聯區(qū)域內的編碼基因進行多個數據庫（NR、GO、KEGG、UniProt、EggNOG）注釋，結果顯示：70個編碼基因中有69個注釋到NR數據庫中，69個注釋到UniProt數據庫中，58個注釋到GO數據庫中，19個注釋到KEGG通路，57個注釋到EggNOG數據庫。

通過對GO數據庫中候選區(qū)域基因參與的生物過程（Biological process）、分子功能（Molecular function）和細胞組分（Cellular component）進行分類分析，結果表明，在生物過程中候選基因主要富集在上皮細胞發(fā)育、腦室系統(tǒng)發(fā)育和氣管發(fā)育等；在細胞組分中，候選基因多富集于細胞核、細胞質、細胞膜；在分子功能中，候選基因富集最多的是組蛋白乙酰轉移酶（圖3）。

為進一步了解候選區(qū)域內基因的生物學功能，對富集到KEGG數據庫中的基因進行分析，結果顯示候選區(qū)域中19個基因分布于34個信號通路中，涉及新陳代謝（Metabolism）、環(huán)境信息加工（Environment information processing）、 細 胞 進 程（Cellular processes）、 遺 傳 信 息 加 工 （Genetic information processing）、有機體系統(tǒng)（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

續(xù)上表

3 討論

家蠶的部分經濟性狀具有性別差異，主要表現在雄蠶比雌蠶體質強健、繭層率和出絲率高，雌蠶比雄蠶全繭量和繭層量高、蛹體重大等，因此學者們推測家蠶性染色體上可能存在控制繭層率、全繭量、繭層量、蛹體重等繭質性狀的基因[19-20]。Zhan等[21]以菁松和蘭10為親本構建回交一代群體（BC1M），利用SSR分子標記對家蠶部分繭質性狀進行了QTL初步定位，結果在性染色體上檢測到了與全繭量、繭層量和蛹體重相關的QTLs，但未檢測到與繭層率相關的基因位點。侯成香等[8]在此基礎上，同樣以菁松和蘭10為親本配置BC1M群體，利用SSR和SNP標記構建家蠶性染色體分子遺傳圖譜，對繭、絲有關的QTLs進行精細掃描，但同樣未掃描到與繭層率相關的QTLs。本研究以菁松和芙蓉為親本構建BC1M群體，利用BSA-Seq分析方法對控制繭層率的基因進行定位，結果將繭層率相關的QTLs定位在第2染色體（4 430～4 930 kb）、第4染色體（12 350～12 920 kb）和第13染色體（3 230～3 730 kb），也未在性染色體上檢測到與繭層率相關的基因位點。鑒于繭層率是繭層量對全繭量的比率，主要由繭層量和蛹體重決定，而前人研究表明性染色體上存在控制繭層量與蛹體重的基因，因此分析認為，繭層率與性別關聯可能是由于其決定因素（繭層量、蛹體重）與性別關聯。

有關繭層率的數量性狀位點定位研究早有報道，魯成等[5]利用大造和C100雜交的BC1群體，構建了AFLP分子標記連鎖圖譜，通過復合區(qū)間作圖將繭層率相關的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期的主要調節(jié)因子，與其他核內周期調節(jié)者構成復雜的調控網絡，精密協(xié)作保證核內周期的順利完成[24]。前人研究發(fā)現，果蠅Ras信號（Ras85D）可能通過調控CycE、CDK2等核內周期調節(jié)因子，參與調控細胞的存活、生長等生理過程。Caldwell等[25]通過激活果蠅前胸腺細胞（由核內復制細胞組成）中的Ras信號，使前胸腺細胞體積明顯增大。Ma等[26]利用GAL4/UAS技術在家蠶后部絲腺中特異性過表達Ras1CA激活了Ras信號，促進了后部絲腺細胞和細胞核體積增大，全繭量有所增加。馬倩等[23]基于RNASeq技術比較了野生型與Ras1CA過表達轉基因家蠶后部絲腺組織中的差異基因，結果發(fā)現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 染色體上定位到與繭層率關聯的區(qū)域；通過KEGG功能注釋及相關文獻檢索，篩選到10個可能與繭層率密切相關的基因，參與了Ras信號通路、細胞周期通路、MAPK信號通路、胰島素信號通路、鈣離子信號通路、核糖體信號通路和蛋白質在內質網上的加工過程。研究結果為高繭層率相關調控基因的精細定位奠定了基礎。