肠道病毒D68型(enterovirus D68,EV-D68)是小RNA病毒科肠病毒属的无包膜正链RNA病毒,主要通过呼吸道飞沫及密切接触途径传播,具有引发区域性乃至全球性大流行的潜在能力。该病毒自1962年首次分离后在全球范围内长期呈现零星散发状态
[1],但2014年起,全球多国报告了EV-D68的暴发流行
[2-4],自此,EV-D68逐渐在全球多地形成了间歇性暴发与长期散发并存的流行特征
[5-6]。国内首例感染病例于2006年被报道
[7],2014年后,受全球流行影响,国内EV-D68的检出率也显著上升
[8]。近年来,EV-D68仍在我国东部地区持续存在散发病例
[9-10]。EV-D68感染以呼吸道症状为主,轻症通常表现为上呼吸道感染,重症可发展为肺炎并诱发哮喘急性发作,甚至累及神经系统
[11-12]。目前尚无特效抗病毒药物或疫苗,临床以对症治疗为主,深入解析其致病机制是开发防控手段的前提。转录组学分析是对细胞内所有基因的表达水平进行全面检测与分析的高通量技术,其中多时间点的动态转录组学分析,能够直接反映宿主基因随感染时间的表达变化
[13]。本研究选用对EV-D68易感的A549细胞作为体外感染模型
[14],通过动态转录组学技术系统解析病毒早期复制、复制高峰期以及晚期宿主稳态失衡等关键阶段的差异基因与信号通路
[15],为深入理解EV-D68的分子致病特征、开发潜在的靶向治疗策略提供了理论依据。
1 材料与方法
1.1 细胞与病毒
本研究所用A549细胞由本实验室保存,采用含10%胎牛血清(公司:ThermoFisher Scientific)和1%青链霉素(公司:Coolaber)的DMEM高糖培养基(公司:ThermoFisher Scientific)进行培养。EV-D68毒株由中国医学科学院提供,并在本实验室扩增保存。
1.2 转录组样品RNA的提取
本实验采用感染复数(multiplicity of infection,MOI)为0.1的EV-D68病毒液接种A549细胞,同时设置未感染的阴性对照组,分别在感染后6 h、12 h和24 h收集2组细胞,所有实验组与对照组均设置3次独立生物学重复。通过Trizol法提取细胞总RNA,利用NanoDrop ND-2000分光光度计(公司:NanoDrop Technologies)检测RNA纯度,并结合琼脂糖凝胶电泳实验评估其完整性,鉴定合格样品用于后续实验。
1.3 文库构建与测序
使用Illumina TruSeq™ RNA文库构建试剂盒,对质检合格的总RNA样品进行文库构建。首先通过Oligo(dT)磁珠富集mRNA,随后进行片段化并反转录合成双链cDNA。经末端修复、加A尾、连接接头及PCR扩增后,利用琼脂糖凝胶电泳与PCR方法对文库的片段大小分布及浓度进行质控检测,确保其片段大小与浓度符合上机要求,质控达标的文库最终在Illumina NovaSeq™ X Plus测序平台上进行双端测序。
1.4 测序数据质量评估
原始下机数据经由Fastp v0.20.0软件进行质控过滤,得到高质量有效读段,过滤标准为:剔除切除接头及两端低质量碱基后长度<100碱基对的读段;去除N碱基数量超过10个的读段;且低质量碱基(Q≤30)占比超过40%的读段也一并剔除。随后,利用HISAT2(v2.2.1)将高质量有效读段比对至人类参考基因组GRCh38.112,获得序列比对映射(Binary Alignment Map,BAM)文件,使用SAMtools(v1.10)对BAM文件进行排序、索引等处理,并使用RSeQC(v2.6.4)软件包对比对结果进行质量评估
[16]。
1.5 基因差异表达分析
采用StringTie(v2.0.4)软件计算基因表达量,并以每千碱基每百万映射片段值(fragments per kilobase per million,FPKM)进行标准化。使用R软件包DESeq2(v1.26.0)对原始计数数据进行差异表达分析,筛选标准为校正后
P值(
Padj)<0.05且|log
2 FoldChange|≥0.58(即基因表达量倍数变化≥1.5倍),将符合条件的基因定义为差异表达基因(differentially expressed genes,DEGs)
[17]。
1.6 GO、KEGG富集分析
对DEGs进行基因本体(Gene ontology,GO)功能注释与京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)通路富集分析。采用Fisher精确检验法检测在DEGs中显著富集的生物学功能条目与信号通路,富集结果的显著性判定标准为Padj <0.05。
1.7 时序分析
使用R软件包Mfuzz(v2.46.0)对DEGs进行时序趋势聚类分析。基于6 h、12 h和24 h各时间点的各基因标准化表达水平FPKM值,通过模糊C均值聚类算法对具有相似表达谱的基因进行归类。根据累计方差变化及生物学意义,划分DEGs到多个表达趋势簇。
1.8 实时荧光定量反转录聚合酶链式反应(Real-time quantitative reverse transcription polymerase chain reaction,RT-qPCR)
使用总RNA提取试剂盒(公司:CWBIO)提取细胞总RNA,并按照反转录试剂盒(公司:TaKaRa)说明书将其反转录为cDNA。以此cDNA为模板,使用SYBR Green PCR MasterMix试剂盒(公司:TaKaRa)进行RT-qPCR实验。引物序列见
表1,以β-actin作为内参基因,每个样本设置3个生物学重复,基因相对表达水平采用2
-ΔΔCt方法进行计算。
1.9 统计学方法
应用GraphPad Prism 10软件分析数据。正态分布的计量资料以±s表示,比较采用独立样本t检验。P<0.05为差异有统计学意义。
2 结 果
2.1 测序数据及其质量分析
本研究对EV-D68病毒感染组和未感染对照组的A549细胞进行了转录组测序,共产出109.46 Gb数据。每个样本获得的原始读段数38 948 266~53 654 150。经过严格质控过滤后,各样本的高质量有效读段数36 326 778~48 423 282,读段有效率均不低于89.74%。所有样本的Q30碱基百分比均高于93.91%,其中最高可达95.75%,GC含量为43.99%~49.38%。质控结果表明,本次测序数据质量可靠,为后续的差异表达基因分析奠定了坚实基础。见
表2。
2.2 EV-D68感染细胞中DEGs分析
差异表达分析结果显示,EV-D68感染A549细胞后DEGs数量在感染过程存在差异。6 h共检测到158个DEGs,其中139个基因表达上调,19个基因表达下调(
图1A);12 h共检测到122个DEGs,其中78个基因表达上调,44个基因表达下调(
图1B);至24 h,DEGs数量增至2 746个,其中1 484个基因表达上调,1 262个基因表达下调(
图1C)。
2.3 EV-D68感染细胞中DEGs的GO富集分析
为系统解析DEGs的生物学功能,对不同时间点的DEGs分别进行GO富集分析。GO注释从生物过程(biological process,BP)、分子功能(molecular function,MF)和细胞组分(cellular component,CC)3个层面进行分类,可见:在6 h组中,共获得227个显著富集的GO条目,其中173个涉及BP、38个涉及MF、16个涉及CC。按显著性排序的前20个条目中,有17个属于BP、3个属于MF,包括与基因转录程序相关的RNA聚合酶Ⅱ介导的转录正向调控、转录顺式调控区域结合;与细胞凋亡相关的细胞凋亡过程、细胞增殖负调控等生物学过程;与免疫和炎症相关的对肿瘤坏死因子(tumor necrosis factor,TNF)的细胞应答、对白细胞介素(interleukin,IL)-1的细胞应答等生物学过程;以及与代谢相关的细胞脂质代谢过程等(
图2A)。
在12 h组,共获得552个显著富集的GO条目,其中453个涉及BP、69个涉及MF、30个涉及CC。显著富集的前20个条目中,17个涉及BP、1个涉及MF、2个涉及CC,与6 h组相比,12 h组中也包括与细胞凋亡相关的细胞增殖负调控、凋亡过程正调控;但与免疫和炎症相关的条目增多,如脂多糖应答、IL-17细胞应答;以及与基因转录程序相关的miRNA转录正调控、H4组蛋白乙酰转移酶复合物、转录因子AP-1复合物等(
图2B)。
在24 h组,共获得178个显著富集的GO条目,其中120个涉及BP、32个涉及MF、26个涉及CC。显著富集的前20个条目中,6个涉及BP、10个涉及MF、4个涉及CC。与6 h组和12 h组早中期感染过程相比,24 h组富集特征发生明显改变,除了包括与基因转录调控相关的RNA聚合酶Ⅱ转录的负调控、DNA模板转录的负调控等条目外,还包括与细胞组分相关的细胞质、细胞核等条目;以及与基因产物本身的分子活性相关的DNA结合、RNA聚合酶Ⅱ顺式调控区域序列特异性DNA结合等(
图2C)。
2.4 EV-D68感染细胞中DEGs的KEGG通路富集
在6 h组,共富集到26条KEGG通路,按
P值排序后选取显著性排名前20的通路,主要涉及免疫与炎症相关通路,包括IL-17信号通路、TNF信号通路、Toll样受体信号通路及NOD样受体信号通路等(
图3A)。随着感染进程推进,在12 h共富集到29条通路,相较于6 h,丝裂原活化蛋白激酶(mitogen-activated protein kinase,MAPK)通路被激活,炎症反应持续存在并通过核因子κB(nuclear factor kappa-B,NF-κB)信号通路放大,同时宿主细胞启动Janus激酶-信号转导与转录激活因子(Janus kinase-signal transducer and activator of transcription,JAK-STAT)信号通路建立干扰素介导的抗病毒状态(
图3B)。在24 h,显著富集17条通路,结果显示,炎症相关通路(如MAPK和TNF信号通路)仍持续活跃,并与补体和凝血级联通路形成协同调控网络,提示炎症反应在该阶段被进一步强化(
图3C)。
2.5 时序分析
为系统解析EV-D68感染过程中A549细胞转录应答的动态变化,本研究对6 h、12 h和24 h的DEGs进行时序聚类分析。结果显示,所有DEGs可归为6个具有显著时序表达特征的基因簇(基因簇1~6,
图4)。结合GO功能与KEGG通路富集分析,对各簇的潜在生物学功能进行解析。
基因簇1的基因表达量在12 h达到峰值,主要富集于脂质代谢相关通路,例如胆固醇代谢过程(Padj=4.3×10―3)、甾醇生物合成(Padj=0.030)及磷脂酶D信号通路(Padj=0.034)等,HMGCS1、HMGCR、FASN和SREBF2等胆固醇与脂肪酸合成的关键调控因子均包含在此簇中。
基因簇2在3个时间点呈现持续上调的表达模式,该簇基因主要参与信号转导与细胞结构调控,包括蛋白丝氨酸/苏氨酸激酶活性(Padj=0.019)以及细胞―基质黏附(Padj=2.2×10―3)等,涉及HSP90AA1、HSP90AB1、PAK1/PAK2及PTK2等基因。
基因簇3的基因表达量在6 h达到顶峰后迅速下降,显著富集于药物代谢及氧化应激相关通路,如细胞色素P450-药物代谢(Padj=0.041)及烟酸与烟酰胺代谢(Padj=2.3×10―3),包括GSTM1和GSTO2、CYP2C8和CYP2B6等基因。
基因簇4的基因表达量在12 h维持相对稳定,而在24 h显著下调,主要涉及小分子分解代谢(Padj=4.6×10―12)及脂质分解代谢(Padj=1.1×10―5)等生物学过程,并富集于脂肪酸降解(Padj=2.2×10―13)等代谢相关通路,代表基因包括ACADVL、PRDX1、IDH1、RRM2等。
基因簇5呈现出早期高表达、中期抑制、晚期表达上调的动态模式,该簇基因主要定位于染色体区(Padj=0.014)、转录调控复合物(Padj=0.020)等细胞组分,并参与染色体结构及细胞周期调控相关过程,如有丝分裂姐妹染色单体分离(Padj=5.6×10―3)及蛋白酶体介导的蛋白质分解代谢(Padj=6.0×10―3)等,该簇包含RHOA、CAV1、SUMO1/2及XPO1等参与细胞周期调控及核质运输的基因。
基因簇6自12 h起持续上调,并在24 h达到高峰,该簇显著富集于免疫应答与炎症反应相关过程,如TNF(Padj=2.7×10―12)、IL-17(Padj=8.7×10―5)和NF-κB信号通路(Padj=1.0×10―4)等。代表基因包括多种关键趋化因子(如CXCL8、CXCL2、CCL2、CXCL1),提示感染后期可能存在炎症反应的级联放大。
2.6 转录组测序结果验证
为了全面验证动态转录组测序结果的可靠性,本研究综合考量了基因差异表达倍数大小及其潜在的生物学功能,选取了4个上调DEGs和4个下调DEGs进行RT-qPCR验证。所选基因转录水平趋势与转录组测序结果一致,提示本研究的转录组测序结果具有较好的可靠性与重复性。见表
3,
4。
3 讨 论
自2014年以来,EV-D68在全球范围内多次引发流行,对儿童健康构成持续威胁。该病毒临床表现多样,不仅能够引发严重的呼吸道感染及神经系统并发症,还与哮喘的发生和加重密切相关
[18]。然而,宿主在EV-D68感染过程中的动态分子应答机制仍缺乏系统性认识。
基于此,本研究选用呼吸道病毒研究中常用的A549细胞模型,开展EV-D68感染后不同时间点(6 h、12 h、24 h)的转录组测序分析,从不同感染时间点来看,在EV-D68感染早期阶段(6 h),差异表达基因聚集在Toll样受体、NOD样受体通路、TNF以及IL-17等炎症信号通路。其中Toll样和NOD样受体通路可分别活化NF-κB和MAPK级联反应,激活早期免疫应答并诱导TNF等促炎因子产生
[19-20];IL-17家族细胞因子能够与TNF共同作用促进中性粒细胞向肺部迁移,上述过程可引起气道高反应性及支气管收缩等病理改变,为后续中性粒细胞性气道炎症及哮喘的发生创造条件
[21-22]。
感染至12 h,宿主的炎症反应信号进一步整合与放大,MAPK、NF-κB和JAK-STAT信号通路共同构成了核心调控网络。其中p38 MAPK信号通路的激活与肺水肿的形成密切相关
[23],JAK-STAT通路可介导干扰素信号转导,强化宿主天然免疫应答
[24]。这种多通路调控的炎症放大效应与呼吸道合胞病毒感染模式相似
[25],可能是多种呼吸道病毒感染引发严重呼吸道炎症反应的共同机制。
感染到了晚期阶段(24 h),宿主应答开始主要涉及组织损伤、重塑及相关病理生理改变等过程,在这一阶段,补体与凝血级联反应通路显著富集,其活化产物可进一步促进炎症细胞浸润、细胞因子释放,该过程与哮喘和急性呼吸窘迫综合征等呼吸道疾病的发病机制密切相关
[26]。补体系统与凝血系统还能够形成膜攻击复合物,这种复合物可直接损伤内皮细胞,引发微血管血栓形成和血管通透性增高,从而加重呼吸衰竭风险
[27-28]。
同时本研究通过时序分析,把不同变化趋势的基因聚类成簇,反映出了病毒感染后宿主细胞的动态应答变化。其中基因簇3主要涉及氧化应激通路,它在感染早期短暂激活后迅速回落,表明宿主可能一开始试图通过清除活性氧来抑制病毒复制
[29]。到感染中期,富集胆固醇代谢等脂质合成通路的基因簇1达到峰值,这种波动很可能为病毒复制所需的膜结构创造了有利条件
[30]。同样呈波动表达的基因簇5,其包含的
RHOA基因可参与细胞骨架调控及气道平滑肌收缩过程,与病毒感染相关的喘息症状密切相关
[31-32],这一簇基因的变化反映了宿主与病毒之间的动态调控平衡。
在整个感染过程中,包含
HSP90等应激相关基因的基因簇2表达则呈现持续上升的趋势,这类基因可能既有助于宿主细胞抵御应激损伤,同时也被病毒利用来提升自身感染效率
[33-34]。随着感染进入晚期,宿主细胞稳态逐渐被破坏,涉及脂肪酸氧化、溶酶体功能的基因簇4基因表达广泛下调,说明宿主细胞的代谢屏障和稳态维持机制可能在此阶段受损,这在一定程度上为病毒的复制和最终释放创造了有利条件
[35]。与此同时,基因簇6基因则在感染中晚期急剧上调,并且显著富集于TNF、IL-17等炎症通路,伴随着趋化因子的大量表达,这一变化印证了感染后期炎症放大甚至炎症风暴的特征。
本研究通过差异表达基因筛选及时间序列分析,揭示了宿主细胞在EV-D68感染过程中的动态转录变化特征。结果显示TNF、MAPK等炎症通路在感染过程中持续激活,同时也展现了早中期的脂质代谢改变,以及晚期补体凝血系统的激活与细胞稳态失衡的过程,对加深EV-D68致病机制的理解,解析其诱导呼吸系统疾病的分子基础提供了重要线索。然而,本研究仅基于体外细胞模型开展分析,缺乏动物模型及临床样本的进一步验证,未来仍需进一步探索。