苜蓿miRs表達(dá)譜分析及跨界潛力miRs初步篩選

賈晶瑩， 李雅輝， 伏兵哲， 馬云， 蔡小艷

（寧夏大學(xué)農(nóng)學(xué)院，寧夏回族自治區(qū)反芻動物分子細(xì)胞育種重點實驗室，銀川 750021）

microRNA（miR）大量存在于動植物體內(nèi)，是由21~23 個核苷酸組成的非編碼小分子RNA[1?2]。miR 通過堿基互補(bǔ)配對的方式與靶基因（mRNA）結(jié)合，是哺乳動物和植物中通過調(diào)控轉(zhuǎn)錄后水平基因的表達(dá)而發(fā)揮作用的主要分子[3]。隨著miRs研究的深入，其在不同生命體間進(jìn)行傳導(dǎo)的現(xiàn)象也被發(fā)現(xiàn)，植物miRs跨物種調(diào)控動物靶基因的重要作用逐漸被人們所認(rèn)識。Zhang 等[4]率先發(fā)現(xiàn)植物miRs可以通過飲食進(jìn)入動物體內(nèi)，并通過與動物內(nèi)源性靶mRNA 結(jié)合而發(fā)揮跨界調(diào)控作用。研究表明，來自大豆和西蘭花的miR159具有抑制癌細(xì)胞生長的能力[5]；金銀花水煎液中的miR2911可被感染甲型流感病毒（IAV）的小鼠通過胃腸道吸收，并通過血液循環(huán)進(jìn)入感染IAV 小鼠的肺部發(fā)揮調(diào)控作用[6]；并且，金銀花中的miR2911 能顯著抑制SARS-CoV-2的復(fù)制，對新型冠狀病毒肺炎具有良好的治療效果[7]。

植物miRs 可通過跨界調(diào)控動物生理功能，而玉米、大豆和苜蓿等植物性飼料是牛羊豬等家畜采食的主要飼草料，相關(guān)學(xué)者綜合這2 個方面聯(lián)系也開展了一些飼料作物跨界調(diào)控的研究。Luo等[8]研究發(fā)現(xiàn)，給豬飼喂一段時間的玉米后，可在豬的組織和血清中檢測到玉米源的miRs，并且這些玉米源miRs 具有調(diào)控豬基因表達(dá)的潛力。Marzano 等[9]研究表明，苜蓿miR-5754和大豆miR-4995 作用于轉(zhuǎn)移相關(guān)肺腺癌轉(zhuǎn)錄物1 (MALAT1)和核副癌細(xì)胞組裝轉(zhuǎn)錄物1(NET1)，從而降低這些致癌靶轉(zhuǎn)錄物的穩(wěn)定性，而轉(zhuǎn)錄物的產(chǎn)物可促進(jìn)細(xì)胞增殖。Zhang 等[4]在以青貯飼料為主食的動物血液中也檢測到植物miRs,且其表達(dá)譜更像是與食物有關(guān)；還在小牛、胎牛和馬等動物血液中檢測到miR-168、miR-156 和miR-166 的表達(dá)。苜蓿是保證奶牛高產(chǎn)奶量和高營養(yǎng)品質(zhì)形成的主要飼草。研究表明，食物及飼料來源的miRs可以改變奶牛組織和循環(huán)中的miRs表達(dá)譜，從而對人體和牛奶乳蛋白產(chǎn)量產(chǎn)生影響[10?11]。由此引發(fā)2 個科學(xué)問題：一是在牛、馬等動物體內(nèi)檢測到的植物源miR-168、miR-156 和miR-166 是否主要來自苜蓿？二是不同苜蓿品種間miRs 是否具有較大差異，這種差異是否會導(dǎo)致動物體內(nèi)miRs表達(dá)譜發(fā)生改變？

基于對上述問題的思考，本研究分別選取‘新鹽52 號’（Xinyan 52，XY52）和‘ 中苜1 號’（Zhongmu 1，ZM1）進(jìn)行RNA-Seq，構(gòu)建2個苜蓿品種的miRs 組學(xué)表達(dá)譜，篩選差異表達(dá)miRs，預(yù)測miRs 靶基因及作用信號通路，并運用RT-qPCR 技術(shù)對篩選出的5個可跨界調(diào)控的miRs和5個差異表達(dá)的miRs進(jìn)行定量分析。以期篩選出‘新鹽52號’和‘中苜1號’苜蓿中具有研究潛力的miRs，為后續(xù)研究苜蓿中可能對奶牛產(chǎn)生跨物種調(diào)控作用的miRs提供篩選基礎(chǔ)。

1 材料與方法

1.1 試驗材料

試驗材料為‘中苜1號’（ZM1）和‘新鹽52號’（XY52）苜蓿，均于2020 年10 月采自寧夏大學(xué)平吉堡科技園區(qū)。‘新鹽52 號’苜蓿是在寧夏培育的苜蓿耐鹽耐旱新品系，在寧夏具有較大的開發(fā)利用價值；‘中苜1號’是審定登記的國家牧草品種，兼具耐鹽和耐旱的特性[12]，在寧夏表現(xiàn)較佳，因此選擇上述2 個苜蓿品種作為研究對象。采集初花期的苜蓿地上部分，并快速用手術(shù)剪將樣品剪碎放入凍存管，然后投入到液氮罐中，帶回實驗室- 80 ℃冰箱凍存?zhèn)溆谩?/p>

Trizol 試劑、反轉(zhuǎn)錄試劑盒、熒光定量試劑盒購自Takara公司（大連）。

1.2 試驗方法

1.2.1苜蓿中總RNA 提取與反轉(zhuǎn)錄 樣品用液氮預(yù)冷過的研缽進(jìn)行研磨，利用Trizol 試劑對樣品進(jìn)行總RNA 提取。使用多功能酶標(biāo)儀（BioTeK SynergyLX）測定總RNA 質(zhì)量濃度，利用2%瓊脂糖凝膠電泳檢測總RNA 的質(zhì)量及完整性。按照反轉(zhuǎn)錄試劑盒說明書將提取的總RNA 反轉(zhuǎn)成cDNA，保存于?20 ℃冰箱，備用。每個品種3次重復(fù)，將測序樣品分別編號為XY52-1、XY52-2、XY52-3和ZM1-1、ZM1-2和ZM1-3。

1.2.2測序數(shù)據(jù)質(zhì)量控制 樣品轉(zhuǎn)錄組測序委托北京百邁克生物科技有限公司完成。過濾測序原始數(shù)據(jù)，去掉低質(zhì)量及較短（≤18個）或較長（≥30個）核苷酸的序列；去除未知堿基N含量≥10%的讀長等，得到高質(zhì)量序列。

1.2.3miRs 鑒定 將高質(zhì)量序列分別與Rfam、Silva 等數(shù)據(jù)庫進(jìn)行比對，過濾核糖體RNA（rRNA）、轉(zhuǎn)運RNA（tRNA）等非編碼RNA（ncRNA）和重復(fù)序列，獲得含有miRs的未注釋的讀長。

比對到參考基因組的讀長與miRBase（v22）（http://www.mirbase.org/）數(shù)據(jù)庫中已知miRs 的成熟序列及其上游2 nt 與下游5 nt 的范圍進(jìn)行比對，鑒定得到已知miRs；利用miRDeep2 軟件，通過調(diào)整參數(shù)及打分系統(tǒng)進(jìn)行新miRs的預(yù)測。

1.2.4篩選差異表達(dá)miRs 使用TPM（transcripts per kilobase of exon model per million mapped reads，每千個堿基的轉(zhuǎn)錄每百萬映射讀取的轉(zhuǎn)錄本）算法對2 個品種中miRs 的表達(dá)量進(jìn)行歸一化處理，運用DESeq2（https://www.bioconductor.org/packages/release/bioc/html/DESeq2.html）將樣品組間的TPM 計算結(jié)果進(jìn)行差異表達(dá)分析，根據(jù)差異倍數(shù)變化（fold change，F(xiàn)C）和錯誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）將|log2(FC)|≥1.00、FDR≤0.01作為篩選差異表達(dá)miRs 的標(biāo)準(zhǔn)，得到差異表達(dá)miRs集。

1.2.5miRs 靶基因預(yù)測及功能注釋 利用

Target Finder（http://targetfinder.org/）軟件根據(jù)已知miRs、新預(yù)測miRs和對應(yīng)物種的基因序列信息進(jìn)行靶基因預(yù)測。使用BLAST 軟件將預(yù)測到的靶基因序列與KEGG（https://www.kegg.jp/）、GO（http://www.geneontology.org）、Pfam（http://pfam.xfam.org/）等數(shù)據(jù)庫進(jìn)行比對，得到靶基因的注釋信息。使用GOseq R packages 進(jìn)行GO 富集分析，KOBAS 軟件用于KEGG 富集分析，富集分析中，按照q<0.05 進(jìn)行篩選（q為多重假設(shè)檢驗校正之后的P值，q值越小，富集顯著性越可靠）。并對差異表達(dá)miRs 靶基因進(jìn)行Pathway 富集分析，利用富集因子分析通路的富集程度。富集因子計算公式如下。

1.2.6實時熒光定量PCR（RT-qPCR）測定miRs表達(dá)量 選擇miR156b-5p[13]、miR166a[14]、miR167a[15]、miR168a[4,15]和miR319a-3p[16]5 個可跨界調(diào)控動物生理功能的植物miRs，以及5個在2個苜蓿品種中差異表達(dá)的miRs進(jìn)行實時熒光定量檢測。

根據(jù)測序得到的miRs 序列設(shè)計引物，試驗所需引物由Primer 5.0 軟件進(jìn)行設(shè)計（表1），由通用生物系統(tǒng)（安徽）有限公司合成。將經(jīng)過質(zhì)量檢驗的miRs 特異性反轉(zhuǎn)錄成cDNA，以5s 作為內(nèi)參進(jìn)行qPCR，反應(yīng)體系詳見表2；反應(yīng)程序如下：95 ℃3 min；95 ℃ 5 s，60 ℃ 30 s，39 個循環(huán)；95 ℃ 10 s。實時熒光定量PCR 結(jié)果用2?ΔΔCt[5]進(jìn)行處理，使用GraphPad Prism 7 軟件進(jìn)行作圖并進(jìn)行差異顯著性分析。

表1 引物信息Table 1 Primer information

表2 熒光定量反應(yīng)體系Table 2 Fluorescence quantitative reaction system

2 結(jié)果與分析

2.1 miRs種類數(shù)量和長度分析鑒定

分析miRs 測序組學(xué)數(shù)據(jù)發(fā)現(xiàn)，‘中苜1 號’和‘新鹽52 號’中分別檢測到656 和703 個miRs，其中新預(yù)測到的miRs 數(shù)量一致，均為223 個（表3）。比較2個苜蓿品種各miRs的TPM值發(fā)現(xiàn)，2個苜蓿品種中優(yōu)勢表達(dá)的前4個已知miRs均為miR5213-5p、miR159a、miR396a-5p和miR166e-3p。

表3 中苜1號和新鹽52號中優(yōu)勢表達(dá)的miRsTable 3 Advantage expression miRs in Zhongmu 1 and Xinyan 52

對2 個苜蓿品種中檢測到的miRs 長度進(jìn)行統(tǒng)計，發(fā)現(xiàn)已知miRs 的長度主要為21 nt（圖1A），新預(yù)測miRs長度主要以24 nt為主（圖1B）。并且2 種苜蓿新預(yù)測的miRs 中均有1 個miR 長度達(dá)到25 nt。

圖1 miRs的長度分布Fig. 1 Length distribution of miRs

2.2 2個苜蓿品種中差異表達(dá)miRs的對比分析

統(tǒng)計分析發(fā)現(xiàn)2個苜蓿品種中有21個miRs差異表達(dá)。其中已知miRs 10個，新預(yù)測miRs 11個,僅novel-miR54 在‘新鹽52 號’中表達(dá)量較‘中苜1號’顯著上調(diào)，其余20 個miRs 表達(dá)量均下調(diào)（圖2）。將在2 個苜蓿品種中差異表達(dá)且TPM 值豐度排在前5 的miRs 進(jìn)行序列比對，剔除重復(fù)序列，依次補(bǔ)位，最終得到novel-miR54、novel miR158、miR5754、miR156f（miR156e、miR156h-5p）、miR5743a（miR5743b）5個miRs（表4）。

圖2 苜蓿差異表達(dá)miRs熱圖Fig. 2 Heat-map of differentially expressed miRs in alfalfa

表4 差異表達(dá)miRs序列Table 4 miRs sequence of differentially expressed

2.3 差異表達(dá)miRs靶基因GO富集分析

在‘中苜1 號’檢測到的656 個miRs 中，有608 個miRs 共預(yù)測到7 536 個靶基因。在‘新鹽52 號’檢測到的703 個miRs 中，有654 個miRs 共預(yù)測到7 933個靶基因。在‘中苜1號’和‘新鹽52號’全部miRs預(yù)測到的靶基因中，有6 783個靶基因與GO數(shù)據(jù)庫比對成功。

在差異表達(dá)miRs 預(yù)測得到的623 個靶基因中，有503個靶基因與GO數(shù)據(jù)庫比對成功，獲得注釋信息。其中在‘新鹽52 號’中表達(dá)量上調(diào)的novel-miR54 預(yù)測得到116 個靶基因；在‘新鹽52號’表達(dá)量下調(diào)的miRs 中，10 個已知miRs 預(yù)測得到255個靶基因，10個新預(yù)測miRs預(yù)測得到252個靶基因。

差異表達(dá)miRs 在分子功能、細(xì)胞組分和生物過程3 個層面分別注釋到444、372 和349 個靶基因。注釋到生物過程層面的靶基因主要涉及的功能或代謝路徑為蛋白質(zhì)磷酸化、氧化還原過程和防御反應(yīng)。注釋到細(xì)胞組分層面的靶基因主要涉及的功能或代謝路徑為膜的組成成分、核和質(zhì)膜。注釋到分子功能層面的靶基因主要涉及的功能條目為ATP結(jié)合、DNA結(jié)合和鋅離子結(jié)合（表5）。差異表達(dá)miRs 靶基因與全部miRs 靶基因顯著富集的功能條目趨于一致，但在生物過程層面的防御反應(yīng)和分子功能層面的DNA 結(jié)合、ADP 結(jié)合等方面富集的靶基因有差異，這些差異通路的靶基因可為揭示2個苜蓿品種功能差異提供基礎(chǔ)素材。

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

2.4 差異表達(dá)miRs靶基因KEGG通路富集分析

差異表達(dá)的21 個miRs 共預(yù)測得到623 個靶基因，其中有215 個靶基因與KEGG 數(shù)據(jù)庫比對成功，獲得注釋，KEGG 的注釋結(jié)果按照KEGG 通路類型進(jìn)行分類，得到靶基因KEGG 通路類型分類。

差異表達(dá)miRs 的靶基因主要富集在5 個KEGG 通路類型，分別是細(xì)胞過程、環(huán)境信息處理、遺傳信息處理、新陳代謝和有機(jī)系統(tǒng)（表6）。對差異表達(dá)miRs的靶基因進(jìn)行通路富集分析，結(jié)果（圖3）發(fā)現(xiàn)，在靶基因注釋到的全部KEGG通路中，內(nèi)吞作用（6.74%）、RNA 轉(zhuǎn)運（6.74%）、ABC 轉(zhuǎn)運蛋白（5.62%）和泛素介導(dǎo)的蛋白水解（5.62%）為最顯著富集通路。

圖3 差異表達(dá)miRs靶基因KEGG通路富集散點圖Fig. 3 Scatter diagram of KEGG pathway enrichment of differentially expressed miRs target gene

表6 差異表達(dá)miRs靶基因KEGG通路富集分析Table 6 Target gene KEGG pathway enrichment of differential expressed miRs

2.5 5 個高豐度差異表達(dá)miRs 的表達(dá)量的檢測驗證

對 novel-miR54、novel-miR158、miR5754、miR156f （miR156e、miR156h-5p） 、miR5743a（miR5743b）這5 個差異表達(dá)的miRs 進(jìn)行RTqPCR，結(jié)果（圖4）表明，5個miRs的表達(dá)量變化趨勢與RNA-seq 結(jié)果一致。就表達(dá)量而言，上述miRs 僅有novel-miR54 在‘新鹽52 號’中的表達(dá)量較‘中苜1 號’極顯著上調(diào)（P<0.01），上調(diào)倍數(shù)為7.09 倍；其余4 個miRs 在‘新鹽52 號’中表達(dá)量均表現(xiàn)為下調(diào)。在品種方面，‘中苜1號’中miR156f的相對表達(dá)量最高，較‘新鹽52 號’極顯著上調(diào)（P<0.01）；‘新鹽52 號’中miR5743a 的表達(dá)量最高。

圖4 差異表達(dá)miRs的相對表達(dá)量Fig. 4 Relative expression level of differentially expressed miRs

2.6 5個已知跨界功能miRs表達(dá)量的驗證

對已知有跨界調(diào)控作用的miR156b-5p、miR166a、miR167a、miR168c-3p 和miR319a-3p 進(jìn)行定量分析，結(jié)果（圖5）表明，同一品種，miR156b-5p、miR167a、miR168c-3p 和miR319a-3p的表達(dá)量均低于miR166a，并且RT-qPCR 結(jié)果與RNA-Seq 結(jié)果具有相同的變化趨勢，說明RNASeq 結(jié)果可信度高。比較不同品種之間的表達(dá)量，僅miR166a 在2 個苜蓿品種中差異顯著，其在‘新鹽52號’中的表達(dá)量極顯著上調(diào)（P<0.01）。

3 討論

本研究通過miRs 組學(xué)分析，篩選出在2 個不同苜蓿品種中表達(dá)量高的miRs，為后續(xù)選定可以跨界調(diào)控奶牛性狀的苜蓿miRs提供了研究基礎(chǔ)。

3.1 關(guān)于2 種苜蓿miRs 表達(dá)譜、差異miRs 及功能代謝通路分析

測序結(jié)果表明，‘中苜1 號’和‘新鹽52 號’苜蓿中分別檢測到656 和703 個miRs，其中已知miRs中TPM 值最高的前4個miRs種類一致，說明2 個苜蓿品種中的高豐度miRs 相似，由此推測在研究miRs 跨界功能時苜蓿品種可以不作為主要因素予以考慮。另外測序和定量分析發(fā)現(xiàn)的novel-miRs為研究苜蓿屬miRs提供了新素材。

2 個苜蓿品種中僅有21 個差異表達(dá)的miRs，這些可能是導(dǎo)致二者功能差異和‘新鹽52 號’苜蓿獨有特性的miRs。通過比較新品系miRs 庫與現(xiàn)有品種的異同，可為鑒別和預(yù)判新品系的利用價值奠定分子基礎(chǔ)，例如差異表達(dá)的miR156。研究表明，miR156過表達(dá)基因型在干旱脅迫期間有利于保持植株更高的氣孔導(dǎo)度；同時，miR156 可以靶向SPL13，而SPL13 表達(dá)水平的降低減少了苜蓿的水分流失[17]，這表明miR156 可顯著改善苜蓿植株的耐旱性。從苜蓿miRs 跨界潛力方面分析，2 個苜蓿品種飼喂奶牛時，品系間高豐度表達(dá)miRs 可能是導(dǎo)致奶牛體況或乳品質(zhì)差異的關(guān)鍵miRs，如‘中苜1 號’中的miR5213-5p 及‘新鹽52號’中的miR159a 和miR166e。其中miR166e 在2個苜蓿品種中表達(dá)豐度均處在前4位，并且其在2個苜蓿品種中的相對表達(dá)量差異顯著，因此，應(yīng)用‘中苜1 號’和‘新鹽52 號’飼喂奶牛時，miR166e可作為發(fā)掘飼喂苜蓿品種與奶牛體內(nèi)基因表達(dá)及生產(chǎn)性能相關(guān)的miRs。

2個苜蓿品種差異表達(dá)miRs靶基因均顯著富集到RNA 轉(zhuǎn)運和ABC 轉(zhuǎn)運蛋白通路。差異表達(dá)miRs具有轉(zhuǎn)運RNA的能力，如果其可以跨界發(fā)揮作用，則可以通過調(diào)節(jié)奶牛體內(nèi)RNA 的轉(zhuǎn)運影響奶牛miRs組織表達(dá)譜。ABC 轉(zhuǎn)運蛋白在人、植物和動物等生物體內(nèi)對脂類的轉(zhuǎn)運及代謝具有重要意義，可從外界轉(zhuǎn)入長鏈脂肪酸及參與高密度脂蛋白的合成[18?19]。不同物種間的ABC 轉(zhuǎn)運蛋白結(jié)構(gòu)具有相似性[20]，miRs 序列也具有保守性，因此推測苜蓿源miRs進(jìn)入奶牛體內(nèi)，也能夠通過靶向奶牛體內(nèi)ABC 轉(zhuǎn)運蛋白進(jìn)而調(diào)節(jié)奶牛體內(nèi)脂質(zhì)代謝。但是奶牛攝入苜蓿源miRs 的量是否能夠達(dá)到發(fā)揮作用的量，以及苜蓿miRs靶向的ABC轉(zhuǎn)運蛋白的調(diào)節(jié)對象是奶牛乳脂還是牛肉脂肪都還需要深入研究。

由于miRs 的TPM 值越高，則豐度越高，被奶牛攝入后留存的可能性越大，因此通過序列比對篩選出novel-miR54、novel-miR158、miR5754、miR156f 和miR5743a 共5 個高豐度差異表達(dá)的miRs，RT-qPCR 分析的表達(dá)量變化趨勢與RNAseq測序結(jié)果一致，證明RNA-seq可信度高。對南方型紫花苜蓿的測序結(jié)果表明，miR156家族占總miRs 數(shù)量的9%，具有較高的表達(dá)豐度[21]，與本研究結(jié)果一致。另外，本研究測序新發(fā)現(xiàn)的miRs中novel-miR54 表達(dá)水平最高，值得深入探索其靶基因和功能作用。

3.2 關(guān)于2 種苜蓿中已知跨界miRs 的表達(dá)量及初步篩選結(jié)果

Kammes 等[22]研究發(fā)現(xiàn)，飼喂苜蓿可以改善奶牛乳品質(zhì)，提高奶牛乳脂含量和牛奶中尿素氮含量；飼喂苜蓿會改變奶牛瘤胃、十二指腸、空腸、肝臟和乳腺組織中與乳效率相關(guān)miRs的表達(dá)，調(diào)節(jié)乳脂含量[11,23?24]。這證明采食苜蓿影響奶牛體內(nèi)miRs 的表達(dá)譜，到底是由于苜蓿營養(yǎng)導(dǎo)致的奶牛乳腺miRs變化還是苜蓿中miRs的作用，還需要進(jìn)一步將本研究測定及篩選出的miRs在奶牛體內(nèi)血液及組織中的表達(dá)量進(jìn)行檢測和分析才能明確。

miR156b-5p、miR166a、miR167a、miR168c-3p和miR319a-3p這5個具有跨界功能的植物miRs在測序和定量分析中均被檢測到，植物miRs序列具有保守性[25]，因此推測上述miRs 在苜蓿中也具有跨界的潛力。由于高表達(dá)量的miRs被奶牛攝入后在奶牛體內(nèi)具有重要作用的可能性更大，因此通過RNA-seq 和RT-qPCR 最終篩選出在2 個苜蓿品種中均高表達(dá)的miR166a，為進(jìn)一步研究苜蓿中miRs在奶牛體內(nèi)發(fā)揮調(diào)控作用提供研究對象。

植物miR166 在其他物種的跨界調(diào)控作用已有相關(guān)報道，王鵬俊[26]研究發(fā)現(xiàn)，給豬飼喂玉米后，在豬的血清中檢測到玉米源的miR166a-3p；賈凌[27]研究發(fā)現(xiàn)，給蠶飼喂桑葉后，在蠶的血淋巴和組織中發(fā)現(xiàn)了桑樹源的miR166b 和miR166c；Lukasik等[28]在健康志愿者的母乳中檢測到5種植物源的miRs，其中包含miR166a。ath-miR166a 在人類樣品和豬母乳外泌體中豐度均為最高[15]；人參水煎劑中的miR166 對小鼠的免疫基因有一定影響[14]。由此可見，苜蓿中miR166a 對奶牛的分子作用機(jī)制有重要的研究價值和應(yīng)用潛力，值得進(jìn)一步挖掘。