牛白血病病毒与牛泡沫病毒miRNAs在自然感染牛中的共表达及其对肿瘤发生与免疫调控的协同作用机制研究

本研究发现自然感染牛白血病病毒（BLV）的牛体内共表达7种BLV来源和3种BFV来源的miRNAs，并鉴定出7个差异表达的内源性miRNAs。通过蛋白互作网络分析揭示病毒与宿主miRNAs靶基因在免疫应答（如细胞因子信号通路）、肿瘤发生（如MAPK信号通路）和染色质重塑等关键通路上存在功能协同，为BLV/BFV共感染致病机制提供了新视角，对开发诊断标志物和干预策略具有重要价值。

ABSTRACT

病毒编码的microRNAs（miRNAs）在病毒-宿主相互作用中扮演关键角色。在感染牛群的病毒中，牛白血病病毒（BLV）和牛泡沫病毒（BFV）均能编码miRNAs，其中BLV miRNAs因在诱导淋巴瘤中的作用备受关注。BLV在美洲流行率高，对奶牛业造成严重经济损失。BLV与BFV共感染在牛群中常见，且泡沫病毒具有跨物种传播潜力，防范其进入食物链至关重要。本研究在自然感染BLV的牛中鉴定出7种BLV来源miRNAs（blv-miR-B1-3p、blv-miR-B2-5p、blv-miR-B2-3p、blv-miR-B3-5p、blv-miR-B3-3p、blv-miR-B4-3p和blv-miR-B5-5p）和3种BFV来源miRNAs（bfv-miR-BF1-5p、bfv-miR-BF1-3p和bfv-miR-BF2-5p）的共表达。此外，发现7个牛源miRNAs（bta-miR-375、bta-miR-133a、bta-miR-677、bta-miR-1、bta-miR-3613a、bta-miR-9-5p和bta-miR-95）在高BLV前病毒载量（PVL）牛与未感染牛间差异表达（|倍数变化| > 1.5且q值 < 0.05）。对病毒和宿主miRNAs靶基因的蛋白互作网络（PPI）分析表明，这些基因功能汇聚于免疫应答、肿瘤发生和染色质重塑等关键通路。尽管BLV和BFV miRNAs靶向不同基因，但这些靶标参与共享生物学过程，这种汇聚可能反映病毒对宿主功能的协同调控，影响疾病进展和BFV传播。这些发现为开发诊断和治疗策略以控制牛病毒持久性和肿瘤发生提供了新视角。

IMPORTANCE

BLV和BFV是感染牛的反转录病毒，均编码miRNAs。BFV miRNAs的功能尚不明确，而BLV miRNAs因潜在参与致癌过程受到关注。BLV-BFV共感染常见，且泡沫病毒有跨物种屏障潜力，防范其进入食物链至关重要。本研究首次报道自然感染牛中3种BFV miRNAs和7种BLV miRNAs的共表达。病毒和宿主miRNAs靶基因的PPI网络分析揭示了与肿瘤发生、免疫调控和染色质重塑相关的关键代谢通路。虽然BLV和BFV miRNAs靶向不同基因，但这些靶标参与共同生物学过程，提示功能汇聚可能影响疾病进展和BFV传播。这些发现为开发控制牛病毒持久性和肿瘤发展的诊断和治疗策略提供了机遇。

INTRODUCTION

病毒-宿主相互作用研究传统上聚焦于病毒蛋白调控宿主蛋白以劫持细胞过程促进病毒复制和致病性。然而，与宿主细胞类似，许多病毒编码调控性非编码RNA，包括小RNA如miRNAs，这些miRNAs为调控宿主基因表达提供了新能力以确保病毒感染成功。miRNAs是酶加工的小RNA，通常长19-24核苷酸（nt），在转录后水平精细调控基因表达。miRNAs的调控功能主要由RNA诱导沉默复合体（RISC）介导，该复合体促进成熟miRNA的2-7核苷酸种子区域与靶mRNA的3'非翻译区（3'UTR）碱基配对。miRNA-mRNA相互作用诱导细胞核酸酶介导的mRNA降解或翻译抑制。估计miRNAs可能调控至少60%的人类编码基因组。

宿主来源miRNAs通过展现既可益于病毒也可益于宿主的双重性质调控病毒感染。例如，直接靶向病毒RNA的宿主miRNAs常增强病毒生命周期。反之，宿主miRNAs可通过靶向编码病毒周期所需必需宿主因子或建立免疫应答和防御机制的mRNAs施加间接负面影响。鉴于其精细调控基因表达的高多样性和小编码尺寸，病毒进化出表达自身miRNAs的能力不足为奇。自首例爱泼斯坦-巴尔病毒编码miRNAs报道以来，DNA病毒编码的miRNAs已被证明是最丰富的，尤其是疱疹病毒科成员表达的miRNAs。相比之下，表达miRNAs的RNA病毒远不常见，通常限于那些具有DNA中间体的病毒。此外，在正链和负链RNA病毒中均鉴定出若干miRNA样分子。

通常，病毒miRNAs遵循与宿主miRNAs相同的生物发生和效应通路。然而，某些反转录病毒，如泡沫病毒科成员（如猿泡沫病毒[SFV]和牛泡沫病毒[BFV]）和正反转录病毒亚科成员（如牛白血病病毒[BLV]），通过采用RNA聚合酶III而非RNA聚合酶II解耦miRNA转录与基因组转录。该机制可能防止miRNA生物发生过程中核酸酶介导的病毒基因组自降解。

一般而言，病毒miRNAs的功能作用涉及下调病毒和宿主蛋白；然而，上调也会发生，要么通过调控作为转录抑制因子的基因间接实现，要么直接实现，因为某些miRNAs可通过结合特定DNA序列或与增强基因表达的因子相互作用而发挥转录激活因子功能。总体而言，这些行动调控基因表达以调节病毒生命周期（如促进或抑制细胞增殖和凋亡）、逃避宿主免疫应答并建立潜伏或持久感染。这在反转录病毒和疱疹病毒中尤为明显，其中病毒蛋白表达在潜伏期高度受限。低抗原性调控分子的产生为病毒持久性同时逃避宿主免疫检测提供了显著优势。

感染牛且编码miRNAs的病毒包括BLV、BFV、牛疱疹病毒1型（BoHV-1）和牛疱疹病毒5型（BoHV-5）。BoHV-1编码至少10种miRNAs，其中两种在体内和体外均涉及维持潜伏期。相比之下，三种BFV来源miRNAs的功能作用和靶基因仍知之甚少。在BLV编码的miRNAs中，blv-miR-B4-3p是研究最深入的。它与miR-29a共享种子序列，后者是一种oncomiR，靶向肿瘤抑制因子如过氧化物酶同源物（PXDN）和HMG盒转录因子1（HBP1），这在鼠类肿瘤诱导研究中已有报道。体外报告基因检测证明blv-miR-B4-3p对HBP1和PXDN的表达具有抑制效应。然而，在绵羊原发性肿瘤B细胞中，blv-miR-B4-3p与HBP1表达呈负相关，而PXDN水平保持不变。此外，近期发现表明，在自然感染BLV并表达blv-miR-B4-3p的牛中，与未感染对照相比，PXDN表达显著下调。值得注意的是，两组中miR-29a表达均未受影响。

BLV流行率在北美农场水平高达85%，在阿根廷主要乳区泌乳牛中超过80%。来自受影响畜群的畜产品国际贸易限制对经济产生负面影响，但更大损失源于产奶量减少和与无BLV畜群相比无症状BLV携带牛的较早淘汰。

涉及BLV和BFV的混合感染在牛中似乎常见。鉴于泡沫病毒跨物种屏障并引起人类人畜共患感染的能力，其在食物链中的存在可能对人类和动物健康构成潜在风险。此外，感染BLV的牛更易共感染BoHV-1。

本研究通过小RNA测序分析在自然感染BLV的牛中鉴定出三种BFV来源miRNAs和七种BLV来源miRNAs的表达。未发现BoHV-1/5的miRNAs。七种牛miRNAs在感染BLV/BFV的牛与未感染BLV的牛之间显示差异表达。对病毒和宿主来源miRNAs潜在靶基因的蛋白-蛋白相互作用（PPI）网络功能分析揭示了与肿瘤发生和免疫应答调控相关的关键代谢通路。

MATERIALS AND METHODS

Selection of animals

BLV感染和非感染样本来自先前表型分析的群体。简言之，从阿根廷中部地区一奶牛场（该地BLV个体平均流行率超过80%）初筛129头成年荷斯坦奶牛（3岁以上，共享相同泌乳期）。由于据报道抗体水平反映体内前病毒载量（PVL），在最终采样前10个月（T1）和5个月（T2）评估抗BLV酶联免疫吸附测定（抗BLV ELISA）（见下文）；平均反应百分比（PR）在T1为122.7 ± 34.8，在T2为146.3 ± 55.6。选择PR最高四分位数（Q；Q4：T1 = 148.6–178.6%，T2 = 194.4–239.6%）和两次检测均阴性的动物，通过定量PCR（qPCR）（后详）在3个月（T3）和采样时（T4）进一步进行PVL量化。基于PVL结果的一致性，最终选择四头持续高PVL（HPVL）奶牛和三头持续qPCR阴性奶牛进行小RNA测序。

Isolation of peripheral blood mononuclear cells

通过颈静脉穿刺采集动物新鲜血样并补充EDTA（225 μM）。当天使用Ficoll-Paque Plus（GE Healthcare, Uppsala, Sweden）密度梯度离心分离外周血单核细胞（PBMCs），遵循制造商协议。保留血浆部分用于抗BLV ELISA血清学。分离后，PBMCs保存在RNAlater溶液（Ambion, Austin, TX）中并于?80°C储存直至使用。

Serology

采用Trono等描述的抗BLV ELISA测定鉴定BLV感染动物。全BLV病毒颗粒作为抗原检测每个样本血浆抗BLV抗体。简言之，使用公式计算样本阳性（S/P）比，称为PR：S/P = [(OD_样本 ? OD_WS?)/(OD_WS+1/7 ? OD_WS?)] × 100，其中OD指光密度，WS?指示阴性血清，WS+1/7是国际标准弱阳性对照血清（以1:7在阴性血清中稀释）。PR大于25%的样本视为阳性（BLV[+]），低于则视为阴性[BLV(?)]。

BLV PVL quantification

使用High Pure PCR Template Preparation Kit（Roche, Penzberg, Germany）从PBMCs提取基因组DNA，遵循制造商说明。使用Blood Genomic DNA AxyPrep kit（Axygen Biosciences, Union City, CA, USA）从全血样本提取的基因组DNA的质量和浓度，使用微量分光光度计（NanoDrop Technologies, Inc., Wilmington, DE, USA）评估。

基于SYBR Green染料检测系统的BLV POL基因PVL qPCR测定按照Petersen等描述进行。每个25 μL qPCR反应包含Fast Start Universal SYBR Green Master Mix（2×；Roche）、正向和反向引物（800 nM；BLVpol_5 f:5′-CCTCAATTCCCTTTAAACTA-3′和BLVpol_3 r:5′-GTACCGGGAAGACTGGATTA-3′；Thermo Fisher Scientific）和200 ng基因组DNA模板。使用Step One Plus系统（Applied Biosystems, Foster City, CA）进行扩增和检测。每个BLV阳性反应的特异性通过熔解温度解离曲线（Tm）分析确认。基于低自然感染率（外周血中1% BLV感染细胞）的假设，PVL值低于1,500拷贝/μg总DNA分类为低，高于该阈值分类为高。

Small RNA sequencing

使用miRNeasy kit（Qiagen, Hilden, Germany）从PBMCs提取总RNA，根据制造商协议。RNA质量（完整性）和浓度通过Agilent 2200 TapeStation系统（Agilent Technologies, Santa Clara, CA, USA）上的数字电泳评估。

使用QIAseq miRNA Library Kit（Qiagen, Hilden, Germany）进行miRNA测序文库构建，该试剂盒在cDNA合成步骤中将独特分子标识符（UMIs）标签整合到适配器中。这些UMI标签通过允许识别和合并具有相同扩增起源的读段，有助于减轻PCR扩增偏倚和测序伪影。为每个样本分配特定8-nt条形码用于多重测序。将多重文库池化至等摩尔浓度，并在NovaSeq运行（Illumina, San Diego, CA, USA）中以75 × 50 bp配置（每车道300?400百万读段对）处理。

Small RNA sequencing reads processing and deduplication

依次进行测序读段的质量控制以确保可靠读段用于无偏miRNome识别和量化。使用umi_tools v1.6软件从每个读段提取12-nt UMI模式。在合并前独立处理正向（Fw）和反向（Rv）读段。使用Cutadapt v5.0移除低质量（质量截断 < 30）、最小长度 <16 nt、存在模糊碱基（N）或适配器序列的读段。

丢弃QC处理产生的单例读段，仅考虑成对Fw和Rv读段。使用BBMap v25.85将正向读段序列与Rv互补序列比对并合并为单一一致读段，最小插入大小16 nt，最小重叠18 nt。合并后，使用Cutadapt v5.0移除Qiagen特异性适配器序列。

采用严格映射策略使用Bowtie2 v5.4进行读段去重。将读段比对到牛参考基因组ARS-UCD1.2（acc. GCF_002263795.1）和编码感染牛miRNAs的病毒基因组，包括BLV（acc. NC_001414.1）、BFV（acc. NC_001831.1）、BoHV-1（acc. NC_063268.1）和BoHV-5（acc. NC_005261.3）。不允许读段前18个碱基（左端）错配。比对到参考序列后，基于其映射坐标和先前提取的UMI标签使用UMICollapse v0.0进行读段去重。

miRNAs identification and quantification

将去重读段映射到从miRBase v22.1获得的牛、BLV、BFV、BoHV-1和BoHV-5的参考前体miRNA序列。使用默认参数通过miRDeep2 v0.3量化成熟miRNA表达，为每个样本生成miRNA表达计数矩阵。

Differential miRNA expression analysis

使用DESeq2 v0.0分析所有样本的原始计数矩阵以评估BLV感染和非感染动物间的差异miRNA表达。基于miRNA表达谱，使用主成分分析（PCA）对感染和未感染动物间的变异模式进行探索性分析。为确保数据方差齐性，使用DESeq2 v0.0中实施的方差稳定变换方法转换原始miRNA计数。基于miRNA表达水平的均方根偏差计算样本间成对距离。

为检验差异miRNA表达，DESeq2模型通过负二项模型调整miRNA计数以考虑样本特定差异，如文库大小、测序组成和miRNA特定偏倚。在零假设下，模型假定miRNA在两种条件[BLV(?)和BLV(+)组]间具有相同平均表达。使用Wald检验评估条件间log₂倍数变化（Log₂FC）的统计显著性。使用Benjamini-Hochberg错误发现率（BH-FDR）方法进行多重比较校正。BH-FDR调整P值（q） < 0.05且|Log₂FC| > 1.5的miRNAs视为统计显著。在R环境中执行PCA和火山图。

miRNAs gene target prediction

使用miRanda v3.3、PITA v0.0和RNAhybrid v1.1，以默认参数，对BLV(+)和未感染牛间显著差异表达（DE）的牛源和病毒来源miRNAs进行潜在牛基因靶标（miRNA:mRNA相互作用）计算机预测。分析使用来自ARS-UCD1.2注释的21,400个牛基因的3'UTR序列。每个基因靶标预测算法采用不同策略：miRanda依赖种子位点比对和RNA-RNA双链的热力学稳定性；PITA整合靶位点可及性；RNAhybrid使用miRNA-mRNA序列杂交的最小自由能。使用以下阈值过滤预测靶标：得分≥140且自由能≤?20 kcal/mol（miRanda）和自由能≤?20 kcal/mol（PITA和RNAhybrid）。仅保留所有三种工具预测的基因靶标（共识基因靶标）用于下游分析。

Gene ontology and graph-based pathway analysis of miRNA target genes

使用STRINGdb对潜在靶基因进行功能注释，采用基因本体（GO）术语和京都基因与基因组百科全书（KEGG）通路。接下来，使用STRINGdb构建已识别靶基因的PPI网络。该平台整合来自多个来源的策划和预测相互作用，包括实验证据、计算预测、科学文献、共表达数据和主要公共数据库。使用默认参数（膨胀参数 = 5）的马尔可夫聚类算法识别网络内高度互连的蛋白质簇。仅考虑超过15个节点（基因）的簇。使用基于超几何分布的Fisher精确检验对识别簇进行GO术语和KEGG通路过度表征分析，术语相似性阈值>0.7。使用R图形环境生成GO术语类别-基因靶标图以可视化显著富集GO术语与其相关基因间的连接。

RESULTS

在本研究中，四头动物在所有采样时间均表现高抗体反应性（IND_5671 = 178.6 ± 28.4；IND_5841 = 191.3 ± 32.8；IND_6021 = 169.8 ± 42.7；和IND_6097 = 183.2 ± 25.4）以及在T-3和T0高前病毒载量（平均前病毒载量：IND_5671 = 63,148.0；IND_5841 = 44,234.3；IND_6021 = 52,902.4；和IND_6097 = 91,255.4）。这些奶牛分类为BLV(+)。相比之下，三头抗BLV ELISA和qPCR持续阴性的动物分类为BLV(?)。

确认总RNA质量（RIN 8?9.6），并使用UMI标签适配器制备小RNA文库。通过多步处理流程获得高质量、非冗余读段。每个样本的描述性特征、质量控制前后测序读段数量、合并和去重以及映射统计汇总于表1。

去多重读段映射到从miRBase v22.1获得的牛、BLV、BFV、BoHV-1和BoHV-5的参考前体miRNA序列（n = 1,060）以识别和量化每个样本中存在的miRNAs。平均在所有样本中识别和量化351种来自牛和病毒来源的miRNAs（映射读段≥5）（表S1）。基于标准化miRNA表达矩阵的样本距离PCA显示BLV(?)样本聚集在一起，形成与BLV(+)样本分开的 distinct 组。PC1解释表达矩阵中62%的方差，而PC2占17%（图1）。

在BLV(+)组中，总共持续表达10种病毒来源miRNAs，包括来自BLV的七种：blv-miR-B1-3p（平均表达 = 27,681.1）、blv-miR-B2-5p（平均表达 = 59,058.0）、blv-miR-B2-3p（平均表达 = 112.5）、blv-miR-B3-5p（平均表达 = 309.6）、blv-miR-B3-3p（平均表达 = 49,312.2）、blv-miR-B4-3p（平均表达 = 65.1）和blv-miR-B5-5p（平均表达 = 951.7）。此外，识别出三种BFV来源miRNAs：bfv-miR-BF1-5p（平均表达 = 6,844.7）、bfv-miR-BF1-3p（平均表达 = 8,257.3）和bfv-miR-BF2-5p（平均表达 = 14,920.6）。相比之下，病毒来源miRNA表达在BLV(?)组中低（平均表达 <30.0）或不可检测，除blv-miR-B2-5p（平均表达 = 67.4）和blv-miR-B3-3p（平均表达 = 49.8；表2）。

将所有来自BLV(?)样本的读段映射到牛参考基因组（acc. ARS-UCD1.2）显示，14,378,553读段中仅716未能比对，表明几乎所有读段源自牛转录本。此外，来自BLV(?)牛的读段比对到BLV（acc. NC_001414.1）和BFV（acc. NC_001831.1）参考基因组的数量分别仅为600和141读段（表1）。而且，相同读段数据集与识别BLV和BFV miRNAs的比对导致642读段至少有一次比对。其中，仅六读段也映射到牛参考基因组，且无读段比对到miRBase数据库中任何其他已知miRNA（n = 48,871成熟miRNAs）。图S1显示来自BLV(+)和BLV(?)牛的所有读段在BLV和BFV参考基因组上的覆盖范围。可检测读段覆盖位于两个病毒基因组的已知miRNA区域。

宿主miRNAs的差异表达分析识别出总共七种牛miRNAs（bta）具有|倍数变化|（FC） >1.5且q值 <0.05，作为BLV(+)和BLV(?)组间的差异表达（BTA-miRNAs-DE）：bta-miR-375、bta-miR-133a、bta-miR-677、bta-miR-1、bta-miR-3613a、bta-miR-9-5p和bta-miR-95（图2；表S2）。

为评估7种差异表达BTA-miRNAs和10种病毒来源miRNAs的潜在功能影响，进行了靶基因预测分析。基于三种独立预测工具的共识，识别出1,518个基因作为推定靶标。其中，281、977和260个基因分别由BTA-miRNAs-DE、BLV-miRNAs和BFV-miRNAs独家靶向（表S3）。未观察到预测靶标间的重叠。然后，使用GO术语和代谢通路对1,518个潜在靶基因进行功能注释，揭示91个显著过度表征的GO生物过程（BP）术语、10个GO分子功能术语和30个GO细胞组件术语。KEGG和Reactome代谢通路的过度表征分析揭示四个显著术语：“Ras信号通路”（bta04014, q = 0.04）、“癌症通路”（bta05200, q = 0.02）、“催产素信号通路”（bta04921, q = 0.024）和“库欣综合征”（bta04934, q = 0.04）。此外，术语“信号转导”（BTA-162582, q = 0.0002）和“代谢”（BTA-1430728, q = 0.01）在Reactome通路分析中显著（表S4）。

开发了一个PPI网络，包含1,518个节点（基因）和6,075条边（相互作用；图S2）。从该网络中，识别出三个高度互连的蛋白质簇（C1、C2和C3）（PPI P值 < 1.0e^?16），其可能代表 specialized 生物功能枢纽（图3）。这些簇的功能过度表征分析揭示涉及免疫应答调控、细胞信号传导和癌症发展机制的关键细胞过程（表S4）。有趣的是，尽管BTA-、BLV-和BFV来源miRNAs的靶基因缺乏重叠，但它们的预测功能汇聚于互连枢纽，表明在通路水平存在潜在 interplay。富集图描绘了每簇前15个富集BP GO术语间的互连，通过靶基因突出功能重叠（图4）。例如，C1的显著过度表征GO BP术语（q < 0.05）主要与免疫调控和细胞因子信号传导过程相关，包括“细胞因子介导的信号通路”、“T细胞活化正调控”、“白细胞凋亡过程调控”和“干扰素-γ产生调控”。富集的KEGG和Reactome通路包括“病毒蛋白与细胞因子和细胞因子受体相互作用”（bta04061, q值 = 1.02e^?16）、“免疫系统中的细胞因子信号传导”和“IL-6型细胞因子受体配体相互作用”（BTA-6788467, q值 = 0.0105）等。对于簇C2，主要过度表征的BP术语和代谢通路与细胞生长和增殖相关（如“MAP激酶活性正调控”、“上皮细胞增殖正调控”、“造血”和“PI3K级联”[BTA-109704, q值 = 0.0034]），以及组织发育（如“解剖结构形态发生”、“血管内皮生长因子信号通路”和“受体酪氨酸激酶信号传导”[BTA-9006934, q值 = 3.13e^?16]）。此外，C2富集与细胞存活和凋亡调控相关的过程（“凋亡过程负调控”和“PI5P、PP2A和IER3调控PI3K/AKT信号传导”[BTA-6811558, q值 = 1.62e^?09]）。最后，在簇C3中，富集的功能术语主要与通过染色质修饰的表观遗传控制和转录调控相关，包括“组蛋白去乙酰化”和“染色质重塑”，通路如“通过乙酰化调控TP53活性”（BTA-6804758, q值 = 2.85e^?08）、“基因表达（转录）”（BTA-74160, q值 = 3.05e^?05）和“染色质修饰酶”（BTA-3247509,