基于轉錄組數據挖掘外源褪黑激素影響水貂卵巢發育的分子機制

高婭薇,彭 弟,孫朝陽,晏子越,崔 凱,馬澤芳

(青島農業大學動物科技學院,青島 266109)

水貂(Mustelavison)屬哺乳綱(Mammalia)、食肉目(Carnivoraes)、鼬科(Mustelidae)、鼬屬(Mustelia)[1],是一種季節性繁殖的珍貴毛皮動物,其毛皮生長和繁殖活動明顯受光周期調控,主要與光周期調節松果體合成褪黑激素(MLT)的分泌實現調控作用有關[2]。目前水貂養殖生產中應用MLT埋植劑(主要成分MLT,6～7月埋植)促進水貂毛皮早熟,一般提前30～85 d[3-4],其作用機理是模擬短日照作用,提高體內MLT水平,進而促使冬毛早熟。在應用過程中發現MLT埋植(7月初)后會使雌、雄水貂發情時間分別相對提前至10、11月份[5],其作用機理與MLT影響綿羊[6-7]、梅花鹿[8]等動物生殖活動的機理相近,即MLT通過下丘腦-垂體-性腺(HPGA)軸調節生殖激素分泌,調控水貂生殖活動,但其分子機制尚不明確。轉錄組學可研究基因轉錄情況及轉錄調控規律,揭示基因轉錄、轉錄調控規律過程中的分子機制[9-10],已有研究報道通過對卵巢組織進行轉錄組學分析,探究山羊[11-14]、綿羊[15-17]、牛[18]、豬[19-22]等卵巢發育的分子調控機制,但有關MLT影響水貂卵巢發育機制的轉錄組學研究尚未見報道。

本團隊晏子越等[23]前期對MLT埋植劑埋植172 d和未埋植MLT埋植劑的配種準備期水貂卵巢進行形態學、組織學觀察,但未對其進行轉錄組學分析,因此本研究在確定MLT影響水貂卵巢表型發育的基礎上,進行轉錄組學分析,篩選出MLT影響卵巢發育的候選基因,探究MLT對水貂卵巢發育影響的分子機制。

1 材料與方法

1.1 試驗設計與飼養管理

于山東省諸城市(119°E,35°N)某水貂養殖場隨機選取3月齡健康、體重均勻(870.00±67.50)g的白色雌性水貂200只,隨機均分為試驗組和對照組,每組100個重復,每個重復1只。試驗于2020年7月7日至2020年12月26日進行,共172 d。于試驗開始當天利用MLT皮下埋植器在試驗組水貂頸部皮下埋植1粒MLT埋植劑(含MLT 18 mg·粒-1),對照組埋植1粒不含MLT的空粒。于試驗結束當天隨機選取試驗組和對照組水貂各4只,采集卵巢,進行轉錄組學分析。

水貂飼養于南北走向、東西通風的棚舍內,采用單籠飼養,棚舍與籠箱均保持清潔干凈,專人管理,定時飼喂,日糧組成及營養水平見表1。

表1 日糧組成及營養水平(干物質基礎)

1.2 樣品采集與處理

于試驗結束當天,隨機選取試驗組和對照組水貂各4只,CO2處死后立即剖腹采集右側卵巢(共計8枚)放入做好標記的凍存管內,置于液氮(-196 ℃)中速凍,待轉錄組分析使用。

1.3 卵巢轉錄組測序

取液氮(-196 ℃)中待測卵巢,采用Trizol法[24]提取總RNA后,利用Nano Photometer NP80分析儀和Agilent 2100檢測其RNA純度、濃度和完整性,所有樣品檢測合格后,按照Hyde等[25]的方法進行cDNA文庫的構建,并用qPCR法定量cDNA文庫有效濃度,保證文庫有效濃度>2 nmol·L-1。然后用SBS技術在Illumina測序平臺進行RNA-seq,得到原始測序數據(Raw Reads)。

1.4 生物信息學分析

1.4.1 原始測序數據處理、質量評估與序列比對 通過Faslx-toolkit軟件對Raw Reads去除有接頭及低質量的Reads,過濾后得到Clean Reads,利用HISAT2軟件與美洲水貂(Neovisonvison)參考基因組序列(https:∥asia.ensembl.org/Neovison_vison/Info/Index,基因組版本號:NNQGG.v01(GCA_900 108 605.1))進行比對,得到Mapped Reads,并統計出各樣品測序結果與參考基因組的比對效率和Mapped Reads在參考基因組外顯子、內含子和基因間區的數目。

1.4.2 SNP分析和可變剪切分析 使用GATK軟件對樣本數據進行變異位點分析,并通過ASprofile軟件對各樣品中存在的AS類型及其表達量進行統計,利用rMATS軟件對各樣品的差異AS事件進行統計。

1.4.3 新基因挖掘及功能注釋 基于所選參考基因組序列,使用StringTie軟件對Mapped Reads進行拼接,并與原有的基因組注釋信息進行比較,尋找原來未被注釋的轉錄區,發掘水貂的新轉錄本和新基因,其中位于基因間區,且長度大于50 nt并含有2個以上外顯子的潛在轉錄本即為新轉錄本,通過DIAMOND、InterProScan、HMMER等軟件將挖掘出的新轉錄本與COG、GO、KEGG、KOG、Pfam、Swiss-Prot 、EggNOG、NR、TrEMBL數據庫進行比對,獲得新轉錄本的注釋信息。

1.4.4 差異表達基因篩選與分析 使用String Tie軟件對Mapped Reads進行組裝和定量,采用FPKM計算基因的表達量,通過DEAeq2軟件包分析差異表達基因,以差異倍數(Fold Change)≥1.5且錯誤發生率(FDR)<0.01為標準篩選差異表達基因。最后通過GO和KEGG數據庫對差異表達基因進行功能分類和富集通路注釋。

1.5 實時熒光定量PCR

從差異表達基因中隨機挑選6個基因,采用qPCR法對測序結果進行驗證。將樣品中提取的總RNA按照試劑盒操作反轉錄得到cDNA。利用Primer premier 5軟件設計引物,引物信息見表2。以cDNA為模板,GAPDH為內參基因。擴增體系:cDNA 5 μL,上、下游引物各1 μL,2×SYBR real-time PCR premixture 5 μL ,RNase free water 3 μL。擴增程序:95 ℃預變性10 min;95 ℃變性10 s,60 ℃退火15 s,72 ℃延伸20 s,共40個循環;4 ℃保存。然后根據目的基因和內參基因的 Ct 值,用2-ΔΔCt法計算基因的相對表達量。Ct值表示在PCR反應中熒光信號達到規定的閾值時經歷的循環數,ΔΔCt計算公式如下:ΔΔCt=(Ct目的基因-Ct內參基因)試驗組-(Ct目的基因-Ct內參基因)對照組。并利用SPSS 20.0中的獨立樣本t檢驗進行方差分析,數據以“平均值±標準誤”表示,P<0.01代表差異極顯著,P<0.05代表差異顯著為判斷標準。

表2 qPCR驗證差異表達基因引物信息

1.6 候選基因篩選

將GO功能分類中與生殖過程相關的差異表達基因和KEGG通路分析中與繁殖通路相關的差異表達基因結合分析,篩選得到候選基因。

2 結 果

2.1 總RNA質量與測序數據質量控制

兩組卵巢組織總RNA的RIN值>7.0,28 S/18 S和OD260 nm/OD280 nm均在1.8～2.0之間,說明RNA濃度、純度、完整性均符合后續轉錄組測序和qPCR驗證試驗的需求。

由表3可知,轉錄組測序數據經質控后共獲得52.47 Gb Clean Reads,各樣品Clean Reads均達到5.96 Gb,各樣品測序數據堿基質量值Q30和Q20的堿基百分比分別在91%和96%以上,GC含量均在45%以上;由表4可知,兩組樣品Reads與參考基因組的比對效率均超過89%,說明整體測序數據質量良好,可用于后續分析。

表3 轉錄組測序數據質量統計表

表4 測序數據序列比對結果統計表

2.2 SNP分析結果

由表5可知,試驗組和對照組的SNPs位點均在120 000個以上,且均主要出現在基因區,且各組SNP突變均包括轉換型、顛換型、雜合型3種類型,其中轉換型所占比例最高,說明SNP位點主要與堿基轉換有關。

表5 SNPs統計結果

2.3 可變剪切分析結果

由表6可知,差異表達基因發生外顯子跳躍(SE)事件的數量最多,說明試驗組與對照組卵巢基因的差異可變剪切主要與SE事件有關。

表6 可變剪切事件結構統計表

2.4 新基因分析結果

本試驗共獲得9 631個新基因,其注釋結果見表7。由該表可知,除COG數據庫外,GO、KEGG、KOG、Pfam、Swiss-Prot、TrEMBL、Egg NOG、NR數據庫注釋到的轉錄本均超過300個,其中TrEMBL和NR數據庫注釋到的新基因分別為1 323個、1 317個,注釋比例明顯高于其他數據庫,可為后續水貂新基因挖掘及其功能分析提供參考。

表7 新基因功能注釋結果統計

2.5 差異表達基因分析結果

兩組樣品中共得到21 270個基因,由圖1可知,兩組共206個差異表達基因(其中36個基因未得到注釋),69個為上調基因,137個為下調基因。

圖中的每一個點表示一個基因,橫坐標表示某一個基因在兩樣品中表達量差異倍數的對數值,縱坐標表示基因表達量的統計學顯著性的負對數值Each point in the graph represents a gene, the abscissa represents the logarithm of the differential multiples of a gene expression in the two samples, and the ordinate represents the negative logarithm of the statistical significance of gene expression

對差異表達基因進行GO分類,共計2 167、1 481、5 727個未被注釋的基因分別富集到BP、CC和MF三個分支中。如圖2所示,被注釋的差異表達基因分別參與了BP、CC和MF的23、17和18個功能亞分類。其中與生殖過程相關的差異表達基因有7個,分別為細胞角蛋白19(KRT19)、谷胱甘肽硫轉移酶Theta1(GSTT1)、5′胺基乙酰丙酸合酶1(ALAS1)、核受體亞科4A組成員1(NR4A1)、AP-1轉錄因子亞基(JUNB)、精子自身抗原蛋白17(SPA17)、鋅指蛋白296(ZNF296),推測MLT對水貂卵巢發育的影響可能與生殖過程相關的差異表達基因有關。

共159個差異表達基因顯著富集在131條通路上,對差異表達基因KEGG的注釋結果按照KEGG通路類型進行分類,其中差異最顯著的前20條通路如圖3所示。在131條通路中有5條與繁殖相關的通路,分別為細胞色素P450對外來生物的代謝(metabolism of xenobiotics by cytochrome P450)、甘氨酸、絲氨酸和蘇氨酸的代謝(glycine,serine and threoine metabolism)PI3K-Akt信號通路(PI3K-Akt signaling pathway)、MAPK信號通路(MAPK signaling pathway)、E2信號通路(estrogen signaling pathway)。由表8可知,富集于這些通路的差異表達基因有7個,推測MLT對水貂卵巢發育的影響可能與差異表達基因富集在與繁殖相關的通路有關。

表8 與繁殖過程相關的KEGG通路

2.6 候選基因篩選結果

將GO功能分類中與生殖過程相關的7個差異表達基因和KEGG通路注釋中與繁殖通路相關的7個差異表達基因結合分析,共有4個差異表達基因同時富集于生殖過程和繁殖相關的信號通路,分別為細胞角蛋白19(KRT19)、谷胱甘肽硫轉移酶Theta1(GSTT1)、5′胺基乙酰丙酸合酶1(ALAS1)、核受體亞科4A組成員1(NR4A1),將其作為候選基因,推測這4個基因可能參與MLT促進水貂卵巢發育的調控過程。

2.7 qPCR驗證結果

由圖4可知,在差異表達基因中隨機篩選的6個基因的qPCR表達趨勢與轉錄組分析結果一致,表明轉錄組分析結果可靠。

圖4 基因qPCR表達量統計圖Fig.4 Expression level statistical diagram of genes by qPCR

3 討 論

水貂毛絨品質性能良好,其皮張是世界上珍貴裘皮,已經成為國際裘皮貿易的三大支柱產品之一[26],具有較高的經濟價值,因此水貂養殖已發展成為特種經濟動物養殖的熱點項目,是農民增收、致富的重要經濟來源[27],也是畜牧業的重要組成部分。水貂是一種受季節性繁殖調控的動物,其生殖腺的發育僅在特定時間段內進行,這一現象是長期自然選擇的結果[28],也是導致水貂繁殖周期較長、繁殖效率較低、繁殖技術發展較慢的主要原因。這種季節性繁殖主要與光周期調節松果體合成和分泌MLT的節律有關。宮洪杰等[29]發現,水貂埋植MLT埋植劑對發情期水貂卵巢發育有促進作用。據趙英等[30]報道,MLT可作用于HPGA,調控下丘腦促性腺激素的分泌,進而對生殖系統的生長發育產生影響,同時與卵巢、睪丸等組織上的MLT受體相結合,直接參與促黃體生成素(LH)、促卵泡生長素(FSH)及性腺激素等的合成和分泌,但有關MLT對水貂卵巢發育影響的轉錄組學分析尚未見報道。因此,本研究以MLT埋植劑埋植172 d和未埋植MLT埋植劑的配種準備期水貂卵巢為研究對象,探討MLT影響水貂卵巢發育的分子機制。

SNP為單核苷酸多態性,在基因組中廣泛存在[31],主要包括轉換、顛換、雜合等類型。本研究中,SNP位點突變類型為轉換型占比最高,高于顛換型和雜合型。這與周明夏等[32]測序研究產前和產蛋高峰期泰州鵝卵巢組織轉換型SNP比例高于顛換型和張蕾等[33]測序高郵鴨卵巢組織樣品的SNP轉換高于顛換均一致,說明影響差異SNP位點的主要是堿基轉換。

可變剪切能夠賦予同一基因更多的表達可能性,從而產生多種不同的蛋白質,這種機制使得基因組在表達時能夠更加靈活,以適應生物體在不同環境下的需求[34]。本試驗5類剪切事件中SE事件所占比例最大。朱志明等[35]在統計不同產蛋期金定鴨卵巢組織的AS事件時發現,開產期與產蛋高峰期的卵巢差異剪切基因發生SE事件的頻率最高,Miao等[36]也提出,在寒羊和多賽特綿羊的卵巢組織中SE事件是最常發生的可變剪切事件類型,由此說明動物卵巢基因的差異可變剪切現象主要與SE事件相關。由此推斷,SE事件可能參與了動物卵巢的發育調控過程。

本研究共得到21 270個基因,9 631個新基因,兩組差異表達基因共206個,其中篩選到與卵巢發育相關的差異表達基因4個,分別為KRT19、GSTT1、ALAS1、NR4A1。本研究發現KRT19基因富集于E2信號通路。這與李德河等[5]對水貂埋植MLT后體內雌二醇(E2)濃度明顯提高的研究結果一致,由此推測KRT19基因功能可能與E2分泌相關,MLT增加了KRT19基因的表達量,這一基因表達的蛋白參與E2信號通路,從而促進E2的分泌。本研究發現,GSTT1基因富集于細胞色素P450對外來生物的代謝通路上。Sahmi等[37]、Kandiel等[38]均報道稱該通路參與了固醇和類固醇激素的生物合成,尤其是性激素的生物合成,說明該通路的GSTT1基因也參與了E2的分泌。據Fehrenbach等[39]報道E2是卵泡成活因子之一,通過激活其在顆粒細胞上的受體(ERα和ERβ)發揮生理功能,刺激卵泡發育至早期竇卵泡階段,由此說明MLT對卵泡發育的作用可能與KRT19基因在E2信號通路和GSTT1基因在細胞色素P450對外來生物的代謝通路上表達刺激E2分泌有關。

本研究發現,ALAS1基因富集在甘氨酸、絲氨酸和蘇氨酸的代謝通路上。蘭道亮等[40]稱,甘氨酸、絲氨酸和蘇氨酸的代謝通路為牦牛卵母細胞成熟的關鍵調控通路。前人對ALAS1基因功能的研究中發現,ALAS1基因在卵巢卵泡發育及排卵過程中具有重要作用,吳華莉等[41]對豬在血清促性腺激素處理后的研究發現,ALAS1基因在卵泡發育2 h和排卵過程12 h時表達量高。由此說明,MLT對卵泡的促進作用可能與ALAS1基因在甘氨酸、絲氨酸和蘇氨酸代謝通路中的調控有關。

本研究發現,NR4A1基因富集在PI3K-Akt信號通路上。O′Shea等[42]研究稱,PI3K-Akt信號通路在卵母細胞成熟中起了關鍵的調控作用;在前人對NR4A1功能的研究中,馬紅梅和劉嘉茵[43]和李梅等[44]在人和小鼠卵巢卵泡中均檢測到NR4A1表達,且小鼠成熟期卵泡和人黃體組織中NR4A1水平顯著高于生長期卵泡,提示NR4A1參與了卵泡的成熟與細胞分化。由此推斷,MLT對卵泡的促進作用可能與其調控NR4A1基因在PI13-AKT信號通路中表達從而促進卵泡成熟和發育有關。

本研究發現,NR4A1基因也存在于MAPK信號通路。韓樹標等[45]稱MAPK信號通路在卵母細胞成熟方面發揮重要作用,其在卵母細胞成熟時發生磷酸化而被激活,MLT對卵泡的促進作用可能也與NR4A1基因在MAPK信號通路中的表達有關。

4 結 論

本研究通過對MLT埋植劑埋植172 d和未埋植MLT埋植劑的配種準備期水貂卵巢進行轉錄組學分析,挖掘到MLT可能通過調節KRT19、GSTT1、ALAS1、NR4A1等基因促進卵巢的發育,研究結果為水貂的育種工作提供了理論基礎。