五个蔷薇科绣线菊族物种的长读长参考转录组：研究序列分歧、物种相关表达变异与叶片发育候选基因的资源

《Plant Direct》：Long-Read Reference Transcriptomes for Five Species of Rosaceae Tribe Spiraeeae: A Resource for Investigating Sequence Divergence, Species-Associated Expression Variation, and Leaf-Developmental Candidate Genes

蔷薇科绣线菊族（Spiraeeae）物种占据多样生境，并在生长型和叶形态上表现出显著变异。然而，该族尚无已发表的参考转录组，其转录组水平变异也尚未被表征。研究人员生成了覆盖三个属的五个物种的长读长（PacBio Iso-Seq）参考转录组：Petrophytum caespitosum (Nutt.) Rydb.、Spiraea lucida Douglas ex Greene、Spiraea splendens Baumann ex K. Koch、Aruncus dioicus (Walter) Fernald var. dioicus和Aruncus aethusifolius (H.Lév.) Nakai。随后，研究人员利用3′-标签RNA测序（TagSeq）对20个个体进行表达定量，以表征表达模式、同源序列分歧和基因库。多变量分析（主成分分析[PCA]和层次聚类）将重复样本解析为三个与采样属相对应的组。由于材料混合了野外采集植物和盆栽植物，研究人员将表达比较视为描述性的，不将其归因于种间差异的遗传控制。成对分析揭示了属间广泛的转录分歧，包括如A. aethusifolius vs. P. caespitosum等比较中的数千个差异表达转录本。功能富集突出了运输与分解代谢、质体组织和生物胁迫响应方面的变异，提示代谢和细胞器相关过程均参与其中。转录因子家族，尤其是bHLH和bZIP，在差异表达转录本中变异最大。加权基因共表达网络分析（WGCNA）恢复了与物种身份相关的模块；由于种间变异在该数据集中占主导，研究人员将其视为共表达结构的描述性总结，而不将其视为物种特异性调控程序的独立证据。由于读段是针对由物种特异性contig组装的参考进行定量，种内比对率（70.5%–85.1%）高于跨种比对率（25.1%–77.0%），种间表达比较是在考虑这种不对称性的前提下解释的。对同源转录本对的序列水平分析表明，蛋白编码基因主要受纯化选择。一小部分保留的同源对（2.7%–5.0%）显示K_a/K_s > 1，但这些集中在同义分歧较低、比值估计不稳定的区域；因此研究人员将其视为快速分歧的筛选水平候选，富集于胁迫响应、膜运输和离子转运注释，并不声称正选择。这些结果为绣线菊族提供了比较转录组框架，并为未来研究蔷薇科生态与形态变异鉴定出候选转录本和蛋白编码差异。

研究背景与问题：比较转录组学可表征近缘种间转录本丰度、转录本结构和共表达模式差异，为非模式植物提供候选基因与假设。蔷薇科经济物种的遗传研究较多，但许多野生或非商业物种缺乏深入遗传数据，限制了对该科完整生态与形态多样性的比较。绣线菊族（Spiraeeae）是蔷薇科一个支持良好的分支，含9属65–100种，形态和生态多样，但尚无参考转录组，转录组水平变异未表征，叶形态和生态分化的遗传基础不清。研究人员因此为五个物种三个属生成PacBio Iso-Seq长读长参考转录组，并用TagSeq对20个个体定量，以描述物种相关表达差异、同源序列分歧和基因库，并提名叶片发育候选基因。由于叶形态与分类单元混杂且未做定量形态测量，研究人员不将表达差异与叶形关联；材料混合野外采集和盆栽，环境/来源效应不能完全分离，表达比较为描述性，不归因于种间差异的遗传控制；共享参考导致比对率不对称，跨种表达比较需谨慎解释。

关键技术方法：样本包括A. aethusifolius n=6、A. dioicus n=3、P. caespitosum n=5、S. lucida n=3、S. splendens n=3，共20个体；野外P. caespitosum subsp. caespitosum采自Idaho的Salmon–Challis National Forest，S. splendens采自California的Lassen National Forest，其余为盆栽；每个物种一个代表个体（Aa7、Ad7、P10、Sl7、Ss1）用于Iso-Seq参考转录组。方法包括Iso-Seq组装与基准通用单拷贝直系同源基因评估（BUSCO）/Trinotate注释，Salmon定量、edgeR/limma-voom差异表达，ShinyGO富集、PlantTFDB转录因子分析、加权基因共表达网络分析（WGCNA），互为最佳命中（RBH）+LASTZ+KaKs_Calculator计算K_a/K_s，以及Swiss-Prot筛选叶发育调控因子并用Nei–Gojobori法估计K_a/K_s。

研究结果：
3.1 转录组从头组装与功能注释：通过Iso-Seq组装、BUSCO和功能注释，A. aethusifolius全长读段最多（847,960）并产生59,679个非冗余转录本；S. splendens读段最少（333,890）但仍有31,518个转录本；S. lucida含最长转录本（8754 bp）；各物种中位转录本长度为1607–2315 bp；四个组装BUSCO完备度为81.7%–82.8%，S. splendens为71.4%且缺失26.6%，受测序深度限制；涉及S. splendens的RBH比较平均4107对，较其他比较的6336对减少约35%；合并参考BUSCO恢复97.5%，但89.0%为重复、8.5%为单拷贝；注释率种间不同，A. dioicus在NR、NT和Swiss-Prot中较高，A. aethusifolius在KO、GO和KOG中较高，P. caespitosum最低。
3.2 层次聚类与PCA：通过top500最变异转录本进行层次聚类和主成分分析（PCA），样品分为三组，对应Aruncus、Spiraea和Petrophytum；PC1和PC2分别解释54.94%和33.75%的总方差，生物学重复按种紧密聚类，符合分类和已知系统关系。
3.3 差异表达转录本（DETs）与功能富集：通过Salmon定量、edgeR/limma-voom对10个成对比较分析，阈值设为|log₂FC|≥2且错误发现率（FDR）<0.01；A. aethusifolius vs. P. caespitosum有6836个上调、4442个下调，A. aethusifolius vs. S. splendens有5350个上调、4135个下调，S. lucida vs. P. caespitosum有3811个上调、2755个下调；近缘种比较DET较少；Venn图和UpSet图显示多个比较共享和特有的DET集合；转录因子（TF）家族中bHLH和bZIP变异最大，WRKY和C2H2中等，GRAS和ERF较稳定；该分析仅描述TF家族在DET中的表示，不检验TF结合位点或基序富集。
3.4 DETs功能富集：通过Swiss-Prot注释、京都基因与基因组百科全书（KEGG）通路和基因本体（GO）富集分析上调DETs；aminoacyl-tRNA biosynthesis、fatty acid metabolism和peroxisome通路仅在A. aethusifolius vs. P. caespitosum和A. aethusifolius vs. S. splendens中富集；ribosome和oxidative phosphorylation仅在A. dioicus vs. S. lucida、S. splendens和P. caespitosum中富集；GO分子功能以binding和catalytic activity为主，包括ion binding、metal ion binding和ATP binding，RNA binding在A. aethusifolius vs. S. splendens中特异富集；细胞组分富集膜结合细胞器，线粒体相关术语仅在S. lucida vs. P. caespitosum中富集，叶绿体相关术语仅在S. lucida vs. A. aethusifolius和S. splendens中富集；生物过程以代谢过程为主，establishment of localization和transport仅在A. aethusifolius vs. P. caespitosum中富集。
3.5 WGCNA：通过top5000最变异转录本构建加权基因共表达网络，恢复与物种身份相关的模块；由于种间变异主导该数据集，模块相关性被视为共表达结构的描述性总结，不作为物种特异性调控程序的独立证据。
3.6 相对序列分歧与K_a/K_s升高候选位点：通过RBH、LASTZ和KaKs_Calculator的YN算法计算10个物种比较的K_s、K_a和K_a/K_s；K_s分布分三层，A. aethusifolius vs. A. dioicus最低（中位K_s=0.032，75.2%低于0.1），S. splendens vs. S. lucida次之（0.066，64.4%低于0.1），Spiraea–Petrophytum属间比较居中（0.140和0.142），Aruncus与其他属比较最高（0.226–0.260），与系统关系一致；多数同源对K_a/K_s远低于1，纯化选择为主导且最稳健的信号；K_a/K_s>1的保留对占2.70%–4.97%，集中于低K_s的不稳定估计区，作为快速序列分歧的筛选候选，不声称正选择；进入候选筛选的RBH保留比例为37.9%–96.1%，与分歧程度反相关；候选注释涉及热激蛋白、ENTH/VHS家族蛋白、金属离子转运蛋白、锌指RNA结合蛋白、核苷酸转移酶、cyclase-associated protein 1、种子储藏、脂质转移、囊泡运输、氧依赖代谢、离子转运、金属解毒、基因组维持和细胞色素P450等。
3.7 叶片发育调控因子注释与序列保守性：通过合并非冗余参考的Swiss-Prot翻译核苷酸搜索（E<1×10^?5）对40,609个转录本（78%）进行基因水平注释，并筛选CIN型TCP、GRF及其GIF/AN3共激活因子、YABBY、HD-ZIP III/REVOLUTA、BELL和WOX同源盒家族；REVOLUTA注释置信度最高，约85%氨基酸一致性；由于所有物种均针对单一合并参考定量，研究人员不解释这些调节因子的表达；class I KNOX程序（STM/BREVIPEDICELLUS）几乎无代表（仅一个转录本），因此不提出KNOX机制解释Aruncus复叶；许多分配为家族水平，如若干CIN-TCP转录本以约48%一致性匹配TCP4，更准确称为CIN-TCP样；可可靠测试且高置信（最佳命中一致性≥55%）的发育家族同源基因均受纯化选择（K_a/K_s≈0.10–0.25），与全基因组模式一致，提示叶形差异更可能来自调控或表达水平演化，而非这些调节因子本身的氨基酸改变，需用定量叶形态和深度匹配的同质园表达谱检验。

讨论总结：讨论指出转录组变异主轴对应属，S. lucida和S. splendens尽管分别来自盆栽和野外仍聚在一起，说明生长条件不是该结构的主要驱动；但物种身份与采样条件不能分离，故描述为物种相关模式。复制数不等使比较功效不均；共享参考下种内比对率高于跨种，跨种表达差异幅度预期膨胀，跨属比较受影响更大。DET数量随系统距离增加；功能富集反复出现细胞器、氧化还原/能量、脂质/脂肪酸和过氧化物酶体相关类别，提示能量平衡和胁迫相关代谢是转录组分化的组成部分，但不证明适应或环境诱导机制。bHLH和bZIP位于生长调控与胁迫信号界面。序列水平上，纯化选择是最清晰信号，K_a/K_s>1位点因低K_s不稳定而不能证明正选择；叶发育调节因子在序列上保守，叶形差异更可能源于调控而非编码改变；class I KNOX程序缺失，故不提出KNOX机制。未来需整合生态测量、生理性状和种群水平采样。

研究结论（译）：研究人员呈现绣线菊族五个物种的长读长转录组，并用其描述转录组范围表达模式、同源编码序列分歧以及关键叶片发育调节因子的恢复和序列保守性。种间表达差异普遍，并按种和属分组；然而，由于物种身份和生长条件未通过实验分离，研究人员将其解释为物种相关模式，不作适应性或机制性声明。最清晰的进化信号来自同源编码序列，其主要受纯化选择，少数位点显示K_a/K_s升高。总体而言，这些数据为一个研究不足的蔷薇科分支提供了比较资源，并提供了一组有限的候选基因、通路和快速演化位点，供受控发育和生态采样检验。