301种果蝇物种的比较基因注释与直系同源基因研究

《PLOS Biology》：Comparative gene annotation and orthology assignments across 301 species of Drosophilidae

本研究提供了一个覆盖301种果蝇物种的大规模、高一致性的蛋白质编码基因注释与直系同源基因（HOG）数据集，并通过系统发育分析揭示了基因组特征（如基因数目、CDS长度）与技术因素（如测序类型、组装连续性）的关系，同时展示了该数据集在分析全基因组GC含量、密码子使用偏性及氨基酸组成等进化问题上的应用潜力，为进化与功能基因组学的比较研究提供了宝贵的资源。

　　

引言

比较基因组学中的核心问题，如新基因起源、基因组大小变异的原因以及基因内容与结构的进化驱动力，其研究都依赖于高质量的基因组注释。本文作者团队针对果蝇科（Drosophilidae）的301个物种，利用比较注释工具包（Comparative Annotation Toolkit, CAT）和BRAKER3，结合可用的RNA-seq和蛋白质证据，生成了蛋白质编码基因注释。本研究采用系统发育比较的方法进行注释，旨在提高一致性和准确性，并生成一套稳健的基因注释与直系同源基因（Orthology）推断结果。作为大型协作计划的一部分，此数据将持续更新并公开，旨在为整个果蝇科的系统生物学研究奠定基础。

方法

本研究方法主要包括基因组选择与处理、注释流程、正交群推断和进化分析。

基因组组装

研究人员从RefSeq数据库、达尔文生命之树计划及多个独立实验室收集了截至2024年2月公开的果蝇科基因组数据。筛选标准包括支架N50大于50 Kbp且BUSCO完整性得分高于90%，最终选定了301个基因组。使用RepeatMasker和Dfam重复序列库进行软屏蔽，为所有后续分析提供基础。为了进行比较，还基于更严格的标准（如BUSCO完整性≥97%、contig N50≥2 Mbp且使用长读长测序技术）筛选出一个包含215个物种的高质量子集。

RNA测序与蛋白质数据

为生成基因注释，使用了CAT流程。首先，通过ENA Portal API为每个物种收集了RNA-seq数据，选择最多10个配对端数据集，并优先选择经Poly-A富集mRNA的数据。将选定的RNA-seq读数归一化至100倍覆盖率并与各自基因组比对。同时，从节肢动物OrthoDB v10蛋白数据库中提取预测蛋白质序列，以提供进化保守的证据。使用miniprot将蛋白质序列比对到基因组上，并转换为提示文件。

参考物种与Cactus比对

除了RNA-seq和蛋白质提示外，CAT流程还尝试进行“转移”注释。研究人员根据基因组完整性和RefSeq注释质量选择了37个物种作为“参考”。通过BUSCO位点构建初步物种树，并使用ETE3 Python包确定每个参考物种所在的“转移亚支”。随后，利用ProgressiveCactus生成全基因组比对，并使用CAT在每个亚支内同时注释多个目标基因组。

运行CAT

研究人员将RefSeq注释文件转换为CAT兼容的格式，并使用了AUGUSTUS的三种模式。其中，CGP参数基于“果蝇12基因组计划”中12个已良好注释的物种进行训练。

与BRAKER3互补

为了补充CAT生成的比较注释并减少可能的参考偏好，额外整合了BRAKER3的从头CDS预测。使用EarlGrey移除转座元件，然后比较CAT与BRAKER3的CDS，选择更长的CDS或保留来自CAT的CDS（以反映亲属间的额外支持）。

注释质量评估

为了评估注释质量，将CAT-BRAKER注释与BRAKER3单独生成的注释以及RefSeq注释进行比较，量化CDS一致性。同时，利用BUSCO和OMArk评估注释的完整性。结果显示，CAT-BRAKER相对于RefSeq的精确度通常更高。

正交群分配与CDS比对

利用OrthoFinder v2.5.5，在301个果蝇物种和一个外群（家蝇）中识别CDS同源性，并识别出“分层直系同源群”。对每个选定的HOG，使用MACSE v2和MAFFT v7.520进行序列比对，最终获得了35,642个HOGs和22,355个高质量比对。

系统发育广义线性混合模型分析

为了研究技术因素和进化历史对基因注释的影响，研究人员使用系统发育广义线性混合模型分析了基因数量和CDS长度。模型将物种作为随机效应，将距离参考物种的系统发育距离、基因组大小、组装连续性、测序读长类型以及RNA-seq可用性作为固定效应预测因子。

结果与讨论

301个物种的基因注释

通过整合CAT和BRAKER3，平均在每个基因组中识别出14,549个基因，平均CDS长度为1.60 Kbp，与金标准参考物种黑腹果蝇（Drosophila melanogaster）非常相似。然而，存在显著变异，反映了物种间的进化差异以及基因组组装质量和RNA-seq可用性的潜在差异。其中，三个物种（Drosophila vulcana， Drosophila miranda， Drosophila punjabiensis）的预测基因数超过20,000个。对于D. miranda，基因数升高与其新性染色体系统的近期进化史有关。而对于D. vulcana和D. punjabiensis，高基因数更可能反映了技术而非生物学因素。

通过多重方法评估注释质量：与RefSeq注释的CDS重叠比较显示，CAT-BRAKER的精确度通常高于单独的BRAKER3；BUSCO完整性评估显示，大多数物种在基因组水平和蛋白质水平上的评估高度一致；OMArk评估揭示了少数基因组（如Drosophila recens， D. miranda， D. vulcana）存在高水平的基因重复或污染，这些基因组大多被排除在高置信度子集之外。

直系同源基因推断

使用OrthoFinder推断CDS直系同源关系。将预测的蛋白质分配到直系同源群和分层直系同源群中。研究人员根据HOGs在物种中的分布模式，将其分为两类：广泛保守的HOGs（包括“通用HOGs”和“古老HOGs”）以及物种或支系特异的HOGs。在全数据集中，超过一半的预测蛋白质编码基因属于通用HOGs。然而，有约20,000个HOGs仅包含少数物种（<30个）的基因，这可能反映了近期进化获得或直系同源基因推断问题。为了获得一个置信度更高的正交群集合，研究人员重点关注那些至少包含一个黑腹果蝇基因的HOGs，最终确定了12,151个这样的HOGs，它们广泛分布于绝大多数果蝇物种中。

使用BUSCO和HOG基因进行系统发育推断

为了进行比较分析，研究人员基于251个单拷贝HOGs和1,824个BUSCO基因分别推断物种树。发现HOG树与BUSCO树高度一致，仅在少数物种关系上存在分歧。大多数内部分支在两种树中都得到了良好支持。

影响注释基因数和CDS长度的因素

通过系统发育广义线性混合模型分析发现，基因数或CDS长度与已建立的参考注释相比没有显著差异，表明研究人员的注释具有可比质量。物种距离其参考物种每增加一百万年的分歧时间，平均会损失约16个基因，而平均CDS长度仅增加一个核苷酸。在全数据集中，“读长类型”是基因数量的最强预测因子：仅基于短读长数据的组装比长读长组装多注释约1,000个基因，但CDS平均短60个核苷酸。组装连续性（contig N50）增加与预测基因数轻微减少相关。包含RNA测序数据出人意料地减少了基因数量（约440个），这可能反映了整体假阳性率的降低或基因模型的连接。研究人员发现，基因组组装大小每增加1 Mbp，平均会增加约13个基因，但平均CDS长度略有下降。在考虑了这些固定效应后，研究人员发现果蝇科不同支系间的基因数和CDS长度差异不大。基因数的系统发育遗传率中等（40%），而平均CDS长度的系统发育遗传率较低（9.7%），表明进化历史起一定作用，但物种特异性因素对性状变异贡献更大。

果蝇科的GC组成与密码子使用偏性

为了展示新注释的效用，研究人员分析了整个果蝇科GC含量的变异及其与密码子使用偏性的关系。基因组GC含量范围从Drosophila neohypocausta的21%到Drosophila nannoptera的49%。编码区的GC含量（范围：41%–57%）高于全基因组平均水平。GC₃在亲缘物种间高度相关，其系统发育遗传率估计为1，表明其在进化支系内高度保守。GC₃与非编码GC含量呈正相关（系统发育相关性为0.52），表明全基因组突变偏倚贡献了编码区和非编码区的碱基组成。

为了量化自然选择在决定密码子使用中的作用，研究人员估计了2倍简并密码子的选择强度（S统计量）。S估计值范围从Drosophila pachea的0.24到Drosophila takahashii的1.08。正如预期，melanogaster， montium和ananassae类群表现出较高的GC₃和S，证实了这些类群中有利于GC结尾密码子的选择更强。然而，willistoni和saltans类群——表现出低GC₃——也显示出相对较高的S，确认了Drosophila willistoni中观察到的AT偏倚（至少部分地）是选择的结果。总体来看，果蝇科的GC₃与S呈负相关，这表明GC₃较高的物种，平均而言，实际上经受的密码子使用选择压力更弱。S与基因组大小之间微弱的正相关（系统发育相关性0.15）表明，基因组较大的物种倾向于经历稍强的密码子使用选择。

氨基酸组成

为了研究密码子使用的变异是否与氨基酸组成的变异相关，研究人员分析了果蝇科所有注释蛋白质中20种氨基酸的相对比例。主成分分析（PCA）显示，亲缘关系更近的物种具有更相似的氨基酸使用模式。主成分1主要根据密码子的GC含量分离物种，高GC₃基因组富含GC含量高的氨基酸（如Pro， Gly， Ala， Arg），低GC₃基因组则富含AT含量高的氨基酸（如Asn， Tyr， Ile）。其余主成分的载荷未显示出与氨基酸生化特性（如N/C比或必需性）相关的清晰模式。

结论

本研究为果蝇科301个物种生成了标准化、同步的多物种编码DNA序列注释，作为持续进行的社区合作努力的一部分，旨在推动对整个果蝇科进行全面的基因组学研究。这些新的注释、直系同源基因分配以及多序列比对将为单基因和全基因组水平的进化研究提供宝贵的资源，并将支持未来在这一关键模式类群中进行适应和功能基因组学研究。