基于RNA-seq的等位基因特异性表达揭示DMDmdx大鼠肌肉中候选表达变异差异

《Journal of Applied Genetics》：RNA-seq-based allele-specific expression reveals candidate expressed-variant differences in DMDmdx rat muscle

　　

杜氏肌营养不良症（DMD）是一种严重的X连锁疾病，其特征是肌营养不良蛋白（dystrophin）缺陷导致进行性肌肉退变。尽管DMD基因中的原发变异已被充分表征，但疾病严重程度常存在差异，提示可能涉及额外的遗传修饰因子。本研究利用RNA-Seq对源自同一繁育谱系的DMDmdx与野生型（WT）Sprague-Dawley大鼠进行了高分辨率等位基因特异性表达（ASE）和单核苷酸多态性（SNP）分析。尽管具有共同谱系，研究人员在DMDmdx大鼠中检测到更多显著的RNA-seq支持的ASE相关SNP（3,440个），而WT对照组为1,883个。然而，这些发现应被解释为需要正交验证的候选等位基因失衡。功能注释和计算机蛋白质稳定性分析预测，若干错义变异，包括Myh2、Myh3、Ttn、Obscn和Aldoa中的变异，可能影响蛋白质功能或稳定性。功能富集分析识别出与肌肉收缩、细胞骨架组织和免疫反应相关的通路。总之，这些发现表明DMDmdx与WT肌肉中RNA-seq支持的ASE谱存在差异，并为独立验证及进一步研究DMD中潜在的疾病修饰机制提供了候选位点。

研究背景与问题

杜氏肌营养不良症（DMD）是一种严重的X连锁神经肌肉疾病，由肌营养不良蛋白（dystrophin）基因变异引起，导致肌营养不良蛋白几乎完全缺失，肌肉细胞对收缩诱导的损伤高度敏感，引发肌纤维坏死、慢性炎症和再生失败，最终导致进行性纤维化和肌肉萎缩。尽管DMD的原发缺陷已被确立，但即使携带相同或相似肌营养不良蛋白变异的患者，其临床病程也表现出显著异质性，提示遗传修饰因子（genetic modifiers）对疾病进展和严重程度有重要影响。单核苷酸多态性（SNP）是基因组中最常见的遗传变异形式，位于顺式调控区域时可直接影响等位基因特异性表达（ASE），即杂合个体中两个等位基因的差异表达。已有候选基因研究表明Spp1、Ltbp4、CD40和Thbs1等基因内的SNP与DMD关键临床参数相关。DMD^mdx大鼠模型由TALEN介导靶向Dmd外显子23产生11核苷酸缺失，导致移码和提前终止密码子，其表型更接近人类DMD的某些方面。研究人员旨在利用RNA-seq衍生的SNP调用和ASE分析，表征DMD^mdx与WT大鼠骨骼肌中的表达变异谱，并优先筛选候选位点以供未来验证和功能研究。

关键技术方法

研究采用雄性DMD^mdx与WT Sprague-Dawley大鼠（每组n=4），来源于法国南特大学和南特兽医学校–ONIRIS，在波兰华沙生命科学大学兽医研究所饲养。取股二头肌（musculus biceps femoris）进行RNA提取，使用Illumina NovaSeq X Plus平台进行PE150测序。主要技术包括：FaQCs进行质量控制，STAR比对至Ensembl mRatBN7.2参考基因组，GATK v4.4.0.0进行SNP检测与过滤，GATK ASEReadCounter结合自定义Python脚本进行ASE定量与双侧精确二项检验（Benjamini–Hochberg校正），Ensembl VEP进行功能注释，DAVID v6.8进行通路富集分析，I-Mutant v2.0预测蛋白质稳定性变化（ΔΔG），SIFT预测错义变异功能影响。

研究结果

RNA-seq统计、SNP与ASE鉴定：八个样本共产生1,389,713,495条双端 reads，去除重复后平均每样本保留103,488,999条 reads。过滤后获得DMD^mdx样本173,175个样本水平SNP观察值和WT样本855,702个。DMD^mdx中3,440个显著ASE SNP有3,356个（97.6%）显示强失衡，WT中1,883个有893个（47.4%）为强失衡。两组后果类别分布相似，提示基因型差异主要是定量的而非定性的。组内个体间ASE计数存在变异，WT2高于其他WT重复，DMD2和DMD4低于其他DMD^mdx重复。

ASE分布：DMD^mdx组1 Mb窗口内SNP计数范围为0至297，WT为0至225，染色体区域间SNP密度存在差异。

DMD^mdx与WT大鼠的表达变异/ASE差异：DMD^mdx组有58个SNP位点在全部四个样本中共享，WT组仅1个（3:61685147，位于Ttn）。这些SNP关联25个基因，分布于7条染色体，Chr 3 SNP密度最高（21个），主要富集于Ttn。最常见核苷酸替换为C>T和G>A各15次。基因水平分析中，DMD^mdx组变异涉及27个独特基因、33个SNP，分布于5条染色体，其中19个SNP为预先存在的变异；WT组仅4个基因（Myh2、Obscn、Ryr1、Ttn）、18个SNP，位于3条染色体。两组间有4个基因（Myh2、Obscn、Ryr1、Ttn）携带错义ASE相关变异，但具体变异位置 largely 组特异性。

错义变异与蛋白质稳定性：Myh2和Obscn中多数变异SIFT评分低于0.05（如0.01、0.03），预测为有害；Ryr1和Ttn中变异SIFT评分较高（如0.67、1.00），预测为耐受。DMD^mdx组中Aldoa显示多个有害低置信度变异（位置164、78、110，均为D/G替换，ΔΔG分别为-4.46、-4.43和-1.71 kcal/mol），Myh2（Q/P at 164、E/G at 1392）和Myh3（T/I at 547和482）也显示降低稳定性的有害变异；部分Ttn变异（V/F at 11,624、L/V at 7644、E/D at 3852、R/C at 3125）及Myl11（K/N at 112）和Ryr1（M/I at 1656）预测增加稳定性。WT组中Myh2部分变异预测降低稳定性，K/T at 1818和T/I at 1025预测增加稳定性；Obscn中G/V at 3982降低稳定性，S/N at 1921和K/N at 977增加稳定性；Ttn多数变异降低稳定性。

功能分析：通路分析显示富集于肌肉细胞中的细胞骨架和马达蛋白，以及与免疫反应相关的通路（移植物抗宿主病、同种异体移植排斥、自身免疫性甲状腺疾病、抗原加工与呈递）。生物学过程分析富集于骨骼肌纤维发育、横纹肌收缩、骨骼肌收缩和肌肉收缩，以及免疫反应和通过MHC class Ib/ER途径的抗原加工。

讨论与结论

讨论部分指出，DMD^mdx组显著ASE SNP总数高于WT，但组内变异幅度提示该观察应作为趋势解读，除非有正式的组间统计比较支持。相当比例ASE相关SNP定位至内含子区域，可能标记顺式调控元件或影响剪接调控。Myh2变异可能加剧营养不良肌肉的收缩功能障碍；Myh3变异可破坏肌动蛋白-肌球蛋白相互作用和ATP利用；Ttn变异可损害其机械特性，降低弹性，增加营养不良肌纤维的脆弱性；Ryr1在DMD中因氧化应激和S-亚硝基化而"泄漏"，破坏钙稳态；Obscn在肌营养不良蛋白缺陷肌管中下调；Acta1变异与严重先天性肌病相关；Eef1a2可增强utrophin A翻译，提示其变异可能影响内源性代偿机制。营养不良肌肉的慢性病理环境（持续炎症和氧化应激）可能促进等位基因偏向信号。研究局限性包括每组n=4的生物学重复数适中、从bulk肌肉RNA-seq reads推断ASE、组间比较可能受ASE信息性杂合位点数量、测序覆盖度、细胞组成和残余参考映射偏倚影响，因此ASE位点应解释为RNA-seq支持的候选等位基因失衡而非确定性验证的等位基因特异性调控事件。

结论部分翻译：总之，本研究在DMD^mdx大鼠模型中揭示了Dmd基因之外的改变的表达变异景观。在Ttn、Myh2、Myh3、Obscn、Ryr1、Aldoa、Myl11和Acta1等对肌肉结构、功能和代谢至关重要的基因中鉴定出众多SNP和预测的有害变异，提供了营养不良肌肉的表达变异景观。这些发现与遗传修饰因子可能促进DMD转录水平和功能变异的观点一致。这些发现提示营养不良肌肉环境与改变的等位基因特异性转录模式相关，尽管潜在机制和功能后果需要进一步验证。值得注意的是，DMD^mdx和WT组均源自同一繁育谱系，降低了观察到的差异反映无关背景分层的可能性。然而，ASE分析显示DMD^mdx组显著ASE相关SNP数量更高，提示DMD变异可能与更广泛的等位基因特异性转录水平效应相关，这与关于肌营养不良蛋白超越结构支持作用的新兴观点一致。研究每组四只大鼠的样本量可能限制检测细微效应，未来研究将受益于更大队列。