苜蓿miRs表达谱分析及跨界潜力miRs初步筛选

贾晶莹， 李雅辉， 伏兵哲， 马云， 蔡小艳

（宁夏大学农学院，宁夏回族自治区反刍动物分子细胞育种重点实验室，银川 750021）

microRNA（miR）大量存在于动植物体内，是由21~23 个核苷酸组成的非编码小分子RNA[1‑2]。miR 通过碱基互补配对的方式与靶基因（mRNA）结合，是哺乳动物和植物中通过调控转录后水平基因的表达而发挥作用的主要分子[3]。随着miRs研究的深入，其在不同生命体间进行传导的现象也被发现，植物miRs跨物种调控动物靶基因的重要作用逐渐被人们所认识。Zhang 等[4]率先发现植物miRs可以通过饮食进入动物体内，并通过与动物内源性靶mRNA 结合而发挥跨界调控作用。研究表明，来自大豆和西兰花的miR159具有抑制癌细胞生长的能力[5]；金银花水煎液中的miR2911可被感染甲型流感病毒（IAV）的小鼠通过胃肠道吸收，并通过血液循环进入感染IAV 小鼠的肺部发挥调控作用[6]；并且，金银花中的miR2911 能显著抑制SARS-CoV-2的复制，对新型冠状病毒肺炎具有良好的治疗效果[7]。

植物miRs 可通过跨界调控动物生理功能，而玉米、大豆和苜蓿等植物性饲料是牛羊猪等家畜采食的主要饲草料，相关学者综合这2 个方面联系也开展了一些饲料作物跨界调控的研究。Luo等[8]研究发现，给猪饲喂一段时间的玉米后，可在猪的组织和血清中检测到玉米源的miRs，并且这些玉米源miRs 具有调控猪基因表达的潜力。Marzano 等[9]研究表明，苜蓿miR-5754和大豆miR-4995 作用于转移相关肺腺癌转录物1 (MALAT1)和核副癌细胞组装转录物1(NET1)，从而降低这些致癌靶转录物的稳定性，而转录物的产物可促进细胞增殖。Zhang 等[4]在以青贮饲料为主食的动物血液中也检测到植物miRs,且其表达谱更像是与食物有关；还在小牛、胎牛和马等动物血液中检测到miR-168、miR-156 和miR-166 的表达。苜蓿是保证奶牛高产奶量和高营养品质形成的主要饲草。研究表明，食物及饲料来源的miRs可以改变奶牛组织和循环中的miRs表达谱，从而对人体和牛奶乳蛋白产量产生影响[10‑11]。由此引发2 个科学问题：一是在牛、马等动物体内检测到的植物源miR-168、miR-156 和miR-166 是否主要来自苜蓿？二是不同苜蓿品种间miRs 是否具有较大差异，这种差异是否会导致动物体内miRs表达谱发生改变？

基于对上述问题的思考，本研究分别选取‘新盐52 号’（Xinyan 52，XY52）和‘ 中苜1 号’（Zhongmu 1，ZM1）进行RNA-Seq，构建2个苜蓿品种的miRs 组学表达谱，筛选差异表达miRs，预测miRs 靶基因及作用信号通路，并运用RT-qPCR 技术对筛选出的5个可跨界调控的miRs和5个差异表达的miRs进行定量分析。以期筛选出‘新盐52号’和‘中苜1号’苜蓿中具有研究潜力的miRs，为后续研究苜蓿中可能对奶牛产生跨物种调控作用的miRs提供筛选基础。

1 材料与方法

1.1 试验材料

试验材料为‘中苜1号’（ZM1）和‘新盐52号’（XY52）苜蓿，均于2020 年10 月采自宁夏大学平吉堡科技园区。‘新盐52 号’苜蓿是在宁夏培育的苜蓿耐盐耐旱新品系，在宁夏具有较大的开发利用价值；‘中苜1号’是审定登记的国家牧草品种，兼具耐盐和耐旱的特性[12]，在宁夏表现较佳，因此选择上述2 个苜蓿品种作为研究对象。采集初花期的苜蓿地上部分，并快速用手术剪将样品剪碎放入冻存管，然后投入到液氮罐中，带回实验室- 80 ℃冰箱冻存备用。

Trizol 试剂、反转录试剂盒、荧光定量试剂盒购自Takara公司（大连）。

1.2 试验方法

1.2.1苜蓿中总RNA 提取与反转录 样品用液氮预冷过的研钵进行研磨，利用Trizol 试剂对样品进行总RNA 提取。使用多功能酶标仪（BioTeK SynergyLX）测定总RNA 质量浓度，利用2%琼脂糖凝胶电泳检测总RNA 的质量及完整性。按照反转录试剂盒说明书将提取的总RNA 反转成cDNA，保存于−20 ℃冰箱，备用。每个品种3次重复，将测序样品分别编号为XY52-1、XY52-2、XY52-3和ZM1-1、ZM1-2和ZM1-3。

1.2.2测序数据质量控制 样品转录组测序委托北京百迈克生物科技有限公司完成。过滤测序原始数据，去掉低质量及较短（≤18个）或较长（≥30个）核苷酸的序列；去除未知碱基N含量≥10%的读长等，得到高质量序列。

1.2.3miRs 鉴定 将高质量序列分别与Rfam、Silva 等数据库进行比对，过滤核糖体RNA（rRNA）、转运RNA（tRNA）等非编码RNA（ncRNA）和重复序列，获得含有miRs的未注释的读长。

比对到参考基因组的读长与miRBase（v22）（http://www.mirbase.org/）数据库中已知miRs 的成熟序列及其上游2 nt 与下游5 nt 的范围进行比对，鉴定得到已知miRs；利用miRDeep2 软件，通过调整参数及打分系统进行新miRs的预测。

1.2.4筛选差异表达miRs 使用TPM（transcripts per kilobase of exon model per million mapped reads，每千个碱基的转录每百万映射读取的转录本）算法对2 个品种中miRs 的表达量进行归一化处理，运用DESeq2（https://www.bioconductor.org/packages/release/bioc/html/DESeq2.html）将样品组间的TPM 计算结果进行差异表达分析，根据差异倍数变化（fold change，FC）和错误发现率（false discovery rate，FDR）将|log2(FC)|≥1.00、FDR≤0.01作为筛选差异表达miRs 的标准，得到差异表达miRs集。

1.2.5miRs 靶基因预测及功能注释 利用

Target Finder（http://targetfinder.org/）软件根据已知miRs、新预测miRs和对应物种的基因序列信息进行靶基因预测。使用BLAST 软件将预测到的靶基因序列与KEGG（https://www.kegg.jp/）、GO（http://www.geneontology.org）、Pfam（http://pfam.xfam.org/）等数据库进行比对，得到靶基因的注释信息。使用GOseq R packages 进行GO 富集分析，KOBAS 软件用于KEGG 富集分析，富集分析中，按照q<0.05 进行筛选（q为多重假设检验校正之后的P值，q值越小，富集显著性越可靠）。并对差异表达miRs 靶基因进行Pathway 富集分析，利用富集因子分析通路的富集程度。富集因子计算公式如下。

1.2.6实时荧光定量PCR（RT-qPCR）测定miRs表达量 选择miR156b-5p[13]、miR166a[14]、miR167a[15]、miR168a[4,15]和miR319a-3p[16]5 个可跨界调控动物生理功能的植物miRs，以及5个在2个苜蓿品种中差异表达的miRs进行实时荧光定量检测。

根据测序得到的miRs 序列设计引物，试验所需引物由Primer 5.0 软件进行设计（表1），由通用生物系统（安徽）有限公司合成。将经过质量检验的miRs 特异性反转录成cDNA，以5s 作为内参进行qPCR，反应体系详见表2；反应程序如下：95 ℃3 min；95 ℃ 5 s，60 ℃ 30 s，39 个循环；95 ℃ 10 s。实时荧光定量PCR 结果用2−ΔΔCt[5]进行处理，使用GraphPad Prism 7 软件进行作图并进行差异显著性分析。

表1 引物信息Table 1 Primer information

表2 荧光定量反应体系Table 2 Fluorescence quantitative reaction system

2 结果与分析

2.1 miRs种类数量和长度分析鉴定

分析miRs 测序组学数据发现，‘中苜1 号’和‘新盐52 号’中分别检测到656 和703 个miRs，其中新预测到的miRs 数量一致，均为223 个（表3）。比较2个苜蓿品种各miRs的TPM值发现，2个苜蓿品种中优势表达的前4个已知miRs均为miR5213-5p、miR159a、miR396a-5p和miR166e-3p。

表3 中苜1号和新盐52号中优势表达的miRsTable 3 Advantage expression miRs in Zhongmu 1 and Xinyan 52

对2 个苜蓿品种中检测到的miRs 长度进行统计，发现已知miRs 的长度主要为21 nt（图1A），新预测miRs长度主要以24 nt为主（图1B）。并且2 种苜蓿新预测的miRs 中均有1 个miR 长度达到25 nt。

图1 miRs的长度分布Fig. 1 Length distribution of miRs

2.2 2个苜蓿品种中差异表达miRs的对比分析

统计分析发现2个苜蓿品种中有21个miRs差异表达。其中已知miRs 10个，新预测miRs 11个,仅novel-miR54 在‘新盐52 号’中表达量较‘中苜1号’显著上调，其余20 个miRs 表达量均下调（图2）。将在2 个苜蓿品种中差异表达且TPM 值丰度排在前5 的miRs 进行序列比对，剔除重复序列，依次补位，最终得到novel-miR54、novel miR158、miR5754、miR156f（miR156e、miR156h-5p）、miR5743a（miR5743b）5个miRs（表4）。

图2 苜蓿差异表达miRs热图Fig. 2 Heat-map of differentially expressed miRs in alfalfa

表4 差异表达miRs序列Table 4 miRs sequence of differentially expressed

2.3 差异表达miRs靶基因GO富集分析

在‘中苜1 号’检测到的656 个miRs 中，有608 个miRs 共预测到7 536 个靶基因。在‘新盐52 号’检测到的703 个miRs 中，有654 个miRs 共预测到7 933个靶基因。在‘中苜1号’和‘新盐52号’全部miRs预测到的靶基因中，有6 783个靶基因与GO数据库比对成功。

在差异表达miRs 预测得到的623 个靶基因中，有503个靶基因与GO数据库比对成功，获得注释信息。其中在‘新盐52 号’中表达量上调的novel-miR54 预测得到116 个靶基因；在‘新盐52号’表达量下调的miRs 中，10 个已知miRs 预测得到255个靶基因，10个新预测miRs预测得到252个靶基因。

差异表达miRs 在分子功能、细胞组分和生物过程3 个层面分别注释到444、372 和349 个靶基因。注释到生物过程层面的靶基因主要涉及的功能或代谢路径为蛋白质磷酸化、氧化还原过程和防御反应。注释到细胞组分层面的靶基因主要涉及的功能或代谢路径为膜的组成成分、核和质膜。注释到分子功能层面的靶基因主要涉及的功能条目为ATP结合、DNA结合和锌离子结合（表5）。差异表达miRs 靶基因与全部miRs 靶基因显著富集的功能条目趋于一致，但在生物过程层面的防御反应和分子功能层面的DNA 结合、ADP 结合等方面富集的靶基因有差异，这些差异通路的靶基因可为揭示2个苜蓿品种功能差异提供基础素材。

表5 miRs靶基因GO功能富集分析Table 5 GO functional enrichment analysis of miRs target genes

2.4 差异表达miRs靶基因KEGG通路富集分析

差异表达的21 个miRs 共预测得到623 个靶基因，其中有215 个靶基因与KEGG 数据库比对成功，获得注释，KEGG 的注释结果按照KEGG 通路类型进行分类，得到靶基因KEGG 通路类型分类。

差异表达miRs 的靶基因主要富集在5 个KEGG 通路类型，分别是细胞过程、环境信息处理、遗传信息处理、新陈代谢和有机系统（表6）。对差异表达miRs的靶基因进行通路富集分析，结果（图3）发现，在靶基因注释到的全部KEGG通路中，内吞作用（6.74%）、RNA 转运（6.74%）、ABC 转运蛋白（5.62%）和泛素介导的蛋白水解（5.62%）为最显著富集通路。

图3 差异表达miRs靶基因KEGG通路富集散点图Fig. 3 Scatter diagram of KEGG pathway enrichment of differentially expressed miRs target gene

表6 差异表达miRs靶基因KEGG通路富集分析Table 6 Target gene KEGG pathway enrichment of differential expressed miRs

2.5 5 个高丰度差异表达miRs 的表达量的检测验证

对 novel-miR54、novel-miR158、miR5754、miR156f （miR156e、miR156h-5p） 、miR5743a（miR5743b）这5 个差异表达的miRs 进行RTqPCR，结果（图4）表明，5个miRs的表达量变化趋势与RNA-seq 结果一致。就表达量而言，上述miRs 仅有novel-miR54 在‘新盐52 号’中的表达量较‘中苜1 号’极显著上调（P<0.01），上调倍数为7.09 倍；其余4 个miRs 在‘新盐52 号’中表达量均表现为下调。在品种方面，‘中苜1号’中miR156f的相对表达量最高，较‘新盐52 号’极显著上调（P<0.01）；‘新盐52 号’中miR5743a 的表达量最高。

图4 差异表达miRs的相对表达量Fig. 4 Relative expression level of differentially expressed miRs

2.6 5个已知跨界功能miRs表达量的验证

对已知有跨界调控作用的miR156b-5p、miR166a、miR167a、miR168c-3p 和miR319a-3p 进行定量分析，结果（图5）表明，同一品种，miR156b-5p、miR167a、miR168c-3p 和miR319a-3p的表达量均低于miR166a，并且RT-qPCR 结果与RNA-Seq 结果具有相同的变化趋势，说明RNASeq 结果可信度高。比较不同品种之间的表达量，仅miR166a 在2 个苜蓿品种中差异显著，其在‘新盐52号’中的表达量极显著上调（P<0.01）。

3 讨论

本研究通过miRs 组学分析，筛选出在2 个不同苜蓿品种中表达量高的miRs，为后续选定可以跨界调控奶牛性状的苜蓿miRs提供了研究基础。

3.1 关于2 种苜蓿miRs 表达谱、差异miRs 及功能代谢通路分析

测序结果表明，‘中苜1 号’和‘新盐52 号’苜蓿中分别检测到656 和703 个miRs，其中已知miRs中TPM 值最高的前4个miRs种类一致，说明2 个苜蓿品种中的高丰度miRs 相似，由此推测在研究miRs 跨界功能时苜蓿品种可以不作为主要因素予以考虑。另外测序和定量分析发现的novel-miRs为研究苜蓿属miRs提供了新素材。

2 个苜蓿品种中仅有21 个差异表达的miRs，这些可能是导致二者功能差异和‘新盐52 号’苜蓿独有特性的miRs。通过比较新品系miRs 库与现有品种的异同，可为鉴别和预判新品系的利用价值奠定分子基础，例如差异表达的miR156。研究表明，miR156过表达基因型在干旱胁迫期间有利于保持植株更高的气孔导度；同时，miR156 可以靶向SPL13，而SPL13 表达水平的降低减少了苜蓿的水分流失[17]，这表明miR156 可显著改善苜蓿植株的耐旱性。从苜蓿miRs 跨界潜力方面分析，2 个苜蓿品种饲喂奶牛时，品系间高丰度表达miRs 可能是导致奶牛体况或乳品质差异的关键miRs，如‘中苜1 号’中的miR5213-5p 及‘新盐52号’中的miR159a 和miR166e。其中miR166e 在2个苜蓿品种中表达丰度均处在前4位，并且其在2个苜蓿品种中的相对表达量差异显著，因此，应用‘中苜1 号’和‘新盐52 号’饲喂奶牛时，miR166e可作为发掘饲喂苜蓿品种与奶牛体内基因表达及生产性能相关的miRs。

2个苜蓿品种差异表达miRs靶基因均显著富集到RNA 转运和ABC 转运蛋白通路。差异表达miRs具有转运RNA的能力，如果其可以跨界发挥作用，则可以通过调节奶牛体内RNA 的转运影响奶牛miRs组织表达谱。ABC 转运蛋白在人、植物和动物等生物体内对脂类的转运及代谢具有重要意义，可从外界转入长链脂肪酸及参与高密度脂蛋白的合成[18‑19]。不同物种间的ABC 转运蛋白结构具有相似性[20]，miRs 序列也具有保守性，因此推测苜蓿源miRs进入奶牛体内，也能够通过靶向奶牛体内ABC 转运蛋白进而调节奶牛体内脂质代谢。但是奶牛摄入苜蓿源miRs 的量是否能够达到发挥作用的量，以及苜蓿miRs靶向的ABC转运蛋白的调节对象是奶牛乳脂还是牛肉脂肪都还需要深入研究。

由于miRs 的TPM 值越高，则丰度越高，被奶牛摄入后留存的可能性越大，因此通过序列比对筛选出novel-miR54、novel-miR158、miR5754、miR156f 和miR5743a 共5 个高丰度差异表达的miRs，RT-qPCR 分析的表达量变化趋势与RNAseq测序结果一致，证明RNA-seq可信度高。对南方型紫花苜蓿的测序结果表明，miR156家族占总miRs 数量的9%，具有较高的表达丰度[21]，与本研究结果一致。另外，本研究测序新发现的miRs中novel-miR54 表达水平最高，值得深入探索其靶基因和功能作用。

3.2 关于2 种苜蓿中已知跨界miRs 的表达量及初步筛选结果

Kammes 等[22]研究发现，饲喂苜蓿可以改善奶牛乳品质，提高奶牛乳脂含量和牛奶中尿素氮含量；饲喂苜蓿会改变奶牛瘤胃、十二指肠、空肠、肝脏和乳腺组织中与乳效率相关miRs的表达，调节乳脂含量[11,23‑24]。这证明采食苜蓿影响奶牛体内miRs 的表达谱，到底是由于苜蓿营养导致的奶牛乳腺miRs变化还是苜蓿中miRs的作用，还需要进一步将本研究测定及筛选出的miRs在奶牛体内血液及组织中的表达量进行检测和分析才能明确。

miR156b-5p、miR166a、miR167a、miR168c-3p和miR319a-3p这5个具有跨界功能的植物miRs在测序和定量分析中均被检测到，植物miRs序列具有保守性[25]，因此推测上述miRs 在苜蓿中也具有跨界的潜力。由于高表达量的miRs被奶牛摄入后在奶牛体内具有重要作用的可能性更大，因此通过RNA-seq 和RT-qPCR 最终筛选出在2 个苜蓿品种中均高表达的miR166a，为进一步研究苜蓿中miRs在奶牛体内发挥调控作用提供研究对象。

植物miR166 在其他物种的跨界调控作用已有相关报道，王鹏俊[26]研究发现，给猪饲喂玉米后，在猪的血清中检测到玉米源的miR166a-3p；贾凌[27]研究发现，给蚕饲喂桑叶后，在蚕的血淋巴和组织中发现了桑树源的miR166b 和miR166c；Lukasik等[28]在健康志愿者的母乳中检测到5种植物源的miRs，其中包含miR166a。ath-miR166a 在人类样品和猪母乳外泌体中丰度均为最高[15]；人参水煎剂中的miR166 对小鼠的免疫基因有一定影响[14]。由此可见，苜蓿中miR166a 对奶牛的分子作用机制有重要的研究价值和应用潜力，值得进一步挖掘。