基于转录组学研究1型糖尿病神经痛大鼠腰膨大处脊髓背角基因表达的改变

李 洁，马力天，金兆丰，郑伟杰，王敬博，毛培军#，叶 欣*

(1空军第986医院内分泌科，西安 710054；2空军军医大学唐都医院消化内科；3潍坊医学院心理学院；468003部队卫生队；*通讯作者,E-mail:joannaxye@126.com；#共同通讯作者，E-mail：myhao2002@163.com)

研究显示，全球约有4.25亿人患有糖尿病，我国约有1.15亿人患有糖尿病[1]。糖尿病不仅给患者本人带来极大痛苦，同时也增加了社会的经济负担。在糖尿病的并发症中，由神经系统损伤引起的临床综合征是最普遍的，具有高发病率和治疗困难的特点。糖尿病神经病变可累及全身周围神经系统任何部分，是一种独特的外周神经系统神经退行性疾病，优先靶向感觉轴突、自主神经轴突，较小程度上靶向运动轴突[1]。英国的一项大型研究表明，约30%～50%的糖尿病神经病变患者会出现神经性疼痛[2]。这种周围神经病变引起的糖尿病神经痛，迄今为止仍是疼痛治疗领域的一大难题，其机制尚不明确，现有的治疗方法也不能有效缓解糖尿病神经痛，因此有必要进一步探讨其发病机制。

脊髓背角是各种感觉的初级整合中枢，其中含有丰富的神经递质、神经肽及受体等[3]。现有的研究表明，脊髓背角的星形胶质细胞和小胶质细胞的活化，及细胞因子的释放可能是维持糖尿病神经痛的重要机制[4]；不仅如此，脊髓背角中的γ-氨基丁酸(γ-aminobutyric acid，GABA)释放减少和/或其受体表达及功能降低导致的抑制性突触传递减弱也促进了该病的发生[5]。本研究利用转录组测序技术研究了1型糖尿病病理性神经痛大鼠脊髓背角中基因表达的变化，为后续研究提供参考。

1 材料与方法

1.1 动物、主要仪器和试剂

清洁级健康雄性Sprague-Dawley(SD)大鼠40只，体质量200～250 g(购自空军军医大学实验动物中心)，4只/笼，光照12 h/12 h明暗交替，室温(23±2)℃，湿度40%～60%，自由饮食。本研究实验动物的使用及处置方法符合3R原则，实验开展前经过空军军医大学动物伦理委员会的批准。动物的生产许可证号为：No. SCXK(Shaan)2019-001；动物实验伦理审批号：IACUC-20201022。

Qubit 2.0荧光计(Invitrogen，Q32866)，微型旋涡混合仪(上海沪西，WH-3)，台式高速低温离心机(Thermo，Thermo Scientific Sorvall Legend Micro 21R)，电泳仪(北京六一，DYY-11)，生物电泳图像分析系统(复日科技，FR-980A)，高通量转录组测序仪(Illumina)。Qubit RNA检测试剂盒(Life，Q32855)，Total RNA Extractor(Trizol)(上海生工，B511311)。

1.2 1型糖尿病病理性神经痛大鼠模型的建立

参考文献所用方法[6,7]先建立1型糖尿病大鼠。40只清洁级健康雄性SD大鼠随机分成两组:对照组和实验组，每组20只。领回动物后，适应性喂养7 d后进行实验。实验组大鼠腹腔单次注射1%链脲佐菌素(STZ)65 mg/kg，注射前需禁食12 h，对照组注射等体积生理盐水。注射后第4天用鱼跃血糖仪检测大鼠空腹血糖浓度，若空腹血糖大于16.7 mmol/L，则认为1型糖尿病大鼠模型构建成功；之后采用von Frey细丝(详见1.3)分别测定两组大鼠后爪机械刺激缩足反射阈值(mechanical withdrawal threshold，MWT)，MWT低于自身基础水平85%则认为1型糖尿病病理性神经痛大鼠模型构建成功[6,8]。

1.3 机械刺激缩足反射阈值测定

将1.2中成模的1型糖尿病大鼠及对照组大鼠分别置于15 cm×15 cm×20 cm的有机玻璃罩内，底部为金属丝制成的网格垫，将不同粗细的von Frey细丝[9]通过底部金属网格刺激大鼠足底。8根von Frey纤维丝对应的机械力值依次为1.4，2.0，4.0，6.0，8.0，10.0，15.0，26.0 g。将大鼠置于金属筛网架上，测试前先将大鼠放于鼠箱内适应30 min，待其安静后开始实验。动物反应包括缩足、舔脚等。每一批大鼠检测完成后，将金属筛网和鼠箱清理干净，鼠箱需喷洒75%酒精消毒，待酒精全部挥发后进行测试。从最小强度开始，每个强度测试5次，每次均记录出现缩足反射的阈值，间隔3～5 s，出现缩足反射3次以上的刺激强度即为大鼠对机械刺激的反应阈值，计算每组的平均值。1型糖尿病大鼠模型构建成功后，分别在第30，60天进行机械痛阈的测定。在第60天痛阈检测结束后，选取MWT低于自身基础水平85%的1型糖尿病大鼠进行后续实验。

1.4 测序步骤、数据评估及质控

转录组测序(RNA-sequencing，RNA-Seq)是利用新一代测序技术进行转录组分析的技术，可全面获得样本几乎所有转录本的序列信息和表达信息，包括编码蛋白质的mRNA和各种非编码RNA。收取大鼠脊髓背角的组织样本时，使用焦碳酸二乙酯(diethyl pyrocarbonate，DEPC)处理过的磷酸盐缓冲液经心脏进行灌注，并使用无RNA酶的剪刀小心剪取脊髓背角组织，经液氮速冻后放入-80 ℃进行保存，之后使用Illumina公司开发的Hiseq高通量测序平台进行转录组测序，测序获得的原始图像数据文件经碱基识别后分析转化为原始测序序列。使用数据质控软件Trimmomatic进行数据处理，主要步骤为：①去除带N碱基的序列；②去除序列中的接头序列；③从序列3′到5′方向开始去除低质碱基(Q<20)；④从序列5′到3′方向开始去除低质碱基(Q<20)；⑤使用滑窗法去除序列尾部质量值在20以下的碱基(窗口大小为5 bp)；⑥去除序列长度小于35 nt的序列本身及其配对序列。质控结束后，使用序列相似性比对工具Blast从质控后的数据中抽取10 000条序列与美国国立生物技术信息中心(National Center for Biotechnology Information，NCBI)的核酸序列资料库(Nucleotide Sequence Database，NT)进行比对，并进行污染检测。

1.5 样本间相关性分析

生物学重复是任何生物学实验所必需的，样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标，计算出Pearson相关系数后绘制样本间相关性分析热图，相关系数越接近1，表示样品之间表达模式的相似度越高。

1.6 PcoA主坐标分析和层次聚类树分析

为了使数据相似性或差异性可视化，使用主坐标分析(principal co-ordinates analysis，PcoA)观察数据间的差异。PcoA是一种非约束性的数据降维分析方法，可以找到距离矩阵中最主要的坐标，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。

此外，还使用层次聚类分析构建树状结构，得到树状关系形式用于可视化分析。层次聚类分析反映的是使用统计算法bray curtis计算所得的样本间距离，该距离反映样本间总体特征分布上的差异。

1.7 组间表达差异分析

为了使不同样本间估计的基因表达水平具有可比性，引入了TPM的概念。TPM(transcripts per million)是计量在RNA池中某个转录本的比例，使不同组、不同样品间估计的基因表达水平具有可比性。TPM同时考虑了测序深度和基因长度，是常用的基因表达水平估算方法[10]。据TPM值分析数据，使用得到的P值以及基因在两组样本中的FoldChange值绘制两组样本数据的显著性差异的火山图，可宏观展示两组样本间基因上下调情况。进一步，根据P值的大小，列出上下调最显著的15个基因。

1.8 差异基因功能富集分析

使用GO富集分析进行差异基因功能的富集分析，能够识别出与生物现象最相关的生物学途径。除了GO富集分析之外，还使用KEGG数据库对差异基因所参与的通路进行富集。采用超几何分布(不放回抽样)对基因进行富集分析，使用富集分析软件(clusterProfiler)进行功能富集分析。校正后的P<0.05时，该功能存在显著富集的情况。富集结果用散点图表示。

1.9 统计学分析

2 结果

2.1 采用腹腔单次注射STZ可使血糖显著升高

与对照组相比，到第4天，20只实验组大鼠中有19只血糖升高，差异有统计意义(P<0.05，见图1)；到第30天时，19只大鼠中有2只死亡，还剩17只大鼠血糖升高；这17只大鼠顺利存活到第60天，说明实验组中共有17只大鼠1型糖尿病模型构建成功。

与对照组相比，*P<0.05图1 造模后不同时间的血糖值Figure 1 Blood glucose values at different time points after modeling

2.2 机械痛阈的测定

分别在注射STZ后的第30，60天检测对照组与2.1中1型糖尿病大鼠的机械痛阈，其阈值低于自身基础水平85%的1型糖尿病大鼠被认为是1型糖尿病神经痛大鼠。结果表明，第30天，在17只1型糖尿病大鼠模型中，有8只大鼠的痛阈降低，差异有统计意义(P<0.05，见图2)；第60天时，这8只大鼠的痛阈依然降低，差异有统计意义(P<0.05，见图2)。

与对照组相比，*P<0.05图2 1型糖尿病神经痛大鼠的机械痛阈Figure 2 MWT of type 1 diabetic neuropathy pain rats

2.3 各样本质量控制后数据信息

通过对原始数据进行质控后，去除其中带接头的、低质量的数据，保证数据质量。其中Q20以上代表了碱基错误率在1%以下，Q30以上代表了碱基错误率在1‰以下。质控后的数据中，各样本的碱基错误率和N碱基的比例极低(见表1)，确保本研究采集到的数据可以用于后续分析。

表1 各样本质量控制后数据信息

2.4 样本间相关性分析

样品间的相关性反映了样品间的相似情况，即不同样品在表达水平的相似度。相关系数越接近1，样品间的相似度越高，样品间的差异基因越少。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。结果表明，同一组之间的样本相关性较好(图中颜色偏黄的区域)；不同组之间的样本相关性指数较低(图中颜色灰白的区域)(见图3)。说明分组可靠，样本选择合理。

颜色块代表相关性指数值，颜色越灰表示样本间相关性指数越低，颜色越黄则相关性指数越高图3 重复相关性分析热图Figure 3 Heatmap of repeated correlation analysis

2.5 PcoA主坐标分析和层次聚类树分析

通过PCoA可以研究数据的相似性或差异性。结果表明，对照组和1型糖尿病神经痛模型组的样本间的差异较大，距离较远(见图4)。层次聚类树分析表明，1型糖尿病神经痛组和对照组大鼠可以明显分开，且每组组内样本重复性很好(见图5)。结合PcoA主坐标分析和样本层次聚类树分析可以说明1型糖尿病神经痛组和对照组的确存在基因转录水平上的差异。

红色代表对照组，蓝色代表1型糖尿病神经痛组图4 PcoA主坐标分析组内数据的相似性与组间数据的差异性Figure 4 PcoA of the similarity within the group and the difference between two groups

树枝的长度代表样本间距离，越相似的样本会越靠近，样本按分组以不同颜色区分，红色代表对照组，蓝色代表1型糖尿病神经痛组图5 层次聚类树分析样本间的相似性Figure 5 Hierarchical clustering tree analysis of the similarity between samples

2.6 1型糖尿病神经痛组与对照组之间差异表达分析

与对照组相比，1型糖尿病神经痛组共有42个基因在转录水平发生了上调，有51个基因在转录水平发生了下调。火山图显示两组样本数据有显著性差异(见图6)，上下调最为显著的15个基因见表2、表3。

图中每个点代表一个基因，红色代表上调的基因，绿色代表下调的基因，黑色代表非差异基因图6 1型糖尿病神经痛组与对照组间差异基因表达火山图Figure 6 The volcano plot of differential gene expression between two groups

表2 1型糖尿病神经痛组中发生上调的基因

表3 1型糖尿病神经痛组中发生下调的基因

2.7 差异基因功能富集分析

GO富集分析结果表明，差异表达的基因富集于酰基磷酸酶活性、主要组织相容性复合物(major histocompatibility complex，MHC)和MHC Ⅰ类蛋白复合物3个GO分类条目(见图7)，提示在1型糖尿病神经痛大鼠模型中，差异表达的基因可能与酰基磷酸酶活性、MHC蛋白复合物和MHCⅠ类蛋白复合物三个功能类群有关。除了对基因本身功能的注释，我们使用KEGG数据库对基因所参与的通路进行富集。结果表明，在1型糖尿病神经痛大鼠模型中，差异表达的基因主要富集在5条通路上，分别是：抗原加工和提呈、氨基苯甲酸酯降解、自然杀伤细胞介导的细胞毒性、细胞黏附分子(CAM)和吞噬体上(见图8)。

P值的大小用点的颜色表示，P值越小则颜色越接近红色，每个功能下包含的差异基因的多少用点的大小来表示图7 1型糖尿病神经痛组和对照组间差异基因的GO富集分析Figure 7 GO enrichment analysis of differential genes between two groups

P值的大小用点的颜色表示，P值越小则颜色越接近红色，每个功能下包含的差异基因的多少用点的大小来表示图8 1型糖尿病神经痛组和对照组间差异基因的KEGG富集分析Figure 8 KEGG enrichment analysis of differential genes between two groups

3 讨论

糖尿病神经痛易发生于糖尿病后期，高血糖会引起周围神经损伤并诱发疼痛，长期高血糖是主要致病因素[11]。我们成功建立了基于STZ导致胰腺β细胞损害的1型糖尿病模型，该模型会使缩足阈值显著下降[12]，与正常大鼠相比，1型糖尿病神经痛组大鼠的机械痛阈显著降低，并渐趋平稳。

既往的研究表明，高血糖会增强神经纤维的氧化应激，随之而来的炎性反应会影响雪旺氏细胞的结构并导致促炎性细胞因子和伤害性感觉神经递质的释放[13]。因此，在某种程度上，糖尿病所致的神经疼痛和神经损伤所致的神经痛具有一些潜在相似的机制和特性。在1型糖尿病神经痛大鼠中，和坐骨神经损伤模型(SNI)相似，脊髓背角的小胶质细胞的激活也促进了中枢神经系统的敏化并参与了疼痛的发生[13,14]；在SNI模型中，脊髓背角星形胶质细胞的活化同样也参与了糖尿病神经痛[15-17]。因此，有必要进一步研究糖尿病神经痛中脊髓背角基因表达的差异，为进一步明确其机制提供参考。

在此次研究中，通过转录组学测序明确了1型糖尿病神经痛大鼠腰膨大处的脊髓背角中共有42个基因在转录水平发生了上调，有51个基因在转录水平发生了下调。在上调的基因中，S100钙结合蛋白A11(S100 calcium binding protein A11，S100a11)是S100家族重要成员之一，可以调节酶活性、调控炎症反应、调节细胞生长和凋亡等[18]。Yuan等[19]通过蛋白质组学的研究提示，S100a11参与了肩袖损伤和肩痛的发生，其具体机制需要进一步研究。3-羟基-3-甲基戊二酰辅酶A合成酶2(Hmgcs2)为酮体生成的关键酶，主要存在于肝脏、肾脏及神经细胞等的线粒体中[20]。现有的研究表明[21]，背根神经节中的Hmgcs2参与了阿片类药物诱导的痛觉过敏，干预Hmgcs2的表达可以在一定水平上缓解芬太尼诱导的痛觉过敏的发生。其在DNP中的具体作用机制尚需进一步研究。

为了检验差异表达的基因是否富集于某些功能中，我们进行了GO富集分析，发现在DNP中，差异基因的改变可能与酰基磷酸酶活性、MHC蛋白复合物和MHC Ⅰ类蛋白复合物的功能相关。近年来，多篇文献已经报道了MHC在神经病理性痛[22]、骨癌痛[23,24]中的重要作用，包括：在小鼠SNI模型的脊髓背角中MHC-I的表达显著增高，进而增强神经元的兴奋性而降低其异质性，加剧了疼痛刺激[22]；骨癌通大鼠脊髓中的MHC-I的表达增加且与GABA能神经元双标，抑制MHC-I的表达能缓解骨癌痛大鼠的痛觉过敏[23]。酰基磷酸酶活性与疼痛之间的关系尚缺乏研究。在KEGG富集分析中，差异表达的基因主要富集在5条通路上，分别是：抗原加工和提呈、氨基苯甲酸酯降解、自然杀伤细胞介导的细胞毒性、细胞黏附分子(CAM)和吞噬体。脊髓背角中这些通路所发生的变化与疼痛之间的关系尚需进一步研究。

综上所述，脊髓背角在1型糖尿病神经痛的发生发展中也发挥重要作用，但其机制有待进一步研究。本研究从整体上揭示了其中的差异基因以及功能，为后续研究提供参考和部分证据，但由于疼痛的发生机制错综复杂，因此需要更多的实验研究来揭示其机制。