富营养饮食重塑代谢性疾病易感猪关键组织mRNA剪接模式研究

夏博策，张凯艺，苗佳坤，彭焕祺，杨 宇，陶 聪，吴添文，王彦芳，杨述林

(中国农业科学院北京畜牧兽医研究所 动物营养学国家重点实验室，北京 100193)

富营养饮食和缺乏锻炼引起的肥胖及相关代谢性疾病，如：2型糖尿病(type 2 diabetes mellitus，T2DM)、代谢性脂肪性肝病、高血压和心血管疾病等严重危害人类健康。动物模型是研究人类疾病的重要手段，相比研究较多的啮齿类动物，猪与人类的生理特征更相似，更适宜于人类代谢相关疾病的转化应用研究。已有研究表明，可通过高脂高糖饮食诱导猪模拟人类富营养饮食引起的代谢性疾病发生发展过程，进而探究其分子病理机制。

可变剪接(alternative splicing，AS)是指一个基因的pre-mRNA通过不同剪接方式形成一个或多个成熟mRNA的过程。AS类型主要包括5种，即外显子跳跃(skipped exon，SE)、内含子保留(retained intron，RI)、外显子互斥(mutually exclusive exon，MXE)、5′端可变剪接位点(alternative 5′ splice site，A5SS)、3′端可变剪接位点(alternative 3′ splice site，A3SS)。AS普遍发生在基因转录过程中，研究表明哺乳动物90%以上的基因可通过AS形成多个可变剪接体，大多数发生AS的基因存在主要可变剪接体，但次要可变剪接体可能会存在不同作用，甚至会引起疾病。已有研究表明，代谢性疾病中的重点基因存在AS，其中瘦素受体(leptin receptor，Lep-R)存在5种可变剪接体，其较长可变剪接体-能够结合瘦素参与能量动态平衡，而较短可变剪接体-、-和-降解瘦素，减少瘦素信号；胰岛素受体(insulin receptor，IR)在人类中存在两种可变剪接体，可变剪接体-通常存在于胎儿阶段，其增加与糖尿病和高血压相关，可变剪接体-通常存在于成年阶段，其在体重减轻时表达量增加。剪接因子(splicing factors，SFs)是一类RNA结合蛋白，参与AS的调控。代谢性疾病相关SFs可被用作疾病的生物标志物和治疗靶点，二甲双胍是一种用于治疗T2DM的药物，可以通过激活和下调来修饰mRNA的AS；剪接因子1可以调控外显子跳跃，在其小分子抑制剂ABX300作用下可改变三种可变剪接体水平，从而抑制肥胖的发生。综上表明，AS在代谢性疾病研究中具有重要意义。

目前，对于代谢性疾病的转录组大多数分析只停留在mRNA表达水平，然而越来越多的研究提示，AS可能是T2DM等代谢性疾病中的一个重要驱动因素。对于代谢性疾病AS研究，单个基因AS研究较多，而整体AS特征研究少有报道。本实验室前期通过CRISPR/Cas9 技术将3个人源T2DM和NAFLD疾病风险基因、和3转入巴马小型猪内，转三基因猪表现出糖耐量受损、脂肪组织炎症、肝炎症、胰腺脂质浸润以及胰岛肥大等早期代谢性疾病特征，并且与人类代谢性疾病具有高度相似性，成功构建了代谢性疾病易感猪模型。

本研究将基于TG-HFHSD和TG-CD组的脂肪、肝、肌肉和下丘脑组织的RNA-seq数据，鉴定AS事件和分析差异剪接基因(differential splicing genes，DSGs)富集到的功能和信号通路，并着重寻找潜在调控肝代谢相关差异AS事件的SFs，旨在为富营养饮食条件下代谢性疾病早期发生机制解析提供一定的参考。

1 材料与方法

1.1 试验材料

本研究选用健康状态、月龄和体重相同的12头代谢性疾病易感公猪(转、和3基因猪)作为试验猪，正常饲料饲喂9个月后，随机将试验猪分为两组，其中8头试验猪饲喂高脂高能量饲料(37%蔗糖、53%基础日粮和10%猪油)2个月作为富营养饮食组(TG-HFHSD)，4头试验猪饲喂正常饲料2个月作为对照饮食组(TG-CD)。试验结束后，试验组平均体重较对照组平均体重高，血清学中TC等部分代谢性疾病相关指标存在显著差异。动物安乐死后，收集脂肪、肝、肌肉和下丘脑组织样本，并立即将其置于液氮中保存。本研究中的巴马小型猪试验均经中国农业科学院北京畜牧兽医研究所实验动物伦理委员会许可(许可证号：IAS2019-12)，单栏饲养，每日定时饲喂2次， 自由饮水，温度16～28 ℃，相对湿度40%～70%。

1.2 试验猪4种代谢相关组织RNA-seq测序

TRIzol试剂(Invitrogen，美国)法提取TG-HFHSD和TG-CD组试验猪脂肪、肝、肌肉和下丘脑组织总RNA后，利用NanoPhotometer spectrophotometer仪器(IMPLEN，美国)和Agilent 2100 bioanalyzer仪器(Agilent，美国)检测RNA纯度(OD/OD及OD/OD比值)，利用琼脂糖凝胶电泳检测RNA完整性(28S、18S、5S三条带存在，并且28S是18S的两倍左右)。NEBNextUltraTM RNA Library Prep Kit试剂盒(NEB，美国)进行cDNA文库构建，Agilent 2100 bioanalyzer(Agilent，美国)对cDNA文库进行质检，在北京诺禾致源有限公司Illumina HiSeq X 平台测序，测序序列读数(reads)为150 bp。

1.3 试验猪4种代谢相关组织测序数据质控比对和可变剪接分析

TG-HFHSD和TG-CD组试验猪脂肪、肝、肌肉和下丘脑组织的原始测序数据(raw data)经过数据过滤、测序错误率检查、GC含量分布检查后得到高质量序列数据(clean data)。HISAT2 v2.0.5作为基因组比对工具，将clean data比对到Ensemble数据库猪参考基因组Sscrofa11.1(网址ftp://ftp.ensembl.org/pub/release-99/fasta/sus_scrofa/)。StringTie(1.3.3b)进行转录本拼接，FPKM值作为基因表达水平单位。

rMATS v3.2.3是一种AS分析软件，主要检测SE、RI、MXE、A5SS和A3SS五种AS事件，可分析各组织样品的AS事件以及对应PSI值，并计算组间剪接差异值和矫正值。使用Draw Venn Diagram网站(http://bioinformatics.psb.ugent.be/webtools/Venn/)对组织间和组织内AS基因进行Venn分析，并绘制相关Venn图。

1.4 试验猪4种代谢相关组织DSGs的GO和KEGG富集分析

g:Profiler网站(https://biit.cs.ut.ee/gprofiler/gost)基于Ensemble数据库，使用网站中的g:GOSt工具可对输入TG-HFHSD和TG-CD组脂肪、肝、肌肉和下丘脑组织组间DSGs集进行GO功能富集和KEGG通路富集分析，并计算矫正值,筛选显著富集GO条目和KEGG通路，采用矫正值小于0.05为筛选标准。使用Hipolt网站(https://hiplot.com.cn)中的气泡图工具绘制GO和KEGG气泡图。STRING网站(https://cn.string-db.org/cgi/input.pl)分析蛋白互作网络。Cytoscape绘制蛋白互作网络图，使用其中的Network-Analyzer功能对网络进行分析，连接度(degree)大于14的节点被认为是网络中的枢纽节点。

1.5 试验猪肝组织SFs获取与差异分析

通过AmiGO2网站(http://amigo.geneontology.org/amigo/landing)和NCBI网站中人和猪基因组功能注释搜索与AS调控相关基因来获取SFs。使用DESeq2(1.16.1)软件分析TG-HFHSD和TG-CD组试验猪肝组织组间差异表达SFs，计算组间基因表达差异值和矫正值。使用Hipolt网站(https://hiplot.com.cn)中的火山图工具绘制基因水平SFs火山图。

1.6 试验猪肝组织SFs与AS的皮尔逊相关性分析

使用R语言中皮尔逊相关性检验函数cor.test计算差异TG-HFHSD和TG-CD组肝组织组间差异SFs mRNA表达水平FPKM值与组间差异AS事件PSI值的相关性，以<0.001和|r|>0.95为相关性筛选条件。Cytoscape软件绘制相关性网络图，使用NetworkAnalyzer功能对网络进行分析，连接度(degree)大于14的节点被认为是网络中的枢纽节点。

1.7 试验猪肝组织AS事件的半定量PCR验证试验

根据Ensemble上的22序列(序列号：ENSSSCT00000063346)和1序列(序列号：ENSSSCT00000029403)，利用NCBI Primer-BLAST针对其发生AS的位置设计引物U2AF2-P和THOC1-P(两个基因引物均可产生两种长度的PCR产物)，引物序列见表1。以试验猪肝组织cDNA 为模板，利用高保真酶2×Phanta@ Max Master Mix(Vazyme，中国)的50 μL体系扩增，高保真酶反应体系(Vazyme，中国)：2×Phanta Mix 25 μL,cDNA模板2 μL，正、反向引物各2 μL,ddHO 19 μL；PCR反应程序：95 ℃ 5 min, 95 ℃ 15 s,60 ℃ 15 s，72 ℃ 30 s/kb，34个循环；72 ℃ 5 min。 产物使用1.5%的琼脂糖凝胶电泳检测，将不同长度片段分别切胶回收送至北京擎科生物公司进行双向测序。

表1 PCR引物序列

2 结 果

2.1 试验猪4种代谢相关组织的AS总体分布特征

原始测序数据(raw data)经质检后得到高质量序列数据(clean data)， 平均每个组织样品clean data 可达8 ～11.45 GB。所有样品平均GC含量为51.71%，平均Q30大于94%，数据过滤后的clean data可用于后续分析。利用rMATS软件鉴定AS事件，结果在试验猪脂肪、肝、肌肉和下丘脑中分别鉴定到48 279、39 536、47 888和52 210个AS事件，分别发生在10 757、9 421、10 023和11 323个AS基因中，平均每个基因发生的AS事件数分别为4.49、4.20、4.78和4.61个，可以发现下丘脑的AS事件和AS基因最多，肝的AS事件和AS基因最少，并且每个基因平均可产生4个以上的AS事件。4种组织鉴定到的5种AS类型中，SE类型事件占比最高，达到77.09%～78.29%，RI、MXE、A5SS和A3SS类型分别为2.28%～2.77%、14.34%～15.26%、1.67%～1.98%、2.97%～3.40%。

通过对AS基因进行venn分析，发现AS基因大部分是4个组织共有的，可达66.77%～80.25%。而各组织特有的AS基因占比较少，仅有3.41%～9.75%(图1e)。各组织中只发生1种AS类型的基因占发生AS总基因比例为54.59%～56.97%，其中只发生SE类型的基因最多；发生2～4种AS类型的基因占比为42.61%～44.87%，其中SE和MXE类型共发生的基因最多；发生5种AS类型的基因占比为0.42%～0.56%(图1a-d)。

a.b.c.d.肝、肌肉、脂肪和下丘脑中5种AS类型的基因Venn图；e. 4种组织总AS基因Venn图a.b,c,d. The Venn diagrams of the 5 AS types of genes in liver, muscle, adipose and hypothalamus, respectively; e. The Venn diagram of the AS genes of the 4 tissues图1 AS基因的Venn分析Fig.1 The Venn analysis of AS genes

2.2 TG-HFHSD和TG-CD组4种代谢相关组织间DSGs的分布特征

以<0.05筛选TG-HFHSD和TG-CD的组间DSGs，脂肪、肝、肌肉和下丘脑中分别得到528、1 070、 570和600个DSGs。脂肪组织中鉴定到的SE、RI、MXE、A5SS和A3SS事件分别为346(65.53%)、 28(5.30%)、110(20.83%)、13(2.46%)和31(5.87%)个；肝组织中鉴定到的SE、RI、MXE、A5SS和A3SS事件分别为329(30.75%)、177(16.54%)、471(44.02%)、34(3.18%)和59(5.51%) 个；肌肉组织中鉴定到的SE、RI、MXE、A5SS和A3SS事件分别为390(68.42%)、20(3.51%)、 102(17.89%)、23(4.04%)和35(6.14%) 个；下丘脑组织中鉴定到的SE、RI、MXE、A5SS和A3SS事件分别为448(74.67%)、12(2.00%)、99(16.50%)17(2.83%)和24(4.00%)个(图2)。各组织综合比较显示，除肝以MXE类型为主，其余组织都是以SE类型为主，肝中大幅增加的DSGs主要分布在RI和MXE两种类型。

图2 4种组织5种类型DSGsFig.2 5 types of DSGs in the 4 tissues

2.3 TG-HFHSD和TG-CD组4种代谢相关组织间DSGs的富集分析

对TG-HFHSD和TG-CD组的脂肪、肝、肌肉和下丘脑中的DSGs进行富集分析，得出脂肪、肌肉和下丘脑的DSGs主要显著富集在细胞内解剖结构和细胞器等与细胞组成相关的GO条目。肝DSGs除富集到上述细胞组成相关的GO条目外，还显著富集到代谢、可变剪接和DNA损伤修复3个方面的GO条目中。代谢相关的GO条目包括，代谢过程、基础代谢过程、脂质代谢、脂肪酸代谢过程、线粒体、过氧化物酶体、α-氨基酸代谢过程、核糖核蛋白复合物；可变剪接相关的GO条目包括，通过剪接体调控mRNA剪接、mRNA加工的调控、mRNA结合、mRNA代谢过程调控、磷酸核苷代谢过程；DNA损伤修复相关的GO条目包括，DNA损伤刺激的反应、DNA修复调控(图3a)。KEGG通路显著富集到代谢、免疫及可变剪接相关的3个方面，其中代谢相关通路包括，代谢途径、丙酮酸代谢、抗坏血酸和醛糖酸代谢、乙醛酸和二羟酸代谢、脂肪酸降解、过氧化物酶体；免疫炎症相关通路主要是补体与凝血级联等；可变剪接相关通路主要是剪接体等(图3b)。

图3 肝DSGs的GO和KEGG富集分析Fig.3 GO and KEGG enrichment analysis of DSGs in liver

肝DSGs富集分析中，GO分析富集到代谢过程(GO:0008152)的DSGs有553个，KEGG分析中富集到代谢途径信号通路(KEGG:01100)的DSGs有152个，两者进行venn分析获得579个与代谢相关的DSGs。将579个代谢相关的DSGs进行STRING蛋白互作网络分析，获得一个有327个节点和1 114个边的互作网络，以连接度(degree)大于14筛选网络枢纽(hub)基因，主要包括3个方面，一是与剪接调控有关的基因25、3B1、7、17、、；二是与代谢相关的基因、1、1、、2、1、17B4；三是与凝血和纤溶相关的基因2、2、、(图4)。

图4 肝代谢相关DSGs的STRING蛋白互作网络分析Fig.4 STRING protein interaction network analysis of DSGs related to liver metabolism

2.4 TG-HFHSD和TG-CD组肝组织组间差异SFs的分析

SFs是调控AS的重要基因。通过AmiGO2网站和NCBI网站中人和猪基因功能富集和注释查询到355个SFs。在TG-HFHSD和TG-CD组肝RNA-seq数据中，通过<0.05和|log(FC)|>1筛选得到56个组间mRNA水平差异表达的SFs，其中上调SFs有46个，下调SFs有10个，其中40B上调倍数最高和3下调倍数最低(图5a)；通过<0.05筛选得到49个组间AS水平差异的SFs。将56个mRNA水平差异的SFs与49个AS水平差异的SFs进行Venn分析，存在11个mRNA水平和AS水平都差异的SFs(图5b)，这些两个水平均差异的SFs可能既可以对其他基因进行剪接调控，又可以对其自身进行剪接调控。

a.SFs的火山图，图中已标注出|log2(Fold Change)|>2的基因名称；b.SFs基本特征:DSG_SF指AS水平差异SFs，DEG_SF指基因水平差异SFsa. The volcano map of SFs, and the gene names with |log2(Fold Change)| > 2 are marked in the figure; b. Basic characteristics of SFs: DSG_SF refers to the differential SFs at AS level, and DEG_SF refers to the differential SFs at mRNA level图5 肝中SFs的mRNA水平和AS水平分析Fig.5 Analysis of mRNA level and AS level of SFs in liver

2.5 TG-HFHSD和TG-CD组肝组织组间差异SFs与组间差异AS的皮尔逊相关性分析

将TG-HFHSD和TG-CD组肝组织56个组间mRNA水平表达差异SFs对应的FPKM值与579个代谢相关组间DSGs对应AS事件的PSI值进行皮尔逊相关性分析，以|r|>0.95和<0.001筛选极显著相关关系，构建了12个重要SFs和222个差异AS事件的相关性网络图，其中剪接因子11和是相关性网络枢纽基因，11和能够分别与222个AS事件中的119和115个存在相关性，在222个AS事件中，49个事件上调，173个事件下调(图6a)。12个SFs中的、1、2AF2、25、5、3H13、2A、1、4、1发生了AS，其中2AF2、25、2A、1发生差异剪接。SFs与AS事件的相关性示例中，2、7和2是图4中STRING蛋白互作网络中的hub基因(图6b)，说明SFs可能会通过影响关键基因来改变肝整体代谢。

a.SFs与AS事件相关性网络图:中间菱形图代表SFs，周围圆形图代表AS事件;SFs与AS事件之间的连接线，红色代表正相关，绿色代表负相关。b.SRSF11和SREK1与AS事件C4BPB-RI、CPB2-MXE、SRSF7-SE和AGXT2-MXE的相关性散点图:一个AS事件由基因名称和AS事件类型组成，例如“C4BPB-RI”中，“C4BPB”为基因名称，“RI”为AS类型a. Network diagram of the correlation between SFs and AS events: the middle diamonds graph represent SFs, and the surrounding circles represent AS events; the connecting line between SFs and AS event, red represents positive correlation, green represents negative correlation. b. Scatter plot of the correlation between SRSF11 and SREK1 and AS events C4BPB-RI, CPB2-MXE, SRSF7-SE and AGXT2-MXE: An AS event consists of gene name and AS event type, for example, in "C4BPB-RI", "C4BPB" is the gene name, and "RI" is the AS type图6 肝中SFs与代谢相关AS的相关性分析Fig.6 Correlation analysis of SFs and metabolism-related AS in liver

2.6 试验猪肝组织AS事件的验证

通过半定量PCR验证代谢性疾病易感猪肝组织中2AF2和1的AS事件(图7a)，DNA双向测序结果显示，2AF2存在外显子2跳跃事件(图7b)，1存在内含子13保留事件(图7c)，分别存在两个可变剪接体。

a.扩增U2AF2和THOC1两个可变剪接体以及内参基因18S的琼脂糖凝胶图；b.U2AF2两个可变剪接体反向测序结果；c.THOC1两个可变剪接体正向测序结果a. Agarose gel image of amplified U2AF2 and THOC1 alternative isoforms and the reference gene 18S; b. Reverse sequencing results of two alternative isoforms of U2AF2; c. The results of forward sequencing of the two alternative isoforms of THOC1图7 肝AS的PCR验证Fig.7 PCR verification of AS in liver

3 讨 论

目前为止，AS已经被发现40多年，基因通过AS产生多个可变剪接体，使蛋白呈现多样性，这个过程几乎参与了机体所有的生理活动。近几年AS相关研究表明，AS与免疫、肥胖和代谢性疾病之间也存在重要联系。为探讨富营养饮食在代谢性疾病发生过程中对AS和SFs的影响，本研究以代谢性疾病易感猪为材料，进行为期12周的高脂高能量饮食诱导，分析了脂肪、肝、肌肉和下丘脑组织的AS事件，结果显示4种组织AS总体特征是相似的，SE在总AS事件中占比最高，这与之前猪和鸡的研究结果一致；也发现每个基因平均可发生4个以上的AS事件，要高于之前猪睾丸AS研究中的3.67个，可能是因为易感猪体内代谢性疾病高风险基因效应增加了单个基因的AS事件。在发生AS的基因中，大多数基因是4种组织共有的，而各组织特有的AS基因占比较少，也与之前研究一致；发生多种AS事件的基因占比接近一半，并且在多AS事件的基因中SE和MXE共发生的基因占比较多，这可能预示着SE和MXE的发生机制存在某种关联性。AS差异分析发现，肝在DSGs特征上与其他组织明显不同，肝组织总AS事件最少，但DSGs最多，并且是以MXE为主，不同于其他3种组织以SE为主的特征，RI也较其它组织显著增加。

对DSGs进行功能富集分析发现，脂肪、肌肉和下丘脑中DSGs主要显著富集在细胞结构相关的GO条目，肝DSGs除富集到与细胞结构相关的GO条目外，GO和KEGG还富集到可变剪接及mRNA代谢、脂质代谢、DNA损伤及修复、补体与凝血级联等过程。有研究发现，高糖饮食能够引起肝DNA损伤，而DNA损伤反应可调节AS，一些错误AS产生的错误可变剪接体可被无义介导降解而未被降解,可变剪接体可能被翻译从而引起疾病。补体与凝血级联可参与免疫炎症过程。将肝中代谢相关DSGs进行蛋白互作网络分析，获得枢纽基因。这些枢纽基因中，25、3B1、7、17、、是与AS调控有关的基因。2、2、、与凝血和纤溶相关，2是一种凝血酶激活的抗纤溶羧肽酶，2是一种凝血因子，是一种纤溶酶原蛋白，是一种纤维蛋白原。、1、1、、2与线粒体相关，其中可参与催化线粒体脂肪酸β氧化的第一步，也有研究表明其参与肝损伤机制；1是一种线粒体定位的酶，与脂质代谢、T2DM和体重减轻恢复相关。1和17B4与脂肪酸代谢相关，其中1是脂肪酸β氧化途径的第一个酶，失调会导致小鼠严重的微囊性肝脂肪变性和的持续激活，其可变剪接体1b比1a能更有效的调节的稳态，与甘油三酯在肝细胞中积累也有一定关系。结果表明，12周富营养饮食诱导的代谢性疾病易感猪肝组织AS发生了显著变化，其中DSGs参与了剪接调控、脂代谢、DNA损伤、代谢性疾病以及补体与凝血级联相关过程。

SFs在AS调控中起着重要作用，在人类胰岛素抵抗肥胖患者肝和轻度或晚期代谢性脂肪性肝病患者中都检测到SFs表达量的改变。本研究在富营养饮食诱导下的代谢性疾病易感猪肝中也发现了mRNA水平表达差异的SFs。参照文献中SFs与AS事件皮尔逊相关性分析的方法，构建出mRNA 水平差异SFs与代谢相关DSGs相关性网络。筛选到的重要SFs中，1是一种可自我调节自身AS的基因，并且与糖酵解有关；1与胆固醇代谢调节相关；1表达下降会导致线粒体功能和DNA损伤增加，也可参与促炎过程；3在非酒精性脂肪肝炎(non-alcoholic steatohepatitis，NASH)患者发病机制中可能起重要作用；13参与脂质代谢和高血脂症相关的代谢综合征发病；2与肝细胞发育成熟有关，并且与酒精性脂肪肝炎(SAH)相关；-1ɑ与脂质代谢、肥胖、胰岛素抵抗和NASH相关；4缺失加剧了高脂诱导的代谢性脂肪性肝病小鼠的肝脂肪变性、胰岛素抵抗和炎症反应，相反，肝4过表达显著减轻了两种不同代谢性脂肪性肝病模型的病理改变。可以看出，这些SFs与糖脂代谢、代谢性疾病之间有一定关系，本研究在相关性网络中筛选到的SFs可能会是调控肝AS并影响肝早期异常代谢的潜在基因。

4 结 论

本研究对富营养饮食条件下代谢性疾病易感猪脂肪、肝、肌肉和下丘脑的AS分布特征进行分析，发现4种组织都可发生AS事件的改变，其中肝改变最显著，主要是RI型和MXE型的改变；富集分析显示，肝DSGs主要与剪接调控、脂代谢、DNA损伤以及补体与凝血级联过程相关；肝差异表达的SFs也会影响糖脂代谢和代谢性疾病发生。综上表明，AS在富营养饮食引起的肝代谢性疾病中具有重要作用，转录组层面的研究应该关注重要SFs通过调控肝AS导致代谢异常的机制，尤其是重要DSGs不同可变剪接体在肝代谢性疾病发生过程中的作用，将为肝代谢性疾病早期发生机制研究提供参考。