利用转录组测序方法研究獭兔毛色相关基因

牛晓艳，任克良，曹 亮，李燕平，郑建婷，冯国亮，黄淑芳

（山西省农业科学院畜牧兽医研究所，太原030032）

遗传育种

利用转录组测序方法研究獭兔毛色相关基因

牛晓艳，任克良，曹 亮，李燕平，郑建婷，冯国亮，黄淑芳

（山西省农业科学院畜牧兽医研究所，太原030032）

利用转录组测序方法分析了不同毛色（白色和海狸色）獭兔皮肤组织中基因的表达差异，旨在找到影响獭兔毛色的候选基因。采用Illumina HiSeqTM2500测序平台对不同毛色獭兔耳组织中所有mRNA进行高通量测序，所得序列经质控、组装后比对到GO、COG、SWISS-PROT数据库中注释，并进行差异表达基因聚类分析和KEGG通路分析。结果表明：测序共获得Unigene 275 625个，差异表达基因12 408个，其中上调表达基因4 904个，下调表达基因7 504个。KEGG通路分析发现，这些差异基因显著地富集在13个代谢通路中，并在细胞色素代谢通路（metabolism of xenobiotics by cytochrome P450）中有8个差异表达的基因，其中CP2F1和CP2BB属于细胞色素P450家族成员，而GSTA4与黑素细胞分化有关，可能在獭兔毛色中发挥作用。

转录组测序；獭兔；毛色；差异表达；GO和KEGG通路分析

獭兔，学名力克斯兔（Rexrabbit），是世界著名的皮用兔品种，因其皮毛具有“短、平、密、细、美、牢”等特点，制作的服饰色彩斑斓、轻柔飘逸而受到人们的喜爱。獭兔毛色众多，已发展出海狸色、蓝色、青紫蓝、红色等二十多种色型，九十多种颜色［1］。中国是养兔大国，年产獭兔皮张和出口量均居世界首位。从世界獭兔目前养殖的色型种类看，主要为白色獭兔。白色皮张可染色制成许多色泽不同的皮毛制品，因此，市场前景十分广阔。目前，国内外市场对獭兔毛皮质量的要求日益提高，而生产中存在的品种退化以及毛色分离现象在一定程度上影响了毛皮的品质。因此，研究影响毛色的相关基因具有一定意义。

哺乳动物的毛色是由体内色素所决定的，主要是酪氨酸源性色素，即黑色素及其衍生物。黑色素主要位于表皮和毛囊中，包括真黑色素和褐黑色素。真黑色素可使皮肤或毛发表现为褐色和黑色，褐黑色素可使皮肤或毛发表现为黄色和红色［2］。目前，已有300多个遗传座位和150多个毛色相关基因被识别，这些基因以不同的路径影响着色素的形成、转运等过程［3］。其中主效基因包括MITF、KIT及KIT配体、TYR、TYRP1、TYRP2、MC1R、ASIP、MLPH、Agouti等。在家兔中，现已证实至少有8个基因座位（或系统）控制其毛色［2］。除基因间的显隐性及互作关系外，环境因素也对动物毛色具有一定的影响，如环境中微量元素的含量、光照、温度、海拔甚至饲养方式等都会影响黑色素的表达，使毛色产生变化。

转录组测序是利用大规模测序技术直接对cDNA序列进行测序，研究动物特定组织或特定时期的所有mRNA的组分，可找到生物体不同时期、不同组织或不同个体间mRNA的表达差异，从而得到基因表达、可变剪切、基因结构优化、新基因发现等分析结果［4-5］。RNA-seq具有通量高、成本低、灵敏度高、无需基因组序列信息等优势，因此在动植物研究中应用较广。

本研究利用Illumina HiSeqTM2500测序平台对不同毛色獭兔的基因表达情况进行分析，找到差异表达的基因，并对其进行功能注释和KEGG通路分析，以期揭示与獭兔毛色相关的候选基因。

1 材料与方法

1.1 材料

动物：獭兔来自山西省农科院畜牧所实验兔场，选择全同胞的白色和海狸色獭兔各1只，日龄均120 d。取大小为0.5 cm3的耳组织样，分装后迅速放入液氮中速冻。

主要试剂：提取总RNA所用Trizol购自Invitrogen公司。

1.2 方法

1.2.1 文库构建与高通量测序 采用Trizol法提取獭兔耳皮肤组织样的总RNA，再用1%的琼脂糖电泳检测总RNA是否降解，最后用安捷伦2100 RNA Nano 6000 AssayKit（安捷伦公司）检测RNA样品的完整性和浓度。

总RNA样本检测合格后，经过mRNA的富集、片段化、反转录、末端加接头后构建成cDNA文库，文库插入片段大小平均250 bp。随后用Illumina HiSeqTM2500进行测序，测序读长为双端125 bp。

1.2.2 Unigene组装与功能注释 HiSeqTM测序所得的原始序列通过去除低质量序列、去接头等过程得到高质量序列。质控后，通过Trinity（20140717）对无参转录组进行组装［6］。采用TransDeorder（20140717）对组装序列的开放阅读框（ORF）进行预测。采用Trinotate（20140717）对组装序列进行注释，参考数据库为Uniprot、Swiss-prot、GO、KEGG等。利用Blast2GO（Gene Ontology）（http：//www. blast2go.com）计算不同term的基因数目，并根据GO聚类的统计结果绘制相应的柱状图。通过与COG数据库（cluster of orthologous group）比对，对未知基因进行直系同源基因的功能注释［7］。

2 结果与分析

2.1 Unigene的组装

深度测序共获得潜在的Unigene 275 625个，潜在的转录本313 343个，转录本长度在201～1 6191 nt之间。其中长度为200～400 nt和400～600 nt的转录本分别为212 889个（占67.94%）和39 109个（占12.48%），长度在4 000 nt以上的转录本为4 351个（占1.39%），N50和N90值分别为958 nt和234 nt。转录本长度的分布情况见图1。

将用于组装的转录本与组装后的转录本用Bowtie2（2.2.3版）进行比对［10］，其中比对上序列的比例占到87.02%，均一性分析结果良好，符合质量要求。

图1 Illumina HiSeqTM测序获得的转录本长度分布

2.2 Unigene的功能注释

2.2.1 Unigene的GO注释 对所获得的275 625个Unigene进行功能注释，对注释到SWISS-PROT、GO、KEGG数据库的基因进行统计，结果如表1所示。通过GO（gene ontology，GO）注释，可确定这些基因的功能分类，结果如图2所示。从比对结果来看，仍有91.06%的基因未比对到任何数据库，这可能与目前公布的家兔基因组信息尚不完全有关。

图2中，Blast2GO软件分析得到每个Unigene蛋白GO功能注释，将所有Unigene比对到GO的3个数据库：生物过程（biological process，BP）、分子功能（molecular function，MF）和细胞组成（cellular component，CC），并做功能分类统计。其中绝大多数基因都比对到细胞过程（cellular process）、代谢过程（metabolic process）、细胞组分（cell part）以及（分子）联接（binding）等类别中。

表1 Unigene的比对结果

图2 GO统计柱状图

2.2.2 Unigene的COG注释 用BLAST软件将组装序列与COG数据库比对，对结果进行过滤后去掉e＞10-5的比对序列，再进行注释后结果如图3所示。共对10 288个Unigene进行注释，分为25个类别。其中有1 951个基因比对到一般机能（general function prediction）类，有1 084个基因比对到转录后修饰，蛋白折叠和分子伴侣 （posttranslational modification，protein turnover and chaperones）类，只有14个基因比对到细胞动能（cellular motility）类，5个基因比对到细胞核结构（nuclear structure）类，另外还有477个基因功能未知（function unknown），可能是测序发现的新基因。

2.3 不同毛色獭兔差异表达基因的分析

2.3.1 差异表达基因的GO注释 将不同毛色的2个样品所获得的转录本用DEGseq软件进行比较，最终获得差异表达基因12 408个，其中上调表达的基因4 904个，下调表达的基因7 504个。对所获得的差异表达基因进行统计，并利用Uniprot、Pfam、GO和KEGG等数据库对差异基因进行功能注释，获得详细描述信息。根据GO数据库中的相关条目，统计得到差异表达基因，如图4所示。

GO统计结果显示：所有差异基因比对到BP数据库中的有2 711个，其中显著富集在免疫反应（immune response）的基因为154个，占比对到该功能类别总数的5.68%；显著富集在免疫过程（immune system process）的基因为196个，占比对到该功能类别总数的7.22%；显著富集在防御反应（defense response）中的基因127个，占总数的4.68%。所有差异基因比对到MF数据库中的有622个，其中显著富集在RNA介导的DNA酶活性（RNA-directed DNA polymerase activity）中的121个，占比对到MF数据库总数的19.45%；显著富集在DNA酶活性（DNA polymerase activity）中的121个，占总数的19.45%。所有差异表达基因比对到CC数据库中的有338个，其中显著富集在胞外区域（extracellular region）的基因为202个，占比对到该功能类别总数的59.76%。

图3 GOG注释统计图

图4 差异表达基因的GO统计柱状图

由图4可以看出，在BP数据库所有上调表达的基因中，共有443个基因归入细胞过程（cellular process）中，占总数的 20.18%，331个基因归入代谢过程（metabolic process）中，占总数的15.08%；所有下调表达的基因中，共有541个基因归入细胞过程中，占总数的20.76 %，440个基因归入代谢过程中，占总数的16.88%。在MF数据库所有上调表达的基因中，共有410个基因归入细胞连接（binding）中，占总数的44.46%；所有下调表达的基因中，共有554个基因归入细胞连接（binding）中，占总数的46.59%；在CC数据库所有上调表达的基因中，共有519个基因归入细胞组分（cell part）中，占总数的31.86%，共有656个基因归入细胞组分（cell part）中，占总数的31.91%。

2.3.2 差异表达基因的KEGG通路分析 用KEGG注释系统将测序所发现的所有转录本注释到KEGG数据库中，共对该数据库中315个代谢通路进行比对，差异表达基因显著富集的通路13个。其中3个与脂类代谢有关，1个与RNA转运有关，2个与糖类代谢有关，1个与氨基酸代谢有关，1个与腮腺炎有关，1个是细胞色素P450合成的通路。表2是差异表达基因显著富集通路的汇总。

表2 不同颜色獭兔耳组织中差异表达基因显著性富集的通路

对细胞色素代谢通路（Map00980）中表达发生变化的基因进行进一步分析，共发现8个表达发生变化的基因，其中上调表达的基因有：GSTM4，ARK72；下调表达的基因有：CP2F1，UD11，ADH6，GSTA4，CP2BB，ST2A1。其中GSTM4（谷胱甘肽转移酶）在多个细胞色素代谢通路中都出现上调表达，而CP2F1和CP2BB都是细胞色素P450家族成员，但它们的具体生物学功能尚未有文献报道。

3 讨论

3.1 獭兔转录组测序

高通量测序的样本可选用组织，也可选用细胞或细胞系，这主要是由不同的试验目的所决定的。黑色素主要位于表皮和毛囊中，所以选择獭兔的耳皮肤组织。已有研究者选用cDNA芯片分析了不同被毛密度獭兔皮肤组织中差异表达的基因［11］。Pan L等［12］采用Solexa测序技术研究了獭兔不同表型的皮肤组织（plaice skin and un-plaice skin）基因表达的差异，并对一个影响表型的关键基因Lamb3内的SNPs进行分析。翁巧琴等［13］研究了MC1R基因与獭兔毛色的关系，最终发现在MC1R基因第40、50位点存在T＞C和A＞G的突变，在第45位点缺失一段长度为23bp的碱基序列，第507位点存在T＞C的碱基突变，这些突变可能与獭兔毛色有关。牛晓艳等［14］研究了獭兔MITF基因部分外显子的多态性，并在外显子1、5、7中发现一些多态与毛色存在强关联。本研究拟通过獭兔耳组织的RNA-Seq，初步筛选一些与毛色相关的新基因。

3.2 样本含量

RNA-Seq的生物学重复和技术重复的成本都比较高，如果采用多个样本的RNA等量混池进行转录组测序，可能由于样本间的个体差异和遗传背景的差异等导致最后结果的误差较大。又因为RNA测序价格的限制，综合考虑后，决定选用全同胞的海狸色獭兔和白色獭兔各1只，并充分考虑到遗传背景的一致性，在数据分析过程中提高了差异显著性阈值，缩小差异基因的范围，再结合差异基因的生物学功能和KEGG通路信息，最大程度地保证了结果的准确性和可参考性。

3.3 细胞色素代谢通路中差异表达基因

由表2可以看出，在13个差异基因显著富集的通路中，只有Map编号为00980的细胞色素P450代谢通路与细胞色素代谢密切相关。该通路中有8个基因表达发生了变化，其中CP2F1（细胞色素P450家族成员2F1）和CP2BB（细胞色素P450家族成员2B11）属于细胞色素家族成员；而UD11（UDP-葡糖甘酸转移酶）、ADH6（乙醇脱氢酶-6）、ARK72（黄曲霉毒素乙醛还原酶家族成员2）、GSTM4（谷胱甘肽S-转移酶）、ST2A1（类固醇磺基转移酶）与毛色的关系尚未见文献报道。Uehara等［15］2009年报道了GSTA4（谷胱甘肽S-转移酶）基因与小鼠耳蜗皮肤组织黑素细胞迁移有关。Christopher等［16］2013年总结了影响猫和狗毛色的基因座，其中包括：MITF、PMEL、TYRP1、MLPH、MC1R、ASIP、KIT等经典基因以及PSMB7、CBD103、TAQPEP等，但未包括本试验中发现的基因。这可能与本试验设计的局限以及样本含量较少有关，但也可能是发现了一些影响獭兔毛色变化的新基因，这些基因的功能以及对毛色影响的具体方式还需要进一步的试验研究。

4 结论

本研究通过转录组测序（RNA-Seq）方法对不同毛色獭兔耳组织中差异表达的基因进行了研究。共发现12 408个差异表达基因，这些基因主要与代谢过程、细胞组分、催化反应有关。通过KEGG通路分析，最终在与细胞色素代谢有关的通路中发现8个差异表达的基因，对这些基因的功能进行进一步分析发现，只有GSTA4与细胞色素迁移有关，可能最终影响到獭兔的毛色。这些候选基因与毛色的关系还需要进一步进行功能验证，最终找到相关分子标记，用于兔群的筛选以及毛色机理的研究。

［1］ 庞有志.为长毛兔和獭兔品种论不平——兼评家兔品种与品系的分类［J］.中国养兔杂志，2001（4）：30-34.

［2］ 庞有志.兔的毛色遗传与家兔育种［J］.中国养兔杂志，2008（10）：30-35.

［3］ 杨广礼.动物毛色形成的调控机制研究［J］.黑龙江畜牧兽医，2014，3（上）：45-48.

［4］ 周华，张新，刘腾云，等.高通量转录组测序的数据分析与基因挖掘［J］.江西科学，2012，30（5）：607-611.

［5］ 岳桂东，高强，罗龙海，等.高通量测序技术在动植物研究领域中的应用［J］.中国科学：生命科学，2012，42（2）：107-124.

［6］ Manfred G.Full-length transcriptome assembly from RNA-Seqdata without a reference genome［J］.Nature Biotechnology，2011，29：644-652.

［7］ Tatusov.A genomic perspective on protein families［J］.Nucleic Acids Research，2000，28：33-36.

［8］ Mortazavi A，Williams B A，McCue K，et al.Mapping and quantifying mammalian transcriptomes by RNA-Seq［J］.Nature Methods，2000，5（7）：621-628.

［9］ WangL，FengZ，WangX，et al.DEGseq：an R package for identifying differentially expressed genes from RNA-Seq data［J］.Bioinformatics，2010，26（1）：136-138.

［10］Langmead B，Trapnell C，Pop M，et al.Ultrafast and memory-efficient alignment of short DNA sequences to the human genome［J］.Genome Biology，2009，10（3）：25-34.

［11］Chen S J，Liu Y J，Gu Z L，et al.Microarray analysis gene expression profiles in Rex rabbit skin［C］//Proceedings of the 10th World Rabbit Congress，Egypt，2012：43-38.

［12］Pan L，Liu Y，Wu X S，et al.Solexa-sequencing based transcriptome study of plaice skin phenotype in Rex rabbits（Oryctolagus cuniculus）［J］.Plos One，2015，10（5）：1-14.

［13］翁巧琴，王进波，王国朋，等.特种獭兔毛色相关基因多态性的研究［J］.中国养兔杂志，2015（3）：4-9.

［14］牛晓艳，任克良，曹亮，等.獭兔MITF基因部分外显子多态性研究［J］.山西农业科学，2013，41（11）：1231-1238.

［15］Uehare S，Izumi Y，KuboY，et al.Specific expression of Gsta4 in mousecochlear melanocytes：novel role for hearing and melanocyte differentiation［J］.PigmentCellMelanomaResearch，2009，22（1）：111-119.

［16］Christopher BK，GregorySB.Genetics ofpigmentation in dogs and cats［J］.Annual ReviewofAnimal Bioscience，2013，1：125-156.

Transcriptome Analysis on Coat Color Related Genes in Rex Rabbits

Niu Xiaoyan，Ren Keliang，CaoLiang，et al
（Institute ofAnimal Husbandry&Veterinary，Shanxi AcademyofAgricultural Science，Taiyuan 030032，China）

The genes differentially expressed in Rex rabbit skin with different colors（white and beaver）were analyzed via transcriptome sequencing in order to reveal candidate genes related to the coat color.The mRNA in the ear tissue was sequenced and analyzed using Illumina HiSeqTM2500 sequencing platform.After quality control and assemble，the sequences acquired were blasted against GO，COG，SWISS-PROT databases，and then cluster analysis and KEGG pathway analysis were performed.In total，275 625 unigenes were acquired，and 12 408 genes were differentially expressed，among which 4 904 were up-regulated and 7 504 down-regulated.KEGGpathwayanalysis revealed that the differentially expressed genes were significantly enriched in 13 signal pathways（P＜0.05）.Eight genes were identified in the pathway of xenobiotics by cytochrome P450，among which CP2F1 and CP2BB belonged tothe cytochrome P450 family，and GSTA4 was related with melanin differentiation.Theymight play important roles in coat color ofRexrabbit.

transcriptome sequencing；Rexrabbit；coat color；differentiallyexpression；GOand KEGGanalysis

S829.2

A

2095-3887（2016）02-0001-07

10.3969/j.issn.2095-3887.2016.02.001

2016-01-26

国家兔产业技术体系项目（CARS-44-B-6）；山西省农科院博士基金项目（YBSJJ1401）

牛晓艳（1984-），女，助理研究员，博士，主要从事海狸色獭兔育种工作。

任克良，Email：keliangren@sohu.com