禾花鲤与建鲤肌间骨miRNAs测序与分析比较

何苹萍 王 卉 韦嫔媛 林 勇 黄 姻 王韶韶 彭金霞 吴铁军

(广西水产科学研究院, 广西水产遗传育种与健康养殖重点实验室, 南宁 530021)

鲤(Cyprinus carpio)是目前世界上养殖历史最悠久, 养殖范围最广的淡水鱼类, 是许多国家的主要养殖品种[1], 在漫长的自然演化和人工驯养背景下, 形成了许多著名的鲤品种, 如黄河鲤、建鲤、黑龙江鲤等。广西全州禾花鲤(又名禾花鱼), 是最具广西地方特色的稻鱼共养品种, 也是鲤的一个独特的地方土著品种。禾花鲤因经稻田长期驯化, 食禾花而得名, 它体型粗短、鳞细皮薄、肉质细嫩、味道鲜美, 特别是肌间骨柔软细小且可食用, 清代乾隆年间曾成为宫廷贡品。建鲤是以荷包红鲤和沅江鲤为亲本选育的良种, 生长速度快、抗病能力强、含肉多。建鲤是我国养殖鲤鱼的代表性品种之一, 在我国广泛养殖[2,3]。

鱼类的骨骼肌间骨又称为肌间刺, 是位于脊椎骨两侧肌间隔中的硬骨小刺, 是由肌膈间的结缔组织直接骨化而形成的[4]。肌间骨在低等真骨鱼类中普遍存在, 随着鱼类的演化, 其形态逐渐复杂化[5]。众多细小的肌间骨一直以来影响多种鱼类的食用品质, 因此, 鱼类肌间骨的形态和发育机制一直是水产研究的热点。其中, 关于真骨鱼类肌间骨的骨化模式研究主要集中在黄鳝(Monopterus albus)[6]、鲢(Hypophthalmichthys molitrix)[7]、鲫(Carassius auratus)[8]、鲤(Cyprinus carpio)[9]和团头鲂(Megalobrama amblycephala)[10], 这些研究主要集中在肌间骨的形态、数量和分布等规律。而关于肌间骨发育分子调控研究的仅仅有罗非鱼(Tilapia)[11]、鲤(Common carp)[11]和团头鲂(Megalobrama amblycephala)[12]。

小分子RNA(microRNA, miRNA)是一类进化保守的、内源性的、非编码的小RNA。miRNA通过与靶基因的3′非翻译区(Untranslated regions, UTR)结合从而抑制靶基因的翻译或降解靶基因[13]。1993年在秀丽隐杆线虫(Caenorhabditis elegans)时序性发育调控过程中发现了首个miRNA(lin-4)[14];2000年又在线虫中找到了第2个调控时序性发育的miRNA基因let-7[15]。miRNA参与许多重要的生物学过程, 如在个体发育、组织分化、病毒感染等逆境适应过程中都具有重要作用。miRNA在鱼类的研究比较多地集中于免疫[16]方面。近年来, miRNA在鱼类发育过程中的作用也越来越受到重视, 鱼类骨骼要经历软骨化和硬骨化的过程, 受多种转录因子的调控, 如Sox9[17]和Runx2[18]等。而miR-140[19]、miR-204[20]和miR-211[20]则调控上述斑马鱼骨化转录因子的表达。以上结果表明, miRNA是调控鱼类骨骼形成的重要因素。研究发现, miRNA一般是通过多个信号通路来调控骨骼发育[21]。但目前为止, miRNA在水产动物肌间骨发育中的作用研究仅仅在团头鲂中开展过[22,23]。

作为广西稻田养鱼的主要本地特色品种禾花鲤, 其相关的基础研究十分缺乏。本研究针对禾花鲤肌间骨细小且柔软可食用的特点, 开展其分子调控机制的研究。采集禾花鲤及建鲤的肌间骨进行miRNAs的高通量测序, 进行生物信息学比较分析,为禾花鲤肌间骨发育的分子调控机制研究奠定基础。

1 材料与方法

1.1 实验动物

本研究所用的实验动物为来自全州的禾花鲤及来自南宁的建鲤(均为24月龄), 均暂养在温度为25℃的水中。

1.2 组织样本采集

分别采集禾花鲤及建鲤的肌间刺数根, 于液氮中速冻。

1.3 小RNA文库构建和测序

用EASYspin Plus骨组织RNA快速提取试剂盒分别提取禾花鲤及建鲤骨头的总RNA, 随后通过15%变性聚丙烯酰胺凝胶电泳分离纯化16—30 nt的小分子RNA, 进行Solexa高通量测序。小RNA文库的Solexa测序由武汉华大生物公司完成。

1.4 生物信息学分析

一般来说, 小RNA的长度区间为18—30 nt。在测序完成后, 首先进行序列去杂, 筛选一定长度范围内的sRNA, 用软件bowtie (http://computing.bio.cam.ac.uk/local/doc/bowtie2.html#)将长度筛选后的sRNA定位到基因组参考序列上, 分析small RNA在参考序列上的分布情况, 然后将上述mapped到参考序列上的reads, 在Rfam (http://Rfam.sanger.ac.uk/)、miRBase (http://www.mirbase.org/)中进行比对, 得到各样品匹配上的sRNA的详细情况, 包括匹配上的已知miRNAs的二级结构, 各样本中miRNAs的序列、长度、出现的次数等信息。进一步通过miRNAs前体的标志性发夹结构, 使用miRDeep2[24]预测新的miRNAs。对各样本中已知和新miRNAs进行表达量的统计, 并用TPM进行表达量归一化处理。参照Audic等[25]发表在Genome Research上的基于测序的差异基因检测方法开发ExpDiff进行样品间表达差异分析。根据差异miRNAs检测结果, 使用R软件中的pheatmap函数进行层次聚类分析。我们用miRanda[26]来对分析得到的已知、novel miRNAs进行靶基因预测。得到各组比较间的差异表达miRNAs后, 根据miRNAs与其靶基因间的对应关系, 对每组差异表达miRNAs的靶基因的集合分别进行GO (Gene Ontology)和KEGG (Kyoto Encyclopedia of Genes and Genomes)富集分析。

2 结果

2.1 小RNAs的基本特征

为了比较禾花鲤及建鲤肌间刺miRNAs的表达差异, 本研究分别从禾花鲤及建鲤肌间刺提取小RNA, 构建了2个小RNA文库。对这2个文库进行Solexa测序, 分别产生约30179426和29750080条序列。去除低质量序列、接头序列和短序列后, 分别从禾花鲤和建鲤的肌间刺文库中筛选出25474895和24625715条高质量序列进行进一步的分析。结果表明, 禾花鲤以及建鲤小RNA长度主要分布在17—30个核苷酸, 其中21个核苷酸的小RNA最丰富, 其次是20个核苷酸、22个核苷酸(图 1)。

2.2 MiRNAs的鉴定

在上述高质量小RNA序列中, 禾花鲤与建鲤分别有20765041 (84.32%)和21326842 (83.72%)条序列与参考基因组比对成功。比对上的小RNA通过与Rfam、miRBase数据库比对, 进行各种小RNA分类, 其中, 禾花鲤miRNAs占其小RNA总量的9.80%, 远低于建鲤(40.2%, 图 2)。从禾花鲤和建鲤中分别鉴定出已知的成熟miRNAs 595和570种。

2.3 MiRNAs的差异表达分析

通过比较样品之间的表达量, 筛选差异表达miRNAs。我们采用ExpDiff进行差异miRNAs筛选。禾花鲤与建鲤相比, 一共有351个差异miRNAs,其中84个上调表达, 267个下调表达。部分显著差异表达的miRNAs见表 1。

图 1 小RNA长度分布Fig. 1 Length distribution of small RNA sequencesA. 建鲤; B. 禾花鲤A. Jian carp; B. Rice flower carp

图 2 小RNA分类注释Fig. 2 Classification annotation of small RNA sequencesA. 建鲤; B. 禾花鲤A. Jian carp; B. Rice flower carp

表 1 部分显著差异表达及成骨发育相关的miRNAsTab. 1 Partial significantly differential expressed or osteogenesisrelated miRNAs

2.4 差异miRNAs的靶基因GO显著性富集分析

我们用WEGO软件对差异miRNA的靶基因做GO功能分类统计, 从宏观上认识差异基因的功能分布特征, 部分结果展示如表 2。GO总共分为三大功能类, 分别描述基因的分子功能(Molecular function)、所处的细胞位置(Cellular component)和参与的生物过程(Biological process)。差异miRNA的靶基因注释最多的为Biological process, 其中以细胞过程(Cellular process)为最多, 其次为单一组织过程(Single-organism process)、新陈代谢过程(Metabolic process)及生物调节(Biological regulation); Molecular function中以结合(Binding)及催化活性(Catalytic activity)为最多; Cellular component中以细胞(Cell)、细胞器(Organelle)为最多。

2.5 差异miRNAs的靶基因Pathway显著性富集分析

在生物体内, 不同基因相互协调行使其生物学功能, 基于Pathway的分析有助于更进一步了解基因的生物学功能。我们以KEGG公共数据库对差异miRNAs作Pathway富集分析。图 3显示了排名前20的差异miRNAs的靶基因富集的Pathway。其中包括signaling pathways regulating pluriotency of stem cells (调节干细胞多能性的信号通路)、axon guidance(轴突导向)等。

表 2 部分差异miRNAs靶基因的GO分类Tab. 2 Partial GO of differentially expressed miRNAs

图 3 前20的差异miRNAs pathway富集统计散点图Fig. 3 The top 20 pathways enriched by the putative target genes of the differentially expressed miRNAs

3 讨论

Li等[27]发现在MC3T3细胞成骨分化过程中,miR-29b可以通过降解抑制成骨的靶基因HDAC4、TGF-β3、ACVR2A的表达促进成骨。团头鲂肌间骨中miR-29的表达量较结缔组织中的低, 推测miR-29在团头鲂中可能起到抑制成骨的作用[23]。miR-29虽然与miR-29b为一个家族的成员, 但功能却不同, 也可能是在不同物种中作用的差异。He等[28]发现miR-20b通过降低PPARγ的表达提高了成骨的核心转录因子Runx2转录, 促进成骨分化。有研究发现miR-27可以抑制靶基因APC, 从而反激活WNT促进成骨分化[29]。miR-210在ST2成骨分化过程中可能通过抑制成骨分化的重要基因AcvR1b来促进成骨[30]。Oskowitz等[31]则从由美国国立卫生研究中心提供的MSCs其诱导分化的成骨细胞中,利用芯片筛选出在成骨细胞中高表达的miRNAs,如hsa-miR-130a、has-miR-199a、hsa-miR-346、hsa-miR-21及hsa-miR-10a等。在本研究中禾花鲤肌间骨中miR-29b、20b-5p、miR-27-3p、miR-210-3p、miR-130a、miR-199a及miR-21的表达与建鲤相比均为下调表达。这些均为上述已报道的促进人类成骨的miRNAs, 因此我们推断禾花鲤通过下调表达促进成骨的miRNAs来抑制成骨。

Eskildscn等[32]发现miR-138通过抑制FAKERK1/2信号通路降低下游成骨的核心结合因子Runx2和OSX的表达来降低成骨作用。Li等[33]发现在BMP2诱导的C2C12细胞成骨分化中, 对其中显著下调的miR-133和miR-135进行研究后发现, miR-133的靶基因为Runx2, 而miR-135的靶基因为mad5,两者共同作用来抑制骨的形成。团头鲂肌间骨中miR-133b的表达量随着鱼龄的增长呈增加的趋势,推测miR-133b在团头鲂肌间骨的形成过程中起到促进的作用[28]。miR-133b虽然与miR-133为一个家族的成员, 但功能有所不同, 也可能是在不同物种中作用的差异。Wei等[34]发现, miR-34b通过两条途径来抑制成骨, 一方面通过阻碍CDK4、CDK6和Cyclin D1蛋白表达来减少成骨细胞的数量; 另一方面通过降低SATB2蛋白的表达来抑制成骨细胞的分化。Gao[35]发现hsa-miR-31、hsa-miR-106a、hsamiR-148a通过抑制BMSCs的分化来抑制成骨的过程。在本研究中禾花鲤肌间骨中的miR-133、miR-135、miR-34b、hsa-miR-31、hsa-miR-148a以及miR-138的表达与建鲤相比均为上调表达。这些均为上述已报道的抑制人类成骨的miRNAs, 因此我们推断禾花鲤通过上调表达抑制成骨的miRNAs来达到抑制成骨的目的。

总之, 我们推测禾花鲤可能是通过下调表达促进成骨miRNAs及上调表达抑制成骨miRNAs来抑制成骨过程, 从而维持其肌间骨细小及柔软的特性。该研究为禾花鲤肌间骨发育的分子调控机制研究奠定了基础。