利用轉(zhuǎn)錄組測序方法研究獺兔毛色相關(guān)基因

牛曉艷，任克良，曹 亮，李燕平，鄭建婷，馮國亮，黃淑芳

（山西省農(nóng)業(yè)科學院畜牧獸醫(yī)研究所，太原030032）

遺傳育種

利用轉(zhuǎn)錄組測序方法研究獺兔毛色相關(guān)基因

牛曉艷，任克良，曹 亮，李燕平，鄭建婷，馮國亮，黃淑芳

（山西省農(nóng)業(yè)科學院畜牧獸醫(yī)研究所，太原030032）

利用轉(zhuǎn)錄組測序方法分析了不同毛色（白色和海貍色）獺兔皮膚組織中基因的表達差異，旨在找到影響獺兔毛色的候選基因。采用Illumina HiSeqTM2500測序平臺對不同毛色獺兔耳組織中所有mRNA進行高通量測序，所得序列經(jīng)質(zhì)控、組裝后比對到GO、COG、SWISS-PROT數(shù)據(jù)庫中注釋，并進行差異表達基因聚類分析和KEGG通路分析。結(jié)果表明：測序共獲得Unigene 275 625個，差異表達基因12 408個，其中上調(diào)表達基因4 904個，下調(diào)表達基因7 504個。KEGG通路分析發(fā)現(xiàn)，這些差異基因顯著地富集在13個代謝通路中，并在細胞色素代謝通路（metabolism of xenobiotics by cytochrome P450）中有8個差異表達的基因，其中CP2F1和CP2BB屬于細胞色素P450家族成員，而GSTA4與黑素細胞分化有關(guān)，可能在獺兔毛色中發(fā)揮作用。

轉(zhuǎn)錄組測序；獺兔；毛色；差異表達；GO和KEGG通路分析

獺兔，學名力克斯兔（Rexrabbit），是世界著名的皮用兔品種，因其皮毛具有“短、平、密、細、美、牢”等特點，制作的服飾色彩斑斕、輕柔飄逸而受到人們的喜愛。獺兔毛色眾多，已發(fā)展出海貍色、藍色、青紫藍、紅色等二十多種色型，九十多種顏色［1］。中國是養(yǎng)兔大國，年產(chǎn)獺兔皮張和出口量均居世界首位。從世界獺兔目前養(yǎng)殖的色型種類看，主要為白色獺兔。白色皮張可染色制成許多色澤不同的皮毛制品，因此，市場前景十分廣闊。目前，國內(nèi)外市場對獺兔毛皮質(zhì)量的要求日益提高，而生產(chǎn)中存在的品種退化以及毛色分離現(xiàn)象在一定程度上影響了毛皮的品質(zhì)。因此，研究影響毛色的相關(guān)基因具有一定意義。

哺乳動物的毛色是由體內(nèi)色素所決定的，主要是酪氨酸源性色素，即黑色素及其衍生物。黑色素主要位于表皮和毛囊中，包括真黑色素和褐黑色素。真黑色素可使皮膚或毛發(fā)表現(xiàn)為褐色和黑色，褐黑色素可使皮膚或毛發(fā)表現(xiàn)為黃色和紅色［2］。目前，已有300多個遺傳座位和150多個毛色相關(guān)基因被識別，這些基因以不同的路徑影響著色素的形成、轉(zhuǎn)運等過程［3］。其中主效基因包括MITF、KIT及KIT配體、TYR、TYRP1、TYRP2、MC1R、ASIP、MLPH、Agouti等。在家兔中，現(xiàn)已證實至少有8個基因座位（或系統(tǒng)）控制其毛色［2］。除基因間的顯隱性及互作關(guān)系外，環(huán)境因素也對動物毛色具有一定的影響，如環(huán)境中微量元素的含量、光照、溫度、海拔甚至飼養(yǎng)方式等都會影響黑色素的表達，使毛色產(chǎn)生變化。

轉(zhuǎn)錄組測序是利用大規(guī)模測序技術(shù)直接對cDNA序列進行測序，研究動物特定組織或特定時期的所有mRNA的組分，可找到生物體不同時期、不同組織或不同個體間mRNA的表達差異，從而得到基因表達、可變剪切、基因結(jié)構(gòu)優(yōu)化、新基因發(fā)現(xiàn)等分析結(jié)果［4-5］。RNA-seq具有通量高、成本低、靈敏度高、無需基因組序列信息等優(yōu)勢，因此在動植物研究中應(yīng)用較廣。

本研究利用Illumina HiSeqTM2500測序平臺對不同毛色獺兔的基因表達情況進行分析，找到差異表達的基因，并對其進行功能注釋和KEGG通路分析，以期揭示與獺兔毛色相關(guān)的候選基因。

1 材料與方法

1.1 材料

動物：獺兔來自山西省農(nóng)科院畜牧所實驗兔場，選擇全同胞的白色和海貍色獺兔各1只，日齡均120 d。取大小為0.5 cm3的耳組織樣，分裝后迅速放入液氮中速凍。

主要試劑：提取總RNA所用Trizol購自Invitrogen公司。

1.2 方法

1.2.1 文庫構(gòu)建與高通量測序 采用Trizol法提取獺兔耳皮膚組織樣的總RNA，再用1%的瓊脂糖電泳檢測總RNA是否降解，最后用安捷倫2100 RNA Nano 6000 AssayKit（安捷倫公司）檢測RNA樣品的完整性和濃度。

總RNA樣本檢測合格后，經(jīng)過mRNA的富集、片段化、反轉(zhuǎn)錄、末端加接頭后構(gòu)建成cDNA文庫，文庫插入片段大小平均250 bp。隨后用Illumina HiSeqTM2500進行測序，測序讀長為雙端125 bp。

1.2.2 Unigene組裝與功能注釋 HiSeqTM測序所得的原始序列通過去除低質(zhì)量序列、去接頭等過程得到高質(zhì)量序列。質(zhì)控后，通過Trinity（20140717）對無參轉(zhuǎn)錄組進行組裝［6］。采用TransDeorder（20140717）對組裝序列的開放閱讀框（ORF）進行預(yù)測。采用Trinotate（20140717）對組裝序列進行注釋，參考數(shù)據(jù)庫為Uniprot、Swiss-prot、GO、KEGG等。利用Blast2GO（Gene Ontology）（http：//www. blast2go.com）計算不同term的基因數(shù)目，并根據(jù)GO聚類的統(tǒng)計結(jié)果繪制相應(yīng)的柱狀圖。通過與COG數(shù)據(jù)庫（cluster of orthologous group）比對，對未知基因進行直系同源基因的功能注釋［7］。

2 結(jié)果與分析

2.1 Unigene的組裝

深度測序共獲得潛在的Unigene 275 625個，潛在的轉(zhuǎn)錄本313 343個，轉(zhuǎn)錄本長度在201～1 6191 nt之間。其中長度為200～400 nt和400～600 nt的轉(zhuǎn)錄本分別為212 889個（占67.94%）和39 109個（占12.48%），長度在4 000 nt以上的轉(zhuǎn)錄本為4 351個（占1.39%），N50和N90值分別為958 nt和234 nt。轉(zhuǎn)錄本長度的分布情況見圖1。

將用于組裝的轉(zhuǎn)錄本與組裝后的轉(zhuǎn)錄本用Bowtie2（2.2.3版）進行比對［10］，其中比對上序列的比例占到87.02%，均一性分析結(jié)果良好，符合質(zhì)量要求。

圖1 Illumina HiSeqTM測序獲得的轉(zhuǎn)錄本長度分布

2.2 Unigene的功能注釋

2.2.1 Unigene的GO注釋 對所獲得的275 625個Unigene進行功能注釋，對注釋到SWISS-PROT、GO、KEGG數(shù)據(jù)庫的基因進行統(tǒng)計，結(jié)果如表1所示。通過GO（gene ontology，GO）注釋，可確定這些基因的功能分類，結(jié)果如圖2所示。從比對結(jié)果來看，仍有91.06%的基因未比對到任何數(shù)據(jù)庫，這可能與目前公布的家兔基因組信息尚不完全有關(guān)。

圖2中，Blast2GO軟件分析得到每個Unigene蛋白GO功能注釋，將所有Unigene比對到GO的3個數(shù)據(jù)庫：生物過程（biological process，BP）、分子功能（molecular function，MF）和細胞組成（cellular component，CC），并做功能分類統(tǒng)計。其中絕大多數(shù)基因都比對到細胞過程（cellular process）、代謝過程（metabolic process）、細胞組分（cell part）以及（分子）聯(lián)接（binding）等類別中。

表1 Unigene的比對結(jié)果

圖2 GO統(tǒng)計柱狀圖

2.2.2 Unigene的COG注釋 用BLAST軟件將組裝序列與COG數(shù)據(jù)庫比對，對結(jié)果進行過濾后去掉e＞10-5的比對序列，再進行注釋后結(jié)果如圖3所示。共對10 288個Unigene進行注釋，分為25個類別。其中有1 951個基因比對到一般機能（general function prediction）類，有1 084個基因比對到轉(zhuǎn)錄后修飾，蛋白折疊和分子伴侶 （posttranslational modification，protein turnover and chaperones）類，只有14個基因比對到細胞動能（cellular motility）類，5個基因比對到細胞核結(jié)構(gòu)（nuclear structure）類，另外還有477個基因功能未知（function unknown），可能是測序發(fā)現(xiàn)的新基因。

2.3 不同毛色獺兔差異表達基因的分析

2.3.1 差異表達基因的GO注釋 將不同毛色的2個樣品所獲得的轉(zhuǎn)錄本用DEGseq軟件進行比較，最終獲得差異表達基因12 408個，其中上調(diào)表達的基因4 904個，下調(diào)表達的基因7 504個。對所獲得的差異表達基因進行統(tǒng)計，并利用Uniprot、Pfam、GO和KEGG等數(shù)據(jù)庫對差異基因進行功能注釋，獲得詳細描述信息。根據(jù)GO數(shù)據(jù)庫中的相關(guān)條目，統(tǒng)計得到差異表達基因，如圖4所示。

GO統(tǒng)計結(jié)果顯示：所有差異基因比對到BP數(shù)據(jù)庫中的有2 711個，其中顯著富集在免疫反應(yīng)（immune response）的基因為154個，占比對到該功能類別總數(shù)的5.68%；顯著富集在免疫過程（immune system process）的基因為196個，占比對到該功能類別總數(shù)的7.22%；顯著富集在防御反應(yīng)（defense response）中的基因127個，占總數(shù)的4.68%。所有差異基因比對到MF數(shù)據(jù)庫中的有622個，其中顯著富集在RNA介導的DNA酶活性（RNA-directed DNA polymerase activity）中的121個，占比對到MF數(shù)據(jù)庫總數(shù)的19.45%；顯著富集在DNA酶活性（DNA polymerase activity）中的121個，占總數(shù)的19.45%。所有差異表達基因比對到CC數(shù)據(jù)庫中的有338個，其中顯著富集在胞外區(qū)域（extracellular region）的基因為202個，占比對到該功能類別總數(shù)的59.76%。

圖3 GOG注釋統(tǒng)計圖

圖4 差異表達基因的GO統(tǒng)計柱狀圖

由圖4可以看出，在BP數(shù)據(jù)庫所有上調(diào)表達的基因中，共有443個基因歸入細胞過程（cellular process）中，占總數(shù)的 20.18%，331個基因歸入代謝過程（metabolic process）中，占總數(shù)的15.08%；所有下調(diào)表達的基因中，共有541個基因歸入細胞過程中，占總數(shù)的20.76 %，440個基因歸入代謝過程中，占總數(shù)的16.88%。在MF數(shù)據(jù)庫所有上調(diào)表達的基因中，共有410個基因歸入細胞連接（binding）中，占總數(shù)的44.46%；所有下調(diào)表達的基因中，共有554個基因歸入細胞連接（binding）中，占總數(shù)的46.59%；在CC數(shù)據(jù)庫所有上調(diào)表達的基因中，共有519個基因歸入細胞組分（cell part）中，占總數(shù)的31.86%，共有656個基因歸入細胞組分（cell part）中，占總數(shù)的31.91%。

2.3.2 差異表達基因的KEGG通路分析 用KEGG注釋系統(tǒng)將測序所發(fā)現(xiàn)的所有轉(zhuǎn)錄本注釋到KEGG數(shù)據(jù)庫中，共對該數(shù)據(jù)庫中315個代謝通路進行比對，差異表達基因顯著富集的通路13個。其中3個與脂類代謝有關(guān)，1個與RNA轉(zhuǎn)運有關(guān)，2個與糖類代謝有關(guān)，1個與氨基酸代謝有關(guān)，1個與腮腺炎有關(guān)，1個是細胞色素P450合成的通路。表2是差異表達基因顯著富集通路的匯總。

表2 不同顏色獺兔耳組織中差異表達基因顯著性富集的通路

對細胞色素代謝通路（Map00980）中表達發(fā)生變化的基因進行進一步分析，共發(fā)現(xiàn)8個表達發(fā)生變化的基因，其中上調(diào)表達的基因有：GSTM4，ARK72；下調(diào)表達的基因有：CP2F1，UD11，ADH6，GSTA4，CP2BB，ST2A1。其中GSTM4（谷胱甘肽轉(zhuǎn)移酶）在多個細胞色素代謝通路中都出現(xiàn)上調(diào)表達，而CP2F1和CP2BB都是細胞色素P450家族成員，但它們的具體生物學功能尚未有文獻報道。

3 討論

3.1 獺兔轉(zhuǎn)錄組測序

高通量測序的樣本可選用組織，也可選用細胞或細胞系，這主要是由不同的試驗?zāi)康乃鶝Q定的。黑色素主要位于表皮和毛囊中，所以選擇獺兔的耳皮膚組織。已有研究者選用cDNA芯片分析了不同被毛密度獺兔皮膚組織中差異表達的基因［11］。Pan L等［12］采用Solexa測序技術(shù)研究了獺兔不同表型的皮膚組織（plaice skin and un-plaice skin）基因表達的差異，并對一個影響表型的關(guān)鍵基因Lamb3內(nèi)的SNPs進行分析。翁巧琴等［13］研究了MC1R基因與獺兔毛色的關(guān)系，最終發(fā)現(xiàn)在MC1R基因第40、50位點存在T＞C和A＞G的突變，在第45位點缺失一段長度為23bp的堿基序列，第507位點存在T＞C的堿基突變，這些突變可能與獺兔毛色有關(guān)。牛曉艷等［14］研究了獺兔MITF基因部分外顯子的多態(tài)性，并在外顯子1、5、7中發(fā)現(xiàn)一些多態(tài)與毛色存在強關(guān)聯(lián)。本研究擬通過獺兔耳組織的RNA-Seq，初步篩選一些與毛色相關(guān)的新基因。

3.2 樣本含量

RNA-Seq的生物學重復(fù)和技術(shù)重復(fù)的成本都比較高，如果采用多個樣本的RNA等量混池進行轉(zhuǎn)錄組測序，可能由于樣本間的個體差異和遺傳背景的差異等導致最后結(jié)果的誤差較大。又因為RNA測序價格的限制，綜合考慮后，決定選用全同胞的海貍色獺兔和白色獺兔各1只，并充分考慮到遺傳背景的一致性，在數(shù)據(jù)分析過程中提高了差異顯著性閾值，縮小差異基因的范圍，再結(jié)合差異基因的生物學功能和KEGG通路信息，最大程度地保證了結(jié)果的準確性和可參考性。

3.3 細胞色素代謝通路中差異表達基因

由表2可以看出，在13個差異基因顯著富集的通路中，只有Map編號為00980的細胞色素P450代謝通路與細胞色素代謝密切相關(guān)。該通路中有8個基因表達發(fā)生了變化，其中CP2F1（細胞色素P450家族成員2F1）和CP2BB（細胞色素P450家族成員2B11）屬于細胞色素家族成員；而UD11（UDP-葡糖甘酸轉(zhuǎn)移酶）、ADH6（乙醇脫氫酶-6）、ARK72（黃曲霉毒素乙醛還原酶家族成員2）、GSTM4（谷胱甘肽S-轉(zhuǎn)移酶）、ST2A1（類固醇磺基轉(zhuǎn)移酶）與毛色的關(guān)系尚未見文獻報道。Uehara等［15］2009年報道了GSTA4（谷胱甘肽S-轉(zhuǎn)移酶）基因與小鼠耳蝸皮膚組織黑素細胞遷移有關(guān)。Christopher等［16］2013年總結(jié)了影響貓和狗毛色的基因座，其中包括：MITF、PMEL、TYRP1、MLPH、MC1R、ASIP、KIT等經(jīng)典基因以及PSMB7、CBD103、TAQPEP等，但未包括本試驗中發(fā)現(xiàn)的基因。這可能與本試驗設(shè)計的局限以及樣本含量較少有關(guān)，但也可能是發(fā)現(xiàn)了一些影響獺兔毛色變化的新基因，這些基因的功能以及對毛色影響的具體方式還需要進一步的試驗研究。

4 結(jié)論

本研究通過轉(zhuǎn)錄組測序（RNA-Seq）方法對不同毛色獺兔耳組織中差異表達的基因進行了研究。共發(fā)現(xiàn)12 408個差異表達基因，這些基因主要與代謝過程、細胞組分、催化反應(yīng)有關(guān)。通過KEGG通路分析，最終在與細胞色素代謝有關(guān)的通路中發(fā)現(xiàn)8個差異表達的基因，對這些基因的功能進行進一步分析發(fā)現(xiàn)，只有GSTA4與細胞色素遷移有關(guān)，可能最終影響到獺兔的毛色。這些候選基因與毛色的關(guān)系還需要進一步進行功能驗證，最終找到相關(guān)分子標記，用于兔群的篩選以及毛色機理的研究。

［1］ 龐有志.為長毛兔和獺兔品種論不平——兼評家兔品種與品系的分類［J］.中國養(yǎng)兔雜志，2001（4）：30-34.

［2］ 龐有志.兔的毛色遺傳與家兔育種［J］.中國養(yǎng)兔雜志，2008（10）：30-35.

［3］ 楊廣禮.動物毛色形成的調(diào)控機制研究［J］.黑龍江畜牧獸醫(yī)，2014，3（上）：45-48.

［4］ 周華，張新，劉騰云，等.高通量轉(zhuǎn)錄組測序的數(shù)據(jù)分析與基因挖掘［J］.江西科學，2012，30（5）：607-611.

［5］ 岳桂東，高強，羅龍海，等.高通量測序技術(shù)在動植物研究領(lǐng)域中的應(yīng)用［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，F(xiàn)engZ，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］翁巧琴，王進波，王國朋，等.特種獺兔毛色相關(guān)基因多態(tài)性的研究［J］.中國養(yǎng)兔雜志，2015（3）：4-9.

［14］牛曉艷，任克良，曹亮，等.獺兔MITF基因部分外顯子多態(tài)性研究［J］.山西農(nóng)業(yè)科學，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

國家兔產(chǎn)業(yè)技術(shù)體系項目（CARS-44-B-6）；山西省農(nóng)科院博士基金項目（YBSJJ1401）

牛曉艷（1984-），女，助理研究員，博士，主要從事海貍色獺兔育種工作。

任克良，Email：keliangren@sohu.com