基于转录组测序的脓毒症小鼠脾脏差异表达基因分析

卢礼悦，雍 玥，高 浩，李利利，宋建钢

脓毒症是机体对感染的反应失调而导致危及生命的器官功能障碍[1]。全球每年约有3 150万例脓毒症患者，其死亡率约为17%[2]。至今尚未发现能成功治疗此病的药物。其原因与脓毒症致病原因多变，发病机制复杂，几乎涉及免疫系统的各个方面有关。因此，全面地探索脓毒症期间病理生理变化对治疗脓毒症有重要意义。

脾脏包含大量免疫细胞，是人体内最大的外周免疫器官，也是细胞与体液免疫的发生中心。感染可诱发脓毒症，强烈激活人体免疫系统，导致局部及全身炎症反应。若由此引发宿主对病原体相关分子模式(pathogen-associated molecular patterns, PAMP)的异常免疫反应，将诱导大量免疫细胞凋亡并发生免疫抑制[3]。脾脏中则出现淋巴细胞减少、吞噬细胞聚集等现象。脾切除患者遭遇感染后，有更大风险罹患脓毒症，说明脾脏在控制全身感染，调节脓毒症发生发展中起着重要作用[4]，但脾脏控制感染的分子机制尚未阐明。

转录组测序技术（RNA-Seq）能全面快速地获取某一物种特定器官或组织在某一状态下的几乎所有转录信息[5]。本研究通过分析脓毒症小鼠脾脏中mRNA表达变化，以了解脓毒症小鼠脾脏mRNA变化规律，为寻找脓毒症潜在治疗靶点提供理论依据。

1 材料与方法

1.1 实验动物 50只8周龄雄性C57BL/6J小鼠购自国家蛋白质科学中心[动物生产许可号：SYXK（沪）2020-0041]。饲养于SPF级动物房，室温（22±2）℃，12 h昼夜交替，自由摄食和水。

1.2 主要试剂与仪器 脂多糖（LPS，Sigma公司）；生理盐水、Trizol等常用试剂（Servicebio公司）；肿瘤坏死因子（TNF）-α（批号：70-EK282/4-96）、白细胞介素（IL）-6（批号：70-EK206/3-96）、IL-1β（批号：70-EK201B/3-96）ELISA试剂盒（联科生物有限公司）；建库试剂盒（New England Biolabs公司）；Oligo (dT)15（北京天根生化科技有限公司）；Reverse Transcriptase M-MLV（Takara公司）；涡旋混合器、全自动研磨仪（Servicebio公司）；水浴锅（天力医疗器械厂有限公司）；酶标检测仪（BioTeK公司）；Real-Time PCR仪（Roche公司）

1.3 分组与模型构建 适应性饲养1周后，按随机数字表法将小鼠分成LPS组（25只）和Sham组（25只）。LPS组小鼠腹腔注射8 mg/kg LPS[6]，Sham组小鼠予等体积生理盐水。

1.4 ELISA检测 造模后3 h采集小鼠心脏血液，离心取上清液，按照ELISA试剂盒步骤检测血清TNF-α、IL-6和IL-1β的含量。依标准品浓度和光密度（OD）值作标准曲线并计算样本浓度。

1.5 HE染色 造模后3 h收集小鼠脾脏，经4%多聚甲醛固定后石蜡包埋，3 μm厚度切片，二甲苯及梯度酒精脱蜡水化后，苏木精-伊红染色，中性树胶封片，于200倍显微镜下观察脾脏组织的形态变化。

1.6 RNA提取、测序及RT-qPCR验证 造模后3 h在小鼠脾脏样本中加入Trizol（批号：A33248），提取总RNA。Nanodrop测定RNA质量，保证A260/280在1.8～2.2之间。NEB试剂盒将RNA分割成30 0bp左右的短片段，加入随机引物进行反转录，制备成测序文库后在illumina X10平台行双端测序。提取RNA后配制反应液，将混合得到的cDNA溶液行Real-Time PCR检测。用ΔΔCt法分析RT-qPCR检测结果。

1.7 测序数据质控和表达量分析 应用cutadapt（https://cutadapt.readthedocs.io/）对原始测序数据进行质控，通过质检后，利用高通量测序比对工具STAR（https://github.com/alexdobin/STAR）将测序结果与小鼠参考基因组（http://asia.ensembl.org/Mus_musculus/Info/Index）进行比对。根据比对结果使用htseq-count（http://htseq.readthedocs.io/en/master/count.html）工具计算基因表达水平。

1.8 表达差异分析及GO和KEGG富集分析 利用 DESeq2（http://bioconductor.org/packages/release/bioc/html/DESeq2.html）获得组间差异表达结果。以P＜ 0.05且 |log2FoldChange|≥ 1作为筛选标准，得到差异表达基因（DEGs）。利用Metascape工具，以P＜ 0.01为标准，将DEGs富集到不同生物过程、细胞组分、分子功能及代谢通路。

1.9 GSEA富集分析 据样本全基因组表达水平排序，由GSEA V4.0.1软件（http://software.broadinstitute.org/gsea/index.jsp）以P＜ 0.05，Q ＜ 0.25为阈值对结果进行过滤，统计不同功能基因集合在排序的基因列表中排名的富集情况。

2.0 统计学处理 采用SPSS 19.0软件进行数据分析，计量资料用±s表示，满足正态性和方差齐性时采用t检验行两两比较，不满足用非参数秩和检验。P＜ 0.05为差异有统计学意义。

2 结果

2.1 脓毒症小鼠存活率及促炎细胞因子检测 与Sham组相比，LPS组小鼠存活率显著降低，差异有统计学意义（P＜0.01）；同时，LPS组小鼠血清中TNF-α、IL-6和IL-1β的表达水平均高于Sham组，差异有统计学意义（P＜0.01，图1）。

图1 两组小鼠存活率和血清中TNF-α、IL-6和IL-1β表达水平比较

2.2 脾脏织形态HE染色 Sham组小鼠脾组织结构正常，红髓与白髓间界限清晰，脾小结结构正常，细胞分布致密，见图2A。LPS组小鼠脾组织结构散乱，红髓与白髓间界限不清，白髓结构破碎，存在大小不等的缺损，脾小结肿胀，淋巴细胞排列稀疏，见图2B。

图2 小鼠脾组织形态改变（×200）

2.3 测序结果和DEGs筛选 各样本测序获得的总读数在3 330 906～4 580 960个之间，质控后剩余读数在561 541～6 520 634个之间，各样本Q 20 ≥98.64，Q 30 ≥ 95.43，且GC碱基占比≥ 51.17%，表明测序结果良好，数据可用于后续分析。两组间已鉴定的11 009个差异表达基因中，筛选出1 616个在脾脏中显著改变的DEGs。相比于Sham组，LPS组有919个下调基因，697个上调基因，见图3A火山图。此外，根据前100个DEGs的丰度信息，从基因和样品两个层面进行聚类，标准化处理得到Z值后绘制热图，见图3B。

图3 小鼠脾组织RNA-Seq结果

2.4 DEGs及GO和KEGG富集分析 通过分析脾组织中的DEGs，发现腹腔注射LPS主要激活了炎症反应信号。例如，显著上调的DEGs有：CC趋化因子（Ccl25、Ccl3、Ccl4、Ccl5）；CXC趋化因子（Cxcl10、Cxcl11、Cxcl2、Cxcl5、Cxcl9）；集落刺激因子3(Csf3)；干扰素激活基因（Ifi204、Ifi205、Ifi207、Ifi208等）；炎症因子及受体（Il6、Il27、Il11ra1、Il15ra、Il21r）；炎性介质（S100a8）；TNF-α诱导蛋白基因（Tnfaip2、Tnfaip3）；凋亡调节因子（Bcl2l1、Bcl2l11）；以及基质金属蛋白酶（MMP13）等。下调的DEGs有Plekho1、蛋白激酶Map2k6和Mapk14，细胞因子IL-16，以及Snx1和Wdr45b。

为防止分析单个基因得出有偏差的结论，利用Metascape软件，对1 616个DEGs进行了GO和KEGG分析。GO分析显示，这些DEGs主要参与了对病毒的反应、T细胞活化、产生TNF超家族细胞因子、白细胞分化、细胞数目的稳态、细胞因子介导的信号通路等生物过程。其分子功能主要包括参与转录因子结合、蛋白激酶结合、激酶活性、GTP酶结合等，见图4A。

KEGG分析结果同样显示了LPS影响免疫系统，强烈激活炎症途径，例如，Toll样受体信号通路、TNF信号通路、NOD-like受体信号通路、IL-17信号通路、胞质DNA传感途径、RIG-I-like受体信号通路、核因子-κB信号通路、趋化因子信号通路、FcγR介导的吞噬作用、IgA产生的肠道免疫网络等；还参与调控遗传信息，例如，RNA转运、剪接体组成；以及参与JAK-STAT信号转导通路和细胞因子及其受体的相互作用；此外，LPS还与自噬、合成糖鞘脂和调节肌动蛋白细胞骨架等相关，见图4 B。

图4 GO和KEGG富集分析

2.5 GSEA富集分析 为避免遗漏部分差异表达不显著却有重要生物学意义的基因，采用GSEA法二次分析样本基因表达水平，结果显著富集到354条GO条目和22条KEGG通路。大部分GO和KEGG与Metascape软件的分析结果重合，GSEA分析支持和验证了Metascape软件的富集结果。

除Metascape分析中提到的信号通路外，GSEA富集到的前三条关键GO通路是对干扰素β的反应（GO:0035456）、核糖体生物发生（GO:0042254）、病毒防御（GO:0051607）；前三条关键KEGG通路是真核生物核糖体生物发生（mmu03008）、胞质DNA传感途径（mmu04623）以及RNA转运（mmu03013），见图5。

图5 GSEA富集分析

2.6 RT-qPCR验证 为验证RNA-seq测序数据的可靠性，随机选择显著表达的6个上调基因（Mmp13、Zbp1、Ddx58、Cxcl10、Lilrb3、Nr1h3） 和4个下调基因（Plekho1、Polr1d、Wdr45b、Snx1）进行RT-qPCR验证（图6）。结果显示基因表达趋势与测序结果一致（表1）。

图6 差异表达基因RT-qPCR验证

表1 RNA-Seq测序DEGs

3 讨论

腹腔注射LPS诱导小鼠全身炎症以模拟脓毒症发生。两组间改变的DEGs有1 616个。经Metascape富集分析发现，LPS主要影响免疫系统，可强烈激活许多与炎症、遗传信息、信号转导和细胞自噬相关的通路。利用GSEA富集方法，对全基因进行的2次分析显著富集到354条GO条目和22条KEGG通路，多数富集结果与炎症反应、遗传信息处理、细胞凋亡有关。RT-qPCR结果显示10个基因表达趋势与测序结果一致。3个基因Nr1h3、Wdr45b和Snx1表达趋势虽然与测序结果一致，但无统计学差异。可能是由于在LPS刺激的早期，基因在脾脏中表达含量不高，没有产生剧烈变动，导致下调倍数不多[7]，具体原因有待进一步验证。

本研究重点关注到胞质DNA传感途径。胞质DNA传感途径是由模式识别受体的特定家族负责检测来自入侵微生物或宿主细胞的外源DNA，从而诱导机体产生先天免疫应答。在本研究的测序结果中，富集到此途径中有4个基因有显著表达差异，分别是Zbp1、Ddx58以及Cxcl10、Polr1d。

传感器Zbp1是坏死驱动因子，参与激活坏死相关的炎症和免疫反应[8-9]。细胞坏死导致胞膜破裂和细胞内容物释放，可加剧炎症反应[10-11]。Zbp1可激活干扰素调节因子（IRF）和NF-κB转录因子，诱导干扰素（IFN）和其他细胞因子产生。本研究发现Zbp1参与Toll样受体信号通路和胞质DNA传感途径，表明LPS可能通过激活这2条通路中Zbp1介导的细胞坏死进而诱发机体炎症和损伤。传感器Ddx58属于胞质模式识别受体家族，参与如炎症和炎症相关疾病、细胞增殖和凋亡等生物学过程[12]。Ddx58可启动NF-κB、STAT和干扰素三种促炎途径[13]。抑制内皮炎症反应可减轻脓毒症多器官衰竭程度[14]。Ddx58通过调节NF-κB途径，介导脓毒症时内皮炎症反应[15]。本研究发现Ddx58参与胞质DNA传感途径、NF-κB信号通路和RIG-I-like受体信号通路，Ddx58可能通过这三条途径诱导内皮炎症，使脓毒症预后不良。LPS组小鼠脾脏中Cxcl10的表达显著增高。Cxcl10既是单核细胞和T细胞的趋化剂，又是一种促炎因子，可将活化的Th1和NK细胞募集到炎症部位[16]。本研究发现Cxcl10与Toll样受体信号通路、TNF信号通路、IL-17信号通路、胞质DNA传感途径、RIG-I-like受体信号通路、趋化因子信号通路和细胞因子及其受体的相互作用有关，Cxcl10可能在早期炎症中起着重要作用。Polr1d是RNA聚合酶I肽D，参与编码RNA聚合酶亚基，与核糖体前RNA以及小RNA的转录有关。目前，Polr1d的研究主要与特雷彻·柯林斯综合征、结直肠癌、多发性硬化有关，在脓毒症中尚无研究，值得进一步探索。

本结果提示胞质DNA传感途径有可能在脓毒症期间发挥加速炎症信号传导的作用，Zbp1、Ddx58、Cxcl10、Polr1d基因可能是其发挥功能的重要靶点。