中国美利奴羊羊毛弯曲相关microRNA的筛选

2021-04-01于丽娟张艳花拉扎特艾尼瓦尔徐新明玛尔孜娅亚森

新疆农业科学 2021年3期

于丽娟, 张艳花，拉扎特·艾尼瓦尔，徐新明，玛尔孜娅·亚森，狄江

(新疆畜牧科学院畜牧研究所，乌鲁木齐 830011)

0 引言

【研究意义】弯曲是羊毛最主要的物理性质之一，与毛纺工艺有着密切的关系。研究表明，降低羊毛纤维弯曲度能减小纤维抗弯曲刚性和纱线厚度并进而增加织品柔软度，其对纺织品柔软度的影响与纤维直径有同等重要的作用[1]。了解羊毛弯曲形成的分子机理，将有助于利用分子育种培育符合市场需求羊毛弯曲类型细毛羊新群体。【前人研究进展】羊毛弯曲与毛囊中细胞的不对称分裂、毛纤维中正副皮质的双边排列、不同角蛋白的不对称组成及角质化过程等相关[2]。对小鼠、人等哺乳动物毛发纤维弯曲的形成机制方面的研究取得了很大进展，Wnt、EDA/EDAR、EGFR、BMP、NF-κB、Shh等经典信号通路与毛发弯曲有着直接或间接的联系[3]。【本研究切入点】MicroRNA(miRNA)为一类大小在18～24个核苷酸的非编码单链RNA，通过与靶mRNA结合位点的特异性结合，抑制靶mRNA翻译或降解靶mRNA来发挥其重要的生理调节作用[4]。研究表明，miRNA几乎参与了机体所有的代谢调节活动，包括肌肉生长发育、细胞周期变化、细胞增殖与凋亡、脂肪代谢、器官形成以及机体的免疫应答等反应[5]。随着miRNA对皮肤和毛囊发育形成的研究的不断深入，miRNAs是毛囊发育过程中重要的转录后调控因子，参与毛囊的发生和周期性生长[6-9]。目前，miRNA对羊毛弯曲性状形成的相关调控机制仍不清楚，需要研究miRNAs在羊毛弯曲形成过程中发挥的作用。【拟解决的关键问题】通过采用高通量测序技术对羊毛弯曲频率有极端差异表型值细毛羊个体皮肤的miRNA组进行测序并比较分析，筛选出与细毛羊羊毛弯曲性状相关的差异表达miRNA，利用miRNAda与TargetScan软件对差异表达的miRNA进行靶基因预测，并对靶基因进行简单的Gene Ontology与KEGG Pathway富集分析，分析差异 miRNA的预测靶基因参与的生物学功能及调控途径，为进一步解析羊毛弯曲形成分子机理、以及借助分子育种手段培育符合市场需求羊毛弯曲类型细毛羊新群体提供理论指导。

1 材料与方法

1.1 材料

绵羊为新疆拜城种羊场的中国美利奴羊，羊群在基本一致的条件下饲养。采集样品主要以羊毛弯曲频率性状为依据。在羊毛品质鉴定中，羊毛弯曲频率一般分为大、中、小3个弯曲类型，其中大弯为1 cm 羊毛纤维长度内有小于3 个弯曲、中弯为1 cm 羊毛纤维长度内有3 ～ 5 个弯曲、小弯为1 cm 羊毛纤维长度内大于5 个弯曲。个体鉴定时，选取大弯周岁母羊个体与小弯周岁母羊个体各10只。体侧肩胛部皮肤剪毛消毒后，用皮肤采样器分别采集其肩胛部皮肤组织1 cm2左右，投入装有RNA保存液(RNAlater)的冻存管中。冷藏条件下带回实验室，-80℃冰箱保存备用。对其毛样进行细度与羊毛弯曲频率测定，选取羊毛弯曲频率有极端差异的大弯、小弯个体各4只的皮肤样品用于后续实验。

1.2 方法

1.2.1 RNA提取、小RNA文库构建及Solexa测序

按照Trizol法提取细毛羊皮肤组织总RNA，采用Qubit 2.0 (美国)试剂盒进行RNA的精确定量；利用NanoPhotometer®分光光度计进行纯度检测，同时利用Bioanalyzer 2100 系统的 RNA Nano 6000 Assay Kit进行质量完整性的检测。样品检测合格后，使用Small RNA Sample Pre Kit构建文库，利用Small RNA的3’及5’端特殊结构(5’端有完整的磷酸基团，3’端有羟基)，以总RNA为起始样品，直接将Small RNA两端加上接头，然后反转录合成cDNA。随后经过PCR扩增，PAGE胶电泳分离目标DNA片段，切胶回收得到cDNA文库。利用安捷伦2100高灵敏度DNA芯片检测文库质量。库检合格后，使用 Illumina Hiseq 2500 上机测序。

1.2.2 测序数据

对测序得到的原始数据(raw data)进行质控和过滤，包括去除低质量、有5’接头污染、没有3’接头序列、没有插入片段及包含poly A/T/G/C的reads。进行质控后得到的数据为clean reads，用于后续分析。采用bowtie软件与羊的参考基因组(Oar_v4.0)进行比对，分析small RNA在参考序列上的分布情况。且利用miRNA所在的位置对已知(保守)miRNA进行鉴定，而新miRNA则通过miREvo和mirdeep2这2个软件进行预测。利用DESeq软件筛选差异表达的miRNA，利用TPM分别计算样本huACF，huBDE表达量取log2，满足Pvalue≤0.05筛选2组之间的差异miRNA。利用miRNAda、TargetScan靶基因预测软件对miRNA进行靶基因预测。根据miRNA与其靶基因的对应关系，对差异miRNA的预测靶基因进行简单的Gene Ontology与KEGG Pathway富集分析。

1.2.3 实时荧光定量PCR验证

随机选择4个差异表达的miRNA，检测其在2组样品的表达量。以U6作为内参基因，分别根据4个miRNAs的成熟序列设计荧光定量PCR的上游引物，下游引物由试剂盒提供，引物由上海生工生物科技有限公司合成。表1

表1 qRT-PCR引物Table 1 The list of primers of qRT-PCR

采用poly(A)加尾法进行cDNA的反转录，具体操作步骤参照康为世纪miRNA cDNA Synthesis Kit进行。反应体系为20 L：2×miRcute增强型miRNA定量预混试剂(含SYBR和ROX)10 μL，反向引物(10 μmol/L )0．4 μL，正向引物(10 μL )0.4 μL ，cDNA 2 μL，无核酶ddH2O。补至20 μL。PCR反应条件：95℃ 预变性15 min；94℃ 变性20 s，60℃ 退火延伸34 s，共40个循环。每个样品每个miRNA均进行3次技术重复，结果采用2-ΔΔCt法计算miRNA的相对表达量，与测序结果对比，检验测序结果的可靠性。

2 结果与分析

2.1 总RNA质量检测

研究表明，提取的总RNA进行质量检测，其OD260/280值均在1.9～2.0，RIN值均高于8.0，总RNA纯度及完整性均达到构建文库要求。表2

表2 总RNA质量检测Table 2 The quality test results of total RNA

2.2 中国美利奴羊皮肤组织测序

研究表明，经Solexa测序和生物信息学分析,去除接头序列以及低质量序列，8个文库共获得81 820 492条总序列，经过滤后共获得79 455 352条纯净序列。细毛羊皮肤组织小RNA序列长度大部分位于21～24 nt，其中长度为22 nt的序列所占比列最大。表3，图1

表3 数据过滤信息Table 3 The information of raw data filtering

图1 小RNA的长度分布Fig.1 Length distribution of sRNAs

2.3 小RNA的分类注释

利用Bowtie 软件将得到的纯净序列在数据库中比对后进行分类注释，得到各种类型sRNA,包含有已知miRNA、新miRNA、rRN A、snRNA、snoRNA、tRNA和其他RNA。其中，各样品中已知的miRNA所占的比例最大。表4

表4 各样本sRNA分类注释Table 4 sRNA classification annotation results for each sample

2.4 小RNA序列基因组比对

针对每个样本预处理后所有序列和去重复后的唯一序列，分别采用Bowtie软件与羊的参考基因组(Oar_v4.0)进行比对分析。获取其在参考基因组上的位置信息。比对结果表明81%以上的片段能比对到绵羊基因组上。61%以上的种类比对到绵羊基因组上。各样品间的比对率存在一定的差异说明样品间小RNA存在一定差异。表5

表5 与参考序列比对Table 5 Alignment results of the reference sequence

2.5 已知miRNA分析和新miRNA预测

将mapped到参考序列上的reads，与miRBase中指定范围序列进行比对，得到各样品匹配上的sRNA的详细情况，包括匹配上的已知miRNA的二级结构，各样本中miRNA的序列、长度、出现的次数等信息。比对上的已知成熟miRNA143条，已知 miRNA前体有106条。

整合miREvo(Wen et al.，2012)和mirdeep2(Friedlander et al.，2011)这些miRNA预测软件来进行新miRNA的分析。在小RNA文库中总共预测得到306条miRNA前体，成熟miRNA也有282条。

2.6 miRNA差异表达

利用DESeq软件筛选差异表达的miRNA，利用TPM分别计算样本huACF，huBDE表达量取log2，满足Pvalue≤0.05筛选2组之间的差异miRNA。共鉴定出差异表达miRNA25个，其中上调表达8个，下调表达17个。

2.7 差异miRNA靶基因预测与功能富集

利用miRNAda与TargetScan靶基因预测软件对差异表达的miRNA进行靶基因预测。预测结果显示，大弯曲与小弯曲组的25个差异表达的miRNA共预测到了15 761个无重复的靶基因。根据miRNA与其靶基因的对应关系，对差异miRNA的预测靶基因进行简单的Gene Ontology与KEGG Pathway富集分析。GO分析结果表明，其主要参与跨膜信号受体活动、信号转导活动、分子转导活动、膜的组成、细胞交流、生物过程、细胞过程的调节等过程。15 761个靶基因注释到252个信号通路。其中，107个靶基因参与了Wnt信号通路，93个靶基因参与了AMPK信号通路，34个靶基因参与了Notch信号通路，161个靶基因参与了MAPK信号通路，68个靶基因参与了NF-kappa B信号通路。图2

图2 差异表达miRNA靶基因GO

2.8 miRNA qRT-PCR 验证

随机选择4个差异表达miRNA(novel_503、oar-miR-133、oar-let-7g、novel_9)进行定量验证，结果表明：novel_503、oar-miR-133、oar-let-7g、novel_9在大弯曲组中上调表达，而在小弯曲组中下调表达，与高通量测序结果一致，证明高通量测序结果的可靠性。图3，表6

图3 miRNA 的qRT-PCR 验证Fig.3 qRT-PCR verification of miRNA

表6 用于qRT-PCR 验证的miRNA测序表达量数据Table 6 Sequencing data of miRNA verified by qRT-PCR

3 讨论

关于中国美利奴羊羊毛经济性状方面的研究主要集中在对羊毛细度、长度等。在中国美利奴羊羊毛弯曲性状方面的研究仅见狄江等[10]通过全基因组关联分析检测与细羊毛弯曲频率性状显著关联的SNPs 。

高通量测序技术已被广泛应用于动植物转录组测序、小RNA测序等方面，在中国美利奴羊羊毛弯曲性状形成方面的研究还未见报道。研究运用高通量测序技术对羊毛弯曲频率有极端差异表型值(大弯与小弯)细毛羊个体皮肤的miRNA组进行测序并比较，共鉴定出425个miRNAs，其中包括143个已知miRNAs和282个新发现的miRNAs。鉴定出的miRNAs进行表达量差异分析发现，在大弯曲组与小弯曲组共筛选到8个上调和17个下调的miRNAs。其中，差异表达miRNA中的oar-miR-133、oar-miR-379-5p的表达量趋势与Liu的研究结果一致[11]。并且，高水平表达的oar-miR-379-5p在滩羊与藏绵羊皮肤毛囊中也检测到[12,13]。差异表达miRNA中的oar-miR-21可能通过调控CNKSR2、KLF3和TNPO1的表达来影响毛囊发育[14]。差异表达miRNA中的oar-let-7g属于let-7家族，该miRNA的靶基因与皮肤组织中毛囊的生长和毛发的质量相关[15],而且研究表明，let-7a、let-7b、let-7c、let-7i等与滩羊被毛卷曲、湖羊羔皮花纹及生长相关[16,17]。推测let-7家族可能参与羊毛弯曲的形成。对于其他差异表达的miRNA,目前还没有文献报道这些miRNAs在羊毛弯曲形成中的具体作用，其功能尚需进一步研究。