驴乳外泌体中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家族及通路，为开发以驴乳为原料的婴幼儿乳粉及功能性乳制品提供一定参考。