全转录组学在畜牧业中的应用

石田培，张莉



全转录组学在畜牧业中的应用

石田培，张莉

中国农业科学院北京畜牧兽医研究所，北京 100193

RNA作为一种大分子参与基因编码、解码、调控、表达等多种生物学过程。目前，对RNA的功能研究主要通过全转录组测序方法来完成。全转录组研究可以对基因结构与功能进行更深层次地分析和探究，揭示基因表达与生命现象之间的内在联系。现阶段，基于高通量测序技术的转录本结构研究、基因表达水平研究及非编码区域功能研究在模式动物、猪、禽类中已大量开展，但在羊上却鲜有报道。本文介绍了利用RNA-seq及Small RNA-seq 技术研究全转录组的一般流程及常用策略，综述了全转录组学技术在畜牧业领域中的研究进展。

基因表达；全转录组；RNA测序技术；高通量测序

在后基因组学研究中，转录组学是解读基因组功能元件和揭示细胞及组织分子机制的基础，在生物表型和基因表达研究中占据了重要地位[1]。转录组一般指狭义转录组，即能够编码蛋白的RNA， 而生物学中的转录本概念则为生物体特定细胞或组织在某一特定条件下的所有转录产物，涵盖了编码RNA (mRNA)和非编码RNA (non-coding RNA, ncRNA)[2]。同相对稳定的基因组比较，转录组会根据生理状态、生长阶段、生长环境的变化而变化[3]，具有高度动态可变的特性。除此之外，转录组图谱中含有丰富的生物学信息，包括基因表达丰度及差异、基因结构、反义转录本、选择性剪切、单核苷酸多态性和基因融合等。

生物体是一个错综复杂的网络，构成复杂有机体的任何简单分子都不能孤立存在或行使功能，如转录组中的各类RNA。这些不同类型的RNA生物学功能各异，参与众多信号通路的调控。单一的mRNA或ncRNA研究缺乏关联性，对多种RNA 信息进行整合分析，探索潜在的网络调控机制则成为转录组学研究的趋势。在传统转录组学(以mRNA[4]为主)研究中，样品总RNA的提取通常利用Oligo (dT)磁珠富集法，该方法能有效发掘富含poly (A)尾巴的RNA[5]，包括mRNA和一小部分长链非编码RNA (long non-coding RNA, lncRNA)，但过滤掉了其他poly (A) minusRNA组分，即无poly (A)尾RNAs，因此导致circRNA (circular RNA)和部分lncRNA的丢失，不能全面地反映转录本的真实情况。Yang等[6]首次提出了基于poly (A) minus富集的RNA-seq方法，可将sno-lncRNAs (small nucleolar RNAs)、ciRNAs (circular intronic RNA)和circRNAs等这类非poly (A)结尾的新型RNA分子富集并保留下来。至此，包括miRNA[7]、lncRNA[8,9]和circRNA[10,11]等非编码RNA的全转录组测序研究应运而生。本文综述了基于RNA-seq 和Small RNA-seq技术的全转录组研究概况，总结了全转录组测序的一般流程和常用策略，最后对畜牧业未来的发展趋势进行展望。

1 全转录组测序概述

全转录组是指在一定时空状态下有机体特定细胞或组织能够转录出的所有转录本[12]，其中蕴含着重要的生物学调控规律。丰富的RNA组学信息可用于开展差异表达ncRNA的靶基因预测、差异表达ncRNA和mRNA的正负相关性分析、ceRNA (com­peting endogenous RNAs)对靶向mRNA的功能富集分析及关键基因挖掘等研究。全转录测序是通过高通量RNA-seq测序技术对RNA序列进行测序，并通过一些定量平台和技术反应其表达量和表达差异，从而形成表达图谱。全转录组学分析就是以表达谱为基础，对RNA进行鉴定和注释，预测相应靶细胞或编码潜能，并基于GO(Gene Ontology)、COG (Clusters of Orthologous Groups)、KEGG (Kyoto Encyclopedia of Genes and Genomes)等数据库进行功能富集、聚类分析、信息挖掘和通路探究等。不同于传统转录组学，全转录组学在mRNA研究基础上，涵盖了lncRNA、circRNA、miRNA、siRNA(small interfering RNA)和piRNA(piwi-interacting RNA)等多种非编码RNA及其之间的调控网络分析。目前，转录组测序主要集中在单一RNA的定量表达和注释分析上，而全转录组学测序可同时对多种RNA进行鉴定及关联分析(图1)。

2 全转录组测序研究的一般技术流程和方法

2.1 测序样品的准备

在样品采集过程中，组织中的RNA 极易受内源或外源RNA酶作用而降解，同样也容易受到蛋白质、DNA、同源和外源酚类等物质的污染，因此，样品质量和保存条件是决定试验结果质量的关键。不同物种不同组织部位的RNA得率大不相同，在动物中，肝脏、脾脏和胰脏组织的得率较大，而肌肉、脑和心脏等组织得率较小。在皮肤、骨骼和毛发中RNA的提取难度较大，效果也较差。在准备样品时，首选新鲜组织，剔除非研究所需的组织类型。若组织体积较大，多管分组置于-80℃或液氮中长期保存。

2.2 建库测序

2.2.1 cDNA文库构建和RNA测序

总RNA质检后便可建库，由于miRNA与其他3类RNA序列长度差异较大，需要使用SE (single- end)和PE (paired-end)两种测序策略同时建立两个文库。链特异性建库的具体步骤如下(图2)：(1)通过试剂盒在样品总RNA中去除rRNA；(2)利用六碱基随机引物合成第一条cDNA链；(3)第二条cDNA链合成时，用dUTP代替dTTP，使链上布满U位点；(4)在3¢端加A，加接头；(5)使用USER酶在尿嘧啶位置上产生一个单核苷酸缺口，借此消化掉第二条链，只保留第一条链；(6) PCR扩增。小RNA建库的具体步骤如下(图3)：(1)在3¢端加接头；(2)在5¢端加接头；(3)反转录扩增；(4)用凝胶纯化富集扩增的cDNA后测序。在测序时，理论上数据量越大越利于后续低丰度基因的完整组装，但实际上并非数据量越大越好，需要根据物种情况及相关研究决定数据量的大小。

图1 全转录组测序一般流程

图2 链特异性建库流程

2.2.2 原始数据预处理

测序产生的数据是一系列不能直接使用的原始数据，主要是FASTQ格式的读段(reads)。其中除实验所需的碱基质量信息之外，还包括测序仪器名称、上机次数、试剂型号等信息。原始数据在建库过程或者测序过程中会产生大量低质量读段、较短的序列、含N序列甚至是一些污染序列(如细菌基因序列)[13]、接头序列、载体序列等，因此在数据分析之前，应对原始数据进行过滤、剪切和校正，以确保后期读段定位、转录本组装、基因定量等流程的顺利开展。目前可用的质控软件是FastQC和NGSQC[14]、Qualimap2[15]、HTQC[16]、QCchain[17]、almostSigni­fican[18]、fastq-clean[19]及FaQCs[20]等，最常用的是FastQC和HTQC。FastQC可作为一个单独的JAVA程序，在速度上远超其他工具，极短时间内就可以运行数千万的读段，输入文件的选择性也较大，如压缩或未压缩的EASTQ或SAM/BAM文件。FastQC除能列出读段的数目及质量编码以外，还能可视化报告碱基内容及质量、读段长度和重复序列等信息。在实际操作中，可联合使用上述质控软件以达最佳过滤效果。

图3 小RNA建库流程

2.2.3 读段定位及转录本重组装

由原始数据进行深度清理及质量控制后获得待分析数据(clean reads)，需要通过一系列软件将读段比对到参考基因组或者转录本上，并根据实际定位情况进行转录本组装[21]。在全转录组数据分析时，选择比对程序时还需考虑剪接比对(unspliced aligners)情况[22]。当生物体没有内含子或进行miRNA测序时，可以使用Quality (MAQ)[23]、Burrow-Wheeler Aligner (BWA)[24]和Bowtie[25]等连续比对软件，这些方法用于识别已知外显子和接头，不能识别由可变剪切所产生的新型外显子。但当读段定位至外显子-外显子的交界接头处或作图到有内含子的基因组时，则必须使用剪接比对程序，如TopHat[26]、Map­Splice[27]、STAR[28]和GSNAP[29]等。可先将读段分成两个短片段再参与定位，同时还可记录分割数据以供查找后续转录本，该方法能够识别由可变剪切产生的新转录本。一般情况下，先利用Bowtie进行初步匹配，用匹配成功的reads来获取基因组覆盖区域，预测剪接点(splice junction)局域。然后利用Tophat将未成功匹配的reads划分到splice junction序列上。如果研究物种没有参考基因组或转录组时，需要将读段自行组装成可用的参考序列，然后将所有测序读段通过从头组装生成重叠群和单一序列，再进行后续数据分析。目前，有多种组装程序可供分析使用[30,31]，如Newbler[32]、MIRA[33]、CAP3[34]、Seq-Man[35]、TGICL[36]、stackPACK[37]、Velvet[38]、AbySS[39]、ALLPATHS2[40]、Oases[41]、SOAPdenovo[42]、Multiple-k method[43]、Scaffolding using translation mapping (STM)、Trinity[44]和PCAP[45]等。

2.2.4 表达量化及差异分析

经过读段定位后，可根据读段在转录本上的分布情况预测基因丰度。一般来说，通过软件分析获得的注释文件中会含有转录本分布信息，能够通过分析读段的匹配情况来识别新的转录本。在高等动物的生命活动调控过程中，存在可变剪切(alterna­tively splice, AS)等一系列复杂的调控系统，致使全部转录物不能直接对应基因。此外，还存在不同转录本外显子共享、读段不能跨越多个剪接位点等问题。因此，需要先进行可变剪接识别[46]。现阶段，常用转录本识别软件包括Stringtie[47]、cufflinks[48]、CIDANE[49]、GRIT[50]、TransComb[51]、iReckon[52]、SLIDE[53]、Montebello[54]、Augustus[55]、IsoLassocia­tion[56]、Scripture[57]、Traph[58]和MITIE等。人()和斑马鱼()等基因信息注释完整的物种可直接进行基因预测定量分析，但对于一些基因注释信息不完整的物种需要先进行转录本的预测。目前已完成的大量转录组测序数据快速完善了遗传数据库，反复检验校对了基因注释信息，推动了分子生物学的研究进展。

测序过程中测序深度、基因片段大小、运用算法、实验批次等因素极易造成误差，所以在定量时应使用标准化的方法消除差异，最常用的样本内标准化方法包括RPKM[59]、FPKM、TMP[60]和KPKM等。根据比对到基因上的reads或者fragments数目，进行基因长度和测序总量归一化后，即可统计表达量。常用的转录本定量软件包括HTSeq[61]、feat­ureCounts[62]、StringTie、Cufflinks、RSEM[63]、Sailfish[64]、kallisto[65]、IsoLasso和NURD[66]等。为确定在不同条件或平台的表达差异情况，经过量化和标准化表达水平的基因仍需利用测序覆盖度估测表达丰度分布[67,68]。

RNA-seq数据的差异表达(differential expression, DE)分析往往以单变量的方式进行，很难拟合出一个普遍适用的统计模型，因此根据生物统计学原理将DE分析软件进行下列分析对比(表1)。为使表达分布符合统计方法的假设，RNA-seq数据会被某种方式归一化。按照归一化在DE分析中的前后顺序可将常用软件分为两类：第一类是先进行标准化处理再通过统计学原理计算表达差异，包括TMM[69]、DESeq[70]、PoissonSeq[71]和UpperQuartile[72]等；第二类是将归一化作为DE分析的步骤，在处理时直接执行归一化，如FPKM、RPKM和TPM，但可信度较差。

表1 差异表达分析软件

2.3 基因功能分析

生物信息学主要用数学及统计学的方法对生物信息进行存储、分配、检索及分析，是一门集合生命科学、自然科学与信息科学的新学科。目前，生物信息学在生命科学研究中应用十分广泛，在全转录组研究中也占据重要地位。在全转录组学分析中主要以差异基因的功能富集、候选基因的功能预测、调控网络的构建等分析为主。当前以GO、COG及KEGG数据库为基础的分析工具包、软件和网站众多，如基于R语言的GSA、PADOG、SAFE、Globaltest、Sigpathway、GAGE、GSVA、PLAGE、ZSCORE、SSGSEA、MRGSE、ANCOVA、CAMERA、SPIA、TopoGSA、ToPASeq、NetGSA、NEA和GOGANPA等工具包；DAVID、GOstat、GenMAPP、GOMiner、Onto-Express、EnrichNet和NOA等网站；MetaCore、Cytoscape和GSEA等软件。以上分析软件各有千秋，功能也不尽相同，在数据处理与分析时，应根据实验目选择合适的方法。如需更系统地反映分子调控机制，还可结合其他基因组学手段如全基因组测序、RNA甲基化、表观修饰等数据进行联合分析。

3 全转录组学技术在畜牧业中的应用

3.1 全转录组测序在猪中的应用

猪()的饲养对畜牧业发展乃至民生都至关重要。不同品种猪的肌纤维类型对猪肉品质有很大影响，但潜在的分子机制仍不清楚。近年来，有关猪非编码RNA的功能研究开展得很多。为探究miRNA在猪骨骼肌中的作用，Mcdaneld等[75]分别选取了增殖中的卫星细胞(4~代)、胚胎(60、90、105日龄)、出生胎儿和成年猪组织进行了全转录组研究，发现了12个新型miRNA与肌肉生长发育密切相关，其中肌肉特异性miR-206在卫星细胞中近乎缺失，但在其他发育阶段表达量却很高；miR-1在成年猪中的表达丰度最高；miR-133在胚胎期和初生时的丰度极低甚至检测不到，但miR-368、miR-376和miR-423-5p在初生猪中极高；miR-432表达量在胚胎发育早期最高，随后逐渐降低。该研究是对猪胚胎发育期骨骼肌miRNA较为全面的转录组分析，为深入探究猪骨骼肌miRNA调控机制提供了理论基础。沈一飞[76]利用RNA-seq和Small RNA-seq技术对瘦肉型约克夏猪和脂肪型金华猪进行了甲状腺组织的mRNA、lncRNA和miRNA共表达鉴定与比较，结果发现差异表达mRNA 492个，差异表达lncRNA 48个和miRNA 18个。通过功能分析和聚类，共有256个RNA(其中18个miRNA、1个lncRNA)参与到同一个调控网络中。在整个调控网络中，发现ssc- miRNA-221-5p、ssc-miRNA-708-5p、ssc-miRNA-532- 3p和novel_12等发挥重要调控作用。Li等[77]对从胎儿期到成年期期间获得的组织混合物制备的10个小RNA测序文库中的荣昌猪(Rongchang pigs) miRNA进行全面检测，通过哺乳动物miRNA、前体发夹(pre-miRNA)、高覆盖率猪基因组装配(2009年4月)和表达序列标签(EST)的分析，将猪miRNAome的所有组成部分扩展到867个pre- miRNAs (623个基因组坐标)，编码1004个miRNA，其中777个是独特的。对47个组织特异性样品中选定的30个miRNA进行qRT-PCR实验，发现测序数据和试验结果一致。Sun等[78]对长白猪(Landrace)和兰塘猪(Lantang pigs)背最长肌进行了全转录组测序研究，从22 469个编码转录物中筛选出547个差异表达mRNA，通过生物信息学分析挖掘出与肌肉生长发育相关的 17个基因。此外还发现差异表达的5566个lncRNA和4360个circRNA。其中，3376个lncRNA和1401个circRNA在Lantang文库中上调表达，而1590个lncRNA和2959个circRNA下调表达。通过结合匹配的miRNA谱分析测序数据，鉴定出26种参与ceRNA网络的海绵载体，包括19种lncRNA、40种circRNA和9种mRNA。全转录组研究提供了一种全新的分析方法，对解析猪肌肉生长发育规律和疾病发生机制具有重要意义。

3.2 全转录组测序在禽类中的应用

我国家禽遗传资源丰富，地方品种各具特色，但由于总体生产水平较低，许多优良性状利用效率不高。随着分子生物学和生物信息学的发展与融合，通过各种技术和手段，已鉴定出一批与生长、繁殖等重要经济性状相关的分子标记和候选基因。Yu等[79,80]选取了6和10日龄的鸡胚左(L)、右(R)卵巢作为样品，通过全转录组测序技术对卵巢的退化进行了研究。在6R样品中产生31 066 414个序列读数，在6L样品中产生31 900 200个序列读数，在10R样品中产生31 400 070个序列读数，在10L样品中读取35 504746个序列。这些数据中，大约73.33%的序列可以定位到鸡胚卵巢的参考基因组序列上。鸡胚卵巢的发育受许多基因和信号通路的调控，通过对差异表达基因比对和通路功能聚类发现了22个与卵巢发育和退化相关的基因。其中转录水平排名前20的卵巢基因可能与碳代谢、吞噬体及钙信号传导密切相关。Glazov等[81]在已发现的miRNA基础上进行了深度挖掘，将3个小RNA文库同时比较并严格区分真正miRNA前体与结构相似的RNA，共鉴定到361个新的miRNA、88个新的miRNA候选物、18个mirtrons (包括6个新的非典型mirtron候选物)和21个mirtron候选物。为识别潜在的禽类特异性miRNA，同时与人、狗()、斑马鱼、爪蟾()和蜥蜴()进行保守进化分析，结果表明只有6种miRNA在非禽类脊椎动物中具有保守性，剩余miRNA可能具有鸟类和/或鸡系特异性。Li等[80]分别从孵育10天、12天、14天和18天的白来航鸡(White Leghorns)蛋胚中收集骨骼肌(胸大肌)，经过RNA-seq技术测序并与已发布的数据库进行比对，筛选到281个新型基因间lncRNA，对这些lncRNA进行保守性分析，并利用UCSC数据库评分，结果表明以上lncRNA的保守性均高于随机非编码序列，但远低于蛋白质编码基因。该研究是首例有关白来航鸡骨骼肌lncRNA的分析，鉴定出的新型lncRNA极大丰富了鸡ncRNA数据库。

近年来，鹅()以其适应性强、生长快、营养成分丰富和投入要求低而备受关注。此前，Kang等[82]通过抑制性消减杂交(supp­ression subtractive hybridization, SSH)方法鉴定了一些与鹅从产卵阶段到产蛋阶段繁殖相关的差异表达基因，Guo等[83]也使用相同的方法在产蛋阶段和育雏阶段发现了若干差异表达的基因。Xu等[84]采用短读序列技术(Illumina)对10只380日龄的雌性浙东白鹅(Goose)的卵巢样本进行了从头转录组装配，使用Illumina RNA-seq和DGE深度测序并绘制出鹅卵巢组织的转录组图谱，得到67 315 996个100 bp的短读数，组装成130 514个独特序列。基于已知蛋白质的BLAST结果，分析鉴定到52 642个目标序列。该研究分析了鹅产蛋、育雏期间的转录变化情况，鉴定出大量与卵泡发育和生殖相关的基因。

病理性肥胖是鸭养殖业所面临的重要问题之一，其分子机制仍然未知。Chen等[85]对家养鸭与野鸭两个品种腹部脂肪进行转录组测序分析，预测了23 699个未注释基因，确定了753个差异表达基因。在北京鸭()中，一些与脂质代谢的相关基因(、和等)和致癌基因(、和等)上调表达，而与肿瘤抑制和免疫相关的基因(、、和等被下调，这些数据表明鸭的肿瘤发生可能与病理性肥胖密切相关。此外，发现280 576个单核苷酸变异在两个品种之间存在差异，包括8641个异构变异，富含参与脂质和免疫相关途径的基因，表明与鸭的代谢功能和免疫相关功能密切相关。

3.3 全转录组测序在反刍动物中的应用

近年来国内牛羊肉市场需求不断增加，这对牛羊的育种和养殖工作提出了更高的要求，只有充分了解牛羊生长发育与繁殖等性状的发生机制，才能提高生产效率，全转录组学则为其提供了全新的研究手段。Di等[86]在滩羊(Tan)和小尾寒羊(STH)的卵巢中鉴定出483个miRNA (包括97个已知的、369个保守的和17个预测的新miRNA)。基于KEGG分析，一些差异表达miRNA的靶基因参与生殖激素相关途径(如类固醇生物合成、雄激素和雌激素代谢和GnRH信号传导途径)以及卵泡和黄体发育相关途径，这对绵羊的繁育具有重要意义。Chang等[87]使用新一代测序技术(Solexa高通量测序技术)研究了绵羊黄体期卵巢组织，鉴定出267种新的miRNA，并利用qRT-PCR和Northern印迹证实了在绵羊卵巢和睾丸中表达的一种新型miRNA (ovis_aries_ovary m0033_3p)。根据序列和结构的一致性，推测ovis_ aries_ovarym0033_3p具有类似于hsa-miR-214-3p的功能，能够参与细胞存活、胚胎发育、繁育生殖和卵巢癌抗性的精细调节。张世芳等[88]采用Solexa技术对5头特克赛尔羊(Texel)进行miRNA深度测序，获得了16 532 850条原始序列读数。通过与哺乳动物成熟miRNA数据库、miRNA前体序列、绵羊基因组数据库的比对分析，更新miRNA前体序列库至1529条，编码的miRNA成熟体序列增至1999条。Miao等[89]对道赛特绵羊(Dorset)和小尾寒羊(Han)卵巢组织测序，鉴定出可能参与繁殖力调节的候选基因，这些候选基因参与各种细胞活动，如代谢级联、催化功能和信号转导。此外，通过miRNA谱分析鉴定了每组绵羊特有的特异性miRNA，发现若干与生殖力调控相关的miRNA。

Billerey等[90]检测了9头利木赞牛犊(Limousin)的胸肌样本，每个文库约获得14~45百万个配对末端读数，发现418种lincRNA (large intergenic non-coding RNAs)，与已知的10 775种蛋白编码基因存在显著差异。Sun等[91]利用Ribo-Zero RNA-seq技术深度剖析了胚胎、犊牛和成年牛骨骼肌的转录组谱，发现犊牛和成年牛之间的表达水平高度相关。在胚胎期有数百个基因显著表达，但在出生后至少减少了10倍，表明这些基因在肌肉发育中具有潜在作用。此外，该研究首次分析了牛骨骼肌中全部转录异构体，发掘出36 694个新型异构体，检测到185 036个SNP和12 428个短插入缺失(InDel)位点。研究发现可变剪接事件、SNP和InDel的数量在胚胎中比在犊牛和成年牛中更多，这表明基因表达在胚胎中最活跃。Cánovas等[92]通过转录组学测序对7头荷斯坦奶牛(Holstein cow)的乳样品进行了SNP筛选，共检测到19 175个差异表达基因，100 734个SNP位点，其中33 045个与荷斯坦奶牛SNP位点重合，这些SNP位于泌乳期间表达基因的编码区中，可用于荷斯坦奶牛乳用性状的基因变异分析和关联研究。

在梅花鹿()的育种研究中，Yao等[93]对鹿茸进行了转录组测序，组装出89 001个独特序列(平均大小450 bp)，发现了一些高度表达的基因参与鹿茸快速生长的调节，包括转录因子、信号分子和细胞外基质蛋白。这些数据是当前梅花鹿最全面的序列资源，为鹿的分子遗传学和功能基因组学的研究提供了理论基础。

3.4 全转录组测序在马属动物中的应用

全转录组研究在马属动物中开展得较少，目前主要是通过转录组学测序技术进行基因挖掘、注释和功能预测。Xie等[94]从头组装了驴()白细胞的转录组，鉴定出264 714个不同序列，预测了38 949个蛋白质片段。通过比较驴、马()和野马()的蛋白质序列，将驴蛋白片段与哺乳动物表型相关联。通过比较驴和马的外耳性状相关蛋白，鉴定出3种与耳形大小相关的蛋白质HIC1、PRKRA和KMT2A。Scott等[95]通过转录组测序对马的lncRNA进行注释，发现了20 800新型转录本，证明了lncRNA独有的特征，包括低表达、低外显子多样性和低水平的序列保守性。该研究结果所提供的候选基因可作为日后lncRNA注释的基准。

近年来发表了大量关于马运动机能、骨骼发育的文献报道，一部分是通过马组织的RNA-seq数据改进蛋白质编码基因的结构注释，另一部分是对RNA序列的分析。如Park等[96]对6匹纯种马运动前后的血液和肌肉进行全转录测序，通过与前人的研究对比，发现超过19 417个新型单基因簇，鉴定出189 973个单核苷酸位点变异(single nucleotide variants, SNV)。使用差异表达分析，确定了多个运动调节基因：血液中有62个上调基因和80个下调基因，肌肉中有878个上调基因和285个下调基因。结果表明，在差异表达的基因中有91个转录因子编码基因，其中包括56个功能未知的转录因子候选物可能与早期调节运动机制相关；此外，还发现了新型RNA表达模式：同一基因的不同可变剪接形式在运动前后表现出反向表达模式。该研究首次提供了马转录组数据和较为全面的分析结果，包括运动前后表达的基因，以及与运动相关的候选基因：6个运动相关基因和91个早期调节转录因子，3个高SNV密度的基因，以及4个交替表达的剪接体。

3.5 全转录组测序在其他特种经济动物中的应用

特种经济动物养殖已成为调整农村产业结构、发展特色经济的新亮点，为了发挥皮毛的最大经济价值，科研人员开始研究关于被毛颜色的调控机制。牛晓艳等[97]首先对不同毛色的獭兔(Rex rabbit)进行测序和差异基因分析，找到12 408个差异表达基因，然后通过KEGG分析将得到的差异基因聚类富集到相关代谢通路上，结果发现8个与黑素细胞分化相关的差异基因。宋兴超[98]利用RNA-seq技术对水貂()被毛色素沉积机理进行了研究，鉴定出不同时期被毛黑素含量的变化，并根据水貂皮肤组织中成熟黑素细胞的分布特点，开展SNPs检测，将 不同基因突变体与水貂毛色表型进行关联分析。mRNA定量表达验证结果表明，和等基因参与了黑素细胞发育、黑素小体前体形成、黑素小体转运和真黑和褐黑色素合成等生物学过程。

4 结语与展望

全转录组学以其精准、系统、直观的技术优势为畜禽重要经济性状功能基因的挖掘、鉴定与验证提供了新的技术平台和手段，并已广泛地运用在临床医学、药学、生物学、水产学和农林学等多个领域，为人类疾病研究、新药研发和动植物育种等开辟了新的思路。但是，全转录组学在畜牧领域的研究较其他领域而言起步较晚，研究还不够深入，尤其在羊上，转录组的研究目前还主要集中在小RNA测序和基因注释上。本研究团队将开展绵羊全转录组研究，并对其生长性状、肉用性状等重要经济性状进行解析和应用。

[1] Lockhart DJ, Winzeler EA. Genomics, gene expression and DNA arrays.,2000, 405(6788): 827–836.

[2] Lindberg J, Lundeberg J. The plasticity of the mammalian transcriptome.,2010, 95(1): 1–6.

[3] Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics.,2009, 10(1): 57–63.

[4] Wong ML, Medrano JF. Real-time PCR for mRNA quantitation.,2005, 39(1): 75–85.

[5] Liu C, LI X, Chen LL. Methods for genome-wide characterization of long noncoding RNAs.,2016, (6): 745–752.

[6] Yang L, Duff MO, Graveley BR, Carmichael GG, Chen LL. Genomewide characterization of non-polyadenylated RNAs.,2011, 12(2): R16.

[7] Carthew RW, Sontheimer EJ. Origins and mechanisms of miRNAs and siRNAs.,2009, 136(4): 642–655.

[8] Chen LL. Linking long noncoding RNA localization and function.,2016, 41(9): 761–772.

[9] Yin QF, Yang L, Zhang Y, Xiang JF, Wu YW, Carmichael G, Chen LL. Long noncoding RNAs with snoRNA ends.,2012, 48(2): 219–230.

[10] Zhang Y, Zhang X, Chen T, Xiang J, Yin Q, Xing Y, Zhu S, Yang L, Chen L. Circular intronic long noncoding RNAs.,2013, 51(6): 792–806.

[11] Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L. Complementary sequence-mediated exon circularization.,2014, 159(1): 134–147.

[12] Qi YX, Liu YB, Rong WH. RNA-Seq and its applications: a new technology for transcriptomics.,2011, 33(11): 1191–1201.祁云霞, 刘永斌, 荣威恒. 转录组研究新技术: RNA- Seq及其应用. 遗传, 2011, 33(11): 1191–1201.

[13] Hong QY, Bi XJ, Wang DN, Li ZZ, Yu H, Xia NS, Li SW. Research progress on RNA-Seq technology.,2017, 37(6): 443–448.洪奇阳, 毕行建, 王大宁, 李子真, 俞海, 夏宁邵, 李少伟. 转录组测序技术研究进展. 中国生化药物杂志, 2017, 37(6): 443–448.

[14] Patel RK, Jain M. NGS QC toolkit: a toolkit for quality control of next generation sequencing data.,2012, 7(2): e30619.

[15] Okonechnikov K, Conesa A, García-Alcalde F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data.,2016, 32(2): 292–294.

[16] Yang X, Liu D, Liu F, Wu J, Zou J, Xiao X, Zhao F, Zhu B. HTQC: a fast quality control toolkit for Illumina sequencing data.,2013, 14: 33.

[17] Zhou Q, Su X, Wang A, Xu J, Ning K. QC-Chain: fast and holistic quality control method for next-generation sequencing data.,2013, 8(4): e60234.

[18] Ward J, Cole C, Febrer M, Barton GJ. Almost significant: simplifying quality control of high-throughput sequencing data.,2016, 32(24): 3850–3851.

[19] Zhang M, Sun H, Fei Z, Zhan F, Gong X, Gao S. Fastq_clean: an optimized pipeline to clean the Illumina sequencing data with quality control. In: IEEE International Conference on Bioinformatics and Biomedicine. 2015, 44–48.

[20] Lo CC, Chain PSG. Rapid evaluation and quality control of next generation sequencing data with FaQCs.,2014, 15(1): 366.

[21] Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq.,2011, 8(6): 469–477.

[22] Yang IS, Kim S. Analysis of whole transcriptome sequencing data: workflow and software.,2015, 13(4): 119–125.

[23] Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality scores.,2008, 18(11): 1851–1858.

[24] Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform., 2009, 25(4): 1754–1760.

[25] Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.,2009, 10(3): R25.

[26] Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq.,2009, 25(9): 1105–1111.

[27] Wang K, Singh D, Zeng Z, Coleman SJ, Huang Y, Savich GL, He X, Mieczkowski P, Grimm SA, Perou CM, MacLeod JN, Chiang DY, Prins JF, Liu J. MapSplice: accurate mapping of RNA-seq reads for splice junction discovery.,2010, 38(18): e178.

[28] Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner.,2013, 29(1): 15–21.

[29] Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads.,2010, 26(7): 873–881.

[30] Kumar S, Blaxter ML. Comparing de novo assemblers for 454 transcriptome data.,2010, 11: 571.

[31] Garg R, Patel RK, Tyagi AK, Jain M. De novo assembly of chickpea transcriptome using short reads for gene discovery and marker identification., 2011, 18(1): 53–63.

[32] Margulies M, Egholm M, Altman WE, Attiya S, Bader JS, Bemben LA, Berka J, Braverman MS, Chen YJ, Chen Z. Genome sequencing in microfabricated high-density picolitre reactors.,2005, 437(7057): 376–380.

[33] Burlibasa C, Vasiliu D, Vasiliu M. Genome sequence assembly using trace signals and additional sequence information. In: German Conference on Bioinformatics. 1999, 45–56.

[34] Huang X, Madan A. CAP3: a DNA sequence assembly program.,1999, 9(9): 868–877.

[35] Swindell SR, Plasterer TN. SEQMAN. Contig assembly.,1997, 70: 75–89.

[36] Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, Tsai J, Quackenbush J.gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets.,2003, 19(5): 651–652.

[37] Miller RT, Christoffels AG, Gopalakrishnan C, Burke J, Ptitsyn AA. A comprehensive approach to clustering of expressed human gene sequence: the sequence tag alignment and consensus knowledge base.,1999, 9(11): 1143–1155.

[38] Zerbino DR, Birney E. Velvet: Algorithms for de novo short read assembly using de Bruijn graphs.,2008, 18(5): 821–829.

[39] Jackman SD, Birol İ. Assembling genomes using short- read sequencing technology.,2010, 11(1): 202–202.

[40] MacCallum I, Przybylski D, Gnerre S, Burton J, Shlyakhter I, Gnirke A, Malek J, McKernan K, Ranade S, Shea TP, Williams L, Young S, Nusbaum C, Jaffe DB. ALLPATHS 2: small genomes assembled accurately and with high continuity from short paired reads.,2009, 10(10): R103–R103.

[41] Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels.,2012, 28(8): 1086–1092.

[42] Li R, Li Y, Kristiansen K, Wang J, Wang J. SOAP: short oligonucleotide alignment program.,2008, 24(5): 713–714.

[43] Surget-Groba Y, Montoya-Burgos JI. Optimization of de novo transcriptome assembly from next-generation sequencing data.,2010, 20(10): 1432–1440.

[44] Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome.,2011, 29(7): 644–652.

[45] Huang X, Wang J, Aluru S, Yang SP, Hillier L. PCAP: a whole-genome assembly program.,2003, 13(9): 2164–2170.

[46] Huh JW, Kim YH, Park SJ, Kim DS, Lee SR, Kim KM, Jeong KJ, Kim JS, Song BS, Sim BW, Kim SU, Kim SH, Chang KT. Large-scale transcriptome sequencing and gene analyses in the crab-eating macaque () for biomedical research.,2012, 13: 163.

[47] Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads.,2015, 33(3): 290–295.

[48] Roberts A, Pimentel H, Trapnell C, Pachter L. Identification of novel transcripts in annotated genomes using RNA-Seq.,2011, 27(17): 2325–2329.

[49] Canzar S, Andreotti S, Weese D, Reinert K, Klau GW. CIDANE: comprehensive isoform discovery and abundance estimation.,2016, 17: 16.

[50] Boley N, Stoiber MH, Booth BW, Wan KH, Hoskins RA, Bickel PJ, Celniker SE, Brown JB. Genome-guided transcript assembly by integrative analysis of RNA sequence data.,2014, 32(4): 341–346.

[51] Liu JT, Yu T, Tao J, Li GJ. TransComb: genome-guided transcriptome assembly via combing junctions in splicing graphs.,2016, 17: 213.

[52] Mezlini AM, Smith EJ, Fiume M, Buske O, Savich GL, Shah S, Aparicio S, Chiang DY, Goldenberg A, Brudno M. iReckon: simultaneous isoform discovery and abundance estimation from RNA-seq data.,2013, 23(3): 519–529.

[53] Li JJ, Jiang CR, Brown JB, Huang H, Bickel PJ. Sparse linear modeling of next-generation mRNA sequencing (RNA-Seq) data for isoform discovery and abundance estimation.,2011, 108(50): 19867–19872.

[54] Hiller D, Wong WH. Simultaneous isoform discovery and quantification from RNA-seq.,2013, 5(1): 100–118.

[55] Stanke M, Keller O, Gunduz I, Hayes A, Waack S, Morgenstern B. AUGUSTUS: ab initio prediction of alternative transcripts.,2006, 34(web server issue): 435–439.

[56] Li W, Feng J, Jiang T. IsoLasso: A LASSO regression approach to RNA-seq based transcriptome assembly.,2011, 18(11): 1693–1707.

[57] Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Xian A, Fan L, Koziol MJ, Gnirke A, Nusbaum C. Ab initio reconstruction of cell type–specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs.,2010, 28(5): 503–510.

[58] Tomescu, Alexandru I, Kuosmanen, Anna, Makinen, Veli, Rizzi R. A novel min-cost flow method for estimating transcript expression with RNA-Seq.,2013, 14(Suppl.5): S15.

[59] Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq.,2008, 5(7): 621–628.

[60] Pachter L. Models for transcript quantification from RNA-Seq. 2013.

[61] Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data.,2015, 31(2): 166–169.

[62] Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.,2014, 30(7): 923–930.

[63] Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.,,2011, 12: 323.

[64] Rob P, Mount SM, Kingsford C. Sailfish enables alignment-free isoform quantification from RNA-seq reads using lightweight algorithms.,2013, 32(5): 462–464.

[65] Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification.,2016, 34(5): 525–527.

[66] Bullard JH, Purdom E , Hansen KD, Dudoit S. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments.,2010, 11: 94.

[67] Ma X, Zhang X. NURD: an implementation of a new method to estimate isoform expression from non-uniform RNA-seq data.,2013, 14: 220.

[68] Oshlack A, Wakefield MJ. Transcript length bias in RNA-seq data confounds systems biology.2009, 4: 14.

[69] Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data.,2010, 11(3): 1–9.

[70] Anders S, Huber W. Differential expression analysis for sequence count data.,2010, 11: R106.

[71] Li J, Witten DM, Johnstone IM, Tibshirani R. Normalization, testing, and false discovery rate estimation for RNA- sequencing data.,2012, 13(3): 523–538.

[72] Rosenbaum PR, Rubin DB. Reducing bias in observational studies using subclassification on the propensity score.,1984, 79(387): 516–524.

[73] Law CW, Chen Y, Wei S, Smyth GK. voom: precision weights unlock linear model analysis tools for RNA-seq read counts.,2014, 15(2): R29.

[74] Li J, Tibshirani R. Finding consistent patterns: A nonparametric approach for identifying differential expression in RNA-Seq data., 2011, 22(5): 519–536.

[75] Mcdaneld TG, Smith TP, Doumit ME, Miles JR, Coutinho LL, Sonstegard TS, Matukumalli LK, Dan JN, Wiedmann RT. MicroRNA transcriptome profiles during swine skeletal muscle development.,2009, 10: 77.

[76] Shen YF. Sequencing and characterization of mRNA, lncRNA, and miRNA in thytoid gland of Yorkshire and Jinhua Pigs[D]. Zhejiang University, 2016.沈一飞. 约克夏猪和金华猪甲状腺组织mRNA、lncRNA和miRNA测序及功能分析[学位论文]. 浙江大学, 2016.

[77] Li M, Xia Y, Gu Y, Zhang K, Lang Q, Chen L, Guan J, Luo Z, Chen H, Li Y, Li Q, Li X, Jiang AA, Shuai S, Wang J, Zhu Q, Zhou X, Gao X, Li X. MicroRNAome of porcine pre- and postnatal development.,2010, 5(7): e11541.

[78] Sun J, Xie M, Huang Z, Li H, Chen T, Sun R, Wang J, Xi Q, Wu T, Zhang Y. Integrated analysis of non-coding RNA and mRNA expression profiles of 2 pig breeds differing in muscle traits.,2017, 95(3): 1092–1103.

[79] Jianning YU, Yan L, Chen Z, Hui LI, Ying S, Zhu H, Shi Z. Investigating right ovary degeneration in chick embryos by transcriptome sequencing., 2017, 63(3): 295–303.

[80] Li T, Wang S, Wu R, Zhou X, Zhu D, Zhang Y. Identification of long non-protein coding RNAs in chicken skeletal muscle using next generation sequencing.,2012, 99(5): 292–298.

[81] Glazov EA, Cottee PA, Barris WC, Moore RJ, Dalrymple BP, Tizard ML. A microRNA catalog of the developing chicken embryo identified by a deep sequencing approach., 18(6): 957–964.

[82] Kang B, Guo JR, Yang HM, Zhou RJ, Liu JX, Li SZ, Dong CY. Differential expression profiling of ovarian genes in prelaying and laying geese.,2009, 88(9): 1975–1983.

[83] Guo J, Tang QP, Zhang SJ, Ma YH, Lu HL. Identification of broodiness-related geese genes by suppression subtractive hybridization.,2011, 42(10): 1477–1784.郭军, 汤青萍, 章双杰, 马月辉, 陆火林, 苏建东, 邹剑敏, 陈宽维, 李慧芳. 利用抑制消减杂交技术筛选鹅就巢行为相关基因. 畜牧兽医学报,2011, 42(10): 1477– 1484.

[84] Xu Q, Zhao WM, Chen Y, Tong YY, Rong GH, Huang ZY, Zhang Y, Chang GB, Wu XS, Chen GH. Transcriptome profiling of the goose () ovaries identify laying and broodiness phenotypes.,2013, 8(2): e55496.

[85] Chen L, Luo J, Li JX, Li JJ, Wang DQ, Tian Y, Lu LZ. Transcriptome analysis of adiposity in domestic ducks by transcriptomic comparison with their wild counterparts.,2015, 46(3): 299–307.

[86] Di R, He J, Song S, Tian D, Liu Q, Liang X, Ma Q, Sun M, Wang J, Zhao W, Cao G, Wang J, Yang Z, Ge Y, Chu M. Characterization and comparative profiling of ovarian microRNAs during ovine anestrus and the breeding season.,2014, 15: 899.

[87] Chang W, Wang J, Tao D, Zhang Y, Jianzhong HE, Shi C. Identification of a novel miRNA from the ovine ovary by a combinatorial approach of bioinformatics and experiments.,2015, 77(12): 1617–1624.

[88] Zhang SF, Wei CH, Lu J, Zhang XN, Zhou XL, Zhang SZ, Wang GK, Cao JX, Zhao FP, Zhang L, Du LX. Identification of the microRNAome in texel sheep by deep sequencing., 2013, 40(9): 19–22.张世芳, 魏彩虹, 陆健, 张小宁, 周鑫磊, 张淑珍, 王光凯, 曹家雪, 赵福平, 张莉, 杜立新. 深度测序鉴定绵羊microRNA转录组. 中国畜牧兽医, 2013, 40(9): 19– 22.

[89] Miao X, Qin QL. Genome-wide transcriptome analysis of mRNAs and microRNAs in Dorset and Small Tail Han sheep to explore the regulation of fecundity.,2015, 402: 32–42.

[90] Billerey C, Boussaha M, Esquerré D, Rebours E, Djari A, Meersseman C, Klopp C, Gautheret D, Rocha D. Identification of large intergenic non-coding RNAs in bovine muscle using next-generation transcriptomic sequencing.,2014, 15: 499.

[91] Sun X, Li M, Sun Y, Cai H, Li R, Wei X, Lan X, Huang Y, Lei C, Chen H. The developmental transcriptome landscape of bovine skeletal muscle defined by Ribo-Zero ribonucleic acid sequencing.2015, 93(12): 5648–5658.

[92] Cánovas A, Rincon G, IslasTrejo A, Wickramasinghe S, Medrano JF. SNP discovery in the bovine milk transcriptome using RNA-Seq technology.,2010, 21(11–12): 592–598.

[93] Yao B, Yu Z, Mei Z, Liu M, Liu H, Li J. De novo characterization of the antler tip of Chinese Sika deer transcriptome and analysis of gene expression related to rapid growth.,2012, 364(1–2): 93–100.

[94] Xie FY, Feng YL, Wang HH, Ma YF, Yang Y, Wang YC, Shen W, Pan QJ, Yin S, Sun YJ, Ma JY.assembly of the donkey white blood cell transcriptome and a comparative analysis of phenotype-associated genes between donkeys and horses.,2015, 10(7): e0133258.

[95] Scott EY, Mansour T, Bellone RR, Brown CT, Mienaltowski MJ, Penedo MC, Ross PJ, Valberg SJ, Murray JD, Finno CJ. Identification of long non-coding RNA in the horse transcriptome.,2017, 18(1): 511.

[96] Park KD, Park J, Ko J, Kim BC, Kim HS, Ahn K, Do KT, Choi H, Kim HM, Song S, Lee S, Jho S, Kong HS, Yang YM, Jhun BH, Kim C, Kim TH, Hwang S, Bhak J, Lee HK, Cho BW. Whole transcriptome analyses of six thoroughbred horses before and after exercise using RNA-Seq.,2012, 13: 473.

[97] Niu XY, Ren KL, Cao L, Li YP, Zhen JJ, Feng GL, Huang SF. Transcriptome analysis on coat color related genes in rex rabbits.,2016, 36(2): 1–6.牛晓艳, 任克良, 曹亮, 李燕平, 郑建婷, 冯国亮, 黄淑芳. 利用转录组测序方法研究獭兔毛色相关基因. 中国草食动物科学, 2016, 36(2): 1–6.

[98] Song XC. Study on the mechanisms of hair pigmentation and skin transcriptome annotation in mink (neovison vison) based onhigh throughput RNA-sequencing[D]. Chinese Acadamy of Agricultural Sciences, 2016. 宋兴超. 水貂被毛色素沉积机理及基于高通量RNA- seq皮肤转录组注释研究[学位论文]. 中国农业科学院, 2016.

Application of whole transcriptomics in animal husbandry

Tianpei Shi, Li Zhang

RNA is a polymeric molecule which is involved in various biological processes including the coding, decoding, regulation, and expression of genes. Whole transcriptome sequencing is the dominant method for studying RNA functions which assists researchers to deepen the exploration and analysis of gene structure and function and to reveal intrinsic links between gene expression and life phenomena. To date, extensive research has been done in animal husbandry models including swine and poultry by using high-throughput RNA sequencing technology. These studies included transcript structure, gene expression level and non-coding region function. In this review, we briefly introduce the general processes and strategies of RNA-seq and small RNA-seq technologies, and summarize the various achievements of the application of whole transcriptomics in animal husbandry.

gene expression; whole transcriptome; RNA-seq; high-throughput RNA sequencing

2018-07-31;

2018-10-22

国家自然科学基金项目(编号：U1503285) 和中国农业科学院基本科研业务费重大项目(编号：Y2017XM02)资助[Supported by the National Natural Science Foundation of China(No. U1503285) and Chinese Academy of Agricultural Sciences(No.Y2017XM02 )]

石田培，在读硕士研究生，专业方向：动物遗传育种与繁殖。E-mail: 1337684764@qq.com

张莉，博士，研究员，博士生导师，研究方向：动物遗传育种与繁殖。E-mail: zhangli07@caas.cn

10.16288/j.yczz.18-218

2018/11/15 9:28:51

URI: http://kns.cnki.net/kcms/detail/11.1913.R.20181115.0928.002.html

(责任编委: 赵要风)