线虫姐妹种间基因表达可塑性比表达分歧更为普遍：发育调控网络的保守与演化

本文推荐：本研究通过比较秀丽隐杆线虫（Caenorhabditis）近缘姐妹种C. remanei与C. latens的跨组织与跨性别转录组，揭示了在形态和发育表型保守的背景下，基因表达在物种间存在广泛分歧（~49%基因），但组织与性别差异仍是表达变异的主要来源。研究发现雄性偏向（male-biased）基因在物种分歧中贡献突出，而组织偏向（tissue-biased）基因则呈现较高的序列与表达保守性，说明发育可塑性的共享选择压力强于物种特异性选择。研究为理解基因调控网络（GRN）演化与发育系统漂变（DSD）提供了重要见解。

　　

1 引言

动态基因调控对生物发育与稳态至关重要，其演化也与适应性密切相关。基因调控元件——如位于基因附近的顺式调控序列（cis-regulatory sequences）和基因组其他位置编码的反式作用因子（trans-acting factors）——相互作用并与其他基因产物共同塑造基因调控网络（GRNs）的结构与功能。GRN通常由模块化的子网络组成，以同步方式执行不同功能，从而允许对基因表达进行适应性可塑性控制，以协调不同细胞类型、发育阶段和性别之间的选择压力冲突。尽管如此，独立进化的种群积累的不同突变可导致GRN的分歧。物种特异性选择是基因表达谱差异的潜在驱动因素，而发育系统漂变（DSD）也可能导致调控分歧和基因表达谱的分化，即使发育表型保持保守。

基因调控机制的模块化和部分独立性使得不同发育细胞类型（即由于差异调控导致的组织偏向表达）的转录组谱可能存在差异化演化。在动物组织中，脑组织显示出最高的保守性，而睾丸中表达的基因在物种间表现出显著的分歧。此外，具有狭窄组织特异性表达模式的基因往往比具有广泛表达和多效性作用的基因进化更快。然而，对哺乳动物组织偏向表达进化的研究表明，同一物种不同组织之间的差异大于跨物种同源组织之间的差异，尽管存在一些相互矛盾的证据，并且在一定程度上取决于所考虑的具体组织。

类似于基因调控机制可导致组织偏向表达，它们也可导致个体间发育可塑性的性别偏向转录组差异。性别偏向基因表达可能源于定量差异或性别限制表达，导致两性之间不同的遗传网络结构。因此，性别偏向基因表达可能同时受到性选择和性别偏向自然选择的塑造，以及被称为性别拮抗的冲突选择压力，从而性别偏向基因表达有助于解决关于性别最佳表达的适应性冲突。此外，雄性偏向基因往往比性别中性基因具有更高的组织特异性，这意味着突变将导致更少的多效效应，因此可以更有效地响应给定的选择压力。性别偏向基因，特别是雄性偏向基因，通常在编码序列和表达水平上比雌性偏向和性别中性基因进化更快。因此，显示性别偏向表达的基因可能对物种间的转录组差异贡献不成比例。

2 方法

2.1 线虫解剖与RNA提取

我们在25°C下使用大肠杆菌OP50在NGM-琼脂平板上培养年轻成虫C. remanei (VX0003) 和 C. latens (VX0088)。为收集腺体和体细胞组织，每个样本解剖约400只年轻成虫，使用0.3 mm × 30 mm针头。解剖时，首先在无菌水中在线虫咽后部切割，使腺体释放到体外，随后通过第二次切割将腺臂与体 carcass 分离。将约400个分离的腺体pooled together 形成一个“腺体”组织样本，其余约400个carcasses形成一个“体细胞”组织样本。每个物种的雄性和雌性分别收集三份腺体和体细胞组织样本。所有样本在Tri-reagent中于-80°C储存约48–72小时，之后使用Zymo RNA Isolation kit (#R2051) 提取RNA。提取的RNA随后用Turbo DNase (#AM2238) 和RiboLock RNase Inhibitor (#EO0381) 处理以降解样本中的DNA。为获得高质量和高浓度的RNA，使用Zymo RNA Clean and Concentrator kit (#R1015) 根据制造商协议进一步纯化样本。

2.2 测序与转录本分析

我们为总共24个样本生成了高质量的mRNA序列（RNA-seq）数据，对应于八个处理（2物种 × 2性别 × 2组织）各三个重复。每个样本平均产生 > 75 million paired-end reads（读长150 bp），使用NovaSeq S4 flowcell生成。cDNA文库制备和下一代测序由加拿大多伦多病童医院应用基因组学中心进行。我们通过结合现有注释和本研究生成的RNAseq数据，使用了C. remanei PX506和C. latens PX534 v2参考基因组的最新注释。这些更新的注释由多伦多大学基因组进化与功能分析中心生成，使用RepeatModeler2默认参数分析基因组中的重复序列。使用Trimmomatic对序列读段进行质量修剪，使用Nubeam-dedup默认参数去除重复读段。使用STAR默认参数将质量验证后的mRNA读段映射到各自物种的参考基因组和更新注释上。使用FeatureCounts量化所有基因的表达，并将后续读段计数数据分析限制在13,785个1:1直系同源基因或使用Orthofinder鉴定的物种特异性单例基因集合上。

2.3 编码序列分歧和密码子使用偏倚估计

使用GffRead v0.12.3从每个基因组组装中检索所有转录本的DNA编码序列。对于每对1:1直系同源基因，使用PRANK v.170427在密码子比对模式下将C. remanei和C. latens的编码序列相互比对。使用HyPhy v2.5.52中的FitMG94脚本估计每个比对后的dN和dS值，并在将所有终止密码子替换为gap后进行。对于每个1:1直系同源基因比对，还使用coRdon v1.14.0计算有效密码子数（ENC），使用每个比对中的C. remanei序列。然后为每个1:1直系同源基因比对计算校正后的dS度量（dS′），以校正密码子使用偏倚对dS估计的影响。使用线性校正dS′ = dS + 0.006926114 * (61 – ENC)，其中0.006926114是为这些比对计算的dS与ENC之间线性回归的斜率（回归中排除dS > 1的比对，因为这些可能代表比对错误或不正确的直系同源推断），61是