今天给大家分享的文献是柞蚕的基因组文章，文章是由福建农林大学，河南科技大学，南阳师范学院以及河南蚕业科学研究院共同完成的。
这篇文章的highlights:
- 柞蚕雄性成虫个体染色体水平的基因组
- 天蚕蛾科大概在81.5百万年前从蚕蛾科中分化出来
- 染色体分裂促进了柞蚕的染色体核型进化
柞蚕是一种半驯化的养殖性的鳞翅目昆虫，主要取食橡树叶，具有很高的经济和生态价值。本文通过三代PacBio测序以及Hi-C染色体捕获技术组装出最终的基因组大小是720.67 Mb，49条染色体，scaffold N50为13.77Mb，重复序列占60.74%，21,431个蛋白编码基因。比较基因组学分析表明柞蚕在大约30.3百万年前与日本柞蚕从共同祖先中分化出来。染色体的分裂促进了染色体数目的增加。
全基因组复制事件(WGD)一般在植物中较常见，而在昆虫中比较罕见。柞蚕的基因组大小及染色体数目要比家蚕大得多，这可能表明可能存在WGD或者segmental duplication (SD)。
为了避免W染色体重复序列增加基因组的复杂度，选择了雄性幼虫用于基因组测序。
28 ugDNA提取自6龄幼虫的表皮组织，采用标准的酚氯仿抽提法。PacBio Sequel II以及Illumina HiSeq平台用于基因组测序。PacBio产生的clean subreads通过WTDBG v2.1进行组装，使用同一个PacBio数据集进行blasr 比对以及Arrow进行矫正。Illumina的数据用BWA和Pilon来polish and improve genome assembly。基因组组装质量通过BUSCO评估。
同一个体收集的脂肪体用于Hi-C建库测序。10mg的样品在1%的甲醛中交联10min，然后加入0.2M的甘氨酸中断交联反应。3D-DNA pipeline用于矫正contigs。校正后的contigs使用ALLHiC assembler构建二倍体scaffold模型。整个的pipeline可以参考：（https://github.com/tangerzhang/ALLHiC/wiki）。
RepeatModeler用于从头构建重复序列数据集结合Repbase，再使用RepeatMasker以及RepeatProteinMask将重复序列屏蔽掉。
取五个个体的混合组织，去除肠道后，提取polyA-RNA并纯化，反转录成cDNA，用于构建Iso-seq文库。通过PacBio Sequel II平台进行全长转录组测序。
蛋白编码基因通过转录组数据，同源比对以及从头预测三种策略。首先，全长转录组数据使用PASA来比对到基因组中获得基因结构特征。其次，Augustus，SNAP and GeneMark 用于从头预测。
gene sets for Augustus and SNAP training were selected from the complete open reading frames prepared by PASA
用于同源比对的序列来自于日本柞蚕以及家蚕，使用GeneWise识别。最后使用EVM构建非冗余的基因结构数据集。
使用BLAST+将基因序列比对到 SwissProt, NT, and NR databases，获得基因功能信息。 Motifs and domains的识别通过HMMER将】比对到Pfam数据库。再通过Trinotate自动注释流程将基因序列比对到 eggNOG , GOand KEGG。
ncRNAs, rRNAs, miRNAs, and snRNAs 通过INFERNAL 比对到Rfam database。而tRNAs 通过tRNAscan-SE鉴定。
将柞蚕基因组序列与日本柞蚕，家蚕，黑脉金斑蝶，小菜蛾，棉铃虫，果蝇和草地贪夜蛾的基因组序列进行ALL-TO-ALL的BLASTP比对，将结果导入OrthoMCL来寻找物种间的直系和旁系同源基因。单拷贝的直系同源基因根据四倍简并位点通过MUSCLE比对形成每个物种的super-sequence。将这八个物种的super-sequence通过PhyML进行建树，果蝇设为外群。
物种分歧时间的估算使用PAML软件包中的MCMCTREE程序，一些物种的分歧时间可以通过TimeTree数据库中查询得到作为标定化石时间。在基因家族聚类以及物种分歧时间估算之后，使用CAFÉ计算基因家族的收缩和扩张，随后收缩或扩增的基因家族进行基因功能的KEGG和GO功能富集分析。
使用MCScan来计算柞蚕和家蚕染色体的共线性关系。
Whole-genome alignment comparison (WGAC) was additionally carried out to identify SD (length ≥ 1 kb and identity ≥ 90%) using LASTZ alignment to a RepeatMasker masked genome.
柞蚕进行基因组survey后的预估基因组大小为675.84-Mb，杂合度1.0%左右，GC含量37%。组装好的基因组大小要大于预估的基因组，这可能是由于其基因组杂合度过高导致的。通过BUSCO评价，二代测序数据的比对率以及三代全长转录组的比对率可以看出其组装出的基因组质量较高。
下面是具体不同测序平台的数据：
通过High-throughput chromosome conformation capture analysis (Hi-C)计算contigs之间的相互作用频率，最终得到49条染色体。
stronger intrachromosomal contacts than interchromosomal contacts, and the interaction probability of the interchromosome contacts rapidly decreased with increasing linear distance.
基于序列和结构的分析得到60.74%的重复序列，21,431 protein-coding genes，以及non-coding RNAs (ncRNAs)。
- A total of 15,256 gene families were identified in the eight species, and 7,265 core gene families were shared by the three Bombycoidea species; however, A. pernyi contained 438 unique gene families not shared by other species. In addition, 1,622 one-copy orthologs were defined across these eight species.
- MCMCTree further estimated that the Bombycidae lineage diverged from the lineage leading to A. pernyi and A. yamamai ~81.5 mya, and A. pernyi diverged from its sister species A. yamamai approximately ~30.3 mya.
- Every chromosome in A. pernyi only corresponded to one chromosome in B. mori, which was consistent with the situation of dot plot. The limited number of syntenic blocks within A. pernyi genome provides a contrasting signature to the many syntenic blocks within a genome that are detected in a species that has experienced WGD.
