喬嬌靈 朱緒偉 段竹君 李 瑩 楊新峰 段建平

(1南陽師范學(xué)院，河南南陽 473061； 2河南省蠶業(yè)科學(xué)研究院，河南鄭州 450008)

柞蠶(Antheraeapernyi)是我國重要的經(jīng)濟(jì)昆蟲，很長一段時(shí)間，由于缺少可用基因組數(shù)據(jù)資源，限制了柞蠶分子生物學(xué)及遺傳改良等科學(xué)研究的發(fā)展。2014年，由遼寧省蠶業(yè)科學(xué)研究所、遼寧省農(nóng)業(yè)科學(xué)院大連生物技術(shù)研究所、吉林省蠶業(yè)科學(xué)研究院、河南省蠶業(yè)科學(xué)研究院、沈陽農(nóng)業(yè)大學(xué)及黑龍江省蠶業(yè)研究所等單位聯(lián)合完成了柞蠶基因組測序[1]，而真正可使用的柞蠶基因組數(shù)據(jù)，直到2020年才得已公布[2]。該研究利用一化性柞蠶雄性個(gè)體為材料，經(jīng)三代測序組裝及三維基因組輔助染色體掛載，將一化性柞蠶雄性個(gè)體基因組組裝至染色體級(jí)別。一般鱗翅目昆蟲性別決定方式為ZW型，ZW型個(gè)體表現(xiàn)為雌性，ZZ型個(gè)體表現(xiàn)為雄性。目前，柞蠶常染色體和Z染色體已被成功組裝，但是W染色體還沒有被組裝。W染色體是雌性個(gè)體最重要的一條性染色體，只有完成雌性個(gè)體基因組測序研究工作，才能獲得W染色體序列數(shù)據(jù)，最終建立完整的柞蠶基因組圖譜。為此，本研究以一化性柞蠶雌性個(gè)體為研究對象，通過二代高通量測序，探討一化性柞蠶雌性個(gè)體基因組特征，以期為柞蠶基因組研究方案的制定奠定基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料

1.1.1 供試柞蠶 一化性柞蠶雌性個(gè)體取自河南省蠶業(yè)科學(xué)研究院云陽柞蠶繁育基地，將雌性個(gè)體除去中腸內(nèi)容物，剩余組織剪碎、液氮速凍備用。

1.1.2 主要試劑 苯酚(分析純)、氯仿(分析純)、Quibit HS Assay Kit文庫檢測試劑盒，美國Invitrogen公司產(chǎn)品；蛋白酶K，中國TransGen公司產(chǎn)品；TruSeq?Nano DNA LT Library Preparation Kit DNA測序文庫構(gòu)建試劑盒，美國Illumina公司產(chǎn)品。

1.1.3 主要儀器設(shè)備 M220 DNA剪切儀，美國Covaris公司產(chǎn)品；Bioanalyzer 2100生物分析儀，美國Agilent公司產(chǎn)品；Illumina Hiseq X Ten測序儀，美國Illumina公司產(chǎn)品；R740服務(wù)器，美國DELL公司產(chǎn)品。

1.2 建庫及測序

參考分子克隆實(shí)驗(yàn)指南[3]的方法提取柞蠶雌性個(gè)體基因組DNA，將該DNA用超聲波打碎為400 bp左右的片段，并補(bǔ)平粘性末端，經(jīng)磷酸化修飾和接頭連接后，利用接頭引物PCR擴(kuò)增并富集該測序文庫。擴(kuò)增后的文庫，經(jīng)磁珠純化后，用Quibit HS Assay Kit文庫檢測試劑盒及Bioanalyzer 2100生物分析儀檢測，確定該文庫的質(zhì)量。

利用TruSeq?Nano DNA LT Library Preparation Kit測序文庫構(gòu)建試劑盒，構(gòu)建插入片段大小約為400 bp的DNA測序文庫，用Bioanalyzer 2100生物分析儀確認(rèn)該文庫的純度和大小，采用Illumina Hiseq X Ten測序儀，利用雙末端測序收集二代原始數(shù)據(jù)，測序片段讀長為150 bp。

1.3 二代數(shù)據(jù)處理

首先進(jìn)行數(shù)據(jù)質(zhì)控。用HTQC toolkit[4]去除原始數(shù)據(jù)(raw reads)中的接頭和低質(zhì)量reads，如N含量大于10%的reads和低質(zhì)量堿基(≤5)占比大于50%的reads。再用FastUniq去除PCR重復(fù)[5]，獲得質(zhì)控?cái)?shù)據(jù)(clean reads)。調(diào)用FastQC v0.11.6(https：//www.softpedia.com/get/Science-CAD/FastQC.shtml)對質(zhì)控?cái)?shù)據(jù)進(jìn)行質(zhì)量展示。以clean reads為輸入數(shù)據(jù)集，設(shè)k-mer=17，用jellyfish[6]中的count和histo命令統(tǒng)計(jì)k-mer頻數(shù)，生成k-mer頻數(shù)表，結(jié)合軟件gce v1.0.0[7]，估算基因組的大小和雜合度?；蚪M大小(G)滿足公式G=knumber/kdepth，其中knumber和kdepth分別為k-mer個(gè)數(shù)及期望測序深度。

用fq2fa-filter-merge命令將clean reads轉(zhuǎn)換為fasta數(shù)據(jù)，導(dǎo)入軟件IDBA v1.1.3[8]，利用de Bruijn graph算法，初步從頭組裝雌性個(gè)體基因組，預(yù)測柞蠶雌性個(gè)體基因組的大概大小。調(diào)用R v 3.6.1語言(https：//www.r-project.org/)進(jìn)行統(tǒng)計(jì)作圖，以500 bp的窗口大小及250 bp的步長值，滑窗統(tǒng)計(jì)基因組GC含量，同時(shí)以5 000 bp為窗口統(tǒng)計(jì)測序深度分布情況。

2 結(jié)果與分析

2.1 柞蠶雌性個(gè)體基因組DNA測序

通常，只需要收集50×二代測序數(shù)據(jù)，就可以滿足基因組特征評估的需要。有研究顯示，一化性柞蠶雄性個(gè)體基因組大小為721 Mb[2]，估計(jì)雌性個(gè)體基因組不會(huì)超過1 Gb。經(jīng)Illumina測序，我們最終收集了58.2 Gb的raw reads(表1)，預(yù)估測序深度達(dá)75×左右，遠(yuǎn)大于50×。再對raw reads進(jìn)行去接頭、去低質(zhì)量reads、去PCR重復(fù)等質(zhì)控處理，得到56.8 Gb的clean reads，93%的堿基質(zhì)量值達(dá)Q30(表1)。

表1 柞蠶雌性個(gè)體二代測序數(shù)據(jù)統(tǒng)計(jì)

采用FastQC展示clean data的質(zhì)控效果，顯示質(zhì)控后的clean data中A與T的含量、及C與G的含量基本相同，沒有堿基不平衡的現(xiàn)象(圖1-A和圖1-B)，堿基整體質(zhì)量值在Q30以上(圖1-C和圖1-D)，說明經(jīng)質(zhì)控后的clean data質(zhì)量可靠，可以進(jìn)行后續(xù)特征分析。

2.2 柞蠶雌性個(gè)體基因組雜合度評估

根據(jù)已有文獻(xiàn)估計(jì)柞蠶雌性個(gè)體基因組大小不會(huì)超過1 Gb，調(diào)查基因組特征時(shí)選用k-mer=17進(jìn)行計(jì)算，當(dāng)k-mer=17時(shí)，從clean data中可得到14 783 706 360個(gè)k-mer。進(jìn)一步繪制17-mer分布曲線，顯示2個(gè)峰，1個(gè)峰在深度約22×處，另1個(gè)峰在深度約11×處(圖2-A)，22×處應(yīng)該是主峰，11×處應(yīng)該是雜合峰，雜合峰位于主峰二分之一處，且峰形已經(jīng)很明顯，說明雌性個(gè)體基因組具有一定的雜合度；主峰后有拖尾現(xiàn)象，暗示雌性個(gè)體基因組重復(fù)序列含量非常高。僅基于主峰的位置，我們估算一化性柞蠶雌性個(gè)體基因組大小在675 Mb左右，實(shí)際大小應(yīng)該比此估算值略大。因?yàn)橐汛_認(rèn)柞蠶雌性個(gè)體基因組雜合，為評估其雜合度，進(jìn)一步選擇1個(gè)模式物種基因組(擬南芥基因組)，模擬對應(yīng)測序深度的短片段數(shù)據(jù)，設(shè)雜合度梯度變化，通過k-mer曲線擬合估計(jì)該個(gè)體基因組雜合度。結(jié)果顯示，當(dāng)雜合度為1.25%時(shí)，該雌性個(gè)體基因組k-mer曲線擬合較好，說明一化性柞蠶雌性個(gè)體基因組雜合度為1.25%(圖2-B)。

A.雌性個(gè)體clean data R1的堿基分布；B.雌性個(gè)體clean data R2的堿基分布；C.雌性個(gè)體clean data R1的質(zhì)量分布；D.雌性個(gè)體clean data R2的質(zhì)量分布。圖1 柞蠶雌性個(gè)體二代測序數(shù)據(jù)質(zhì)控結(jié)果

A. 17-mer分布曲線；B. k-mer擬合曲線。圖2 柞蠶雌性個(gè)體基因組k-mer分布曲線

2.3 柞蠶雌性個(gè)體基因組GC含量及測序深度分布特征

采用clean reads預(yù)組裝該雌性個(gè)體基因組，結(jié)果顯示初組裝的序列(contig)總長約為676.9 Mb，contig N50長度為520 bp。將所有contig進(jìn)一步拼裝，獲得的序列(scaffold)總長約為764.9 Mb，scaffold N50長度為2 409 bp(表2)。預(yù)組裝的基因組大小，與已經(jīng)發(fā)表的雄性個(gè)體基因組大小相差不大，與上述預(yù)估的結(jié)果基本一致，說明雌性個(gè)體基因組大小應(yīng)該在760 Mb左右。但預(yù)組裝的基因組序列scaffold條數(shù)太多，達(dá)852 519條，組裝質(zhì)量較差。

表2 柞蠶雌性個(gè)體基因組二代數(shù)據(jù)預(yù)組裝

利用上述二代預(yù)組裝結(jié)果，劃窗統(tǒng)計(jì)柞蠶雌性個(gè)體基因組的GC含量，結(jié)果顯示雌性個(gè)體基因組中GC含量約為36%(圖3-A)。將clean data回貼至預(yù)組裝的基因組，劃窗統(tǒng)計(jì)二代測序數(shù)據(jù)的深度，結(jié)果顯示與前述評估一致，有2個(gè)峰，第1個(gè)應(yīng)該是雜合峰，第2個(gè)應(yīng)該是主峰(圖3-B)，進(jìn)一步證實(shí)該雌性個(gè)體基因組雜合度較高。

A.雌性個(gè)體基因組GC分布圖；B.雌性個(gè)體基因組測序深度。圖3 柞蠶雌性個(gè)體基因組GC分布及測序深度

3 討論

從人類基因組計(jì)劃開展至今，測序技術(shù)已經(jīng)更新了三代。目前，以PacBio和Nanopore為代表的第三代測序技術(shù)發(fā)展，結(jié)合染色體構(gòu)象捕獲技術(shù)(HiC)，可快速將基因組組裝至染色體。染色體級(jí)別的基因組組裝，是當(dāng)今基因組計(jì)劃研究的基本要求。一化性柞蠶雄性個(gè)體基因組測序采用了“二代評估、三代組裝、二代糾錯(cuò)及HiC輔助掛載”策略，并獲得了所有48條常染色體序列和1條Z染色體序列[2]。因此，要將柞蠶雌性個(gè)體基因組組裝至染色體水平，需要基于二代測序，評估雌性個(gè)體基因組特征，為組裝策略的制定提供數(shù)據(jù)支撐。本研究結(jié)果顯示雌性個(gè)體基因組GC含量為36%，與雄性個(gè)體相當(dāng)，但雌性個(gè)體基因組雜合度為1.25%，大于雄性個(gè)體基因組雜合度的1.00%[2]，暗示雌性個(gè)體基因組也屬于高復(fù)雜度基因組，雌性個(gè)體基因組中重復(fù)序列含量可能高于雄性個(gè)體，單純二代數(shù)據(jù)組裝，效果可能不理想。我們二代預(yù)組裝，獲得由852 519條scaffold組成的雌性個(gè)體基因組數(shù)據(jù)，也印證了高重復(fù)序列含量會(huì)影響基因組組裝質(zhì)量，雌性個(gè)體基因組組裝需要引入三代測序數(shù)據(jù)。因此，將來開展一化性柞蠶雌性個(gè)體基因組測序研究，需要對其基因組DNA進(jìn)行二代測序和較高深度的三代測序，同時(shí)引入新的HiC算法，如ALLHiC[9]，來保證基因組組裝的質(zhì)量，才可能獲得W染色體數(shù)據(jù)，最終形成完整的柞蠶基因組圖譜。