谷胱甘肽促進牛體外受精胚胎發育的轉錄組初探

于學穎，郭芹芹，郝海生，孫尉峻，趙學明，朱化彬，楊　凌，杜衛華*

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

?

谷胱甘肽促進牛體外受精胚胎發育的轉錄組初探

于學穎1，2，郭芹芹2，郝海生2，孫尉峻2，趙學明2，朱化彬2，楊凌1*，杜衛華2*

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

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

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

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

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

1材料與方法

1.1主要試劑

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

1.2卵母細胞的采集、體外受精及胚胎培養

從當地屠宰場采集牛卵巢，置于保溫瓶中，3 h內運回實驗室，溫度控制在35～37 ℃。使用37 ℃含有100 IU·mL-1青霉素和100 mg·mL-1硫酸鏈霉素的生理鹽水將卵巢洗凈，將18號針頭與真空蠕動泵連接，抽取3～6 mm的竇狀卵泡，體視顯微鏡下挑選卵丘卵母細胞復合體(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培養箱中成熟培養22 h。

精液從液氮中取出后，空氣中平衡10 s，投入37.5 ℃的水浴鍋中進行解凍；用5 mL洗精液(m-BO medium[16]，10 mmol·L-1咖啡因、4 mg·mL-1BSA)進行稀釋，然后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培養箱中共孵育8 h。用吸卵針脫除外周的顆粒細胞后，受精卵移入前期培養液(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-細胞胚胎移入后期培養液(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%非必需氨基酸)繼續培養，每隔48 h半量換液。培養過程中，分別于受精后48和132 h收集8～16-細胞胚胎和桑椹胚用于測序和熒光定量PCR檢測。處理組胚胎的培養液中添加3 mmol·L-1GSH，對照組則不添加GSH。

1.3構建cDNA文庫并測序

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

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

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

Gene Ontology(GO)是基因功能國際標準分類體系，根據試驗目的篩選差異基因后，研究差異基因在Gene Ontology中的分布狀況。采用GOseq[17]分析方法，準確地計算出GO term被差異基因富集的概率，闡明樣本差異在基因功能上的體現。

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

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

1.6新轉錄本的功能預測

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

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

1.7熒光定量PCR

根據測序結果中基因表達值進行初步篩選，結合 KEGG Pathway 和GO分析結果，確定與GSH代謝、氧化應激反應和細胞凋亡相關的10個基因進行熒光定量PCR反應，驗證測序結果，其中GAPDH為內參基因。參考NCBI基因序列，利用primer 5.0設計引物，引物序列見列表1。

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

表1　基因引物信息

1.8測序結果的驗證

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

2結果

2.1測序數據的質量評估

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

2.2有效序列與?；蚪M比對結果

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

2.3有效序列在牛基因組上的分布情況

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

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

表2　測序數據質量指標匯總

ML_T、ML_C、C8C16_T、C8C16_C分別為GSH處理組桑椹胚、對照組桑椹胚、GSH處理組8～16-細胞期胚胎和對照組8～16-細胞期胚胎。Raw reads為原始序列。Clean reads為有效數據。Error rate是測序錯誤率，Q20、Q30分別指測序錯誤率分別≤1%、≤0.1%的堿基數目比例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　測序數據與牛基因組的比對結果

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

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差異表達基因的篩選

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

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

2.5差異表達基因的功能注釋

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

左邊區域為處理組中特異存在的差異基因，右邊區域為對照組中特異存在的差異基因，重疊區域為在兩個組中均存在的差異基因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差異表達基因的代謝通路和信號轉導途徑分析

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

表4　差異基因KEGG富集結果

2.7新轉錄本的生物信息學分析

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

首先對3個基因的新轉錄本進行ORF預測，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分析結果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 程序對基因編碼序列進行功能結構域的預測，發現在GSK3B(NovelIso00290.1)基因的氨基酸序列中29～243位發現絲氨酸/蘇氨酸蛋白激酶的催化結構域(Serine/Threonine protein kinases catalytic domain，S-TKc)(圖4a)。在GSK3B(NovelIso00290.2)基因氨基酸序列29～46位發現低復雜區(Low-complexity region，LCR)(圖4b)。在SMAD3(NovelIso00499.1)基因的氨基酸序列中分別發現了1～26位的MAD 同源(MAD homology 1，MH1)結構域、28～36位LCR結構域 和 125～296位的dwarfins蛋白家族的B結構域(Domain B in dwarfin family proteins，DWB)(圖4c)。SMAD3(NovelIso00499.2)基因的氨基酸序列中同樣位點發現同樣的結構域。BCL6(NovelIso00109.1)編碼的蛋白序列未發現功能結構域。

a.GSK3B(NovelIso00290.1)基因中預測到的絲氨酸/蘇氨酸蛋白激酶的催化結構域(S-TKc)；b.GSK3B(NovelIso00290.2)基因中預測到的低復雜區(LCR，粉色框)；c.SMAD3(NovelIso00499.1)基因中發現的MAD同源(MH1)結構域、LCR結構域和dwarfins蛋白家族B結構域(DWB)。SMAD3(NovelIso00499.2)基因的氨基酸序列中同樣的結構域。BCL6(NovelIso00109.1)基因編碼的蛋白序列中未發現功能結構域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　基因編碼序列功能結構域預測Fig.4　Prediction of the function domains of genes

表5　定量PCR檢測結果與轉錄組比較結果

2.8基因的定量PCR與測序結果比對

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

3討論

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

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

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

BCL6與基因表達調控和B細胞受體信號通路有關。測序中發現BCL6新轉錄本缺失了第4、5、6、7、9、10外顯子序列，在第8個外顯子之后發現長度為360 bp的新外顯子，說明GSH可能影響BCL6基因表達情況。

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

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

綜上表明，本研究應用RNA-seq技術對GSH處理的8～16-細胞和桑椹期的體外胚胎進行測序，在驗證數據準確性和代表性的基礎上，篩選與GSH處理相關的差異表達基因，對其生物學功能和代謝途徑進行富集分析；發現大量的新轉錄本，也對其RNA和蛋白結構進行生物信息學預測；為從分子水平揭示GSH提高胚胎抗氧化能力和促進胚胎發育的機制奠定研究基礎。

參考文獻(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，FURNUS 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)；基本科研業務費重點項目(2013ywf-zd-2)；家畜胚胎工程與繁殖創新團隊(ASTIP-IAS06-2016)

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

中圖分類號：S823

文獻標志碼：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