驢乳外泌體中miRNA的測序與分析

柴玉霞，王新宇，岳喜慶，楊 梅,*

（1.沈陽農業大學食品學院，遼寧 沈陽 110866；2.沈陽農業大學畜牧獸醫學院，遼寧 沈陽 110866）

驢乳是營養豐富、生物學價值較高并具有一定保健功能的天然食品。驢乳中除脂肪含量顯著低于人乳外，其他化學成分都與人乳相近。驢乳中蛋白質含量較低，且主要是乳清蛋白，易于消化吸收，適合于嬰幼兒腎功能發育不完善以及腸道對滲透壓耐受能力較弱的特點；乳糖含量高，可促進嬰幼兒腸道有益微生物增殖以及鈣吸收。驢乳含有多種活性功能性成分，具有抑菌、提高機體免疫力、調節腸道菌群、抗氧化、抗疲勞、抑瘤等作用，可作為膳食補充劑，提供傳統嬰兒配方乳粉所缺乏的生物活性和抗感染成分。此外，驢乳具有高適口性和低致敏性，適用于嬰幼兒配方乳品的研發。

外泌體廣泛存在于哺乳動物體液中，由于其獨特的形成過程，可攜帶多種蛋白質、核酸和脂質等，由細胞主動向外分泌到達受體細胞，從而參與細胞間的信息傳遞和物質交流。乳外泌體中微小RNA（miRNA）受外泌體膜結構的保護作用，可穩定進入嬰兒腸道，免受核酸酶消化，作為信號分子傳遞給其他細胞，調節受體細胞的生物學活性。此前，研究者們采用不同的測序技術，分別對人乳、羊乳、牛乳、豬乳等哺乳動物乳的外泌體miRNA進行了探索。目前對于驢乳的研究大多集中于其化學成分、理化性質及加工特性的研究，對驢乳外泌體miRNA的種類及功能性卻少見報道。

本研究構建驢乳及人乳外泌體的小RNA（sRNA）文庫，應用Illumina測序技術對驢乳及人乳外泌體中的miRNA進行測序分析，并對鑒定出的miRNA進行生物信息學分析，分析其功能及差異性。通過對驢乳及人乳外泌體miRNA對比分析，有助于加深對驢乳外泌體miRNA的種類及功能性的認知，為開發以驢乳為原料的更適合嬰幼兒需求的配方乳粉提供理論依據。

1 材料與方法

1.1 材料與試劑

驢乳采集于大連圣鴻道荙驢業科技發展有限公司養殖場。所選取的德州驢（=30）標準為：年齡2～3 歲、首次分娩、體質量250～350 kg、飲食與飼養條件一致。人乳樣本由30 名健康母親捐贈，捐贈者均為25～30 歲、自然分娩。該研究方案獲得沈陽農業大學和相關機構指南的批準。采集的驢乳與人乳樣品置于干冰中運輸并保存于超低溫冰箱（－80 ℃）。實驗前為排除個體差異，將30 份驢乳樣品隨機分成3 組，每組10 份樣品進行充分混合，人乳樣品采用相同方法處理。

TRIzol試劑、6% Novex TBE凝膠 美國Invitrogen公司；TruSeqsRNA樣品制備試劑盒、HiSeq X試劑盒美國Illumina公司。

1.2 儀器與設備

wonbio-96型高通量研磨儀 上海萬柏生物科技有限公司；5424R型離心機 美國Eppendorf公司；2100型生物分析儀 美國安捷倫科技有限公司；JY600C型通用電泳儀 北京君意東方電泳設備有限公司；ND-2000型紫外-可見光分光光度計 上海基因公司；FR-1000型凝膠成像分析系統 上海復日科技有限公司。

1.3 方法

1.3.1 乳外泌體的提取

取5 mL乳樣離心（2 500h、4 ℃、10 min），重復兩次以清除細胞殘渣和脂肪顆粒，吸取上清液，離心（12 000h、4 ℃、30 min）以清除所有細胞殘渣。吸取上清液，離心（120 000h、4 ℃、4 h），用磷酸鹽緩沖液溶解沉淀物，再次離心（100 000h、4 ℃、1 h），收集底部沉淀。用磷酸鹽緩沖液溶解并稀釋提取的外泌體沉淀后用透射電鏡進行觀察。

1.3.2 總RNA的提取、文庫構建以及測序

取200 μL乳樣在研磨儀中進行液氮研磨（10 s、50 Hz），加入1 000 μL預冷TRIzol試劑，隨后將樣品溶液離心（13 000h、4 ℃、5 min），加入200 μL氯仿，充分振蕩混勻，離心（13 000h、4 ℃、15 min），吸取上清液至新的離心管中，加入1.2 倍預冷無水乙醇，在－20 ℃沉淀2～4 h后再次離心（13 000h、4 ℃、5 min），棄除上清液。加入1 000 μL體積分數75%乙醇溶液清洗沉淀，再次離心（13 000h、4 ℃、5 min），收集沉淀。

用1%瓊脂糖凝膠評估RNA降解和污染情況。使用紫外-可見光分光光度計測量RNA濃度。最后使用生物分析儀評估RNA完整性。僅用高質量RNA（OD/OD=1.8～2.2，OD/OD≥2.0，RNA完整值≥8，28S/18S≥1.0，＞3 μg）構建測序文庫。

每個樣品取3 μg RNA構建文庫。用TruSeq sRNA樣本制備試劑盒連接3’接頭，并將隨機引物與過量的3’端隨機接頭雜交。連接5’末端，加入M-MuLV逆轉錄酶合成第1鏈。隨后進行聚合酶鏈式反應擴增，使用8%聚丙烯酰胺凝膠純化擴增產物，篩選其中長度140～160 bp范圍內的片段進行回收和溶解以便于后續實驗。對構建的RNA文庫質量進行評估，聚類生成后，對驢乳及人乳外泌體miRNA的樣品溶液進行測序分析。

測序完成后，統計所有測序片段的質量波動結果以及堿基分布情況，并對每個樣本的堿基質量、堿基錯誤率進行分析。然后使用Fastp軟件對原始數據進行篩選，去除驢乳和人乳外泌體miRNA樣品溶液中所有質量低的堿基序列，然后對篩選后制備的sRNA序列的數量及種類進行統計分析，分析其堿基分布情況，獲得miRNA的偏好性信息，對不同長度sRNA的首位堿基及不同位置堿基的偏好性進行統計分析。

1.3.3 sRNAs的定位及序列特征分析

使用Rfam數據庫（http://rfam.xfam.org/）進行注釋，去除其中非miRNA序列，同時對比對成功的非miRNA序列進行種類和數目統計。去除rRNA、小核RNA（snRNA）、轉運RNA（tRNA）、核酶等非已知miRNA序列后，利用bowtie（http://bowtie-bio.sourceforge.net/index.shtml）對照參考基因組數據。通過比對mirBaSe（http://www.mirbase.org/），完全匹配上的序列用于已知miRNA表達譜的計數與分析。

1.3.4 已知miRNA鑒定和新miRNA識別

以已有的基因組序列和數據庫中miRNA序列為參考，對注釋后的驢乳及人乳樣品溶液分別進行比較分析，篩選出已知的miRNA。由于驢在miRBase中沒有對應的miRNA信息，根據miRNA在近緣物種間的保守性，選擇與目標物種相近、且在miRBase數據庫中有miRNA信息的近緣物種馬，作為驢乳外泌體中已知miRNA信息的參考。截取不能比對上Rfam和miRBase的sRNA周圍序列，預測其二級結構，并根據二級結構中能量值和Dicer酶切位點等信息鑒定出新的miRNA。

1.3.5 統計與生物信息學分析

對驢乳及人乳外泌體中檢測到的所有miRNA進行靶基因預測，并分別對其預測出的靶基因進行功能描述和通路分析。以基因本體論（Gene Ontology，GO）數據庫為參考，將對比上的miRNA的靶基因進行統計分析，并預測其功能性。以京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）數據庫為參考，采用富集分析方法，對其參與的通路進行統計分析。以Fisher檢驗統計富集顯著程度。使用Holm、Sidak、Bonferroni等方法降低假陽性率。

2 結果與分析

2.1 驢乳和人乳外泌體的表征

如圖1所示，驢乳和人乳樣本中均可見大小均勻、形態一致的圓形囊泡，直徑約為100 nm，且呈經典茶托樣結構，與此前對外泌體形態結構的研究結果相符。

圖1 透射電子顯微鏡下驢乳（A）和人乳外泌體（B）的形態Fig.1 Morphology of exosomes in donkey (A) and human milk (B)under transmission electron microscope

2.2 測序數據質量及序列長度分布特征

在驢乳和人乳外泌體中分別鑒定到10 038 218、15 285 595 條原始序列（表1），質量控制后，驢乳中得到純凈序列4 123 989 條，占原始序列的41.08%。人乳中純凈序列共4 953 715 條，占比32.41%（表2）。后續分析時，長度＜18 nt或＞32 nt、含N序列等低質量序列均已去除。

表1 驢乳及人乳外泌體中原始測序數據統計結果Table 1 Statistics of raw sequencing data of exosomes in donkey milk and human milk

表2 驢乳及人乳外泌體中質量控制后序列信息Table 2 Sequence data of exosomes in donkey milk and human milk after quality control

如圖2所示，驢乳與人乳的純凈序列長度主要分布在18～22 nt，驢乳樣本中sRNA序列長度主要集中在20 nt，占比26.07%，與驢乳相比，人乳中sRNA序列長度分布更均勻。

圖2 驢乳及人乳外泌體中sRNA序列長度分布Fig.2 sRNA sequence length distribution of exosomes in donkey milk and human milk

通過堿基偏好性統計可以了解miRNA的加工和作用機制，同時也能間接反映樣品質量、建庫、測序實驗的準確性。對測得的所有純凈sRNA序列的各位點堿基分布情況進行分析統計，從圖3可以看出，驢乳與人乳外泌體中，不同長度sRNA首位堿基均對G有顯著偏好性，驢乳外泌體sRNA的10、12號位點對堿基U具有顯著偏好性。與驢乳相比，人乳外泌體sRNA各位點堿基占比分布較均勻，其中1、2、3、28號位點對堿基G有顯著偏好性，在25、26、27號位點對堿基C有顯著偏好性。

圖3 驢乳及人乳外泌體中不同長度和不同位點的堿基偏好性分析Fig.3 Base preference analysis of different lengths and sites for exosomes in donkey milk and human milk

2.3 sRNA序列的分類注釋

如表3所示，驢乳外泌體中sRNA序列共有4 123 989 條，其中與數據庫比對成功的有3 161 039 條，占比76.65%。人乳外泌體中sRNA序列共有4 953 715 條，與數據庫比對成功的有3 402 874 條，占比68.69%。驢乳和人乳外泌體sRNA序列中均含有rRNA、miRNA、核仁小RNA（snoRNAs）、snRNA、tRNA等，其中rRNA數量最多，占比分別為73.14%（驢乳）和64.16%（人乳）。

表3 人乳和驢乳外泌體中sRNA序列與Rfam數據庫比對結果Table 3 Comparison of sRNA sequences of exosomes in donkey milk and human milk with the Rfam database

將rRNA、snRNA、tRNA等非已知miRNA序列去除后，將剩余的sRNA序列與參考基因組進行比對。驢乳外泌體中合并相同序列后繪制到染色體的總序列數為382 363 條，比對成功的序列有166 294 條，占比43.49%。人乳外泌體中合并相同序列后繪制到染色體的總序列數為459 184 條，比對成功的序列有215 063 條，占比46.84%。

2.4 已知miRNA的鑒定及新miRNA預測

由表4可知，在人乳和驢乳外泌體miRNA中，分別鑒定到36 種和16 種成熟miRNA，其中家族成員、和在人乳和驢乳外泌體中均被檢測到。家族因其廣泛的功能特性而備受關注。研究發現，可以結合靶mRNA的互補位點，影響蛋白翻譯過程或改變mRNA的穩定性，進而影響基因表達。Chen Ting等通過Solexa測序方法在豬乳外泌體中共檢測到9 種家族成員，這些miRNA在仔豬消化道產生免疫球蛋白中起關鍵作用。Izumi等在牛乳外泌體中檢測到、、、、和，其中表達量最高。本研究在人乳和驢乳外泌體中都鑒定到家族成員，且顯著表達，然而，在人乳與驢乳外泌體中家族成員所具有的生理功能仍有待進一步研究。

表4 驢乳及人乳外泌體中已知miRNA統計結果Table 4 Statistics of miRNA in donkey milk and human milk exosomes

將不能與Rfam和miRBase成功比對的sRNA與參考基因組比對，預測新miRNA，驢乳中最終預測得到196 個新miRNA序列，總表達量為181 274，表達量超過1 000的miRNA序列共27 個。人乳中預測到256 個新miRNA序列，總表達量為1 034 642.56，表達量超過1 000的miRNA序列共36 個。這些預測到的新miRNA能夠極大豐富物種miRNA序列信息。

續表4

2.5 靶基因預測及功能分析

2.5.1 靶基因GO功能性分析

對已知miRNA的靶基因進行GO二級分類注釋和統計分析，如圖4～6所示，驢乳和人乳外泌體miRNA參與的生物過程和細胞組成相似，都參與細胞代謝過程、單一生物體過程、刺激反應、生物過程調節、細胞及細胞器的構成等條目，但不同乳中注釋到同一GO條目的靶基因不同。此外，人乳與驢乳外泌體miRNA參與的分子功能具有顯著差異。驢乳外泌體已知miRNA中，共19 條靶基因被注釋到分子功能類的3 個條目，其中17 條靶基因與結合作用相關，其余2 條靶基因分別參與催化活性和分子功能調節器。人乳外泌體已知的成熟miRNA中，20 366 條靶基因被注釋到分子功能類的18 個條目，其中9 998 條靶基因與結合作用有關，其余靶基因分別參與催化活性、信號傳感器活性、核酸轉錄因子活性、轉運蛋白活性、分子傳感器活動、分子功能調節器、轉錄因子活性，蛋白質結合、結構分子活性、電子載體活性等條目。

圖4 驢乳及人乳外泌體miRNA在生物過程層面的GO注釋Fig.4 GO annotation of miRNAs in exosomes of donkey milk and human milk at the biological process level

圖5 驢乳及人乳外泌體miRNA在分子功能層面的GO注釋Fig.5 GO annotation of miRNAs in exosomes of donkey milk and human milk at the molecular function level

圖6 驢乳及人乳外泌體miRNA在細胞組成層面的GO注釋Fig.6 GO annotation of miRNAs in exosomes of donkey milk and human milk at the cell component level

2.5.2 靶基因通路富集分析

驢乳中對預測到的miRNA靶基因進行KEGG通路富集分析，已知miRNA中共有1 204 個靶基因被注釋到272 條通路中，其中7 條獲得顯著富集（＜0.05），主要為醛固酮的合成與分泌（KEGG編號ko04925）、FcγR介導的吞噬作用（ko04666）、RNA聚合酶（ko03020）、Hedgehog信號通路（ko04340）、光轉導（ko04745）、HIF-1信號通路（ko04066）、肌萎縮側索硬化（ko05014），最顯著富集的通路為醛固酮的合成與分泌。醛固酮是一種在腎上腺皮質（腎小球帶）外層合成和分泌的類固醇激素，醛固酮通過吸收鈉和水在調節系統性血壓中起重要作用。新miRNA序列共有3 139 個靶基因被注釋到300 條通路中，其中有9 條獲得顯著富集（＜0.05），主要參與Hedgehog信號通路（ko04340）、線粒體調節（ko04139）、肌萎縮側索硬化（ko05014）、MAPK信號通路（ko04011）、GnRH信號通路（ko04912）、mRNA監測途徑（ko03015）、Wnt信號通路（ko04310）、神經營養蛋白信號傳導通路（ko04722）、醛固酮的合成與分泌通路（ko04925）。驢乳中新miRNA序列的靶基因富集最多的通路是癌癥相關通路（ko05200）。

人乳外泌體已知miRNA中共19 722 個靶基因被富集到324 條通路，其中103 條通路獲得顯著富集（＜0.05），其中神經活性配體-受體相互作用（ko04080）、細胞因子-細胞因子受體相互作用（ko04060）、胰島素抵抗（ko04931）、cAMP信號通路（ko04024）、癌癥中的miRNA（ko05206）、逆行內源性大麻素信號（ko04723）、膽堿能突觸（ko04725）、谷氨酸能突觸（ko04724）、AGE-RAGE信號通路（ko04933）、趨化因子信號通路（ko04062）獲得極顯著富集（＜0.001）。新miRNA序列共有23 376 個靶基因被注釋到325 條通路中，其中有90 條獲得顯著富集（＜0.05），其中281 個靶基因被富集到神經活性配體-受體相互作用（ko04080）。

磷脂酶D是負責產生脂質第二信使磷脂酸的必需酶，該酶參與基本的細胞過程、括膜運輸、肌動蛋白細胞骨架重塑、細胞增殖和細胞存活。甘油磷脂在機體中發揮著重要作用，除參與生物膜、膽汁的構成外，還能夠識別蛋白質，調控細胞膜的信號傳遞過程。細胞表面蛋白可以通過糖基磷脂酰肌醇的糖脂結構附著在細胞膜上。研究表明，這些與脂質代謝相關的信號通路在牛乳、羊乳和豬乳中均有相關表達。與該結果類似，本實驗在驢乳和人乳外泌體miRNA的靶細胞中均富集到了磷脂酶D信號通路、甘油磷脂代謝和糖基磷脂酰肌醇錨生物合成等通路（表5）。驢乳外泌體miRNA中脂質代謝相關通路的發現，豐富了驢乳外泌體功能性的理論研究，能夠進一步對以驢乳為原料的嬰幼兒乳制品的研發提供理論依據。

表5 驢乳及人乳外泌體miRNA中脂質代謝相關通路Table 5 Numbers of target genes of miRNAs in donkey milk and human milk exosomes involved in lipid metabolism-related pathways

3 結 論

利用Illumina測序技術和生物學信息分析技術對驢乳和人乳外泌體中sRNA進行測序分析，并對其所含miRNA的種類、數量及其功能性進行鑒定與分析。驢乳外泌體中sRNA序列長度主要集中在20 nt，不同長度sRNA首位堿基均對G有顯著偏好性。在驢乳外泌體中鑒定到16 種已知miRNA，預測到196 個新miRNA序列，已知miRNA中、、eca-miR-429和表達量較高。生物信息學分析表明驢乳外泌體miRNA的功能性與人乳外泌體miRNA相似，主要參與細胞代謝過程、結合作用、細胞及細胞器的構成，并參與了磷脂酶D代謝、甘油磷脂代謝和糖基磷脂酰肌醇錨生物合成等脂質代謝相關通路。本研究對驢乳外泌體中miRNA的種類和功能進行初步探索，發現一些共有的miRNA家族及通路，為開發以驢乳為原料的嬰幼兒乳粉及功能性乳制品提供一定參考。