于學(xué)穎，郭芹芹，郝海生，孫尉峻，趙學(xué)明，朱化彬，楊　凌，杜衛(wèi)華*

(1．河北工程大學(xué)，邯鄲 056000； 2.中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193)

?

谷胱甘肽促進(jìn)牛體外受精胚胎發(fā)育的轉(zhuǎn)錄組初探

于學(xué)穎1，2，郭芹芹2，郝海生2，孫尉峻2，趙學(xué)明2，朱化彬2，楊凌1*，杜衛(wèi)華2*

(1．河北工程大學(xué)，邯鄲 056000； 2.中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193)

摘要：旨在探究牛體外受精胚胎的培養(yǎng)液中添加3 mmol·L-1谷胱甘肽(Glutathione，GSH)后，8～16-細(xì)胞期和桑椹期的牛體外胚胎在轉(zhuǎn)錄組水平上的變化。收集8～16-細(xì)胞期和桑椹期的體外受精胚胎，構(gòu)建文庫，通過Illumina平臺進(jìn)行測序；篩選差異基因，并進(jìn)行功能注釋和代謝途徑分析；采用定量PCR對10個基因的表達(dá)水平進(jìn)行檢測，以驗證測序結(jié)果的準(zhǔn)確性；在獲得的新轉(zhuǎn)錄本中挑選3組進(jìn)行RNA和蛋白結(jié)構(gòu)的預(yù)測。通過OOSP1、THAP9等10個基因的定量檢測、測序結(jié)果的質(zhì)量評價(堿基錯誤率、Q20、Q30、GC含量)、reads與?；蚪M的比對及其在基因組的分布統(tǒng)計分析，均說明該測序結(jié)果準(zhǔn)確、可靠。測序獲得差異表達(dá)基因4 100個(其中，處理組8～16-細(xì)胞期胚胎和桑椹胚中3 952個，對照組中884個)，已知基因2 225個，未知基因1 875個。這些差異基因主要富集到與ATP合成、呼吸鏈、DNA合成、翻譯過程和物質(zhì)代謝相關(guān)的84條GO terms上，參與細(xì)胞黏附、檸檬酸循環(huán)、谷胱甘肽代謝、溶酶體以及氨基酸代謝等27條通路。另外，與已知的轉(zhuǎn)錄本相比，5個新轉(zhuǎn)錄本中均發(fā)現(xiàn)不同長度的新外顯子；預(yù)測的蛋白結(jié)構(gòu)中除包含各自的保守結(jié)構(gòu)域(絲氨酸/蘇氨酸蛋白激酶的催化結(jié)構(gòu)域、MAD同源結(jié)構(gòu)域)外，還發(fā)現(xiàn)2個新的蛋白結(jié)構(gòu)域(低復(fù)雜區(qū)、dwarfins蛋白家族B結(jié)構(gòu)域)。綜上表明，在培養(yǎng)液中添加GSH能導(dǎo)致胚胎轉(zhuǎn)錄譜的變化，顯著提高胚胎的體外發(fā)育能力。

關(guān)鍵詞：牛體外胚胎；轉(zhuǎn)錄組測序；定量PCR；新轉(zhuǎn)錄本；功能注釋；代謝通路

胚胎發(fā)育是一個精密且高度協(xié)調(diào)的過程，為確保早期胚胎發(fā)育周期的有序進(jìn)行，維持細(xì)胞內(nèi)氧自由基(Reactive oxygen species，ROS)產(chǎn)生和清除的平衡至關(guān)重要。在正常的細(xì)胞代謝過程中，各個系統(tǒng)都會產(chǎn)生ROS并發(fā)揮一定的生物學(xué)作用，包括信號轉(zhuǎn)導(dǎo)及功能調(diào)控等[1]。但ROS過量積累會對細(xì)胞產(chǎn)生有害影響，如DNA損傷、脂質(zhì)中多不飽和脂肪酸氧化、氨基酸氧化、或氧化輔助因子對特定酶的滅活[2]。F.Calzi等[3]研究表明，體外培養(yǎng)時低氧環(huán)境下培養(yǎng)的胚胎發(fā)育能力要顯著優(yōu)于高氧環(huán)境；輸卵管內(nèi)的氧氣濃度低于子宮內(nèi)，約為大氣壓的40%[4-5]。

谷胱甘肽(Glutamine，GSH)是一種廣泛存在于動物體細(xì)胞和配子中的三肽，主要有兩種存在形式：氧化型(GSSG)和還原型(GSH)。還原型GSH會在谷胱甘肽過氧化物酶的作用下被氧化為GSSG，GSSG可通過谷胱甘肽還原酶再次被還原成GSH[6]。GSH能清除細(xì)胞內(nèi)的ROS，并維持胞內(nèi)的氧化還原平衡[7]；對于卵母細(xì)胞，GSH含量隨其成熟的進(jìn)程而逐漸升高[8]；對胚胎，胞內(nèi)GSH和ROS含量能夠決定合子的發(fā)育潛力[7]。因此，GSH常添加于卵母細(xì)胞成熟液和胚胎培養(yǎng)液中以提高胚胎的體外發(fā)育能力[9-10]，尤其是促進(jìn)胚胎的致密化[11]。另外，GSH可通過表觀遺傳重編程調(diào)控胚胎的發(fā)育潛能[12]，降低胞內(nèi)GSH含量，影響DNA的甲基化修飾[13]。關(guān)于GSH促進(jìn)胚胎體外發(fā)育的研究仍停留在氧化還原平衡的層面，對胚胎內(nèi)部發(fā)生的分子生物學(xué)事件、GSH對胚胎整個基因組在RNA水平的影響還沒有得到完全地揭示。筆者前期的研究表明，體外培養(yǎng)液中添加GSH(3 mmol·L-1)能顯著提高牛體外受精胚胎的囊胚率，但對GSH合成、代謝相關(guān)基因和胚胎發(fā)育相關(guān)基因的表達(dá)水平?jīng)]有顯著影響[10]。因此，對GSH處理胚胎的基因組轉(zhuǎn)錄水平進(jìn)行整體性、系統(tǒng)性研究顯得至關(guān)重要。隨著高通量測序的不斷發(fā)展和單細(xì)胞測序技術(shù)的誕生[14-15]，胚胎轉(zhuǎn)錄組研究成為可能。本研究以GSH處理的牛體外受精胚胎為對象，未處理組胚胎為對照，采用單細(xì)胞測序技術(shù)進(jìn)行轉(zhuǎn)錄組測序，對測序結(jié)果進(jìn)行驗證；篩選差異表達(dá)基因，并進(jìn)行基因功能注釋、信號通路及網(wǎng)絡(luò)分析；同時對新轉(zhuǎn)錄本的功能進(jìn)行生物信息學(xué)預(yù)測。

1材料與方法

1.1主要試劑

本研究所用試劑，如無特殊說明均購自西格瑪奧德里奇(上海)貿(mào)易有限公司。熒光定量PCR試劑盒Power SYBR?Green PCR Master Mix和Power SYBR?Green Cells-to-CtTMKit購自美國Life Technologies公司，細(xì)胞裂解試劑盒SMARTer?UltraTMLow RNA Kit for Illumina?Sequencing和cDNA全長擴增試劑盒Advantage?2 PCR Kit購自美國Clontech公司，純化試劑盒AMPure XP beads購自美國Beckman Coulter公司。精液購自山東奧克斯生物技術(shù)有限公司。

1.2卵母細(xì)胞的采集、體外受精及胚胎培養(yǎng)

從當(dāng)?shù)赝涝讏霾杉Ｂ殉玻糜诒仄恐校? h內(nèi)運回實驗室，溫度控制在35～37 ℃。使用37 ℃含有100 IU·mL-1青霉素和100 mg·mL-1硫酸鏈霉素的生理鹽水將卵巢洗凈，將18號針頭與真空蠕動泵連接，抽取3～6 mm的竇狀卵泡，體視顯微鏡下挑選卵丘卵母細(xì)胞復(fù)合體(Cumulus-oocyte complex，COC)，在成熟液(TCM199 medium、10% FBS、0.01 IU·L-1FSH、10 IU·L-1LH、1 μg·mL-1雌二醇、100 ng·mL-1IGF、50 ng·mL-1EGF、100 IU·mL-1青霉素、100 μg·mL-1鏈霉素)中洗3次，放入含750 μL成熟液的四孔板中，每孔放置50枚COC，在38.5 ℃、5% CO2、100%濕度的CO2培養(yǎng)箱中成熟培養(yǎng)22 h。

精液從液氮中取出后，空氣中平衡10 s，投入37.5 ℃的水浴鍋中進(jìn)行解凍；用5 mL洗精液(m-BO medium[16]，10 mmol·L-1咖啡因、4 mg·mL-1BSA)進(jìn)行稀釋，然后500 g離心8 min，棄上清，再次離心洗滌；最后用受精液(m-BO medium，4 mg·mL-1BSA)將精子密度稀釋到2×107個·mL-1。在含有15枚COC的50 μL受精滴中加入50 μL精液，于38.5 ℃、5% CO2、100%濕度的CO2培養(yǎng)箱中共孵育8 h。用吸卵針脫除外周的顆粒細(xì)胞后，受精卵移入前期培養(yǎng)液(NaCl 109.50 mol·L-1、KCl 3.10 mol·L-1、NaHCO326.20 mol·L-1、MgCl2·6H2O 0.80 mol·L-1、KH2PO31.19 mol·L-1、丙酮酸鈉0.4 mol·L-1、葡萄糖1.5 mol·L-1、半乳糖酸鈣5 mol·L-1、BSA 6 mg·mL-1、L-谷氨酰胺1 mol·L-1、2%必需氨基酸、1%非必需氨基酸)。48 h后，將4～8-細(xì)胞胚胎移入后期培養(yǎng)液(NaCl 109.50 mol·L-1、KCl 3.10 mol·L-1、NaHCO326.20 mol·L-1、MgCl2·6H2O 0.80 mol·L-1、KH2PO31.19 mol·L-1、丙酮酸鈉0.4 mol·L-1、葡萄糖1.5 mol·L-1、半乳糖酸鈣5 mol·L-1、10% FBS、L-谷氨酰胺1 mol·L-1、2%必需氨基酸、1%非必需氨基酸)繼續(xù)培養(yǎng)，每隔48 h半量換液。培養(yǎng)過程中，分別于受精后48和132 h收集8～16-細(xì)胞胚胎和桑椹胚用于測序和熒光定量PCR檢測。處理組胚胎的培養(yǎng)液中添加3 mmol·L-1GSH，對照組則不添加GSH。

1.3構(gòu)建cDNA文庫并測序

以3 mmol·L-1GSH處理的8～16-細(xì)胞期(C8C16_T)和桑椹期(ML_T)的胚胎為處理組，GSH未處理的同一階段胚胎(C8C16_T和 ML_C)為對照組進(jìn)行測序。每組20個胚胎，兩個生物學(xué)重復(fù)，共8個cDNA文庫進(jìn)行測序。

胚胎細(xì)胞裂解后進(jìn)行第一鏈cDNA合成，純化后通過LD-PCR擴增得到雙鏈cDNA。定量檢測合格后，使用Covaris系統(tǒng)超聲打斷雙鏈cDNA，獲得的短片段進(jìn)行末端修復(fù)、加A尾并連接測序接頭，選擇大小在200 bp左右的片段PCR富集得到最終的cDNA文庫。建庫后，使用Qubit2.0進(jìn)行初步定量，再使用Agilent 2100檢測文庫的insert size，最后使用qPCR準(zhǔn)確定量文庫的有效濃度(>2 nmol·L-1)，保證文庫質(zhì)量。文庫檢驗合格后，把不同文庫按照有效濃度及目標(biāo)上機數(shù)據(jù)量的需求混池后進(jìn)行HiSeq/MiSeq測序。得到的原始序列中去除帶接頭的、低質(zhì)量、帶引物的序列，得到有效序列。選擇牛參考基因組(版本號：UMD3.1)，通過TOPHAT(v2.0.6)軟件將所有序列比對到參考基因組上。

1.4差異表達(dá)基因的GO功能注釋和富集分析

Gene Ontology(GO)是基因功能國際標(biāo)準(zhǔn)分類體系，根據(jù)試驗?zāi)康暮Y選差異基因后，研究差異基因在Gene Ontology中的分布狀況。采用GOseq[17]分析方法，準(zhǔn)確地計算出GO term被差異基因富集的概率，闡明樣本差異在基因功能上的體現(xiàn)。

1.5差異表達(dá)基因的代謝通路富集分析

KEGG(Kyoto encyclopedia of genes and genomes)是有關(guān)Pathway的主要公共數(shù)據(jù)庫[18]。通過Pathway顯著性富集能確定差異表達(dá)基因參與的最主要生化代謝途徑和信號轉(zhuǎn)導(dǎo)途徑。使用KOBAS(2.0)，以KEGG Pathway為單位，應(yīng)用超幾何檢驗，找出與整個基因組背景相比，在差異表達(dá)基因中顯著性富集的Pathway。

1.6新轉(zhuǎn)錄本的功能預(yù)測

根據(jù)測序結(jié)果，選取BCL6(NovelIso00109.1)、GSK3B(NovelIso00290.1，NovelIso00290.2)和SMAD3(NovelIso00499.1，NovelIso00499.2)等與早期胚胎發(fā)育相關(guān)基因的新轉(zhuǎn)錄本進(jìn)行簡單的功能預(yù)測分析。

通過NCBI的ORF Finder(http：//www.ncbi.nlm.nih.gov/gorf/gorf.html)程序先對3個基因的新轉(zhuǎn)錄本進(jìn)行開放閱讀框(Open reading frame，ORF)分析，獲得ORF的長度、起始和終止位置；再輸入新轉(zhuǎn)錄本的序列信息，獲得其編碼的蛋白序列；采用SMART(http：//smart.embl-heidelberg.de)程序?qū)Φ鞍仔蛄羞M(jìn)行功能結(jié)構(gòu)域預(yù)測，與PFAM蛋白質(zhì)結(jié)構(gòu)域數(shù)據(jù)庫進(jìn)行多重序列比對。

1.7熒光定量PCR

根據(jù)測序結(jié)果中基因表達(dá)值進(jìn)行初步篩選，結(jié)合 KEGG Pathway 和GO分析結(jié)果，確定與GSH代謝、氧化應(yīng)激反應(yīng)和細(xì)胞凋亡相關(guān)的10個基因進(jìn)行熒光定量PCR反應(yīng)，驗證測序結(jié)果，其中GAPDH為內(nèi)參基因。參考NCBI基因序列，利用primer 5.0設(shè)計引物，引物序列見列表1。

收集8～16-細(xì)胞(受精后48 h)期和桑椹胚(受精后132 h)期的牛胚胎，采用7900HT system系統(tǒng) (Applied Biosystems，美國)進(jìn)行檢測。10 μL PCR反應(yīng)體系：2×Master Mix 5 μL，上下游引物(10 μmol·L-1)各0.2 μL，cDNA模板1 μL，Rnase-Free Water 3.6 μL，混勻樣品。PCR反應(yīng)條件為95 ℃預(yù)變性30 s；95 ℃變性5 s，60 ℃退火15 s，72 ℃延伸15 s，40個循環(huán)；陰性對照用0.2 μL Rnase-Free Water代替模板。試驗對每個樣品檢測進(jìn)行3個技術(shù)重復(fù)。通過定量PCR獲得Ct值，采用2-ΔΔCt計算得到基因的相對表達(dá)量。

表1　基因引物信息

1.8測序結(jié)果的驗證

比較熒光定量PCR方法獲得的8～16-細(xì)胞期與桑椹期胚胎的P-value和FC(Fold Change)，與轉(zhuǎn)錄組測序獲得的相應(yīng)結(jié)果進(jìn)行比較。其中，P-value是通過比較兩個時期的基因相對表達(dá)水平計算得到的，顯著性分析按雙側(cè)P<0.05定義，F(xiàn)C=log2(桑椹胚時期基因的相對表達(dá)水平/8～16-細(xì)胞時期基因的相對表達(dá)水平)。全部數(shù)據(jù)采用SAS 9.2軟件進(jìn)行單因素方差分析。

2結(jié)果

2.1測序數(shù)據(jù)的質(zhì)量評估

根據(jù)表2，以ML_T1_1為例，測序得到原始序列(Raw reads)42 426 942條，有效序列(Clean reads)36 950 922條，總的數(shù)據(jù)量(Clean bases)達(dá)3.65G。測序錯誤率(Error rate)為0.03%，低于0.5%，在可接受范圍內(nèi)。測序錯誤率小于1%的堿基(Q20)在整個數(shù)據(jù)中大于96%，錯誤率小于0.1%的堿基(Q30)占整個數(shù)據(jù)92%以上，即準(zhǔn)確測序的堿基比例較高。GC含量39.36%，處于公認(rèn)的40%～60%?？傊?個測序樣品中錯誤率均控制在0.04%以下，Q20、Q30的比例均在88%以上，GC含量也在40%左右，因此該測序數(shù)據(jù)準(zhǔn)確可靠，可用于后續(xù)的分析。

2.2有效序列與?；蚪M比對結(jié)果

測序獲得的有效序列與?；蚪MAC_000160.1比對結(jié)果見表3。以處理組桑椹胚為例，其中約58 000 000條序列比對到參考基因組上，比對率約92%，其中88% 序列唯一比對到參考基因組上。在允許前6～7個堿基錯配的前提下，8個樣品中至少有86%的序列唯一比對到基因組上，其中部分序列與可變剪接有關(guān)，如處理組桑椹胚中約19%序列比對在剪接結(jié)合區(qū)。所以，8個樣品中均有91%以上的序列比對到?；蚪M上。

2.3有效序列在?；蚪M上的分布情況

將比對到基因組上的有效序列按定位區(qū)域外顯子(Exon)、 內(nèi)含子(Intron)和基因間區(qū)(Intergenic)進(jìn)行統(tǒng)計，見圖1。以處理組桑椹胚為例，其中67.0%的序列分布在Exon上，10.9%的序列分布在Intron上，剩余22.1%分布在Intergenic。

圖1　序列在基因組不同區(qū)域的分布Fig.1　Distribution of reads in different regions of bovine genome

表2　測序數(shù)據(jù)質(zhì)量指標(biāo)匯總

ML_T、ML_C、C8C16_T、C8C16_C分別為GSH處理組桑椹胚、對照組桑椹胚、GSH處理組8～16-細(xì)胞期胚胎和對照組8～16-細(xì)胞期胚胎。Raw reads為原始序列。Clean reads為有效數(shù)據(jù)。Error rate是測序錯誤率，Q20、Q30分別指測序錯誤率分別≤1%、≤0.1%的堿基數(shù)目比例ML_T、ML_C、C8C16_T、C8C16_C are groups of morulae with GSH treatment，morulae without treatments，8-16-cell stages embryos with or without GSH treated，respectively.Raw reads are obtained by conversion from raw sequence reads data.Clean reads refers to the valid data which discard the raw reads of low quality sequence reads and fitting contamination reads.Error rate refers to the sequencing error rate.Q20，Q30 are the proportion of bases with sequencing error rate ≤1% or ≤0.1% respectively

表3　測序數(shù)據(jù)與牛基因組的比對結(jié)果

ML_T、ML_C、C8C16_T、C8C16_C分別為GSH處理組桑椹胚、對照組桑椹胚、GSH處理組8～16-細(xì)胞期胚胎和對照組8～16-細(xì)胞期胚胎

ML_T、ML_C、C8C16_T、C8C16_C are groups of morulae with GSH treatment，morulae without treatments，8-16-cell stages embryos with or without GSH treated，respectively

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

來自4個組別、2次重復(fù)的8個樣品的測序結(jié)果共獲得314 135 947條序列，各個重復(fù)間的皮爾遜相關(guān)系數(shù)保證樣品制備和測序技術(shù)的重復(fù)性。分別在8～16-細(xì)胞期胚胎和桑椹胚中檢測到13 489和11 527個基因，檢測到的基因約占?；蚪M中基因(22 000個)的50%。

根據(jù)測序結(jié)果，比較處理組的8～16-細(xì)胞期胚胎和桑椹胚，共發(fā)現(xiàn)3 952個差異基因(2 156個已知基因和1 787個未知基因)；而比較對照組的8～16-細(xì)胞期胚胎和桑椹胚，則發(fā)現(xiàn)884個差異基因(363個已知基因和521個未知基因)，其中736個基因在兩個組中均存在(圖2)。在 4 100個差異基因中，已知基因2 225個，未知基因1 875個。已知差異基因中，60個僅在對照組中表達(dá)，1 862個僅在處理組表達(dá)，303個在兩組中均表達(dá)。處理組中，有下調(diào)基因193個，上調(diào)基因170個；對照組中，1 148個基因下調(diào)，1 017個基因上調(diào)。

2.5差異表達(dá)基因的功能注釋

通過GO分析發(fā)現(xiàn)，獲得的差異基因主要富集到與ATP合成、呼吸鏈、DNA合成、翻譯過程和物質(zhì)代謝相關(guān)的84條GO terms上，如質(zhì)子轉(zhuǎn)運ATP合成酶復(fù)合物(Proton-transporting ATP synthase complex)、NADH脫氫酶復(fù)合物(NADH dehydrogenase complex)、呼吸鏈復(fù)合體I(Respiratory chain complex I)、DNA聚合酶活化(DNA polymerase activity)、線粒體功能(Mitochondrial part)、核糖體代謝過程(Ribosome metabolic process)等，而對照組中的差異基因只富集到與呼吸鏈和ATP相關(guān)的16條terms上，如呼吸鏈(Respiratory chain)，ATP生物合成過程(ATP biosynthetic process)，ATP合成與質(zhì)子轉(zhuǎn)運(ATP synthesis coupled proton transport)等。

左邊區(qū)域為處理組中特異存在的差異基因，右邊區(qū)域為對照組中特異存在的差異基因，重疊區(qū)域為在兩個組中均存在的差異基因The left area represents the differentially expression genes in treated group，the right area represents the differentially expressed genes in control group，and the overlapping area represents common genes of differentially expressed in two groups圖2　差異基因維恩圖Fig.2　Venn diagram of differentially expression genes

2.6差異表達(dá)基因的代謝通路和信號轉(zhuǎn)導(dǎo)途徑分析

通過KEGG Pathway信號通路分析，共發(fā)現(xiàn)27條通路，對照組中有礦物質(zhì)吸收、細(xì)胞黏附、葉酸合成等途徑，而處理組則發(fā)現(xiàn)了檸檬酸循環(huán)(TCA循環(huán))、谷胱甘肽代謝、溶酶體以及氨基酸代謝等多條通路(表4)。

表4　差異基因KEGG富集結(jié)果

2.7新轉(zhuǎn)錄本的生物信息學(xué)分析

轉(zhuǎn)錄組測序共獲得325 585條新轉(zhuǎn)錄本，根據(jù)新轉(zhuǎn)錄本信息和文獻(xiàn)查閱，確定3組與早期胚胎發(fā)育相關(guān)基因的新轉(zhuǎn)錄本進(jìn)行功能預(yù)測分析，包括BCL6(NovelIso00109.1)、GSK3B(NovelIso00290.1，NovelIso00290.2)和SMAD3(NovelIso00499.1，NovelIso00499.2)。

首先對3個基因的新轉(zhuǎn)錄本進(jìn)行ORF預(yù)測，BCL6(NovelIso00109.1)ORF長246 bp(圖3a)，起始密碼子位于 1 bp處，終止密碼子位于246 bp處，編碼81個氨基酸殘基。GSK3B(NovelIso00290.1)ORF長735 bp(圖3b)，起始密碼子位于2 bp處，終止密碼子位于736 bp 處，編碼244個氨基酸殘基。GSK3B(NovelIso00290.2)ORF的起始密碼子位于2 394 bp 處，終止密碼子位于2 870 bp 處，長477 bp 的 ORF(圖 3c)，編碼158個氨基酸殘基。SMAD3(NovelIso00499.1)ORF開始于1 806 bp 處，終止于2 768 bp 處，長963 bp，編碼320個氨基酸殘基(圖 3d)。SMAD3(NovelIso00499.2)也找出一條同樣長度的ORF(圖 3e)，但起始位點不同，起始于662 bp處，終止于1 624 bp處，編碼320個氨基酸殘基。

a～e.BCL6(NovelIso00109.1)、GSK3B(NovelIso00290.1)、GSK3B(NovelIso00290.2)、SMAD3(NovelIso00499.1)和SMAD3(NovelIso00499.2)序列的ORF分析結(jié)果a-e.The ORF analysis of BCL6 (NovelIso00109.1)，GSK3B (NovelIso00290.1)，GSK3B (NovelIso00290.2)，SMAD3 (NovelIso00499.1)，SMAD3(NovelIso00499.2) sequences,respectively圖3　基因序列的 ORF 分析Fig.3　ORF analysis of genes sequences

采用SMART 程序?qū)蚓幋a序列進(jìn)行功能結(jié)構(gòu)域的預(yù)測，發(fā)現(xiàn)在GSK3B(NovelIso00290.1)基因的氨基酸序列中29～243位發(fā)現(xiàn)絲氨酸/蘇氨酸蛋白激酶的催化結(jié)構(gòu)域(Serine/Threonine protein kinases catalytic domain，S-TKc)(圖4a)。在GSK3B(NovelIso00290.2)基因氨基酸序列29～46位發(fā)現(xiàn)低復(fù)雜區(qū)(Low-complexity region，LCR)(圖4b)。在SMAD3(NovelIso00499.1)基因的氨基酸序列中分別發(fā)現(xiàn)了1～26位的MAD 同源(MAD homology 1，MH1)結(jié)構(gòu)域、28～36位LCR結(jié)構(gòu)域 和 125～296位的dwarfins蛋白家族的B結(jié)構(gòu)域(Domain B in dwarfin family proteins，DWB)(圖4c)。SMAD3(NovelIso00499.2)基因的氨基酸序列中同樣位點發(fā)現(xiàn)同樣的結(jié)構(gòu)域。BCL6(NovelIso00109.1)編碼的蛋白序列未發(fā)現(xiàn)功能結(jié)構(gòu)域。

a.GSK3B(NovelIso00290.1)基因中預(yù)測到的絲氨酸/蘇氨酸蛋白激酶的催化結(jié)構(gòu)域(S-TKc)；b.GSK3B(NovelIso00290.2)基因中預(yù)測到的低復(fù)雜區(qū)(LCR，粉色框)；c.SMAD3(NovelIso00499.1)基因中發(fā)現(xiàn)的MAD同源(MH1)結(jié)構(gòu)域、LCR結(jié)構(gòu)域和dwarfins蛋白家族B結(jié)構(gòu)域(DWB)。SMAD3(NovelIso00499.2)基因的氨基酸序列中同樣的結(jié)構(gòu)域。BCL6(NovelIso00109.1)基因編碼的蛋白序列中未發(fā)現(xiàn)功能結(jié)構(gòu)域a.The S-TKc in GSK3B (NovelIso00290.1) predicted；b.The LCR in GSK3B (NovelIso00290.2)；c.The MH1，LCR and DWB in SMAD3 (NovelIso00499.1).The same domains were found in SMAD3 (NovelIso00499.2).However there was no domains found in BCL6 (NovelIso00109.1) sequence圖4　基因編碼序列功能結(jié)構(gòu)域預(yù)測Fig.4　Prediction of the function domains of genes

表5　定量PCR檢測結(jié)果與轉(zhuǎn)錄組比較結(jié)果

2.8基因的定量PCR與測序結(jié)果比對

通過比較整體的變化趨勢和基因表達(dá)差異倍數(shù)來判斷定量PCR與轉(zhuǎn)錄組測序結(jié)果的一致性。當(dāng)基因的 FC值在定量PCR和測序結(jié)果中都>1或<1時，說明該基因在兩種檢測中均為上調(diào)或下調(diào)，表達(dá)量增加或降低，這兩種情況均表明定量PCR檢測結(jié)果與轉(zhuǎn)錄組測序是一致的。經(jīng)過比對，所選10個基因的表達(dá)倍數(shù)及整體變化趨勢與測序結(jié)果均保持一致，表明轉(zhuǎn)錄組測序結(jié)果的可靠性。

3討論

在胚胎培養(yǎng)液中添加GSH能提高牛體外受精胚胎的囊胚率和發(fā)育能力[10，19]，然而其分子機制并沒有完全揭示。本研究采用單細(xì)胞測序技術(shù)，揭示GSH添加后牛胚胎整個基因組的基因表達(dá)模式變化，獲得大量的差異表達(dá)基因，其中14個差異基因富集到谷胱甘肽的代謝通路上，此外差異基因還富集到26條其他相關(guān)信號傳導(dǎo)和代謝通路上，如檸檬酸循環(huán)(TCA循環(huán))、溶酶體以及氨基酸代謝等?？梢娫谂囵B(yǎng)液中添加GSH后，除自身的代謝途徑外還引發(fā)了其他的對8～16-細(xì)胞期胚胎發(fā)育至桑椹胚有重要作用的生物學(xué)過程，促進(jìn)胚胎的體外發(fā)育，顯著提高體外胚胎的囊胚率。

轉(zhuǎn)錄組測序因其具有全方位、高通量等優(yōu)點而廣泛應(yīng)用于動植物研究中，尤其在基因功能及其調(diào)控網(wǎng)絡(luò)的揭示方面有著重要地位，然而單細(xì)胞測序在動物胚胎發(fā)育調(diào)控方面的研究還在起步階段。2011年果蠅單囊胚的轉(zhuǎn)錄組測序揭示基因mRNA的表達(dá)豐度與性別決定機制的關(guān)系[20]。人類早期胚胎的轉(zhuǎn)錄組測序闡明了外胚層細(xì)胞和胚胎干細(xì)胞分子調(diào)控機制的差異，證實胚胎干細(xì)胞的多能性，為人類早期胚胎發(fā)育和干細(xì)胞研究提供借鑒[21]。對于家畜牛單胚胎的轉(zhuǎn)錄組測序首次發(fā)現(xiàn)等位基因的序列存在差異，同時獲得大量的新SNP 位點[22]；而且，附植前的胚胎基因組在轉(zhuǎn)錄水平上發(fā)生巨大變化[23]。綜上表明，轉(zhuǎn)錄組測序技術(shù)在基因表達(dá)、基因變異和早期胚胎發(fā)育的轉(zhuǎn)錄調(diào)控中有著廣泛的應(yīng)用前景。

本研究中所有測序結(jié)果的各項質(zhì)量指標(biāo)均在正常范圍內(nèi)。測序樣品中堿基錯誤率均控制在0.04%以下，Q20、Q30的比例均在88%以上，GC含量也在40%左右。測序結(jié)果與牛基因組比對發(fā)現(xiàn)，8個樣品中均有91%以上的序列比對到?；蚪M上，其中唯一比對的序列占總數(shù)的86%，說明測序數(shù)據(jù)覆蓋整個?；蚪M。另外，至少有67%以上的序列分布到外顯子區(qū)，少數(shù)的則定位到內(nèi)含子和基因間區(qū)，可能是出現(xiàn)的新轉(zhuǎn)錄本或非編碼RNA，或者是目前已知基因的內(nèi)含子中尚未發(fā)現(xiàn)的新外顯子。然而，轉(zhuǎn)錄組測序數(shù)據(jù)量巨大，每一個生物信息學(xué)分析步驟都可能出現(xiàn)假陽性結(jié)果，所以有必要采用定量PCR對其進(jìn)行驗證，這也是評定測序結(jié)果的重要標(biāo)志之一。本研究對10個基因的RNA表達(dá)水平進(jìn)行定量檢測，得到的結(jié)果均與測序結(jié)果相吻合，充分說明測序的可靠性和準(zhǔn)確性。

BCL6與基因表達(dá)調(diào)控和B細(xì)胞受體信號通路有關(guān)。測序中發(fā)現(xiàn)BCL6新轉(zhuǎn)錄本缺失了第4、5、6、7、9、10外顯子序列，在第8個外顯子之后發(fā)現(xiàn)長度為360 bp的新外顯子，說明GSH可能影響B(tài)CL6基因表達(dá)情況。

GSK3B可調(diào)節(jié)牛卵母細(xì)胞減數(shù)分裂進(jìn)程，參與卵母細(xì)胞和卵丘細(xì)胞中MAPK信號通路[24]。本研究發(fā)現(xiàn)GSK3B(NovelIso00290.1)轉(zhuǎn)錄本比已知序列多出一個外顯子(長114 bp)，GSK3B(NovelIso00290.2)轉(zhuǎn)錄本中也發(fā)現(xiàn)兩個新外顯子(長度分別為2 187和689 bp)。在蛋白水平，兩個新轉(zhuǎn)錄本都存在S-TKc元件，它是GSK3蛋白家族中存在的一個保守結(jié)構(gòu)域，可作為分子結(jié)合位點(如ATP的結(jié)合位點)，或氨基酸多肽的結(jié)合位點參與蛋白磷酸化過程，與多數(shù)細(xì)胞活性有著密切的聯(lián)系。LCR也是在新轉(zhuǎn)錄本中發(fā)現(xiàn)的由簡單氨基酸組成的蛋白質(zhì)序列[25]，雖然其功能模型尚未清楚，但有研究表明它在生物學(xué)功能中發(fā)揮關(guān)鍵作用[26-28]。A.Coletta研究表明，含有LCR結(jié)構(gòu)的蛋白在生物總庫的交互數(shù)據(jù)集(The biological general repository for interaction datasets，the BioGrid)中分布更加廣泛[25]。

SMAD3是SMAD蛋白家族成員，該家族蛋白根據(jù)功能分為受體調(diào)節(jié)型、通用型和抑制型3類，SMAD3則是受體調(diào)節(jié)型。測序結(jié)果中SMAD3(NovelIso00499.1)新轉(zhuǎn)錄本中發(fā)現(xiàn)長度為1 695和193 bp的新外顯子，而SMAD3(NovelIso00499.2)轉(zhuǎn)錄本中發(fā)現(xiàn)長度為551和193 bp的新外顯子。新轉(zhuǎn)錄本的蛋白結(jié)構(gòu)預(yù)測后發(fā)現(xiàn)MH1和DWB結(jié)構(gòu)域，前者是SMAD3蛋白中存在的一個保守結(jié)構(gòu)域，主要參與TGF-β信號通路的傳導(dǎo)[29-30]；后者是dwarfins蛋白家族共有結(jié)構(gòu)域，與哺乳動物中TGF-β磷酸化和細(xì)胞生長的調(diào)控有關(guān)[31]。在核因子I(NF-1)或CCAAT盒結(jié)合的轉(zhuǎn)錄因子(CTF)中也發(fā)現(xiàn)DWB結(jié)構(gòu)域的存在。

綜上表明，本研究應(yīng)用RNA-seq技術(shù)對GSH處理的8～16-細(xì)胞和桑椹期的體外胚胎進(jìn)行測序，在驗證數(shù)據(jù)準(zhǔn)確性和代表性的基礎(chǔ)上，篩選與GSH處理相關(guān)的差異表達(dá)基因，對其生物學(xué)功能和代謝途徑進(jìn)行富集分析；發(fā)現(xiàn)大量的新轉(zhuǎn)錄本，也對其RNA和蛋白結(jié)構(gòu)進(jìn)行生物信息學(xué)預(yù)測；為從分子水平揭示GSH提高胚胎抗氧化能力和促進(jìn)胚胎發(fā)育的機制奠定研究基礎(chǔ)。

參考文獻(xiàn)(References)：

[1]DEVASAGAYAM T P，TILAK J C，BOLOOR K K，et al.Free radicals and antioxidants in human health：current status and future prospects[J].JAssocPhysiciansIndia，2004，52：794-804.

[2]BROOKER R J.Genetics：analysis and principles[M].New York：McGraw-Hill Science，2011.

[3]CALZI F，PAPALEO E，RABELLOTTI E，et al.Exposure of embryos to oxygen at low concentration in a cleavage stage transfer program：reproductive outcomes in a time-series analysis[J].ClinLab，2012，58(9-10)：997-1003.

[4]FISCHER B，BAVISTER B D.Oxygen tension in the oviduct and uterus of rhesus monkeys，hamsters and rabbits[J].JReprodFertil，1993，99(2)：673-679.

[5]LEESE H J.Metabolic control during preimplantation mammalian development[J].HumReprodUpdate，1995，1(1)：63-72.

[6]LUBERDA Z.The role of glutathione in mammalian gametes[J].ReprodBiol，2005，5(1)：5-17.

[7]TAKAHASHI M，NAGAI T，OKAMURA N，et al.Promoting effect of beta-mercaptoethanol oninvitrodevelopment under oxidative stress and cystine uptake of bovine embryos[J].BiolReprod，2002，66(3)：562-567.

[8]DE MATOS D G，F(xiàn)URNUS C C.The importance of having high glutathione (GSH) level after bovineinvitromaturation on embryo development effect of beta-mercaptoethanol，cysteine and cystine[J].Theriogenology，2000，53(3)：761-771.

[9]GARDINER C S，SALMEN J J，BRANDT C J，et al.Glutathione is present in reproductive tract secretions and improves development of mouse embryos after chemically induced glutathione depletion[J].BiolReprod，1998，59(2)：431-436.

[10]SUN W J，PANG Y W，LIU Y，et al.Exogenous glutathione supplementation in culture medium improves the bovine embryo development afterinvitrofertilization[J].Theriogenology，2015，84(5)：716-723.

[11]PERREAULT S D，BARBEE R R，SLOTT V L.Importance of glutathione in the acquisition and maintenance of sperm nuclear decondensing activity in maturing hamster oocytes[J].DevBiol，1988，125(1)：181-186.

[12]TIMME-LARAGY A R，GOLDSTONE J V，IMHOFF B R，et al.Glutathione redox dynamics and expression of glutathione-related genes in the developing embryo[J].FreeRadicBiolMed，2013，65：89-101.

[13]LERTRATANANGKOON K，WU C J，SAVARAJ N，et al.Alterations of DNA methylation by glutathione depletion[J].CancerLett，1997，120(2)：149-156.

[14]PIRAS V，TOMITA M，SELVARAJOO K.Transcriptome-wide variability in single embryonic development cells[J].SciRep，2014，4：7137.

[15]DE DOMENICO E，OWENS N D，GRANT I M，et al.Molecular asymmetry in the 8-cell stage Xenopus tropicalis embryo described by single blastomere transcript sequencing[J].DevBiol，2015，408(2)：252-268.

[16]BRACKETT B G，OLIPHANT G.Capacitation of rabbit spermatozoainvitro[J].BiolReprod，1975，12(2)：260-274.

[17]YOUNG M D，WAKEFIELD M J，SMYTH G K，et al.Gene ontology analysis for RNA-seq：accounting for selection bias[J].GenomeBiol，2010，11(2)：R14.

[18]KANEHISA M，ARAKI M，GOTO S，et al.KEGG for linking genomes to life and the environment[J].NucleicAcidsRes，2008，36：D480-D484.

[19]ALI K，ABD-ALLAH S，KANDEIL M，et al.Impact of exogenous glutathione duringinvitromaturation，fertilization and development of buffalo oocytes[M].Fifth Scientific Conference，Ras Elbar，2007：1-12.

[20]LOTT S E，VILLALTA J E，SCHROTH G P，et al.Noncanonical compensation of zygotic X transcription in early Drosophila melanogaster development revealed through single-embryo RNA-seq[J].PLoSBiol，2011，9(2)：e1000590.

[21]YAN L，YANG M，GUO H，et al.Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells[J].NatStructMolBiol，2013，20(9)：1131-1139.

[22]CHITWOOD J L，RINCON G，KAISER G G，et al.RNA-seq analysis of single bovine blastocysts[J].BMCGenomics，2013，14：350.

[23]JIANG Z，SUN J，DONG H，et al.Transcriptional profiles of bovineinvivopre-implantation development[J].BMCGenomics，2014，15：756.

[24]UZBEKOVA S，SALHAB M，PERREAU C，et al.Glycogen synthase kinase 3B in bovine oocytes and granulosa cells：possible involvement in meiosis duringinvitromaturation[J].Reproduction，2009，138(2)：235-246.

[25]COLETTA A，PINNEY J W，SOLIS D Y，et al.Low-complexity regions within protein sequences have position-dependent roles[J].BMCSystBiol，2010，4：43.[26]FONDON J W III，GARNER H R.Molecular origins of rapid and continuous morphological evolution[J].ProcNatlAcadSciUSA，2004，101(52)：18058-18063.

[27]VERSTREPEN K J，JANSEN A，LEWITTER F，et al.Intragenic tandem repeats generate functional variability[J].NatGenet，2005，37(9)：986-990.

[28]PHATNANI H P，GREENLEAF A L.Phosphorylation and functions of the RNA polymerase II CTD[J].GenesDev，2006，20(21)：2922-2936.

[29]ATTISANO L，LEE-HOEFLICH S T.The Smads[J].GenomeBiol，2001，2(8)：3010.

[30]BROWN K A，PIETENPOL J A，MOSES H L.A tale of two proteins：differential roles and regulation of Smad2 and Smad3 in TGF-beta signaling[J].JCellBiochem，2007，101(1)：9-33.

[31]YINGLING J M，DAS P，SAVAGE C，et al.Mammalian dwarfins are phosphorylated in response to transforming growth factor beta and are implicated in control of cell growth[J].ProcNatlAcadSciUSA，1996，93(17)：8940-8944.

(編輯程金華)

doi:10.11843/j.issn.0366-6964.2016.07.008

收稿日期：2016-01-22

基金項目：國家科技支撐項目(2012BAD12B01-2)；基本科研業(yè)務(wù)費重點項目(2013ywf-zd-2)；家畜胚胎工程與繁殖創(chuàng)新團隊(ASTIP-IAS06-2016)

作者簡介：于學(xué)穎(1989-)，女,河北廊坊人，碩士，主要從事家畜胚胎工程研究，E-mail: xueyingyu2015bj@163.com *通信作者：杜衛(wèi)華,博士,副研究員,E-mail：dwh@iascaas.net.cn；楊凌,博士,副教授,E-mail:yangling@hebeu.edu.cn

中圖分類號：S823

文獻(xiàn)標(biāo)志碼：A

文章編號:0366-6964(2016)07-1363-10

Transcriptome of Bovine IVF Embryos Treated with Glutathione

YU Xue-ying1，2，GUO Qin-qin2，HAO Hai-sheng2，SUN Wei-jun2，ZHAO Xue-ming2，ZHU Hua-bin2，YANG Ling1*，DU Wei-hua2*

(1.CollegeofAgriculture，HebeiUniversityofEngineering，Handan056000，China；2.InstituteofAnimalScience，ChineseAcademyofAgriculturalSciences，Beijing100193，China)

Abstract:To elucidate the effects of exogenous glutathione (GSH) on embryos transcriptome，bovine IVF embryos at 8-16-cell and morula stage were collected for high-throughput RNA sequencing (RNA-seq).Quantitative determination of genes expression was used to verify the RNA-seq data.Gene ontology and KEGG pathway enrichment analysis of differential expression genes were performed.The accuracy and validity of RNA-seq data were confirmed with combination of error rate of base，mapping of reads to bovine genome，distribution region analysis and consistency between qPCR and RNA-seq results of ten genes.A total of 4 100 genes were differentially expressed between embryos at two stages and 3 952 genes were identified as expressed in embryos treated with GSH.A high enrichment of genes involved in 84 GO terms，such as ATP synthesis，respiratory chain，DNA synthesis and translation.Pathway analysis revealed 27 pathways involved in cell adhesion，TCA cycle，GSH metabolism and amino acid metabolism.Additionally，2 known conserved domains (S-TKc，MH1) and 2 new domains (LCR，DWB) were found in 5 new transcripts by prediction，compared with known transcripts.To sum up，GSH treatment can result into the comprehensive change at transcriptional level of bovine IVF embryos and improve their development consequently.

Key words:bovine IVF embryo；RNA-Seq；qRT-PCR；novel isoform；GO analysis；KEGG pathway analysis