NDV 基因Ⅶ型与基因Ⅱ型感染CEF 后microRNA 表达谱分析

薄一丹，庞洪泽，赵 款,2，雷白时,2，蒋文明，袁万哲,2，张武超,2，刘聚祥

（1.河北农业大学 动物医学院，河北 保定 071001；2. 河北省兽医生物技术创新中心, 河北 保定 071001；3.中国动物卫生与流行病学中心，山东 青岛 266032）

新城疫（Newcastle disease, ND）是由新城疫病毒（Newcastle disease virus, NDV）引起的一种烈性、高度致死性禽类传染性疾病［1］。NDV 属于副粘病毒科，禽腮腺炎病毒属，单股、负链、不分节段的RNA 病毒［2］。目前国内以基因Ⅶ型为主要流行毒株，可造成免疫鸡群中出现非典型新城疫症状［3］。基因Ⅶ型NDV 强毒株与早期强毒株F48E9株（Ⅸ型）和Hert33 株（Ⅳ型）相比，其在动物机体内的复制能力更强，可以引起更加严重的免疫系统损伤［4］。

microRNA（miRNA）又称微小RNA，一般由20 ～25 个核苷酸组成内源性的单链RNA，广泛存在于各种动植物细胞中［5］。miRNA 为非编码RNA，本身并不能翻译成蛋白质，可通过与靶基因3’UTR 区域结合起到对目的基因的调控作用［6］。mi RNA 通过上述这种调控方式，在病原体入侵以及机体自身的免疫反应过程中均发挥着重要的调节作用，如miRNA-32 可以通过靶向病毒基因组起到抑制翻译的作用［7］；然而病原体自身也会编码产生miRNA，用来帮助病毒抵抗机体的免疫反应，最终实现免疫逃逸［8］，如人巨细胞病毒（hCMV）可以通过下调miR-UL112 来抑制NK 细胞对病毒的杀伤作用［9］。

Ingle H 等人在2015 年率先证明miR-485 在调节NDV 感染中发挥重要作用［10］，但直到现在关于miRNA 与NDV 之间相互作用的研究仍较少，因此，深入解析NDV 感染鸡胚成纤维细胞（Chicken embryo fibroblasts,CEF）后miRNA 的表达情况，对研发新型防治药物具有重要意义。本研究利用我国现流行的NDV 基因Ⅶ型流行毒株HB 株与疫苗毒株La Sota 株分别感染CEF 后，通过高通量技术对比分析La Sota 株与HB 株感染CEF 36 h 后，宿主细胞miRNA 的差异表达情况，为进一步揭示miRNA在NDV 感染宿主细胞过程中发挥的调控作用提供数据基础。

1 材料与方法

1.1 病毒株、实验动物及主要试剂

NDV 基 因Ⅶ型HB 株（GenBank 登 录 号：MK342603）、La Sota 毒 株（GenBank 登 录 号：AF077761）均由中国动物卫生与流行病学中心保存；SPF 鸡胚购自于山东昊泰实验动物繁育有限公司。总RNA 提取试剂Trizol、RNAiso for Small RNA、Mir-XTmmiRNA First-Strand Synthesis Kit 和Mir-XTmmiRNA qPT-PCR TB Green Kit 均购自宝日医生物技术（北京）有限公司； 1xPBS、青霉素链霉素混合溶液购自于北京索莱宝科技有限公司；FBS 胎牛血清、DMEM 培养基、胰酶均购自于Gibco 公司；TPCK 胰酶购自于Sigma 公司。

1.2 NDV 一步生长曲线的绘制

将HB 株、La Sota 株 分 别 按 照MOI=1 感 染CEF 细胞，在感染后12、24、36、48、60 和72 h收集病毒液，反复冻融3 次后，按照Reed-Muench法测定并计算各时间点TCID50，绘制一步生长曲线。

1.3 测序样品制备

共分为2 组：HB 株（MOI=1）感染CEF 组及La Sota 株（MOI=1）感染CEF 组，每组设置3 个重复。细胞接毒后培养36 h，弃去上清，加入Trizol 收集细胞悬液，干冰保存，送至南京派森诺基因科技有限公司。

1.4 文库构建和sRNA 深度测序

使用 TruSeq sRNA 样品制备试剂盒（Illumina,USA）合成并扩增 cDNA，从 PAGE 凝胶中提取和纯化 145 ～160 bp PCR 扩增片段，根据制造商的说明，合格文库中的 DNA 片段在 Illumina HiSeq 2500仪器（Illumia Inc,USA）上进行测序。

1.5 miRNA-seq 数据处理及分析

1.5.1 已知miRNA 及novel miRNA 的鉴定 测序完成后，对所有的测序数据进行质量筛选，其中碱基长度小于15 nt 大于30 nt 的数据以及未匹配上的数据都将被删除过滤。过滤后数据通过与miBase 数据库参考鸡源基因组(http://ftp.ensembl.org/pub/release103/fasta/gallus _gallus/dna/Gallus_gallus.GRCg6a.dna.toplevel.fa.gz）进行对比来注释sRNA，筛选对比已知miRNA。对已知miRNA 进行首位碱基以及各点位碱基偏好性进行分析。通过miRNA 预 测 软 件miR Evo（miR Evo_vl.1）和Mirdeep2（mirdeep2_0_0_05）预测novel miRNA，并对其二级结构、成熟体以及靶基因进行预测。

1.5.2 差异表达miRNA 的筛选 根据各样品中 miRNA 的表达量数据，采用 DESeq（version 1.18.0，Anders S 和 Huber W，2010） 对 miRNA进行差异表达分析，并按照表达量倍数差异（|log2FoldChange|>1）和表达差异，显著性（P-value＜0.05）筛选出差异的保守 miRNA，通过R 语言Pheatmap 软件对所有miRNA 和样品进行双向聚类分析；R 语言ggplot2 软件绘制差异表达基因火山图。

1.5.3 差异表达miRNA 靶基因GO 功能及KEGG通路分析 分别使用topGO 以及KOBAS v2.0 对差异表达miRNA 靶基因进行GO 以及KEGG 富集分析，分析时利用GO 项目注释的差异miRNA 靶基因对每个项目的miRNA 靶基因列表和数目进行计算，然后通过计算P-value（显著富集的标准为P-value＜0.05），从而确定差异miRNA 靶基因行使的主要生物学功能以及富集通路。

1.6 实时荧光定量 PCR（RT-qPCR）验证

通过实时荧光定量PCR 检测HB 株与La Sota株在感染CEF 后miRNA 表达量，以U6 作为内参基因，随机选取5 个差异表达miRNAs，根据其序列结果设计引物如表1 所示，均由上海生工生物工程有限公司合成。数据处理采用2-△△Ct法计算。

表1 RT-qPCR 扩增引物Table 1 Primers for RT-qPCR

2 结果与分析

2.1 NDV 不同毒株感染CEF 后一步生长曲线

HB 株为基因Ⅶ型流行毒株，La Sota 株为基因Ⅱ型疫苗毒株，按照MOI=1 感染CEF 细胞，绘制HB 株与La Sota 株的生长动力学曲线（图1）。

图1 NDV 不同株感染CEF 细胞的一步生长曲线Fig.1 One-step growth curve of CEF cells infected with different NDV strains

结果显示，2 株病毒在感染初期随感染时间的延长，其滴度不断升高，均于感染后60 h 达到峰值，其中HB 株增殖滴度在60 h 以内均高于La Sota 株，且在24 ～36 h 存在统计学差异，感染后60 h 病毒滴度开始下降，72 h 时La Sota 株滴度高于HB 株但并无显著性差异。由于在病毒感染36 h 时细胞产生明显病变，因此选择在36 h 收集样品进行测序。

2.2 小RNA 文库测序统计

通过Hiseq Single-End 模式测序，对下机后原始数据进行接头去除，根据序列质量进行剪切，过滤序列用于后续分析，HB 感染组最终获得纯净序列占原始数据的51.33%；La Sota 感染组纯净序列占原始数据的49.91%（表2）。

表2 测序质量数据分析Table 2 Quality of sequencing data

2.3 sRNA 长度分布和分类注释

过滤掉长度＜18 nt 以及＞40 nt 的sRNA，发现2 组sRNA 长度均集中在20 ～23 nt 之间（图2）。其中，长度22 nt 的序列占比最高，HB 感染组为18.71%、La Sota 感染组为25.79%。将sRNA 对比鸡源参考基因后，对sRNA 进行分类注释，2 个文库中顺序均为know miRNA>piRNA>unknow miRN A>intron>repeat>exon>snoRNA>snRNA>tRNA>no vel miRNA，know miRNA 占比最高，HB 感染组占49.4%（图3A），La Sota 感染组占58.2%（图3B）。

图2 序列长度分布情况Fig.2 Sequence length distribution

图3 sRNA 分类注释Fig.3 sRNA classification annotation

2.4 sRNA 分类注释

2.4.1 已知miRNA 的鉴定 将长度为18 ～40 nt sRNA与参考基因组已知miRBase 中鸡（G.gallus）miRNA 序列进行对比，其中HB 感染组中共发现554 条miRNA成熟体以及440 条miRNA 前体；La Sota 感染组中共发现613 条miRNA 成熟体以及479 条miRNA 前体。通过统计2 个已知miRNA 文库发现，长度在20 ～24 nt的大部分以首位碱基以U 为主，其次为A（图4）；2 个文库中2 ～8 位种子序列以G 碱基为主（图5）。

图4 miRNA 首位碱基偏好性Fig.4 miRNA first base preference

图5 miRNA 个点位碱基偏好性Fig.5 miRNA first base preference

2.4.2 novel miRNA 预测 根据与基因组的对比结果，运用软件分析未注释到任何信息的序列，进行novel miRNA 序列分析，预测结果表明HB 感染组中共发现45 条novel miRNA；La Sota 感染组中共发现56 条novel miRNA，其中表达量前5 的novel miRNA 如表3 所示，同时通过软件预测表达量前5的novel miRNA 前体二级结构，发现具有典型的发卡结构（图6）。

表3 表达量前5 的novel miRNA 序列Table 3 The top 5 novel miRNA sequences in expression

图6 novo miRNA 前体二级结构预测Fig.6 Prediction of the secondary structure of novo miRNA precursors

2.5 差异表达miRNA 分析

2.5.1 聚类分析 将La Sota 株与HB 株感染CEF 36 h 后的miRNA 表达情况进行聚类分析（图7），横坐标表示基因，每一列为1 个样本，红色表示高表达基因，蓝色表示低表达基因，颜色从红到蓝，表示 log10(FPKM+1)从大到小，可以看出NDV 在感染CEF 36 h 后，存在较多的miRNA 差异表达。

图7 聚类分析Fig.7 Cluster analysis

2.5.2 差异表达miRNA 分析 将La Sota 感染组与HB 感染组miRNA 表达水平进行比较（图8），其中横坐标表示log2FoldChange 值，纵坐标表示-log10（P-value）值，红色表示上调基因，蓝色表示下调基因。按照表达量倍数差异（|log2FoldChange|>1）和表达差异显著性（P-value ＜0.05）从La Sota与HB 表达差异的miRNA，筛选出差异的保守miRNA。共统计到10 个差异表达的miRNA，包括7 个表达上调，以及3 个表达下调（见表4）。

图8 La Sota 株与HB 株感染CEF 细胞的差异表达基因火山图Fig.8 Volcano plot of differentially expressed genes in CEF cells infected with La Sota and HB strains

表4 差异表达miRNA 统计Table 4 Statistics of differentially expressed miRNA

2.5.3 GO 功能注释和KEGG 信号通路富集分析为了解NDV 不同毒株感染CEF 后参与的体内生物学过程，对差异表达的miRNA 靶基因进行GO 功能注释，挑选P-value 最小即富集

将最显著的前20 个GO term 条目进行展示（表5），发现细胞成分CC 主要富集在：细胞器、细胞器膜以及囊泡上；生物过程BP 主要富集在：细胞以及蛋白定位和细胞成分的调节；分子功能MF 主要富集在核苷酸结合与GTP 结合上。KEGG 富集结果是通过 Rich factor、FDR 值和富集到此 pathway 上的miRNA 靶基因个数来衡量富集的程度，挑选富集最显著的前20 个KEGG 通路（表6），发现主要集中在细胞代谢、生物合成以及应激信号转导等通路上。

表5 差异表达miRNA 靶基因G0 功能富集分析Table 5 GO Functional enrichment analysis of differentially expressed miRNA target genes

表6 差异表达miRNA 靶基因KEGG 富集分析Table 6 KEGG enrichment analysis of differentially expressed miRNA target genes

2.6 RT-qPCR 验证

随机选取5 个miRNA，分别为：gga-miR-204、gga-miR-155、gga-miR-3536、gga-miR-223、ggamiR-1680-5p，通过实时荧光定量PCR 方法对随机选取的miRNA 的表达情况进行验证（图9），可以看出RT-qPCR 结果与测序RNA-seq 结果趋势一致，证明了测序结果的可靠性。

图9 RT-qPCR 验证Fig.9 RT-qPCR verification

3 讨论与结论

随着对miRNA 调控功能的不断研究，越来越多的miRNA 被证明参与调控宿主与病原体之间的相互作用，有研究证明miRNA 可以通过靶向M 以及L 蛋白，对NDV 的复制造成抑制；miRNA 也参与调控NDV 在Hela 细胞中的复制［11-14］。

在 HB 感染组共发现554 条miRNA 成熟体、440 条miRNA 前 体 以 及45 条novel miRNA；La Sota 感染组共发现613 条miRNA 成熟体、479 条miRNA 前体以及56 条novel miRNA，丰富了鸡源miRNA 数据库。通过对miRNA 碱基偏好性分析表明，20 ～24 nt 首位碱基主要以U 为主，研究表明当首位碱基为U 时有利于对miRNA 前体进行识别与切割［15］；成熟miRNA2 ～8 位点为种子序列可以靶向miRNA 的3’端进行互补配对并对其发挥降解作用［16］，测序结果表明2 个文库中2 ～8 位碱基均以G 为主。

通过GO 和KEGG 途径分析显示，不同NDV毒株感染在宿主的生物合成以及代谢途径上具有明显富集，因为病毒感染过程中通常会掠夺宿主细胞的营养从而改变宿主细胞的代谢网络［17］。Zhou 等［18］推测gga-miR-3538 可能通过参与病毒-载体相互作用、氧化磷酸化以及细胞生长来发挥作用。Lim等［19］研究表明gga-miR-1764 可以参与葡萄糖、电解质和氨基酸等营养物质的运输过程；La Sota 感染组与HB 感染组在C 型凝集素受体信号通路（CLRs）上具有明显的富集，如C 型凝集素受体CD69 就可以通过调节miR-155 及其靶蛋白来调控T 细胞的发育与激活；MAPKs 通路主要通过Raf/MEK/ERK 信号刺激基因产生特异性变化，从而调节细胞的分化与增殖［20］，也有研究证明了gga-miR-142-3p 通过靶向TAB2 负调控NF-κB 和MAPKs 信号通路减轻炎症，通过促进细胞增殖抑制细胞凋亡从而抵御鸡支原体的感染。根据以上研究对差异表达miRNA功能进行推测：gga-miR-3538 与gga-miR-1764-3p 可能在NDV 感染过程中通过调控宿主的生物合成、代谢以及营养运输过程中途径发挥作用；ggamiR-155 在感染NDV 不同毒株后，均有不同程度的上调，证明gga-miR-155 可能通过CLRs 参与适应性免疫的激活；La Sota 感染组gga-miR-142-3p 表达量上调，推测弱毒La Sota 株可能通过上调ggamiR-142-3p 减轻感染后的炎症反应，抑制细胞凋亡，这也与NDV 弱毒株在感染后造成的炎症反应以及临床症状较轻相符合。

miRNA 可以通过调控细胞凋亡进而影响病毒在细胞内的复制。如Li 等［21］表明gga-miR-233 通过靶向FOCO3抑制鸡支原体在DF-1 细胞中的增殖并诱导细胞凋亡；Wang 等［22］证实了miR-365 通过抑制IGF1 的表达促进细胞凋亡； Xu 等［23］通过试验证明miR-187-5p 通过抑制蛋白的表达，进而抑制细胞凋亡。在本研究中，测序结果表明gga-miR-233表达量上调，gga-miR-365-1-5p、gga-miR187-5p 表达下调。此前研究结果表明NDV V 蛋白可以抑制细胞凋亡，并且V 蛋白抑制细胞凋亡的能力与毒力有关［24］。表明筛选出的差异表达miRNA 在NDV 感染中可能通过与靶基因相互结合来调节细胞凋亡，并且强弱毒之间的凋亡水平并不相同。

综上所述，本研究通过高通量测序技术，系统对比分析La Sota 株与HB 株感染CEF36 h 后miRNA 的表达情况，共筛选出10 个差异表达的miRNA，通过与先前报道的研究相比较推测这些差异表达miRNA 可能在调控代谢以及细胞凋亡等方面发挥重要作用，为进一步揭示NDV 感染CEF 以及不同毒株在感染过程中miRNA 发挥的调控作用提供了研究基础和理论依据。