成年滩羊和小尾寒羊皮肤毛囊差异表达miRNA 的筛选

何玉龙，杨 瑞，庄仁杰，李 勇，周学章，吴月红*

（1.浙江理工大学生命科学与医药学院，浙江杭州 310018；2.宁夏大学西部特色生物资源保护与利用教育部重点实验室，宁夏银川 750021；3.浙江理工大学校医院，浙江杭州 310018）

被毛生长是一个复杂的生理过程，受环境、营养、代谢及基因调控等因素影响，其中基因调控是决定性因素[1]。从基因水平上揭示绵羊被毛生长特点对于改良其品质具有重要价值。MicroRNAs（miRNAs）在调控基因表达方面起重要作用。随着绵羊相关miRNA 的不断发现，其在皮肤毛囊中的表达也逐渐受到重视，可能在皮肤毛囊发育中发挥着重要作用[2-5]。

小尾寒羊具有生长速度快、产羔率高、性能遗传稳定、适应性强等特点，是我国北方地区的优势绵羊品种[6-7]。利用小尾寒羊生产裘皮的历史悠久，但随着市场需求的变化，小尾寒羊逐步向肉用方向选育[8]。在提高小尾寒羊肉产量的同时，加强裘皮性状品系的选育具有重要意义。但由于目前并未将小尾寒羊毛用性状作为主要育种目标，导致其毛用性状未有明显提高[7,9]。滩羊主要分布在宁夏及其毗邻地区，是中国唯一生产“二毛裘皮”的绵羊品种[10]。目前，关于滩羊[11]和小尾寒羊[12]皮肤毛囊组织中miRNA 的研究已有报道，但关于两者皮肤毛囊中miRNA 的比较研究未见报道。为了解滩羊与小尾寒羊在裘皮性能差异的分子机制，本研究通过高通量测序技术比较滩羊与小尾寒羊皮肤毛囊中miRNA 差异，对于从miRNA 水平上解析皮毛性状形成的分子机理具有重要意义。

1 材料与方法

1.1 样品采集 滩羊来源于宁夏回族自治区吴忠市红寺堡区天源良种羊繁育养殖有限公司，小尾寒羊来源于内蒙古察哈尔右翼中旗某养殖户，分别采集成年滩羊（命名为TY_1）和小尾寒羊（命名为XWHY_1）各3只，每组个体的样品混合组成样品池。样品采集时间为2015 年2—3 月。剪毛消毒后用手术剪刀剪取背部相同位置约3～5 cm2皮肤组织，去除皮下脂肪组织等，剪成约0.5～1.0 cm2小块，PBS 溶液清洗后投入RNA 保护剂RNAlater（Thermo Fisher Scientific，美国）中。冷藏条件下尽快带回实验室，-80℃冰箱保存。

1.2 总RNA 提 取、小RNA 文 库 构 建 及Solexa 测 序将样本用液氮充分研磨后，利用Trizol 试剂（Invitrogen，美国）按照产品说明书进行皮肤组织总RNA 的提取，得到的总RNA 分别采用Nanodrop2000（美国）、Qubit 2.0（美国）和Agilent 2100 分析仪（美国）检测其纯度、浓度和完整性。样品检测合格后，以总RNA为起始样品，使用small RNA Sample Pre Kit 试剂盒（Illumina ，美国）构建文库，首先在small RNA 5'端和3' 端连接接头，反转录合成cDNA，PCR 扩增后利用PAGE 胶电泳筛选目的片段，切胶回收得到的片段即为small RNA 文库。质控合格后，利用HiSeq2500 进行高通量测序，测序读长为single-end（SE）50 nt。

1.3 测序数据分析 按照标准的miRNA 测序分析流程进行测序数据统计、小RNA 分类注释、已知miRNA的鉴定、新发现miRNA 的预测、miRNA 表达量及差异表达分析、差异表达miRNA 靶基因预测及注释等。

1.4 实时荧光定量PCR 验证miRNA 测序结果 为验证测序结果的可靠性，对鉴定出的miRNA 进行实时荧光定量PCR 检测验证。以总RNA 为初始模板，随机选取4 个miRNAs，miRNAs 反转录引物为通用颈环序列（GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCA CTGGATACGAC）和目标miRNA 3' 末端6 个碱基的反向互补序列构成。使用Revert Aid First Strand cDNA Synthesis Kit 试剂盒（Thermo Fisher Scientific，美国）按照产品说明书进行cDNA 第一链合成，使用SYBR®Premix Ex Taq 试剂盒（TaKaRa，大连）进行qRT-PCR验证，以U6 为内参基因。miRNA 相对表达量采用2-△△Ct法进行计算。miRNA 荧光定量PCR 引物见表1。所有引物由生工生物工程（上海）有限公司合成。

表1 实时荧光定量PCR 引物

1.5 差异表达miRNA 的筛选、靶基因预测及参与信号通路分析 利用IDEG6 对上述鉴定出的miRNAs 进行表达差异分析，选择差异表达的miRNAs，分别使用miRanda 和RNAhybrid 软件进行靶基因预测。对差异表达miRNA 的靶基因进行基因本体数据库（Gene ontology，GO）富集分析和京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析。

2 结果与分析

2.1 小RNA 文库序列统计 如表2 所示，滩羊（TY_1）和小尾寒羊（XWHY_1）样品分别通过HiSeq 2500 高通量测序，产出的原始数据（Raw reads）数量大于16.35 million（M），对其进行去除低质量数据、5'接头、3'接头及污染等处理后，样品中过滤后序列所占的比例分别为85.68% 和87.38%，质量值≥30 的碱基所占比例（Q30）分别为95.37%和95.40%。

表2 测序数据统计

2.2 sRNA 分类注释 利用Bowtie 软件将上述过滤后数据分别与Silva、GtRNAdb、Rfam 和Repbase 数据库进行序列比对，过滤核糖体RNA（rRNA）、转运RNA（tRNA）、核内小RNA（snRNA）、核仁小RNA（snoRNA）等非编码RNA（non-coding RNA，ncRNA）以及重复序列，在TY_1 和XWHY_1 样品分别获得包含miRNA 的未注释的reads（Unannotated reads）10 852 746 和14 135 887，用于后续分析，结果见表3。

表3 sRNA 分类注释统计

2.3 与参考基因组比对 利用miRDeep2 软件将上述未注释Reads 与参考绵羊基因组（Oar_v3.1）进行序列比对分析，获取其在参考基因组上的位置信息。如表4 所示，TY_1 和XWHY_1 样品比对到绵羊基因组的比例分别为52.14% 和51.22%，说明50% 以上的Reads 能比对到绵羊基因组上。

表4 与参考基因组比对信息统计表

2.4 miRNA 的鉴定 对于上述比对到参考基因组的序列，利用miRDeep2 软件进行已知及新发现miRNA的鉴定。即将18～30 nt 的核苷酸序列比对到miRBase（V21.0）数据库中特定物种上，鉴定其已知的miRNA。对于未比对到参考基因组的序列，通过碱基数目延伸，进行miRNA 结构预测，获得新发现miRNA（novel miRNAs）。经过比对共检测到561 个miRNAs，其中包括138 个已知miRNAs 和423 个新发现miRNAs。

2.5 miRNA 的长度分析 对2 个样本中鉴定到的138个已知miRNA 和423 个新发现miRNA 进行核苷酸长度统计。如图1 所示，大部分miRNA 序列长度在20～24 nt，已知miRNA 和新发现miRNA 中所占比例分别为96.38%和95.98%。其中长度为22 nt 的序列所占比例最高，分别为42.75%和51.54%；其次为21 nt 和23 nt，其长度分布符合Dicer 酶切割的典型特征，说明测序结果质量可靠，可用于后续分析。

图1 miRNA 长度分布

2.6 miRNA 表达量分析对TY_1和XWHY_1样本中鉴定出的miRNA 进行表达量分析，并用TPM 算法对表达量进行归一化处理。TPM 归一化处理公式：TPM=Readcount×1000000/Mapped Reads。其中，Readcount表示比对到某一miRNA 的reads 数目；Mapped Reads表示比对到参考基因组上的reads 数目。部分表达量较高的已知miRNA 如表5 所示。oar-let-7 家族在TY_1和XWHY_1 样本中表达量均较高。

表5 部分高表达已知miRNA 表达量数据（前20 位）

2.7 miRNA qRT-PCR 验证 随机选取4 个miRNA，通过qRT-PCR 的方法验证测序结果的准确性，结果如图2 所示。通过比较qRT-PCR 与高通量测序的结果发现（表6），这4 个miRNA 在TY_1 和XWHY_1 中表达的变化趋势基本一致，证明高通量测序结果的可靠性。

图2 miRNA 的qRT-PCR 验证分析

表6 用于qRT-PCR 验证的miRNA 测序表达量数据

2.8 差异表达miRNA 分析 利用IDEG6 对上述鉴定出的miRNAs 进行表达差异筛选，获得TY_1 与XWHY_1样品间差异表达的miRNA，筛选标准为|log2（FC）|Co且FDRg2 为71，其中FC 为差异倍数（Fold Change）；FDR 为错误发现率（False Discovery Rate）。结果在样品TY_1 与XWHY_1 共筛选到63 个上调miRNAs 和16个下调miRNAs。进一步对在2 个样品中差异表达的79个miRNAs 进行分析，发现43 个为已知miRNA，36 个为新发现miRNA。最后对43 个差异表达的已知miRNA进 行 聚 类 分 析（ 图3），oar-miR-127、oar-miR-379-5p、oar-miR-10a、oar-miR-411a-5p 和oar-miR-493-5p等在TY_1 和XWHY_1 样品中表达量均较高，而oarmiR-376 家族（包括oar-miR-376b、c、d 和e）和oarmiR-3959-3p 在2 个样品中的表达量差异较显著（差异倍数大于4），提示其可能在不同品种绵羊皮肤毛囊发育中起重要的调控作用。

2.9 差异表达miRNA 靶基因预测及富集分析 根据筛选到的差异表达miRNA 与绵羊的基因序列信息，用miRanda 和RNAhybrid 软件进行靶基因预测。使用BLAST 软件将预测靶基因序列与GO 和KEGG 等数据库比对，获得靶基因的注释信息分别为3 886 和4 449个。最后对样品TY_1 和XWHY_1 间差异表达miRNA靶基因进行GO 富集分析和KEGG 中通路类型分类统计，GO 统计发现其主要参与代谢过程、催化活性、细胞进程和细胞组分等过程（图4）；KEGG 通路的结果表明，差异表达miRNA 的4 449 个靶基因富集到113个信号通路上，其中在嘌呤代谢、内吞作用和糖酵解/糖异生等信号通路上富集显著性较高。

3 讨 论

皮肤毛囊的生长发育受到多个信号分子及信号通路的协同调控，是一个复杂的分子调控过程[13]。从miRNA 的角度简析绵羊皮肤毛囊生长发育的机制，对于最终从分子水平上揭示绵羊皮肤毛囊形成的分子机理及其调控机制具有重要作用。

有研究证实，滩羊在1 月龄时可宰取理想的羔皮，而小尾寒羊从胎内到3 或4 月龄都可宰取理想的羔皮[14]。为探讨这种差异形成的分子机制，本研究运用高通量测序技术分析比较了TY_1 和XWHY_1 皮肤毛囊组织中的miRNA 表达差异。测序原始数据经过去低质量、去污染和去接头等一系列处理，在滩羊和小尾寒羊皮肤毛囊组织样品中分别得到16 350 847 个和17 891 155 个高质量序列读数。通过与miRBase 数据库（V21.0）比对，在2 个样品中共鉴定出561 个miRNA，其中138 个为已知miRNA，423 个为新发现miRNA。通过筛选TY_1 与XWHY_1 间差异表达miRNA 发现63 个为上调miRNAs，16 个为下调miRNAs。这些差异表达miRNA中43 个为已知miRNA，36 个为新发现miRNA。

图3 差异表达的已知绵羊miRNA 聚类分析

图4 差异表达miRNA 靶基因的基因本体（GO）富集分析

相对于妊娠45 d 的绒山羊，在妊娠55 d 和65 d 样本中oar-let-7 和oar-miR-200 家族的miRNA 表达量都显著上调，提示其可能是毛囊发育中关键的miRNA[15]。而在本研究中前20 个高表达的已知miRNA 中oar-let-7 家族有5 个，分别为oar-let-7f、oar-let-7a、oar-let-7g、oarlet-7i 和oar-let-7c；而oar-miR-200 家族有2 个，分别为oar-miR-200b 和oar-miR-200c。 另 外，oar-miR-200 家族在羊驼皮肤和毛囊发育中也可能起重要作用[16]。已有研究报道，oar-miR-30a-5p 和oar-miR-23a 等可通过调控沉默信息调节因子2 相关酶1（Silent mating type information regulation 2 homolog 1，SIRT 1）参与羊毛卷曲生长[17]，而在本研究中也检测到oar-miR- 30a-5p和oar-miR-23a 的高表达，提示其可能参与小尾寒羊和滩羊毛发卷曲生长。另外，本研究在滩羊和小尾寒羊皮肤中检测到高水平表达的oar-miR-125b 在羊驼背部皮肤中表达量也较高[16]。与小尾寒羊相比，在滩羊皮肤毛囊中高水平表达的oar-miR-127 和oar-miR-379-5p，在藏绵羊皮肤毛囊中也检测到[18]。而oar-miR-376 家族（包括oar-miR-376b、c、d 和e）和oar-miR-3959-3p在2 个样品中的表达量差异最显著（差异倍数大于4），提示其可能在不同品种绵羊皮肤毛囊发育中起重要的调控作用。GO 统计差异表达miRNA 的靶基因发现，其主要参与代谢过程、催化活性、细胞进程和细胞组分等过程；KEGG 通路的结果表明，差异表达miRNA 的靶基因主要富集到嘌呤代谢、内吞作用和糖酵解/糖异生等信号通路上。

4 结 论

本实验利用高通量测序技术在滩羊和小尾寒羊皮肤毛囊组织中筛选得到了79 个差异表达的miRNAs，并对其靶基因进行预测，推测其可能参与了绵羊皮肤毛囊发育的调控，但其具体机制还需要进一步验证。