张 茜，王 瑞，魏皓月，张 硕，黄延芹*

1.山东中医药大学，山东 济南 250014；2.山东中医药大学附属医院，山东 济南 250011

2 型糖尿病（type 2 diabetes mellitus, T2DM）是一种多因素影响的慢性全身性代谢性疾病，以胰岛素抵抗（insulin resistance, IR）和（或）胰岛素分泌缺陷引起的以血糖水平升高为主要特征表现。 T2DM因其逐年攀升的患病率和死亡率已成为影响人类健康、生活与经济负担的重要原因，相关西医治疗在临床应用的收效上较为局限[1]。 健脾消渴方以中医经典书著《丹溪心法》中的消渴方为基，全方以黄芪为君药，以黄连、天花粉、生地黄、川牛膝、佩兰为佐使，六味药材共奏清热健脾、活血化瘀之功，其临床疗效已被证实[2]。 中药复方治疗疾病具有多组分、多靶点协同作用的特点，具体作用机制有待明确。本研究通过构建三元转录网络和蛋白质相互作用 （proteinprotein interaction, PPI）网络，探究健脾消渴方的关键作用靶点和相关生物学过程。

长链非编码RNA（long noncoding RNA, lncRNA）是一类不具备蛋白质翻译功能且长度超过200 个碱基的RNA，越来越多的研究显示lncRNA 的差异表达在T2DM 的相关机制中至关重要[3]。 lncRNA 虽缺少可读框无法参与蛋白质翻译过程，但却具有微小RNA（microRNA, miRNA）的特异性结合位点。lncRNA充当“微小RNA 海绵”吸附于miRNA，通过减少miRNA 与下游靶向mRNA 的结合，进而影响靶基因的调控表达。这种lncRNA-miRNA-mRNA 三元转录网络又称为竞争性内源RNA（competing endogenous RNA, ceRNA）调控机制[4]。 lncRNA-miRNAmRNA 三元转录网络在T2DM 病理生理变化中扮演着重要角色，其相关药物作用靶点的筛选和作用机制的明确仍是当前研究的热点[5]。本文运用高通量测序技术筛选T2DM 模型大鼠胰腺组织差异靶点，运用生物信息学方法探讨健脾消渴方对T2DM 大鼠胰腺组织lncRNA-miRNA-mRNA 转录网络的影响，以期为中药治疗T2DM 提供新的方向。

1 材料与方法

1.1 动物

选用4 周龄的SPF 级雄性Wistar 大鼠，体质量（130±10） g，购自北京维通利华实验动物技术有限公司，动物许可证号SCXK（京）2016-0006。 大鼠饲养于山东中医药大学附属医院SPF 级动物房。 饲养环境温度为（22±4） ℃，湿度50%～60%，通风，昼夜交替周期12 h，自由进食饮水。 该动物实验由山东中医药大学附属医院实验动物伦理委员会批准通过，伦理批准号AWE-2019-001。

1.2 药材与主要试剂

黄芪、黄连、天花粉、生地黄、川牛膝及佩兰饮片均购自江阴天江药业公司；普通标准饲料（北京科澳协力公司）；高糖高脂饲料（北京科澳协力公司）；链脲佐菌素STZ（美国默克公司）；柠檬酸钠（天津致远公司）；柠檬酸（天津鼎盛鑫公司）；大鼠胰岛素试剂盒（武汉华美公司）；胰高血糖素酶联免疫试剂盒（武汉华美公司）；miRNasy Mini（德国QIAGENG 生物公司）；VAHTS Total RNA-seq（H/M/R）文库制备试剂盒（南京诺唯赞生物科技股份有限公司）；PCR引物（武汉赛维尔生物科技有限公司）。

1.3 主要仪器

稳豪型血糖仪及配套血糖试纸（美国强生公司）；安捷伦生物分析仪2100（美国安捷伦技术公司）；QubitR○3.0 荧光仪 （美国生命技术公司）；Illumina NovaSeq 6000（美国因美娜公司）；Veriti96 孔型梯度PCR 仪（美国ABI 公司）；Nano Drop One 分光光度计、Multiskan GO1510 全波长酶标仪、Arktik Thermal Cycler 96 孔型PCR 仪（美国赛默飞世尔科技公司）。

1.4 主要方法

1.4.1 健脾消渴方的制备 将黄芪、黄连、天花粉、生地黄、 川牛膝和佩兰按30∶9∶12∶12∶9∶12 混合，在600 mL 水中浸泡30 min。用蒸馏水煮沸至100 ℃持续30 min，然后分离汤剂。 再加入400 mL 水，并根据相同的方法制备第二煎剂15 min。 将药液制备为以药物量计1 g/mL 的浓缩液，存储至4 ℃冰箱中备用。

1.4.2 T2DM 动物模型的建立及健脾消渴方的治疗

31 只大鼠随机选取10 只作为对照组（C 组），给予常规饲料正常饲养。 其余21 只大鼠予高糖高脂饲料饲养8 周，而后腹腔注射35 mg/kg 链脲佐菌素构建T2DM 动物模型。待72 h 后行尾静脉血糖检测，空腹血糖连续2 次≥16.7 mmoL/L，代表模型成功建立。 将造模成功的大鼠按随机数字表法分为2组：模型组（M 组，n=10）和健脾消渴方组（JPXK 组，n=11）。 C 组腹腔注射等量柠檬酸缓冲液（0.1 mmol/L，pH 4.5）以作对照。根据人与动物的剂量换算公式与已有研究结果，JPXK 组于第9～16 周灌胃予健脾消渴方[6.46 mg/（kg·d）]，C 组与M 组选用等量的生理盐水进行灌胃。 经8 周治疗后，解剖大鼠，取胰腺作进一步检查。

1.4.3 空腹血糖（fasting blood glucose, FBG）、空腹胰岛素（fasting insulin, FINS）及血清胰高血糖素（glucagon,GC）测定 分别于灌胃开始前后，禁食不禁水12 h 后进行尾部采血，血糖仪测定治疗前后的FBG。 ELISA法测定FINS 及血清GC， 计算胰岛素抵抗指数（home ostasis model assessment of insulin resistance, HOMAIR）。HOMA-IR=FBG×FINS/22.5。上述实验内容均从各组中分别随机抽取10 个样本进行测定。

1.4.4 高通量测序技术 从各组中分别随机抽取3份大鼠胰腺组织，应用TRIzol 进行总RNA 的提取，后续分离和纯化RNA 严格按照试剂盒说明书进行。使用Nanodrop 检测总RNA 纯度后，经DNA 聚合酶I 和核糖核酸酶H 合成cDNA，并在其末端添加“A”碱基进行末端修复。 采用PCR 对扩增产物进行纯化和浓缩，建立cDNA 文库。 QubitR○3.0 荧光分析仪定量文库中的cDNA 浓度，使用Fastx 工具包软件检测原始测序序列并对结果进行统计分析。 同时，采用seqtk 工具读取原始数据进行检测，根据预先设定的标准对其进行质量控制，并对参照基因组的特定读数进行标注，以便进行转录。 差异性（differential expression, DE）-lncRNA 和DE-mRNA 按照|lgFC|＞1，P＜0.05 的标 准 筛 选，DE-miRNA 按 照P＜0.05 的标准筛选[6]。

1.5 JPXK 组大鼠胰腺DE-lncRNAs 靶基因的预测

将JPXK 组与M 组筛选出DE-lncRNAs，顺式预测与反式预测目的基因。顺式预测将DE-lncRNAs上下游10 kb 左右的mRNA 基因选为目的基因。 反式预测基于互补序列配对原理，使用BLAST 比对获得与lncRNAs 互补的mRNA。 然后用RNAplex 软件计算与mRNA 互补的lncRNAs 的热力学参数，并选择e≤30 的序列。 将顺式预测和反式预测所得靶基因与DE-lncRNA 导入Cytoscape 构建靶基因网络，明确核心调控基因与核心lncRNA。

1.6 JPXK 组大鼠胰腺DE-lncRNA 和DE-mRNA的富集分析

为了阐明基因的生物学功能和相关信号通路，我们使用R 包clusterProfiler 根据GO 和KEGG 数据库进行富集分析。 将健脾消渴方干预的JPXK 组的相关DE-lncRNA 靶基因和DE-mRNA 依次输入，使用Fisher's 精确检验进行富集计算，对生物过程、细胞成分和分子功能进行GO 和KEGG 富集分析。 具体原理是对GO 和KEGG 数据库条目中的差异表达基因进行注释映射，计算每个GO 和通路条目中的基因数量，然后使用超几何分布检验进行统计。 在通过多次假设检验校正计算出的P 值后，将P＜0.05 作为阈值，筛选出满足该条件的GO 和KEGG条目。 依据P 值进行排序，选择显著性最高的前20个GO 条目与KEGG 通路使用ggplot2 绘制气泡图，分析可能涉及的主要生物过程与关键通路。

1.7 qRT-PCR 验证

将测序检测后的剩余样本进行qRT-PCR 验证。随机分别从3 组选取3 个DE-lncRNAs（NONRATT-013183.2、NONRATT008447.2 和NONRATT027846.2）和3 个DE-mRNAs（Nnat、Lcat 和trpm6）。 严格按照试剂盒的操作程序进行试验，内参选用GAPDH，比较基因的相对表达量选用2-ΔΔCt法计算。上述六种引物序列如表1 示。

表1 qRT-PCR 引物序列

1.8 健脾消渴方治疗后的DE-lncRNAs、DE-mRNAs 共表达分析与lncRNA-miRNA-mRNA 转录网络的构建

对JPXK 组DE-lncRNAs 与DE-mRNAs 进 行pearson 相关性检验，按r＞0.9、P＜0.05 选取正相关表达的差异DE-lncRNAs 与DE-mRNAs 关系对。依据目标基因的表达值，对多样本重复数据（每组大于3个生物学重复）运用回归模型分析和种子序列匹配的方法学分析，计算健脾消渴方治疗后miRNA与lncRNA、mRNA 相应的两两调控关系。基于miRNAmRNA、miRNA-lncRNA、lncRNA-mRNA 之间的相关性联系构建了具备miRNA 海绵吸附作用特点的调控网络。 将网络中mRNA 导入在线工具Metascape 数据库，对其进行GO 与KEGG 富集分析。

1.9 PPI 网络分析

为了进一步确定调节核心，我们将ceRNA 网络中的mRNA 提交至目前拥有最多生物与蛋白质以及广泛和多样化的STRING（https://stringdb.org/）数据库，物种设置为“Rattus norvegicus”，将0.4 设置为最小相互作用阈值，并隐藏了独立存在的节点。 利用Cytoscape v.3.8.2 进行可视化的节点网络分析，绘制PPI 网络结构图。

1.10 统计分析

所有实验数据的统计分析均使用SPSS 26.0 软件包进行，所有数据为计量资料。 测量数据以“±s”表示。 组间比较先进行数据的正态性和方差齐性检验，符合条件者采用单因素方差分析检验(one-way Anova)，组间比较采用LSD 方法；不满足正态性分布或方差不齐性时，用非参数检验，两组比较采用秩和检验。 当P＜0.05 时，数据被认为具有显著差异性且有统计学意义。

2 结果

2.1 各组大鼠给药前后FBG、FINS、GC 的测定

造模前，C 组、M 组、JPXK 组大鼠的空腹血糖测量数值均属于血糖正常范围以内。 大鼠造模后连续两次测量结果FBG≥16.7 mmol/L，提示M 组、JPXK组两组大鼠造模成功。 给药前M 组和JPXK 组大鼠FBG 较空白组均升高（P＜0.05）；与治疗前相比，治疗后JPXK 组FBG 显著降低，相较于同时间点的M 组明显下降（P＜0.05）。相较于C 组，M 组的GC、HOMAIR 水平升高，FINS 水平下降（P＜0.05）。与M 组比较，JPXK 组GC、HOMA-IR 水平下降明显，FINS 水平升高（P＜0.05）。 详见表2。

表2 各组大鼠FGB、FINS、GC、HOMA-IR 水平比较

2.2 DE-lncRNAs 与富集分析

测序结果发现大鼠造模成功后，筛选的385 个DE-lncRNAs 中有181 个上调、204 个下调；经健脾消渴方干预后，大鼠胰腺组织有354 个DE-lncRNAs，其中177 个上调，177 个下调 （图1A、 图1B）。 对JPXK 组与M 组比较筛选后的DE-lncRNAs 靶基因通过GO 与KEGG 数据库进行富集分析，结果发现GO 功能的分子功能主要为核糖核酸酶活性和作用于RNA 的催化活性；细胞成分以细胞器膜、内质网膜和吞噬泡的组成成分为主；RNA 代谢的各类调控、染色体分离的调控和脂质代谢过程是主要生物过程（图1C）；KEGG 通路富集分析主要涉及鞘脂代谢、Toll 样受体信号通路、脂质与动脉粥样硬化、糖尿病并发症中的AGE-RAGE 信号通路、FoxO 信号通路、糖酵解相关代谢途径等通路（图1D）。

图1 DE-lncRNAs 和功能富集分析

2.3 DE-mRNAs 与富集分析

M 组大鼠相较于C 组大鼠有492 个DE-mRNAs，其中229 个上调，263 个下调；经健脾消渴方干预后，大鼠胰腺组织有590 个DE-mRNAs，其中462个上调，128 个下调。健脾消渴方干预后DE-mRNAs的GO 功能富集的主要生物过程涉及对肽类激素的反应、糖脂代谢和血管发育等；细胞组成主要为锚定连接和脂滴等；分子功能包括糖胺聚糖结合等。 由图可知，脂肪细胞脂解的调节、PPAR 信号通路、AMPK信号通路、蛋白质消化吸收和PI3K-Akt 信号通路等是KEGG 富集通路的主要分析结果。 详见图2。

图2 DE-mRNAs 和功能富集分析

2.4 差异表达基因的验证

采用qRT-PCR 的方法对随机选择的6 个DElncRNAs 和DE-mRNAs 验证测序结果。 与测序结果相比，PCR 结果表明，各组之间6 种lncRNA 和mRNA 的差异表达趋势在图3 中基本一致。因此，将所有差异表达的基因均列入后续分析。

图3 qRT-PCR 结果

2.5 健脾消渴方治疗T2DM 大鼠的ceRNA 网络构建

按照P＜0.05 的标准筛选JPXK 组的miRNA 测序结果得到32 个差异表达的miRNAs。 将上述差异表达的lncRNA、miRNA 与mRNA 建立相关miRNA-mRNA、miRNA-lncRNA 关系对，结合呈正相关表达的lncRNA-mRNA 建立有统计学意义的miRNA 的海绵吸附作用的调控网络，得到由87 个mRNA，16 个miRNA，33 个lncRNA 构成的ceRNA 网络（图4）。其中NONRATT027642.2 和rno-miR-92a-1-5p 分别在lncRNA 类别和mRNA 类别中节点值最高。 ceRNA 网络中87 个mRNA 的GO 功能富集的生物过程主要为脂肪细胞分化、蛋白酪氨酸/苏氨酸磷酸酶活性、棕色脂肪细胞分化和葡萄糖跨膜转运的正调控等。 KEGG 富集通路主要为丙酸、丙酮酸、丙氨酸、天冬氨酸和谷氨酸代谢、蛋白质消化吸收、酪氨酸代谢、IR、AMPK 信号通路和PI3K-Akt 信号通路等，详见图5。其中丙酸代谢信号通路富集程度最高，故选取此通路绘制通路图，见图6。 有色节点代表健脾消渴方治疗T2DM 的差异靶标富集在该条通路上的靶点，靶点对于此通路发挥下调作用选取绿色显示，上调用红色表示。 图6 结果表明，健脾消渴方治疗T2DM 差异靶点在本条通路中的数个位点均显示下调作用，由此可知健脾消渴方可能通过调控多条信号通路、干预一条通路中的多个靶点，进而影响了胰腺功能和T2DM 的病情进展。

图4 lncRNA-miRNA-mRNA 转录网络

图5 ceRNA 网络中mRNAs 功能富集分析

图6 丙酸代谢信号通路图

2.6 PPI 网络分析与构建

为显示各个靶点之间的关联性，将ceRNA 网络中mRNA 基因数据导入STRING 数据库进行PPI分析和调试。分析过程中将互动分数设定为高可信度，借助Cytoscape 3.7.2 软件运用节点网络分析的可视化处理获取PPI 网络图。 mRNA 之间的相互关系如图7 所示。 由图7 可知，PPI 网络图共82 个节点，30 条关系路线边组成，PPI 富集P 值为0.000163，Degree值变大则连接线段也随之增粗。

图7 ceRNA 网络中mRNAs 的PPI 网络图

3 讨论

根据病因病机和临床表现特点，T2DM 应属于中医学“消渴病”，又名为“消中”“热中”“脾瘅”和“消瘅”等，历代医家多持“阴虚燥热”为传统立论依据。在当前社会经济环境与饮食等因素的影响下，T2DM肥胖患者比例在患病人群基数中逐年升高，近代医家又从脏腑和气血津液等方面丰富和发展了消渴病的病机理论和相关辨证要点，尤其强调了“健脾”之法的重要地位[7-8]。 课题组结合临床实践治疗T2DM脾虚血瘀证的诊治，创立健脾消渴方。方中重用黄芪为君，健益脾气，使气充以助运化；臣以黄连、生地黄和天花粉清热解毒与滋养阴液并重，有清脾胃伏火而不伤阴之妙；佩兰芳香化湿，可佐君药黄芪醒脾开胃、分清化浊；辅以川牛膝补益肝肾、活血祛瘀，培补先天肾阴以滋后天之脾，使邪去正复。 临床研究表明，健脾消渴方在控制患者血糖水平、提高FINS、降低IR 等方面疗效显著[9-10]。

lncRNA 在DNA 复制、DNA 转录、RNA 翻译和RNA 剪接等基因调控过程中起着重要作用[11]。lncRNA 通过修饰调节转录过程在T2DM 中发挥作用，还可通过多种信号通路直接上调或下调炎症因子的表达影响胰岛β 细胞的氧化应激或凋亡[12]。 基于高通量测序探讨整个转录组水平的药物靶向作用位点已广泛运用于定位中医证候特征性标志物和中药疗效机制的研究中，有助于针对中药复方药物复杂性特点阐明作用机制。 目前，尚未见健脾消渴方对T2DM 大鼠胰腺组织ceRNA 网络的相关研究报道。本研究结果表明，健脾消渴方具有降低FBG、GC和提升FINS 的作用，显著改善胰岛素分泌。 此外，还发现了354 个差异表达的lncRNAs 和590 个差异表达的mRNAs，可能涉及健脾消渴方治疗T2DM的作用机制。 本研究基于高通量测序结果与生物信息学分析构建健脾消渴方干预的T2DM 大鼠胰腺lncRNA-miRNA-mRNA 三元转录网络，进一步探究健脾消渴方多靶点协同治疗T2DM 的作用机制。

本研究发现T2DM 模型大鼠经健脾消渴方干预后有354 个DE-lncRNAs（上调DE-lncRNAs177个、下调DE-lncRNAs177 个）。 mRNA 作为lncRNA 分子功能的下游效应器，从健脾消渴方干预后的结果中筛选出共590 个DE-mRNAs，其中462 个上调、128 个下调。为证明测序结果的准确性，本研究挑选了6 个与T2DM 病情进展密切相关的lncRNA 与mRNA 进行验证。 通过对DE-lncRNAs 及DE-mRNAs 的潜在功能进一步的GO 和KEGG 通路富集分析，预测结果表明健脾消渴方抗T2DM 的作用机制主要涉及细胞器膜的组成、RNA 代谢的各类调控和糖脂代谢等生物功能以及鞘脂代谢、脂肪细胞脂解的调节、Toll 样受体信号通路、PPAR 信号通路、AMPK 信号通路、AGE-RAGE 信号通路、FoxO 信号通路、PI3K-Akt 信号通路和糖酵解相关代谢途径。

线粒体和内质网是与血糖水平调控相关的两种重要细胞器，二者共享结构和功能通信以维持细胞环境稳态[13-14]。 线粒体相关内质网膜（mitochondriaassociated membranes, MAM）是两种细胞器之间密切的接触区域，在细胞信号传导、代谢物运输和细胞凋亡等各种细胞过程中起着重要作用，其中就包括参与胰岛素相关信号通路[13-14]。 研究表明，MAM 可能是胰岛素信号传导的关键调节剂，参与人体不同组织的IR[14]。 脂质代谢异常也是T2DM 病情进展的关键因素之一，有72%～85%的T2DM 患者存在脂质代谢紊乱并主要表现为脂蛋白异常和甘油三酯（triglyceride, TGs）增加、脂蛋白动力学改变[15-16]。 在大多数T2DM 患者中，极低密度脂蛋白、高密度脂蛋白等脂蛋白代谢异常通常与IR 有关，脂蛋白脂肪酶、胆固醇酯转移蛋白和肝脂肪酶的活性受其影响[17]。上述报道表明细胞器膜的组成和糖脂代谢等生物功能具有改善IR 和恢复胰岛素分泌的作用，也验证了本研究的准确性。 同时KEGG 富集显示鞘脂代谢与脂肪细胞脂解调节等信号通路，或与健脾消渴方治疗T2DM 相关，关于健脾消渴方对T2DM 患者脂质代谢的具体机制作用值得进一步研究。 Toll 样受体（Toll-like receptors, TLR）信号通路是影响胰腺内免疫稳态的经典通路，TLR 是控制炎症和免疫反应的关键调节因子，在炎症介导的IR 的发病机制中具有重要作用。 研究表明单核和巨噬细胞通过该信号通路激活TNF-α、IL-1β、IL-6 等促炎细胞因子不利于葡萄糖摄取和脂肪酸代谢，激活IKB 激酶-β 和c-Jun 氨基末端激酶减少胰岛素受体底物，在过分驱动的代谢变化中导致IR 的发生[18-20]。

课题组前期研究发现，健脾消渴方可通过AMPK/mTOR 信号通路提高胰岛β 细胞功能[21]，这与本研究KEGG 通路富集结果一致。 此外，进一步的富集分析显示，构建的该ceRNA 调控网络与脂肪细胞分化、葡萄糖跨膜转运的正调控、丙酸代谢和丙酮酸代谢等多种生物途径有关。 同时，网络中的mRNA 相关PPI 网络分析结果也显示出较强的相互作用关系。本文研究结果也印证了中药通过体内的整体作用网络实现对T2DM 治疗的学术观点。 丙酸是一种短链脂肪酸，通过诱导饱腹激素减少食物摄入量及抑制与脂肪酸合成相关的基因表达限制肝脏中脂质的积累[22]。 研究表明，T2DM 患者肠道微生物群以及血清中丙酸生产较常人有所增加；体外实验中丙酸处理还可诱导脂肪细胞自噬，并抑制mTOR 信号传导，故丙酸代谢异常引起丙酸积累可能会影响胰岛素敏感性[23]。 健脾消渴方可能通过改善丙酸代谢防止脂肪细胞过分自噬及肝脏脂肪变性，防止T2DM 发生进一步恶化。 本研究亦显示该网络中节点值最高的NONRATT027642.2 和rno-miR-92a-1-5p 可能是健脾消渴方干预ceRNA 海绵靶向调控下游因子的关键，进而参与丙酸代谢和糖脂代谢等信号通路，影响T2DM 的进展过程。 但针对NONRATT027642.2 和rno-miR-92a-1-5p 有关的ceRNA轴仍需进一步的实验验证。

综上所述，本研究以T2DM 大鼠为研究对象，以健脾消渴方对胰腺组织lncRNA 及mRNA 表达谱的影响为基础构建lncRNA-miRNA-mRNA 三元转录网络。 这些差异基因作为潜在的药物作用靶点能够影响多条信号通路发挥治疗T2DM 的作用，为以后的研究提供了理论基础和方向。