2. 中国农业科学院北京畜牧兽医研究所, 北京 100193;
3. 农业部农业基因数据分析重点实验室, 深圳 518124;
4. 岭南现代农业科学与技术广东省实验室, 深圳 518124
2. Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193, China;
3. Genome Analysis Laboratory of the Ministry of Agriculture and Rural Affairs, Shenzhen 518124, China;
4. Lingnan Guangdong Laboratory of Modern Agriculture, Shenzhen 518124, China
环状RNAs(circular RNAs,circRNAs)是一类通过反向剪接(backsplice)方式形成具有闭合环状结构的内源性非编码RNA分子,其不具有5′端帽子和3′端poly(A)尾巴结构,能够稳定存在于各种类型的真核细胞中。1976年,Sanger等[1]首次从植物病毒中发现了circRNAs分子。近年来,随着生物信息学与测序技术的发展,越来越多的circRNAs在古细菌[2]、酵母[3]、拟南芥[4]、水稻[5]、小鼠[6]、大鼠[6]、人[6]和猪[7]等众多物种中被广泛发现。一般而言,circRNAs可以分为外显子型、内含子型和外显子-内含子结合型。此外,circRNAs还具有种类多样、结构稳定、细胞和组织特异、时空表达特异和进化保守等特点[7-9]。由于circRNAs能够抵抗RNase R的降解作用,具有高度的稳定性,因此,在疾病诊断和治疗中可作为潜在的稳定生物分子标志物[10]。随着对circRNAs研究的深入,其越来越多的潜在功能被挖掘。circRNAs能够吸附微小RNAs(miRNAs)发挥海绵(sponge)作用,少量circRNAs还能被翻译为小肽分子发挥功能。此外,circRNAs可以与转录调控元件结合或与蛋白互作,进而调控基因的转录[11],并且circRNAs还可以发挥m6A(N6-methyladenosine)修饰作用,促进蛋白质翻译的有效启动[12]。还有些circRNAs可以增加增强子的活性[13],以及参与调控其他生物学过程,如在细胞增殖、分化、自噬和细胞凋亡中具有重要作用。
Li等[14]在大白和莱芜猪皮下脂肪中鉴定出29 763条circRNAs,差异表达分析发现,上调的circRNAs有70条,下调的circRNAs有205条,通过构建circRNA-miRNA-mRNA互作网络发现,circRNAs能够调控成脂分化和脂质代谢。Ran等[15]在猪30和180天两个发育阶段睾丸中鉴定出14 108条circRNAs,其中3 220条在两个发育阶段中都表达,2 819和8 069条circRNAs分别只在30和180天表达,这些数据为在分子水平上理解猪睾丸发育或精子发生的调控提供了重要信息。Venø等[16]在6个发育阶段的猪胚胎大脑中(E23、E42、E60、E80、E100、E115)共鉴定出9 377条circRNAs,从E23到E60的发育过程中,circRNAs的表达量逐渐增多,并且在E60达到顶峰,从E60到E80的发育过程中,circRNAs的数目急剧下降,从E80到E115 circRNAs的数目缓慢下降。表明circRNAs的表达具有高度的时空动态特征,揭示了circRNAs在哺乳动物大脑发育过程中的重要功能。本课题组前期在贵州小型猪9种组织中共鉴定出5 934条circRNAs,其中149条与骨骼肌发育相关[7]。近年来的研究发现,circRNAs在猪的多种生命过程中具有复杂的功能和调控机制。
本研究通过对3个猪种16个发育阶段的9种组织进行混合建库测序,系统鉴定猪circRNAs,并分析其各自分子特征,为进一步研究circRNAs的功能提供丰富的信息。
1 材料与方法 1.1 样品收集及测序本研究所用到的转录组数据来自于本课题组前期对长白、通城和五指山猪9种不同组织和16个生长发育阶段的样品混合测序[17]。发育阶段包括:妊娠33、40、50、55、60、65、70、75、80、85、90、95、100和105 d的胚胎以及出生后0和10 d。组织包括长白、通城和五指山3个品种猪的背最长肌、心、脾、肺、肝、肾、胃、小肠和腿肌。利用Illumina Genome Analyzer II测序平台进行测序,获得1.18亿条总reads,长度为90 bp的双端链特异性序列。
1.2 测序数据的质量控制首先利用FASTQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)软件对测序数据进行质量控制并去除接头序列。然后,通过FASTX-Toolkit(http://hannonlab.cshl.edu/fastx_toolkit/)软件去除低质量的reads。最后,为了保证所有reads均为双末端,使用BBmap(http://sourceforge.net/projects/bbmap/)软件对剪切后的reads进行修复。
1.3 circRNAs的鉴定与生物信息学分析按照图 1A的流程对circRNAs进行鉴定。首先,对BWA建立索引,接着利用BWA MEM[18]将质控后的数据与Ensembl 100上的猪参考基因组11.1版本进行比对。随后,采用CIRI_v2.0.6[19]中的CIRI2.pl脚本进行circRNAs的鉴定,不少于2个junction reads被鉴定为circRNAs。接下来,对鉴定的circRNAs与circAtlas[20](http://circatlas.biols.ac.cn/)数据库进行比对以及分子特征分析,并基于分析结果筛选出3条与骨骼肌生长发育相关的高表达circRNAs,然后结合circAtlas数据库查询统计其物种多样性分数(MCS):
![]() |
A.circRNA鉴定流程图;B.每条染色体上的circRNAs数目;C.不同类型的circRNAs百分比 A. Pipeline for the identification of circRNAs; B. The number of circRNAs on each chromosome; C. The percentage of circRNAs in different types 图 1 circRNAs的鉴定与分布特征 Fig. 1 Identification and distribution characteristics of circRNAs |
$ {\rm{MCS = Ns + Nt \times Ni}} $ |
其中,Ns为cir cRNA同源的物种数;
本研究利用microtar[21]、RNAhybrid[22]和miranda[23]这3种软件分别预测了sus-MYH2_0025、sus-CDK13_0002和sus-FANCL_0003的miRNA靶标。为了提高结果的准确性,取同时在两种软件中都预测到的miRNAs作为最终的预测结果。共预测到32个miRNAs靶标,与骨骼肌发育相关的有20个。
2 结果 2.1 circRNAs的鉴定与分布特征在1~18号染色体上共鉴定到8 216条circRNAs (X、Y和线粒体上共鉴定到270条circRNAs,在后续分析中不作为研究对象),且主要分布在1、6和13号染色体(图 1B)。另外,还统计了circRNAs的类型,发现91%的circRNAs属于外显子型,5%的circRNAs属于内含子型,仅4%的circRNAs位于基因间区(图 1C)。
2.2 猪circRNAs的序列特征为了研究猪circRNAs的分子特征,本研究分析了circRNAs的序列特征和表达水平,以及其对应的蛋白编码基因的特征。统计发现,本研究所鉴定的circRNAs由3 462个基因转录,共来自于10 951个转录本。接下来,比较分析mRNAs和circRNAs的长度密度分布(图 2A),发现circRNAs的平均长度短于mRNAs(circRNAs和mRNAs的平均长度分别为350 nt和3 322 nt)。另外,还发现绝大多数(88.9%)的基因产生1~5条circRNAs,其中38%的基因只产生1条circRNA,1个基因(ENSSSCG00000004979,MYO9A)最多可产生23条circRNAs(图 2B)。
![]() |
A.circRNAs和mRNAs的长度分布密度;B.每个基因产生的circRNAs数目 A. The length distribution density of circRNAs and mRNAs; B. The number of circRNAs generated by each gene 图 2 circRNAs序列和蛋白编码基因的特征 Fig. 2 The features of circRNAs sequence and protein-coding gene |
由图 3可知,circRNAs的平均GC含量为47.11%,mRNA的平均GC含量为46.66%,两者无显著差异。
![]() |
图 3 circRNAs和mRNAs的GC含量比较 Fig. 3 The comparison of GC content between circRNAs and mRNAs |
为了研究每条染色体上产生circRNAs与mRNAs数目的相关性,利用R语言中的cor.test函数计算circRNAs与mRNAs的相关系数。结果发现,两者的相关系数R=0.84, (P=1.285×10-5)。表明染色体上circRNAs与mRNAs的数目呈显著正相关(图 4)。
![]() |
灰色区域代表 95%的置信区间 Gray area represents 95% confidence interval 图 4 每条染色体上的circRNAs与mRNAs的相关性 Fig. 4 The correlation between circRNAs and mRNAs on each chromosome |
基于对circRNAs表达量的分析,本研究筛选出3条高表达的circRNAs (sus-MYH2_0025、sus-CDK13_0002和sus-FANCL_0003)。在circAtlas[20]数据库中查询发现,sus-MYH2_0025与猪骨骼肌生长发育相关,位于12号染色体上(chr12:55199495|55262682-),为外显子型circRNA,其宿主基因为MYH2,junction reads为334条,长度为91 nt;circAtlas数据库利用MultiMSOAR v2.0软件[24]计算得出,物种多样性分数(MCS)为1.087 3,组织特异性指数为0.638,表明该circRNA具有高度的组织特异性。sus-CDK13_0002与猪骨骼肌生长发育相关,且存在于大脑、睾丸、肺、心、肝、肾和脾等多种组织中,位于18号染色体上(chr18:54393586|54404485-),为外显子型circRNA,其宿主基因为CDK13,junction reads为264条,长度为1 142 nt,物种多样性分数为7.0,组织特异性指数为0.012;值得一提的是,该circRNA在猪、人和小鼠中高度保守,并且在狗和鸡中也存在。sus-FANCL_0003也与猪骨骼肌生长发育相关并且存在于大脑、睾丸、肺、心、肝、肾和脾等多种组织中,位于3号染色体上(chr3:83349428|83384802+),为外显子型circRNA,其宿主基因为FANCL,junction reads为213条,长度为637 nt,物种多样性分数为2.866 67,组织特异性指数为0.35;另外,该circRNA在猪和大鼠中高度保守。
从miRBase[25]数据库中获得457条猪miRNAs序列,接着利用microtar、RNAhybrid和miranda软件对sus-MYH2_0025、sus-CDK13_0002和sus-FANCL_ 0003进行靶标预测。为了降低结果的假阳性,将两种软件都预测到的靶标作为最终的预测结果。结果发现,可以与sus-MYH2_0025靶向结合的miRNAs有11个,分别为miR-381-5p、miR-1839-5p、miR-31、miR-193a-5p、miR-9808-3p、miR-143-3p、miR-9824-5p、miR-17-3p、miR-9784-5p、miR-4331-3p和miR-22-3p(图 5A),可以与sus-CDK13_0002靶向结合的miRNAs有5个,分别为miR-370、miR-148b-3p、miR-9810-3p、miR-6782-3p和miR-2320-5p(图 5B),可以与sus-FANCL_0003靶向结合的miRNAs有4个,分别为miR-9857-5p、miR-9804-5p、miR-769-3p和miR-9820-5p(图 5C)。
![]() |
A.sus-MYH2_0025结合的miRNA靶标; B. sus-CDK13_0002结合的miRNA靶标;C. sus-FANCL_0003结合的miRNA靶标。中心代表circRNA;外周代表与该circRNA结合的miRNAs A. The miRNA target of sus-MYH2_0025; B. The miRNA target of sus-CDK13_0002;C. The miRNA target of sus-FANCL_0003. The center represents the circRNA; the outer periphery represent the miRNAs bound to the circRNA 图 5 与sus-MYH2_0025、sus-CDK13_0002和sus-FANCL_0003结合的靶标miRNAs预测 Fig. 5 miRNAs prediction target to sus-MYH2_0025, sus-CDK13_0002 and sus-FANCL_0003 |
猪不仅是重要的动物蛋白来源,也是一种理想的大动物模型,研究其转录组及其表达特征对动物育种和生物医学应用都具有重要意义。本研究为了系统鉴定猪circRNAs并分析其分子特征,对3个品种猪16个生长发育阶段的9种组织进行了混合测序。通过对circRNAs的系统鉴定与分析,发现91%的circRNAs属于外显子型,5%的circRNAs属于内含子型,仅4%的circRNAs属于基因间型。这些circRNAs主要分布在1、6和13号染色体,平均长度为350 nt。Liang等[7]在贵州小香猪9种组织中发现,68.40%的circRNAs属于外显子型,21.93%的circRNAs属于基因间型。Yan等[26]发现,猪脾上74.31%的circRNAs属于外显子型,20.36%的circRNAs属于内含子型,仅5.33%的circRNAs属于基因间型。Chen等[27]在长白×约克夏猪的垂体中发现,91%的circRNAs属于外显子型,6%的circRNAs属于内含子型,3%的circRNAs属于基因间型等。以上研究表明,circRNAs主要位于外显子区域,并且外显子型circRNAs是目前关注的焦点。
Memczak等[28]研究发现,81条小鼠circRNAs序列在人上保守。23.6%的小鼠circRNAs在大鼠中保守[29]。另一研究发现,20%的小鼠circRNAs在人细胞系中保守[30],这一结果高于小鼠和人中仅有4%的circRNAs保守的结果[28],原因在于样品测序深度不同。小鼠大脑中15%~20%的circRNAs在猪大脑中保守[16],5%~10%的人大脑circRNAs在猪大脑中保守[29]。Liang等[7]研究发现,20.20%的猪circRNAs序列与人的保守,16.96%的猪circRNAs序列与小鼠的保守。
此外,有研究报道,sus-MYH2_0025的靶标miR-22-3p能抑制平滑肌细胞的增殖和迁移[31]。与sus-CDK13_0002靶向结合的其中一个miRNA为miR-2320-5p, 其序列为UGGCACAGGGUCCA-GCUGUCGG,在猪的发育过程中与脂肪沉积密切相关[32]。另外发现,sus-FANCL_0003结合的miR-793参与凋亡以及衰老等信号通路[33]。
研究发现,sus-MYH2_0025与质膜修复和肌动蛋白介导的细胞收缩,如心肌细胞收缩有关,其宿主基因MYH2是MyHC(肌球蛋白重链)基因家族中的一员,具有高度的发育阶段和时空表达特异性。MyHC亚型位于成熟骨骼肌的快速氧化糖酵解纤维中,并且在猪MyHC基因的上游调控区域中表现出表达分化特异性[34]。另外,在大鼠中,IIA和IIX MyHC基因在围产期表达[35]。猪骨骼肌MyHC基因的2号内含子中含有正调控元件,该元件可以结合骨骼肌特异性因子。sus-CDK13_0002的宿主基因为CDK13,由该基因编码的蛋白质是细胞周期蛋白依赖性丝氨酸/苏氨酸蛋白激酶家族中的一员。CDK13在细胞周期调控和mRNA加工中发挥重要作用,并可能参与造血调节。与CDK13存在相互作用的基因为ENO1、ENO2、TGFB1、ELL、COMP、ICK、NAA50、EED、NTS。此外,CDK13还与一些疾病有关,如大肠癌、钱德勒综合征、角膜营养不良、知觉性聋以及肝癌。sus-FANCL_0003的宿主基因为FANCL,该基因编码泛素连接酶,并且与DNA修复、蛋白质单泛素化、细胞对DNA损伤刺激的反应、蛋白质泛素化以及调节细胞增殖有关。GWAS研究发现,FANCL与高山山羊脓包病[36]和范可尼贫血病(一种常染色体或X连锁隐性疾病)有关[37]。在人和小鼠中,FANCL与细胞增殖、细胞分化、细胞凋亡、脂质代谢过程、蛋白质代谢过程有关。
4 结论本研究利用本课题组前期测序的转录组数据,在猪基因组中鉴定了8 486条circRNAs, 并发现sus- MYH2_0025、sus-CDK13_0002和sus-FANCL_ 0003与猪骨骼肌生长发育相关,进而构建了circRNA-miRNA调控网络,为猪circRNAs的功能和机制研究提供基础。
[1] | SANGER H L, KLOTZ G, RIESNER D, et al. Viroids are single-stranded covalently closed circular RNA molecules existing as highly base-paired rod-like structures[J]. Proc Natl Acad Sci U S A, 1976, 73(11): 3852–3856. |
[2] | DANAN M, SCHWARTZ S, EDELHEIT S, et al. Transcriptome-wide discovery of circular RNAs in Archaea[J]. Nucleic Acids Res, 2012, 40(7): 3131–3142. |
[3] | WANG P L, BAO Y, YEE M C, et al. Circular RNA is expressed across the eukaryotic tree of life[J]. PLoS One, 2014, 9(6): e90859. |
[4] | YE C Y, CHEN L, LIU C, et al. Widespread noncoding circular RNAs in plants[J]. New Phytol, 2015, 208(1): 88–95. |
[5] | LU T T, CUI L L, ZHOU Y, et al. Transcriptome-wide investigation of circular RNAs in rice[J]. RNA, 2015, 21(12): 2076–2087. |
[6] | WERFEL S, NOTHJUNGE S, SCHWARZMAYR T, et al. Characterization of circular RNAs in human, mouse and rat hearts[J]. J Mol Cell Cardiol, 2016, 98: 103–107. |
[7] | LIANG G M, YANG Y L, NIU G L, et al. Genome-wide profiling of Sus scrofa circular RNAs across nine organs and three developmental stages[J]. DNA Res, 2017, 24(5): 523–535. |
[8] | JECK W R, SORRENTINO J A, WANG K, et al. Circular RNAs are abundant, conserved, and associated with ALU repeats[J]. RNA, 2013, 19(2): 141–157. |
[9] | SALZMAN J, CHEN R E, OLSEN M N, et al. Cell-type specific features of circular RNA expression[J]. PLoS Genet, 2013, 9(9): e1003777. |
[10] | MENG S J, ZHOU H C, FENG Z Y, et al. CircRNA:functions and properties of a novel potential biomarker for cancer[J]. Mol Cancer, 2017, 16(1): 94. |
[11] | HUANG A Q, ZHENG H X, WU Z Y, et al. Circular RNA-protein interactions:functions, mechanisms, and identification[J]. Theranostics, 2020, 10(8): 3503–3517. |
[12] | YANG Y, FAN X J, MAO M W, et al. Extensive translation of circular RNAs driven by N6-methyladenosine[J]. Cell Res, 2017, 27(5): 626–641. |
[13] | LI Z Y, HUANG C, BAO C, et al. Exon-intron circular RNAs regulate transcription in the nucleus[J]. Nat Struct Mol Biol, 2015, 22(3): 256–264. |
[14] | LI A, HUANG W, ZHANG X, et al. Identification and characterization of CircRNAs of two pig breeds as a new biomarker in metabolism-related diseases[J]. Cell Physiol Biochem, 2018, 47(6): 2458–2470. |
[15] | RAN M L, WENG B, CHEN B, et al. Strand-specific RNA sequencing in pig testes identifies developmentally regulated genes and circular RNAs[J]. Genes Genom, 2017, 39(10): 1083–1094. |
[16] | VENØ M T, HANSEN T B, VENO S T, et al. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development[J]. Genome Biol, 2015, 16: 245. |
[17] | YANG Y L, ZHOU R, ZHU S Y, et al. Systematic identification and molecular characteristics of long noncoding RNAs in pig tissues[J]. Biomed Res Int, 2017, 2017: 6152582. |
[18] | LI H, DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform[J]. Bioinformatics, 2009, 25(14): 1754–1760. |
[19] | GAO Y, WANG J F, ZHAO F Q. CIRI:an efficient and unbiased algorithm for de novo circular RNA identification[J]. Genome Biol, 2015, 16(1): 4. |
[20] | WU W Y, JI P F, ZHAO F Q. CircAtlas:an integrated resource of one million highly accurate circular RNAs from 1070 vertebrate transcriptomes[J]. Genome Biol, 2020, 21(1): 101. |
[21] | THADANI R, TAMMI M T. MicroTar:predicting microRNA targets from RNA duplexes[J]. BMC Bioinformatics, 2006, 7 Suppl 5(Suppl 5): S20. |
[22] | KRVGER J, REHMSMEIER M. RNAhybrid:microRNA target prediction easy, fast and flexible[J]. Nucleic Acids Res, 2006, 34(S2): W451–W454. |
[23] | BETEL D, WILSON M, GABOW A, et al.The microRNA.org resource: targets and expression[J].Nucleic Acids Res, 2008, 36(S1): D149-D153. |
[24] | SHI G Q, PENG M C, JIANG T. MultiMSOAR 2.0:an accurate tool to identify ortholog groups among multiple genomes[J]. PLoS One, 2011, 6(6): e20892. |
[25] | KOZOMARA A, GRIFFITHS-JONES S. miRBase:annotating high confidence microRNAs using deep sequencing data[J]. Nucleic Acids Res, 2014, 42(D1): D68–D73. |
[26] | YAN Z Q, JIANG T T, WANG P F, et al. Circular RNA expression profile of spleen in a Clostridium perfringens type C-induced piglet model of necrotizing enteritis[J]. FEBS Open Bio, 2018, 8(10): 1722–1732. |
[27] | CHEN Z T, PAN X C, KONG Y R, et al. Pituitary-derived circular RNAs expression and regulatory network prediction during the onset of puberty in Landrace×Yorkshire crossbred pigs[J]. Front Genet, 2020, 11: 135. |
[28] | MEMCZAK S, JENS M, ELEFSINIOTI A, et al. Circular RNAs are a large class of animal RNAs with regulatory potency[J]. Nature, 2013, 495(7441): 333–338. |
[29] | YOU X T, VLATKOVIC I, BABIC A, et al. Neural circular RNAs are derived from synaptic genes and regulated by development and plasticity[J]. Nat Neurosci, 2015, 18(4): 603–610. |
[30] | GUO J U, AGARWAL V, GUO H L, et al. Expanded identification and characterization of mammalian circular RNAs[J]. Genome Biol, 2014, 15(7): 409. |
[31] | HUANG S, WANG M, WU W, et al. Mir-22-3p inhibits arterial smooth muscle cell proliferation and migration and neointimal hyperplasia by targeting HMGB1 in arteriosclerosis obliterans[J]. Cell Physiol Biochem, 2017, 42(6): 2492–2506. |
[32] | LI G X, LI Y J, LI X J, et al. MicroRNA identity and abundance in developing swine adipose tissue as determined by Solexa sequencing[J]. J Cell Biochem, 2011, 112(5): 1318–1328. |
[33] | CHEN X, XIE X, XING Y, et al. MicroRNA dysregulation associated with red blood cell storage[J]. Transfus Med Hemother, 2018, 45(6): 397–402. |
[34] | CHANG K C, FERNANDES K, DAUNCEY M J. Molecular characterization of a developmentally regulated porcine skeletal myosin heavy chain gene and its 5' regulatory region[J]. J Cell Sci, 1995, 108(4): 1779–1789. |
[35] | LAFRAMBOISE W A, DAOOD M J, GUTHRIE R D, et al. Emergence of the mature myosin phenotype in the rat diaphragm muscle[J]. Dev Biol, 1991, 144(1): 1–15. |
[36] | MINOZZI G, MATTIELLO S, GROSSO L, et al. First insights in the genetics of caseous lymphadenitis in goats[J]. Ital J Anim Sci, 2017, 16(1): 31–38. |
[37] | HODSON C, COLE A R, LEWIS L P C, et al. Structural analysis of human FANCL, the E3 ligase in the Fanconi anemia pathway[J]. J Biol Chem, 2011, 286(37): 32628–32637. |