蘭道亮，熊顯榮，柴志欣，艾　鷖，黃　偲，李　鍵*

(1.西南民族大學(xué) 青藏高原研究院，成都 610041；2.西南民族大學(xué) 生命科學(xué)與技術(shù)學(xué)院，成都 610041)

?

牦牛發(fā)情期卵巢比較轉(zhuǎn)錄組學(xué)研究

蘭道亮1，2，熊顯榮2，柴志欣1，艾鷖1，黃偲2，李鍵1，2*

(1.西南民族大學(xué) 青藏高原研究院，成都 610041；2.西南民族大學(xué) 生命科學(xué)與技術(shù)學(xué)院，成都 610041)

旨在進一步了解牦牛發(fā)情期卵巢的分子機制，解析牦牛繁殖的特殊性。本研究應(yīng)用RNA-seq技術(shù)對牦牛和平原黃牛發(fā)情期卵巢進行轉(zhuǎn)錄組高通量測序及全基因組差異表達(dá)模式比對分析。通過比較分析牦牛和黃牛卵巢轉(zhuǎn)錄組數(shù)據(jù)，共篩選出1 307個差異表達(dá)基因，其中661個基因表達(dá)量上調(diào)和646個基因表達(dá)量下調(diào)。進一步功能分析表明，這些差異基因涉及多種 GO分類及KEGG通路。其中GO分類注釋顯示，差異基因與細(xì)胞粘附、激素調(diào)控等生物學(xué)過程存在密切關(guān)聯(lián)，同時鈣離子結(jié)合、陽離子跨膜轉(zhuǎn)運等分子事件表現(xiàn)活躍。KEGG通路分析顯示，補體和凝血級聯(lián)通路的富集水平最高，其次為細(xì)胞色素P450相關(guān)通路。晝夜節(jié)律等一些新型通路，盡管與生殖功能沒有明顯關(guān)聯(lián)，但也表現(xiàn)出顯著富集。本研究首次對比牦牛和黃牛發(fā)情期卵巢轉(zhuǎn)錄組數(shù)據(jù)，篩選并分析相關(guān)差異基因。該研究結(jié)果為進一步闡述牦牛卵巢的基本分子機理提供基礎(chǔ)，同時也為全面理解牦牛繁殖特異性提供新的思路。

牦牛；發(fā)情期卵巢；轉(zhuǎn)錄組；比較分析；分子機制

牦牛(Bosgrunniens)素有 “高原之舟”的美譽，是主要分布在中國青藏高原及其毗鄰高山或亞高山地區(qū)的特有物種，也是當(dāng)今世界上已知唯一一種能夠在高海拔地區(qū)(平均海拔高度3 000米)生存的牛屬動物[1-2]。牦牛能夠很好地適應(yīng)高山草原環(huán)境，并在惡劣的高原環(huán)境條件下(如空氣稀薄、低溫和草料短缺)自由生長和繁殖[3]。同時牦?？梢詾楫?dāng)?shù)啬撩裉峁┠?、肉、毛等豐富的生產(chǎn)和生活資料，是高原地區(qū)不可或缺的重要畜種[1，4-5]。與生活在平原地區(qū)的普通牛類相比，牦牛的性成熟發(fā)育得更慢，其生殖能力普遍更低[6-7]。一頭成年牦牛的平均繁殖率僅為48.61%，其中超過一半的牦牛為每兩年生育一胎或每三年生育兩胎。此外，雌性牦牛發(fā)情率也很低，超過90%產(chǎn)后雌性牦牛無法在同年的發(fā)情期發(fā)情[8-9]。

卵巢是雌性哺乳動物的重要生殖器官。它具有眾多功能，包括提供能受精的卵母細(xì)胞、分泌生殖激素以及維持雌性動物發(fā)情周期。卵巢功能直接影響著雌性動物的繁殖能力[10]。每個發(fā)情周期過程中，卵巢都會經(jīng)歷增殖、侵襲、分化和細(xì)胞凋亡；這些正常的生理變化直接影響和決定了雌性動物的排卵、受精率和產(chǎn)仔數(shù)[11]。與其他普通牛類相比，牦牛的卵巢較小，卵巢系膜更短，其位置也相對固定，但與前者的整體結(jié)構(gòu)是相似的[12]。迄今為止，針對牦牛卵巢的研究多集中于其形狀和解剖結(jié)構(gòu)，而對牦牛卵巢的分子基礎(chǔ)及其分子機制仍然尚不明確。卵巢功能的實現(xiàn)是一個復(fù)雜的過程，這個過程涉及大量基因的轉(zhuǎn)錄調(diào)控。另外，基因表達(dá)的差異性是物種進化的重要組成部分，也是產(chǎn)生生物多樣的根本途徑[11，13]。因此，為了解牦牛卵巢的分子機制和牦牛繁殖的特殊性，需要對其轉(zhuǎn)錄組學(xué)進行研究。隨著高通量測序技術(shù)的發(fā)展，RNA高通量測序(RNA sequencing)為大范圍的轉(zhuǎn)錄學(xué)研究提供強有力的工具，與傳統(tǒng)方法相比，其優(yōu)勢十分明顯[11]。而且隨著牦?；蚪M測序的完成也使進一步的轉(zhuǎn)錄組學(xué)研究成為可能[1]。本研究應(yīng)用RNA-seq技術(shù)對牦牛和平原黃牛發(fā)情期卵巢進行轉(zhuǎn)錄組學(xué)測序及比對分析，深入了解牦牛卵巢在發(fā)情周期中的差異表達(dá)基因及其調(diào)控通路，為進一步研究牦牛卵巢的分子機制及全面了解牦牛繁殖的特異性提供基礎(chǔ)。

1　材料與方法

1.1樣本采集

參照文獻[14-15]的方法對牦牛的發(fā)情期進行鑒定。從高原屠宰場(中國四川紅原，緯度31°51′ N～33°19′ N，經(jīng)度101°51′ E～103°23′ E；平均海拔高度為3 600 m)隨機選擇3頭處于自然發(fā)情周期且體型相當(dāng)?shù)慕】?歲雌性麥洼牦牛。同樣從高原屠宰場(中國四川成都，緯度30°05′ N～31°26′ N，經(jīng)度102°54′ E～104°53′ E；平均海拔高度為750米)隨機選擇3頭處于自然發(fā)情期的健康4歲雌性黃牛(Bostaurusdomesticus)。屠宰后立即采集牦牛和黃牛的卵巢組織并進行觀察確認(rèn)后，將組織放液氮中保存?zhèn)溆谩?/p>

1.2RNA提取、cDNA文庫構(gòu)建和Illumina測序

應(yīng)用TRIzol試劑(Life Technologies公司，美國)提取牦牛和黃牛卵巢組織的總RNA。提取后的各RNA樣品均以無RNase的DNase I(TaKaRa公司，大連)進行處理，以消除可能存在基因組DNA污染。然后分別等量混合3個牦牛和黃牛的RNA樣品，組成兩個RNA池(RNA Pool，各30 μg)。使用Oligotex mRNA小量提取試劑盒(Qiagen公司，德國)從RNA池中分離和純化mRNA。整個流程 RNA 的質(zhì)量和含量通過Agilent 2100生物分析儀(Agilent公司，美國)進行測定。使用 TruSeq RNA Sample Prep Kit試劑盒(Illumina公司，美國)，參照其說明書標(biāo)準(zhǔn)流程構(gòu)建隨機片段測序文庫。即用中段試劑在恒溫混勻儀中對純化的mRNA進行片段化(約200 bp)，而后將其作為cDNA第一鏈合成的模板。隨后使用RNase H、dNTP及DNA聚合酶I合成第二鏈cDNA。純化和配對末端修復(fù)后，將cDNA片段連接到測序接頭并通過聚合酶鏈反應(yīng)(PCR)對其進行擴增，以此獲得最終的雙末端(Paired-end，PE)測序文庫。使用Agilent2100生物分析儀和ABI StepOnePlus實時PCR系統(tǒng)執(zhí)行質(zhì)控(QC)測試后，在Illumina HiSeq 2500平臺上對文庫進行測序。

1.3轉(zhuǎn)錄組數(shù)據(jù)分析

對Hiseq 2500測序產(chǎn)生的原始讀數(shù)進行QC測試后，通過去除接頭序列、空序列及低質(zhì)量測序序列(≤Q20質(zhì)量值)，將原始讀數(shù)過濾為凈測序序列(Clean reads)。使用SOAPaligner/SOAP2軟件，將過濾后序列與牦牛基因組(版本1.0)進行比對[16]。然后利用比對數(shù)據(jù)計算這些測序序列在基因組及基因上的分布，并統(tǒng)計覆蓋率。根據(jù)RPKM法 (Reads per kilobase transcriptome per million mapped reads)計算基因表達(dá)水平。采用 “數(shù)字基因表達(dá)譜顯著性”算法(The significance of digital gene expression profiles)界定牦牛和普通牛類文庫之間的差異表達(dá)基因(Differential expressed genes，DEGs)[17]，同時用P值檢測對應(yīng)達(dá)到統(tǒng)計學(xué)顯著水平的差異基因表達(dá)。偽發(fā)現(xiàn)率(FDR)≤0.001和log2Ratio的絕對值≥1作為確定基因表達(dá)差異顯著性的閾值。通過GO和KEGG通路富集分析(P≤0.05)，對DEGs的功能進行注釋。

1.4Real-time PCR (qRT-PCR) 驗證

從測序的差異表達(dá)基因中隨機選擇15個基因進行qRT-PCR驗證(基因、引物信息見附表1)。使用TRIzol試劑(Life Technologies公司，美國)提取同期取樣的牦牛和黃牛卵巢組織總RNA。用無核酸酶ddH2O將各RNA樣品的濃度調(diào)整至1 μg·μL-1，并用SuperScript?Ⅲ第一鏈合成系統(tǒng)(Life Technologies公司，美國)在20 μL的反應(yīng)系統(tǒng)中對2 μg的總RNA進行反轉(zhuǎn)錄。然后使用SYBRH Premix Ex TaqTM Ⅱ(Tli RNaseH Plus)試劑盒(TaKaRa公司，大連)和Bio-Rad CFX96 TouchTM實時PCR系統(tǒng)(Bio-Rad公司，美國)進行qRT-PCR操作。反應(yīng)條件：95 ℃ 1 min，95 ℃ 10 s，60 ℃ 40 s，總共40個循環(huán)。反應(yīng)結(jié)束后進行最終熔解度曲線分析。采用β-Actin作為內(nèi)參基因，根據(jù)2-ΔΔCt法計算各基因的相對表達(dá)量。

2　結(jié)　果

2.1轉(zhuǎn)錄組測序數(shù)據(jù)分析

轉(zhuǎn)錄組測序數(shù)據(jù)已遞交于NCB數(shù)據(jù)庫(接受號：SRX335453、SRX1302580)。去除低質(zhì)量序列(即僅包含接頭、接頭序列的讀數(shù)和空讀序列)后，在牦牛轉(zhuǎn)錄測序文庫中共獲得53 653 032條過濾測序序列，包含4 828 772 880 堿基(bp)。在黃牛轉(zhuǎn)錄測序文庫中共獲得54 855 396條過濾測序序列，包含4 936 985 640堿基(bp)。堿基組成和質(zhì)量分析結(jié)果顯示，原始測序數(shù)據(jù)表現(xiàn)出均衡的堿基組成，且牦牛和普通黃牛文庫的Q20(即堿基質(zhì)量值≥20)的比率分別為94.7%和94.5%，這一結(jié)果表明，文庫構(gòu)建成功且測序質(zhì)量良好。本研究以牦?；蚪M作為參考基因組。比對分析顯示，牦牛文庫中有33 168 751條(61.82%)測序序列、普通黃牛文庫中有31 730 973條測序序列(57.84%)可比對至牦?；蚪M?；蚋采w率統(tǒng)計表明，牦牛和普通黃牛發(fā)情期卵巢文庫中分別有16 992條和16 904條牦?；虻玫接成?。此外，筆者發(fā)現(xiàn)一定比例的、無法映射至現(xiàn)有基因的測序序列(約30%)，但是這些測序序列均位于基因組上，為進一步挖掘牦牛新轉(zhuǎn)錄子或基因提供數(shù)據(jù)基礎(chǔ)。

2.2差異表達(dá)基因篩選

通過比較分析牦牛和黃牛發(fā)情期卵巢轉(zhuǎn)錄組數(shù)據(jù)，共篩選出1 307個差異表達(dá)基因，其中661個基因表達(dá)量上調(diào)和646個基因表達(dá)量下調(diào)。為驗證這些基因表達(dá)數(shù)據(jù)，筆者隨機選擇15個差異表達(dá)基因，以同期取樣的牦牛和黃牛卵巢組織為樣本進行qRT-PCR檢測。結(jié)果顯示，這些差異基因的基本表達(dá)模式表現(xiàn)出與轉(zhuǎn)錄組數(shù)據(jù)相似的變化趨勢(表1)，表明轉(zhuǎn)錄組測序結(jié)果可靠。

2.3差異基因功能分析

為深入了解差異表達(dá)基因的生物學(xué)意義，筆者對這些差異表達(dá)基因進行GO和KEGG富集分析。GO富集分析表明，差異表達(dá)基因在分布于生物學(xué)過程(Biological processes，BP)、細(xì)胞組分(Cellular components，CC) 及 分 子 功 能 (Molecular functions，MF)3大類中的GO類別上均出現(xiàn)富集。按生物過程分類，共有3 067個GO類別出現(xiàn)富集，富集前10的類別(Top 10)已列于表1中。在富集前10的GO類別中，“生物粘附”(Biological adhesion)、“細(xì)胞粘附”(Cell adhesion)和“細(xì)胞間粘附”(Cell-cell adhesion)等3個粘附相關(guān)GO類別的富集水平最高。此外，與激素相關(guān)的兩個GO類別，即“類固醇激素刺激響應(yīng)”(Response to steroid hormone stimulus)和“雌激素刺激響應(yīng)”(Response to estrogen stimulus)，也富集于前10 GO類別中。在分子功能的分類下，共有667個GO類別出現(xiàn)富集，位列前10的顯著富集GO類別已列于表2中。在位列前10的GO類別中，“鈣離子結(jié)合”(Calcium ion binding)的富集水平最高。在前10富集GO類別中，陽離子跨膜轉(zhuǎn)運相關(guān)(Cation transmembrane transport-related)GO類別占據(jù)了絕大部分(6/10)。在細(xì)胞組分的分類下，共有427個GO類別出現(xiàn)富集。在位列前10的富集GO類別中(表3)，所有類別均屬于細(xì)胞外(Extracellular parts)或質(zhì)膜部分(Plasma membrane parts)。KEGG富集分析表明，DEGs中共有243個通路出現(xiàn)富集，在位列前10的通路中(表4)，“補體和凝血級聯(lián)”(Complement and coagulation cascades)通路的富集水平最高，其次 “細(xì)胞色素P450引起的外源性物質(zhì)代謝”(Metabolism of xenobiotics by cytochrome P450)和“藥物代謝-細(xì)胞色素P450”(Drug metabolism-cytochrome P450)等與細(xì)胞色素P450相關(guān)通路也得到顯著富集。一些與卵巢生殖功能無顯著關(guān)聯(lián)的新型通路如“晝夜節(jié)律”通路(Circadian rhythm)也出現(xiàn)在富集水平前10的通路中。

表1生物學(xué)過程前10富集 GO 分類

Table 1Top 10 enrichment GO categories in the biological processes of DGEs

GO分類GOterm集群頻率a/%Clusterfrequency基因組頻率b/%GenomefrequencyP值CorrectedP-valueBiologicaladhesion9.44.68.88e-08Celladhesion9.34.61.54e-07Cell-celladhesion4.62.00.00078Anatomicalstructuredevelopment32.325.50.00139Anatomicalstructuremorphogenesis17.612.30.00142Developmentalprocess35.828.70.00182Responsetosteroidhormonestimulus4.92.30.00355Responsetoestrogenstimulus2.91.10.00577Responsetoxenobioticstimulus2.30.80.03572Glutathionemetabolicprocess1.20.30.03925

a.不同GO分類中的差異基因占所有GO分類總差異基因(965)的比例；b.該GO分類包含的基因占基因組所有基因(15 573)的比例

a.The proportion of differentially expressed genes with biological processes GO categories in total differentially expressed genes(965)with GO annotation；b.The proportion of genes with GO annotation in the entire yak genome(15 573)

表2分子功能過程前10富集 GO 分類

Table 2Top 10 enrichment GO categories in the molecular functions of DGEs

GO分類GOterm集群頻率a/%Clusterfrequency基因組頻率b/%GenomefrequencyP值CorrectedP-valueCalciumionbinding9.44.00.00232Metaliontransmembranetransporteractivity9.32.40.00237Cationchannelactivity4.61.80.00647Transmembranetransporteractivity32.36.50.01108Glycosaminoglycanbinding17.61.10.01124Monovalentinorganiccationtransmembranetransporteractivity35.82.20.01168Cationtransmembranetransporteractivity4.93.70.01386Gatedchannelactivity2.91.90.01587Inorganiccationtransmembranetransporteractivity2.33.10.02024Carbohydratederivativebinding1.21.20.03270

a.不同GO分類中的差異基因占所有GO分類總差異基因(944)的比例；b.該GO分類包含的基因占基因組所有基因(15 811)的比例

a.The proportion of differentially expressed genes with molecular functions GO categories in total differentially expressed genes(944)with GO annotation；b.The proportion of genes with GO annotation in the entire yak genome(15 811)

表3細(xì)胞組分過程前10富集 GO 分類

Table 3Top 10 enrichment GO categories in the cellular components of DGEs

GO分類GOterm集群頻率a/%Clusterfrequency基因組頻率b/%GenomefrequencyP值CorrectedP-valueExtracellularregion19.811.77.08e-13Extracellularregionpart12.06.41.32e-09Plasmamembranepart17.311.23.74e-07Extracellularspace9.05.04.77e-06Cellperiphery31.524.53.42e-05Intrinsictoplasmamembrane11.07.00.00024Plasmamembrane30.223.90.00045Cellsurface5.62.90.00065Integraltoplasmamembrane10.56.80.00078Extracellularmatrix4.72.50.00364

a.不同GO分類中的差異基因占所有GO分類總差異基因(1 054)的比例；b.該GO分類包含的基因占基因組所有基因(17 161)的比例

a.The proportion of differentially expressed genes with cellular components GO categories in total differentially expressed genes(1 054)with GO annotation；b.The proportion of genes with GO annotation in the entire yak genome(17 161)

表4富集前 10 (Top 10) KEGG 通路列表

Table 4Top 10 enrichment KEGG pathways in the DGEs

信號通路Pathway基因數(shù)(通路頻率/%)aGenenumber(Pathwayfrequency)基因數(shù)(基因組頻率/%)bGenenumber(Genomefrequency)P值CorrectedP-valueExtracellularregion31(2.68)206(1.09)2.910213e-06Extracellularregionpart83(7.17)938(4.95)0.00039214Plasmamembranepart20(1.73)145(0.76)0.00051927Extracellularspace66(5.70)720(3.80)0.00060855Cellperiphery14(1.21)89(0.47)0.00095037Intrinsictoplasmamembrane15(1.30)117(0.62)0.00501102Plasmamembrane34(2.94)366(1.93)0.009756936Cellsurface60(5.19)729(3.84)0.01112557Integraltoplasmamembrane12(1.04)94(0.50)0.01174199Extracellularmatrix6(0.52)36(0.19)0.02055234

a.不同通路中的差異基因占所有通路總差異基因(1 157)的比例；b.該通路包含的基因占基因組所有基因(18 965)的比例

a.The proportion of differentially expressed genes with pathway annotation in total differentially expressed genes(1 157)with pathway annotation；b.The proportion of genes with GO annotation in the entire yak genome(18 965)

3　討　論

為了解牦牛卵巢的分子機制和牦牛繁殖的特殊性，本研究應(yīng)用RNA-seq技術(shù)對牦牛和平原黃牛發(fā)情期卵巢進行轉(zhuǎn)錄組學(xué)測序及比對分析，共篩選獲得1 307個差異表達(dá)基因，進一步功能分析這些差異表達(dá)基因涉及多個GO分類和KEGG信號通路。在生物過程GO分類中，發(fā)現(xiàn)“生物粘附”(Biological adhesion)、“細(xì)胞粘附”(Cell adhesion)和“細(xì)胞間粘附”(Cell-cell adhesion)等3個粘附相關(guān)GO類別的富集水平最高。以往的研究表明，卵巢中的卵泡處于無血管的微環(huán)境中，卵母細(xì)胞、卵丘細(xì)胞和其他卵泡細(xì)胞間存在廣泛的間隙連接，從而形成一個功能完備的聯(lián)合體。卵母細(xì)胞主要通過細(xì)胞黏附和連接與其周圍細(xì)胞(卵泡顆粒及膜細(xì)胞)進行通訊[18-19]。本研究結(jié)果表明，處于發(fā)情周期中的牦牛和黃牛卵巢間可能存在不同的細(xì)胞粘附方式。在分子功能GO分類中，在前10富集的GO類別中，“鈣離子結(jié)合”(Calcium ion binding)的富集水平最高。作為第二細(xì)胞內(nèi)信使，鈣離子調(diào)節(jié)著細(xì)胞中許多重要的生理和病理過程。除了調(diào)節(jié)哺乳動物卵母細(xì)胞的增長和減數(shù)分裂、參與激素和生長因子分泌外，鈣離子也與卵母細(xì)胞的減數(shù)分裂成熟調(diào)節(jié)密切相關(guān)[20-21]?！扳}離子結(jié)合”類別在本研究中的富集水平最高，這說明牦牛和黃牛卵巢間的卵母細(xì)胞減數(shù)分裂成熟調(diào)節(jié)可能存在差異。在前10富集GO類別中，陽離子跨膜轉(zhuǎn)運相關(guān)(Cation transmembrane transport-related)GO類別占據(jù)大部分(6/10)。細(xì)胞膜的離子活性與細(xì)胞內(nèi)環(huán)境調(diào)節(jié)密切相關(guān)，如pH、滲透壓、營養(yǎng)吸收和膜電位；它在神經(jīng)遞質(zhì)分泌和激素代謝調(diào)節(jié)中也起著重要的作用[22]。另外，離子跨膜轉(zhuǎn)運建立的細(xì)胞膜濃度梯度為一系列細(xì)胞代謝活動提供能量來源。在許多陽離子跨膜轉(zhuǎn)運活動中，鈣和鋅離子的轉(zhuǎn)運在卵母細(xì)胞成熟調(diào)節(jié)中發(fā)揮著重要作用[23-25]。本研究中，陽離子跨膜轉(zhuǎn)運相關(guān)GO類別的富集證實，與黃牛類卵巢相比，牦牛卵巢的離子跨膜轉(zhuǎn)運調(diào)節(jié)有所不同，但其確切功能仍有待進一步研究。在細(xì)胞組分的分類下，位列前10的富集GO類別中，所有類別均屬于細(xì)胞外或質(zhì)膜部分。這一結(jié)果與生物過程和分子功能分類一致。例如，粘附、激素和陽離子跨膜轉(zhuǎn)運相關(guān)類別均發(fā)生于細(xì)胞外或質(zhì)膜部分，進一步證實富集的準(zhǔn)確性。

在KEGG富集分析中，其中“補體和凝血級聯(lián)”(Complement and coagulation cascades)通路的富集水平最高。補體是在具有酶活性組織液中發(fā)現(xiàn)的一組不耐熱性球蛋白，在先天性和獲得性免疫反應(yīng)中發(fā)揮著重要作用[26]。補體通過經(jīng)典、凝集素和旁路這3種途徑，可以激活一系列級聯(lián)反應(yīng)，進而觸發(fā)補體分子效應(yīng)。激活后的補體可以調(diào)節(jié)吞噬細(xì)胞和免疫細(xì)胞，并增強炎癥反應(yīng)；還可以通過溶解細(xì)胞壁直接殺滅某些細(xì)菌或細(xì)胞[27]。血凝是指血漿中的纖維蛋白通過外在和內(nèi)在凝血級聯(lián)通路轉(zhuǎn)化為不溶性纖維蛋白的過程[28]。“補體和凝血級聯(lián)”通路屬于一種內(nèi)源性代謝級聯(lián)。除了主體生理功能，該通路也與許多其他生理和病理調(diào)節(jié)過程密切相關(guān)。例如，補體系統(tǒng)在許多生物過程中都起著重要作用，其中包括繁殖、發(fā)育、干細(xì)胞分化和組織再生[29-30]。補體系統(tǒng)的Bf和DAF因子在卵細(xì)胞和配子的卵泡發(fā)育和成熟中扮演重要角色[31]。如圖1所示，3種級聯(lián)方式激活補體的關(guān)鍵基因的表達(dá)水平都出現(xiàn)升高。這一結(jié)果說明，包括細(xì)胞裂解、肌肉收縮、趨化、吞噬細(xì)胞復(fù)原和炎癥在內(nèi)的下游生理功能可能得到增強。在凝血級聯(lián)通路中發(fā)現(xiàn)，激活凝血的外在通路中基因僅出現(xiàn)少量表達(dá)。表明，包括血管損傷、纖維蛋白降解產(chǎn)物和一氧化氮生物合成在內(nèi)的主要生理功能可能出現(xiàn)下降。然而，內(nèi)在通路中的基因也出現(xiàn)大量表達(dá)，表明包括炎癥和細(xì)胞粘附在內(nèi)的相關(guān)下游生理功能可能被加強。導(dǎo)致細(xì)胞粘附過程的強化與GO類別生物過程的富集結(jié)果是相一致的，證明細(xì)胞粘附過程是牦牛和黃牛卵巢生理活動的重要差異之一。該結(jié)果表明，“補體和凝血級聯(lián)”通路在牦牛卵巢的生理活動中具有重要作用，除能夠調(diào)節(jié)免疫力、炎癥外還可能與繁殖和發(fā)育密切相關(guān)。該通路可能與牦牛生活的特殊環(huán)境有關(guān)。為了抵御惡劣的高原環(huán)境引起的卵巢損傷和功能障礙，牦?？赡苄纬梢环N比黃牛更強的分子防御機制。

在位列前10的KEGG通路中，“細(xì)胞色素P450引起的外源性物質(zhì)代謝” (Metabolism of xenobiotics by cytochrome P450)和“藥物代謝-細(xì)胞色素P450”(Drug metabolism-cytochrome P450)等與細(xì)胞色素P450相關(guān)通路也出現(xiàn)富集。細(xì)胞色素P450為一類亞鐵血紅素-硫醇鹽蛋白的超家族，主要包括CYP1、CYP2和CYP3亞家族。細(xì)胞色素P450是一種重要的代謝酶系統(tǒng)，可以催化體內(nèi)多種內(nèi)源性和外源性化合物的氧化反應(yīng)，更改化學(xué)物質(zhì)中的活性基團或有毒化合物，以及改變藥物的療效。細(xì)胞色素P450參與外源物質(zhì)的生物轉(zhuǎn)化，并與內(nèi)源性物質(zhì)的代謝存在關(guān)聯(lián)，同時它也是機體代謝調(diào)節(jié)中一種重要的化合物[32]。此外，細(xì)胞色素P450參與了固醇和類固醇激素的生物合成，尤其是性激素的生物合成，在生殖激素調(diào)節(jié)中發(fā)揮著重要的作用[33-34]。這一結(jié)果與生物過程中激素相關(guān)GO類別的富集相一致，證明牦牛和黃牛卵巢之間存在著激素調(diào)節(jié)系統(tǒng)方面的差異。缺氧，包括高海拔缺氧，能夠顯著改變各細(xì)胞色素P450亞家族中的基因表達(dá)水平。缺氧誘導(dǎo)因子1(HIF-1)，一種在低氧刺激下調(diào)節(jié)耐缺氧基因的重要因子，尤其是在高海拔缺氧調(diào)節(jié)扮演重要角色。研究表明作為細(xì)胞色素P450成員的CYP3A6啟動子可以與HIF-1基因的特異性和轉(zhuǎn)錄激活相結(jié)合，從而誘導(dǎo)耐低氧基因的表達(dá)[35]。細(xì)胞色素P450的成員CYP17A1和CYP2E1基因也與藏族人存在密切聯(lián)系，能夠幫助他們適應(yīng)高原的低氧環(huán)境[36]。牦牛長久以來一直生活在青藏高原。鑒于本研究中有兩個細(xì)胞色素P450通路出現(xiàn)顯著富集，可能也與牦牛對高原低氧環(huán)境的適應(yīng)能力有關(guān)。

加粗黑框基因為差異表達(dá)基因。下同Bold boxes indicate the genes that were differentially expressed.The same as below圖1　補體和凝血級聯(lián)通路Fig.1　Differentially expressed genes in the complement and coagulation cascades pathways

一些與卵巢生殖功能無顯著關(guān)聯(lián)的通路也出現(xiàn)在富集水平前10的通路中，比如 “晝夜節(jié)律”(Circadian rhythm)通路。不同生物的生理和代謝活動以及行為均遵循著某種晝夜節(jié)律，這種特殊的機制稱之為“生物鐘”。該機制的分子基礎(chǔ)有賴于相關(guān)生物鐘基因的協(xié)調(diào)表達(dá)；產(chǎn)生的生物鐘蛋白會形成一個負(fù)向轉(zhuǎn)錄-翻譯反饋環(huán)路，這個環(huán)路維持著晝夜節(jié)律并幫助哺乳動物適應(yīng)環(huán)境的變化[37]。哺乳動物中已確定的主要生物鐘蛋白包括CLOCK蛋白、BMALL蛋白、CRY家族以及PER家族，其中CLOCK和BMALL蛋白在晝夜節(jié)律調(diào)節(jié)中處于核心地位[38]。本研究中，“晝夜節(jié)律”通路出現(xiàn)富集，且與黃牛數(shù)據(jù)相比，牦牛的CLOCK基因表達(dá)量下調(diào)(圖2)。筆者推測，這可能與牦牛生活的特殊環(huán)境有關(guān)。由于青藏高原的高海拔、缺氧環(huán)境與強紫外線，牦牛對環(huán)境條件如溫度和光線的微弱波動感受更加敏感。同時牦牛的繁育活動具有典型的季節(jié)性特征，易受到季節(jié)性變化的影響，這可能也與該通路具有相關(guān)聯(lián)系。該通路在牦牛卵巢中的確切功能，仍需進一步研究證實。

圖2　晝夜節(jié)律通路Fig.2　Differentially expressed genes in the circadian rhythm pathway

4　結(jié)　論

應(yīng)用 RNA-Seq 技術(shù)首次對牦牛和黃牛發(fā)情期卵巢進行了轉(zhuǎn)錄組學(xué)比對分析，共篩選出1 307個顯著差異表達(dá)基因，同時對這些差異表達(dá)基因進行相關(guān)功能分析。本研究結(jié)果為進一步了解牦牛卵巢的基本分子機制，以及牦牛繁殖的特殊性提供分子基礎(chǔ)。

[1]QIU Q，ZHANG G，MA T，et al.The yak genome and adaptation to life at high altitude[J].NatGenet，2012，44(8)：946-949.

[2]XIONG X，F(xiàn)U M，LAN D，et al.Yak response to high-altitude hypoxic stress by altering mRNA expression and DNA methylation of hypoxia-inducible factors[J].AnimBiotechnol，2015，26(3)：222-229.

[3]MIAO F，GUO Z，XUE R，et al.Effects of grazing and precipitation on herbage biomass，herbage nutritive value，and yak performance in an alpine meadow on the Qinghai-Tibetan Plateau[J].PLoSOne，2015，10(6)：e0127275.

[4]FU M，CHEN Y，XIONG X，et al.Establishment of mammary gland modelinvitro：culture and evaluation of a yak mammary epithelial cell line[J].PLoSOne，2014，9(12)：e113669.

[5]GUO X，LONG R，KREUZER M，et al.Importance of functional ingredients in yak milk-derived food on health of Tibetan nomads living under high-altitude stress：a review[J].CritRevFoodSciNutr，2014，54(3)：292-302.

[6]LONG R J，ZHANG D G，WANG X，et al.Effect of strategic feed supplementation on productive and reproductive performance in yak cows[J].PrevVetMed，1999，38 (3)：195-206.

[7]PRAKASH B S，SARKAR M，MONDAL M.An update on reproduction in yak and mithun[J].ReprodDomestAnim，2008，43(Suppl 2)：217-223.

[8]FU M，XIONG X R，LAN D L，et al.Molecular characterization and tissue distribution of estrogen receptor genes in domestic yak[J].Asian-AustralasJAnimSci，2014，27(12)：1684-1690.

[9]XIAO X，ZI X D，NIU H R，et al.Effect of addition of FSH，LH and proteasome inhibitor MG132 toinvitromaturation medium on the developmental competence of yak (Bosgrunniens) oocytes[J].ReprodBiolEndocrinol，2014，12：30.

[10]PAN L，GONG W，ZHOU Y，et al.A comprehensive transcriptomic analysis of infant and adult mouse ovary[J].GenomicsProteomicsBioinformatics，2014，12(5)：239-248.

[11]MIAO X，LUO Q.Genome-wide transcriptome analysis between small-tail Han sheep and the Surabaya fur sheep using high-throughput RNA sequencing[J].Reproduction，2013，145(6)：587-596.

[12]CUI Y，YU S J.An anatomical study of the internal genital organs of the yak at different ages[J].VetJ，1999，157(2)：192-196.

[13]BECKER J，HACKL M，RUPP O，et al.Unraveling the Chinese hamster ovary cell line transcriptome by next-generation sequencing[J].JBiotechnol，2011，156(3)：227-235.

[14]權(quán)凱，張兆旺.母牦牛的繁殖特性[J].畜牧與飼料科學(xué)，2005，26(5)：30-32.

QUAN K，ZHANG Z W.Reproductive characteristics of female yak[J].AnimalHusbandryandFoodScience，2005，26(5)：30-32.(in Chinese)

[15]ZI X D，HE S M，LU H，et al.Induction of estrus in suckled female yaks (Bosgrunniens) and synchronization of ovulation in the non-sucklers for timed artificial insemination using progesterone treatments and Co-Synch regimens[J].AnimReprodSci，2006，92(2)：183-192.

[16]HU Q，MA T，WANG K，et al.The yak genome database：an integrative database for studying yak biology and high-altitude adaption[J].BMCGenomics，2012，13：600.

[17]AUDIC S，CLAVERIE J M.The significance of digital gene expression profiles[J].GenomeRes，1997，7(10)：986-995.

[18]EPPIG J J.Intercommunication between mammalian oocytes and companion somatic cells[J].Bioessays，1991，13(11)：569-574.

[19]GILCHRIST R B，RITTER L J，ARMSTRONG D T.Oocyte-somatic cell interactions during follicle development in mammals[J].AnimReprodSci，2004，83：431-446.

[20]MAKKI M，SABOORI E，SABBAGHI M A，et al.Effects of selenium，calcium and calcium ionophore on human oocytesinvitromaturation in a chemically defined medium[J].IranJReprodMed，2012，10(4)：343-348.

[21]ROZINEK J，RAJMON R，PETR J，et al.Ultrastructural localisation of calcium deposits in pig ovarian follicles[J].AnimReprodSci，2006，91(2)：123-132.

[22]HEDIGER M A.Membrane permeability.The diversity of transmembrane transport processes[J].CurrOpinCellBiol，1997，9(4)：543-546.

[23]KIM A M，BERNHARDT M L，KONG B Y，et al.Zinc sparks are triggered by fertilization and facilitate cell cycle resumption in mammalian eggs[J].ACSChemBiol，2011，6(7)：716-723.

[24]ANCHORDOQUY J M，ANCHORDOQUY J P，SIRINI M A，et al.The importance of having zinc duringinvitromaturation of cattle cumulus-oocyte complex：role of cumulus cells[J].ReprodDomestAnim，2014，49(5)：865-874.

[25]AZUMA T，KONDO T，IKEDA S，et al.Effects of EDTA saturated with Ca2+(Ca-EDTA) on pig，bovine and mouse oocytes at the germinal vesicle stage during maturation culture and the involvement of chelation of Zn2+in pronuclear formation induction by Ca-EDTA[J].Reproduction，2002，124(2)：235-240.

[26]SUNYER J O，BOSHRA H，LORENZO G，et al.Evolution of complement as an effector system in innate and adaptive immunity[J].ImmunolRes，2003，27(3)：549-564.

[27]APPLEDORN D M，MCBRIDE A，SEREGIN S，et al.Complex interactions with several arms of the complement system dictate innate and humoral immunity to adenoviral vectors[J].GeneTher，2008，15(24)：1606-1617.

[28]MOISEYEV G，GIVLI S，BAR-YOSEPH P Z.Fibrin polymerization in blood coagulation-a statistical model[J].JBiomech，2013，46(1)：26-30.

[29]HARRIS C L，MIZUNO M，MORGAN B P.Complement and complement regulators in the male reproductive system[J].MolImmunol，2006，43(2)：57-67.

[30]MASTELLOS D，GERMENIS A E，LAMBRIS J D.Complement：an inflammatory pathway fulfilling multiple roles at the interface of innate immunity and development[J].CurrDrugTargetsInflammAllergy，2005，4(1)：125-127.

[31]HASTY L A，LAMBRIS J D，LESSEY B A，et al.Hormonal regulation of complement components and receptors throughout the menstrual cycle[J].AmJObstetGynecol，1994，170(1)：168-175.

[32]DANIELSON P B.The cytochrome P450 superfamily：biochemistry，evolution and drug metabolism in humans[J].CurrDrugMetab，2002，3(6)：561-597.

[33]KANDIEL M M，WATANABE G，TAYA K.Ovarian expression of inhibin-subunits，3beta-hydroxysteroid dehydrogenase，and cytochrome P450 aromatase during the estrous cycle and pregnancy of Shiba goats (Caprahircus)[J].ExpAnim，2010，59(5)：605-614.

[34]SAHMI M，NICOLA E S，PRICE C A.Hormonal regulation of cytochrome P450 aromatase mRNA stability in non-luteinizing bovine granulosa cellsinvitro[J].JEndocrinol，2006，190(1)：107-115.

[35]FRADETTE C，DU SOUICH P.Hypoxia-inducible factor-1 and activator protein-1 modulate the upregulation of CYP3A6 induced by hypoxia[J].BrJPharmacol，2003，140(6)：1146-1154.

[36]SIMONSON T S，YANG Y，HUFF C D，et al.Genetic evidence for high-altitude adaptation in Tibet[J].Science，2010，329(5987)：72-75.

[37]EDERY I.Circadian rhythms in a nutshell[J].PhysiolGenomics，2000，3(2)：59-74.

[38]PIGGINS H D.Human clock genes[J].AnnMed，2002，34(5)：394-400.

(編輯程金華)

Comparative Transcriptome Analysis between Yak and Cattle Estrus Ovary

LAN Dao-liang1，2，XIONG Xian-rong2，CHAI Zhi-xin1，AI Yi1，HUANG Cai2，LI Jian1，2*

(1.InstituteofQinghai-TibetanPlateau，SouthwestUniversityforNationalities，Chengdu610041，China；2.CollegeofLifeScienceandTechnology，SouthwestUniversityforNationalities，Chengdu610041，China)

For elaborate the general physiological molecular features of the estrus ovary of yaks and particularity of yak reproduction.RNA-seq technology was applied to analyze transcriptome data comparatively between the yak and plain cattle estrus ovaries in this study.Compared with cattle，genome-wide divergent expression patterns during the ovary transcriptome data revealed that 1 307 genes were significantly and differentially expressed，in which 661 genes were upregulated and 646 genes were downregulated.Functional analysis showed that the differentially expressed genes were involved in various GO categories and KEGG pathways.GO annotations indicated that the genes were closely related to cell adhesion，hormonal biological processes，whereas the calcium ion binding，cation transmembrane transport molecular events were significantly active.KEGG pathway analysis showed that the complement and coagulation cascade pathway was the most enriched，followed by the cytochrome P450 related pathways.Moreover，several novel pathways，such as circadian rhythm，were significantly enriched despite having no evident associations with the reproductive function.This study is the first to compare the ovary transcriptome data between yak and cattle and to identify and analyze the related differential genes.Our findings provide a theory support for further investigation of the molecular mechanism of yak ovary and offer new insights to understand comprehensively the specificity of yak reproduction.

yak；estrus ovary；transcriptome；comparative analysis；general molecular mechanism

10.11843/j.issn.0366-6964.2016.09.011

2016-03-31

西南民族大學(xué)中央高?；究蒲袠I(yè)務(wù)費專項資金資助(2015NZYTD01)

蘭道亮(1984-)，男，畬族，江西德興人，副研究員，博士，主要從事高原動物遺傳育種研究，E-mail：landaoliang@163.com；Tel：028-85522552

李鍵，教授，E-mail：lijian@swun.cn

S823.8+5.2

A

0366-6964(2016)09-1830-10