一化性柞蚕雌性个体基因组特征分析

乔娇灵 朱绪伟 段竹君 李 莹 杨新峰 段建平

(1南阳师范学院，河南南阳 473061； 2河南省蚕业科学研究院，河南郑州 450008)

柞蚕(Antheraeapernyi)是我国重要的经济昆虫，很长一段时间，由于缺少可用基因组数据资源，限制了柞蚕分子生物学及遗传改良等科学研究的发展。2014年，由辽宁省蚕业科学研究所、辽宁省农业科学院大连生物技术研究所、吉林省蚕业科学研究院、河南省蚕业科学研究院、沈阳农业大学及黑龙江省蚕业研究所等单位联合完成了柞蚕基因组测序[1]，而真正可使用的柞蚕基因组数据，直到2020年才得已公布[2]。该研究利用一化性柞蚕雄性个体为材料，经三代测序组装及三维基因组辅助染色体挂载，将一化性柞蚕雄性个体基因组组装至染色体级别。一般鳞翅目昆虫性别决定方式为ZW型，ZW型个体表现为雌性，ZZ型个体表现为雄性。目前，柞蚕常染色体和Z染色体已被成功组装，但是W染色体还没有被组装。W染色体是雌性个体最重要的一条性染色体，只有完成雌性个体基因组测序研究工作，才能获得W染色体序列数据，最终建立完整的柞蚕基因组图谱。为此，本研究以一化性柞蚕雌性个体为研究对象，通过二代高通量测序，探讨一化性柞蚕雌性个体基因组特征，以期为柞蚕基因组研究方案的制定奠定基础。

1 材料与方法

1.1 试验材料

1.1.1 供试柞蚕 一化性柞蚕雌性个体取自河南省蚕业科学研究院云阳柞蚕繁育基地，将雌性个体除去中肠内容物，剩余组织剪碎、液氮速冻备用。

1.1.2 主要试剂 苯酚(分析纯)、氯仿(分析纯)、Quibit HS Assay Kit文库检测试剂盒，美国Invitrogen公司产品；蛋白酶K，中国TransGen公司产品；TruSeq®Nano DNA LT Library Preparation Kit DNA测序文库构建试剂盒，美国Illumina公司产品。

1.1.3 主要仪器设备 M220 DNA剪切仪，美国Covaris公司产品；Bioanalyzer 2100生物分析仪，美国Agilent公司产品；Illumina Hiseq X Ten测序仪，美国Illumina公司产品；R740服务器，美国DELL公司产品。

1.2 建库及测序

参考分子克隆实验指南[3]的方法提取柞蚕雌性个体基因组DNA，将该DNA用超声波打碎为400 bp左右的片段，并补平粘性末端，经磷酸化修饰和接头连接后，利用接头引物PCR扩增并富集该测序文库。扩增后的文库，经磁珠纯化后，用Quibit HS Assay Kit文库检测试剂盒及Bioanalyzer 2100生物分析仪检测，确定该文库的质量。

利用TruSeq®Nano DNA LT Library Preparation Kit测序文库构建试剂盒，构建插入片段大小约为400 bp的DNA测序文库，用Bioanalyzer 2100生物分析仪确认该文库的纯度和大小，采用Illumina Hiseq X Ten测序仪，利用双末端测序收集二代原始数据，测序片段读长为150 bp。

1.3 二代数据处理

首先进行数据质控。用HTQC toolkit[4]去除原始数据(raw reads)中的接头和低质量reads，如N含量大于10%的reads和低质量碱基(≤5)占比大于50%的reads。再用FastUniq去除PCR重复[5]，获得质控数据(clean reads)。调用FastQC v0.11.6(https：//www.softpedia.com/get/Science-CAD/FastQC.shtml)对质控数据进行质量展示。以clean reads为输入数据集，设k-mer=17，用jellyfish[6]中的count和histo命令统计k-mer频数，生成k-mer频数表，结合软件gce v1.0.0[7]，估算基因组的大小和杂合度。基因组大小(G)满足公式G=knumber/kdepth，其中knumber和kdepth分别为k-mer个数及期望测序深度。

用fq2fa-filter-merge命令将clean reads转换为fasta数据，导入软件IDBA v1.1.3[8]，利用de Bruijn graph算法，初步从头组装雌性个体基因组，预测柞蚕雌性个体基因组的大概大小。调用R v 3.6.1语言(https：//www.r-project.org/)进行统计作图，以500 bp的窗口大小及250 bp的步长值，滑窗统计基因组GC含量，同时以5 000 bp为窗口统计测序深度分布情况。

2 结果与分析

2.1 柞蚕雌性个体基因组DNA测序

通常，只需要收集50×二代测序数据，就可以满足基因组特征评估的需要。有研究显示，一化性柞蚕雄性个体基因组大小为721 Mb[2]，估计雌性个体基因组不会超过1 Gb。经Illumina测序，我们最终收集了58.2 Gb的raw reads(表1)，预估测序深度达75×左右，远大于50×。再对raw reads进行去接头、去低质量reads、去PCR重复等质控处理，得到56.8 Gb的clean reads，93%的碱基质量值达Q30(表1)。

表1 柞蚕雌性个体二代测序数据统计

采用FastQC展示clean data的质控效果，显示质控后的clean data中A与T的含量、及C与G的含量基本相同，没有碱基不平衡的现象(图1-A和图1-B)，碱基整体质量值在Q30以上(图1-C和图1-D)，说明经质控后的clean data质量可靠，可以进行后续特征分析。

2.2 柞蚕雌性个体基因组杂合度评估

根据已有文献估计柞蚕雌性个体基因组大小不会超过1 Gb，调查基因组特征时选用k-mer=17进行计算，当k-mer=17时，从clean data中可得到14 783 706 360个k-mer。进一步绘制17-mer分布曲线，显示2个峰，1个峰在深度约22×处，另1个峰在深度约11×处(图2-A)，22×处应该是主峰，11×处应该是杂合峰，杂合峰位于主峰二分之一处，且峰形已经很明显，说明雌性个体基因组具有一定的杂合度；主峰后有拖尾现象，暗示雌性个体基因组重复序列含量非常高。仅基于主峰的位置，我们估算一化性柞蚕雌性个体基因组大小在675 Mb左右，实际大小应该比此估算值略大。因为已确认柞蚕雌性个体基因组杂合，为评估其杂合度，进一步选择1个模式物种基因组(拟南芥基因组)，模拟对应测序深度的短片段数据，设杂合度梯度变化，通过k-mer曲线拟合估计该个体基因组杂合度。结果显示，当杂合度为1.25%时，该雌性个体基因组k-mer曲线拟合较好，说明一化性柞蚕雌性个体基因组杂合度为1.25%(图2-B)。

A.雌性个体clean data R1的碱基分布；B.雌性个体clean data R2的碱基分布；C.雌性个体clean data R1的质量分布；D.雌性个体clean data R2的质量分布。图1 柞蚕雌性个体二代测序数据质控结果

A. 17-mer分布曲线；B. k-mer拟合曲线。图2 柞蚕雌性个体基因组k-mer分布曲线

2.3 柞蚕雌性个体基因组GC含量及测序深度分布特征

采用clean reads预组装该雌性个体基因组，结果显示初组装的序列(contig)总长约为676.9 Mb，contig N50长度为520 bp。将所有contig进一步拼装，获得的序列(scaffold)总长约为764.9 Mb，scaffold N50长度为2 409 bp(表2)。预组装的基因组大小，与已经发表的雄性个体基因组大小相差不大，与上述预估的结果基本一致，说明雌性个体基因组大小应该在760 Mb左右。但预组装的基因组序列scaffold条数太多，达852 519条，组装质量较差。

表2 柞蚕雌性个体基因组二代数据预组装

利用上述二代预组装结果，划窗统计柞蚕雌性个体基因组的GC含量，结果显示雌性个体基因组中GC含量约为36%(图3-A)。将clean data回贴至预组装的基因组，划窗统计二代测序数据的深度，结果显示与前述评估一致，有2个峰，第1个应该是杂合峰，第2个应该是主峰(图3-B)，进一步证实该雌性个体基因组杂合度较高。

A.雌性个体基因组GC分布图；B.雌性个体基因组测序深度。图3 柞蚕雌性个体基因组GC分布及测序深度

3 讨论

从人类基因组计划开展至今，测序技术已经更新了三代。目前，以PacBio和Nanopore为代表的第三代测序技术发展，结合染色体构象捕获技术(HiC)，可快速将基因组组装至染色体。染色体级别的基因组组装，是当今基因组计划研究的基本要求。一化性柞蚕雄性个体基因组测序采用了“二代评估、三代组装、二代纠错及HiC辅助挂载”策略，并获得了所有48条常染色体序列和1条Z染色体序列[2]。因此，要将柞蚕雌性个体基因组组装至染色体水平，需要基于二代测序，评估雌性个体基因组特征，为组装策略的制定提供数据支撑。本研究结果显示雌性个体基因组GC含量为36%，与雄性个体相当，但雌性个体基因组杂合度为1.25%，大于雄性个体基因组杂合度的1.00%[2]，暗示雌性个体基因组也属于高复杂度基因组，雌性个体基因组中重复序列含量可能高于雄性个体，单纯二代数据组装，效果可能不理想。我们二代预组装，获得由852 519条scaffold组成的雌性个体基因组数据，也印证了高重复序列含量会影响基因组组装质量，雌性个体基因组组装需要引入三代测序数据。因此，将来开展一化性柞蚕雌性个体基因组测序研究，需要对其基因组DNA进行二代测序和较高深度的三代测序，同时引入新的HiC算法，如ALLHiC[9]，来保证基因组组装的质量，才可能获得W染色体数据，最终形成完整的柞蚕基因组图谱。