綿羊睪丸差異基因及蛋白質互作網絡關系研究

劉在霞，段 仕，孫燕勇，呂 琦，付紹印，何小龍，張文廣，劉永斌

(1.內蒙古農業大學動物科學學院，內蒙古自治區動物遺傳育種與繁殖重點實驗室，呼和浩特 010018；2.內蒙古自治區農業基因組大數據工程研究中心，呼和浩特 010018；3.內蒙古自治區農牧業科學院，呼和浩特 010031；4.和田師范專科學校，和田 848000)

綿羊在全球畜牧業中占有很大比例，而繁殖效率是畜牧業收入水平的關鍵指標，大多數綿羊品種是單胎，季節性繁殖，繁殖率低[1]，迫切需要采用兩年三產的繁殖模式。在提高綿羊繁殖效率的實踐中，對母羊的研究較多，最常見的方法是利用生殖輔助技術來改善發情時間[2]、增加母羊的排卵次數[3]和進行人工授精[4]。但綿羊的成功繁殖取決于綿羊兩性，獲能、受精和胚胎發生過程中的失敗可能是公羊精子的原因[5]。研究發現，其他物種精子中有超過14 000個轉錄本在受精時被轉移到卵子中，這些轉錄影響早期胚胎發育和受孕是否成功，Hodge等[6]收集了3個綿羊品種的精液，在不同品種和射精質量對比中共鑒定出754個差異表達基因(DEGs)，這些基因在維持精子發生、受精、受孕、胚胎發育和后代生產性能中發揮著重要作用。睪丸是產生精子和分泌雄激素的重要器官，睪丸發育和精子發生是一個高度復雜和精確的過程，在發育過程中睪丸的內部結構和功能發生了一系列的變化[7-8]。隨著睪丸大小的改變，產生精子的效率也會改變。營養不足影響性成熟雄性睪丸質量，降低生精效率和精子運動速度，并增加性成熟雄性綿羊的精子DNA損傷[9]。此外，睪丸發育和精子發生主要受mRNA調控，這是動態并具有階段性的[10]。蛋白質是一切生物的物質基礎，蛋白質-蛋白質互作(protein-protein interaction，PPI)在任何生長、繁殖和新陳代謝中都起著極其重要的作用，即使在單個細胞也是如此[11]，PPI通過連接相關的蛋白質，在每個單細胞中啟動各種功能或結構蛋白的作用[12]。Jansen等[13]使用貝葉斯網絡方法預測全基因組酵母中的PPI，用串聯親和純化試驗(TAP)驗證了預測，提供了酵母PPI的全面視角。為全面探索綿羊睪丸在不同品種、不同營養水平和不同發育階段差異基因及PPI網絡，本研究使用課題組前期利用機器學習算法預測的綿羊全基因組820 072對PPI關系[14]，構建睪丸的PPI網絡，識別可能調控綿羊睪丸發育和精子發生的關鍵基因，以期進一步提高綿羊的繁殖效率，為今后的綿羊繁殖研究提供理論依據。

1 材料與方法

1.1 綿羊睪丸轉錄組數據獲取

綿羊睪丸轉錄組數據來源于NCBI的高通量測序數據存儲庫SRA(sequence read archive)，獲取碼：PRJEB19199、PRJNA282344和PRJNA421437，共獲得74個轉錄組數據，包括50只特克塞爾羊與蘇格蘭黑頭羊雜交后代(TB)、16只澳洲美利奴綿羊(8只在65 d內增重10%的高飼糧綿羊(M+)和8只在65 d內減重10%的低飼糧綿羊(M-))、9只不同發育階段的小尾寒羊(ST，2、6和12月齡各3只)。

1.2 睪丸中基因表達量的計算

利用SRAtoolkit軟件的fast-dump參數將高通量測序數據轉變為fastq文件，利用FastQC軟件對原始測序數據進行質量評估，使用Trimmomatic軟件過濾接頭、低質量和未知堿基數目過多的reads，獲得高質量的clean data；利用Hisat2軟件將clean data與綿羊參考基因組(Oar_v4.0)比對，利用Stringtie軟件對轉錄本進行組裝，之后采用R軟件的Ballgown包計算基因在每個樣本中的轉錄本表達量；使用FPKM(Fragment Per Kilobase of per Million mapped reads)值來衡量基因的表達水平。

1.3 PPI網絡構建及網絡中關鍵基因的選擇與分析

PPI數據集來自于課題組前期采用隨機森林(random forest,RF)和十折交叉驗證法構建的預測模型，對28 592個綿羊蛋白質最終預測到820 072對蛋白質具有互作關系[14]。將820 072對互作關系中的睪丸表達基因和睪丸DEGs使用Cytoscape軟件(version 3.7.2)構建睪丸PPI網絡。 使用Cytoscape軟件中的MCODE(molecular complex detection)插件[15]對PPI網絡中的節點進行密度聚類，計算網絡中最重要的節點和集合(cluster)，PPI網絡中選擇集合標準如下：Degree Cutoff≥2;Node Score Cutoff=0.2;K-core≥2;Max Depth=100。

1.4 睪丸DEGs分析

綿羊睪丸DEGs數據集來自3類：①不同綿羊品種(8只M+和8只TB綿羊共16個RNA-Seq)，使用R語言limma軟件包對2個品種(M+vsTB)進行DEGs分析，在差異基因檢測過程中，將log2|FoldChange|≥1且P<0.05作為篩選標準；②通過文獻挖掘16只不同飼喂營養水平的澳洲美利奴綿羊(M+vsM-)[16]；③9只不同發育階段(2、6、12月齡)ST，分別代表睪丸發育過程中性成熟前、性成熟、性成熟后3個階段[17]。

1.5 構建加權綿羊睪丸差異基因共表達網絡

構建加權基因共表達網絡(weight gene co-expression network analysis,WGCNA)[18]，為減少因測序批次不同帶來的誤差，選取同一批次的樣本進行WGCNA分析，共表達網絡的2 058個DEGs表達量來自8個M+樣本。采用分層聚類和動態樹切割算法對模塊進行識別，最小模塊數設置為30個基因，并將合并高度相似的模塊的最小高度設置為0.25。

1.6 功能富集分析

為了研究模塊及其基因潛在功能，使用在線生物信息學數據庫Metascape (https:∥metascape.org/)對所選模塊中的基因進行富集分析，并使其服從GO功能和KEGG通路富集分析，P<0.01作為篩選標準(Min Overlap=3，Min Enrichment=1.5)，利用Cytoscape插件AutoAnnotate自動識別基因集簇,并對基因集進行注釋。

2 結 果

2.1 綿羊睪丸轉錄組數據分析

2.1.1 睪丸表達基因分析 通過對74個樣本進行轉錄組分析，共鑒定到18 650個睪丸表達基因。為降低基因表達的假陽性，篩選3%的樣本中FPKM≠0，且至少有一個樣本FPKM≥3，最終獲得11 884個基因用于后續分析，其中425個基因位于X染色體，11 150個基因位于常染色體，309個為“unplaced”。

2.1.2 睪丸PPI網絡構建 對睪丸中最終獲得的11 884個基因在820 072對互作關系中提取PPI，僅有8 291個基因獲得12 621個蛋白質(節點，Node)，產生237 366對PPI(邊，Edge)，占綿羊總PPI的28.9%，并構建237 366對的PPI網絡(圖1A)。對PPI網絡進行MCODE分析，共確定366個集合，其中集合1的分值最大為142.69，且關系對最多，是由200個Nodes和18 535個Edges構成；集合7中Nodes最多為411個，對排名前100的集合繪制柱狀圖(圖1B)。

2.2 綿羊睪丸DEGs分析

2.2.1 綿羊睪丸DEGs的PPI網絡構建 不同綿羊品種(M+vsTB)共獲得735個DEGs，其中上調基因324個，下調基因411個；不同營養飼喂水平(M+vsM-)共獲得2 243個DEGs；不同發育階段共獲得13個DEGs。綜上，去除重復的基因后共獲得2 058個綿羊睪丸DEGs。對2 058個DEGs在237 366對PPI中提取互作關系，構建差異基因PPI網絡，發現僅有1 644個DGEs對應2 020個蛋白質，產生46 169對PPI，對PPI網絡進行MCODE分析，確定了87個集合(圖2)。選擇前4個得分高的集合，集合1得分最高為25.203，有13個基因；集合2得分為8.759，有26個基因；集合3得分為8.605，有45個基因；集合4得分為8.125，有12個基因(圖3)。

A,睪丸DEGs的PPI網絡；B,PPI蛋白集合邊、點統計A,PPI network of testis DEGs;B,Nodes and Edges counts of PPI圖2 睪丸DEGs的PPI網絡分析Fig.2 PPI network analysis of testis DEGs

A～D，集合1～4代表PPI網絡的核心基因A-D，Cluster 1-4 represented the core gene of PPI network圖3 PPI網絡MCODE分析Fig.3 PPI network MCODE analysis

2.2.2 構建加權綿羊睪丸DEGs共表達網絡 對2 058個DEGs進行WGCNA分析，以確定數據集中的共表達基因，參數設置為：beta=18，R2>0.8。根據基因表達量進行相關度的聚類，聚類度較高的基因被分配到一個模塊中，不同顏色代表不同的模塊，而灰色模塊內基因代表其相關性較低，無法被分配到其他模塊，不進行后續分析，本研究共獲得7個共表達模塊(圖4A)，模塊中的基因數目為51～929個(圖4B)。

A，基因平均連鎖層次聚類樹狀圖；B，模塊中的基因數量，模塊由直方圖顏色指定A,Average linkage hierarchical clustering dendogram of the genes；B,The number of genes in the module，modules were designated by color histogram圖4 加權基因共表達網絡分析Fig.4 Weighted gene coexpression network analysis

2.2.3 DEGs功能富集分析 利用Metascape數據庫探索模塊的功能和途徑，Blue模塊共富集到106個功能術語，主要參與雄性生殖系統相關的功能，包括雄性配子產生(gamete generation)、繁殖(fertilization)、精子發生(spermatogenesis)、鞭毛運動(cilium movement)、細胞周期的過渡(cell cycle transition)、AMPK信號通路(AMPK signaling pathway)等。使用Cytoscape插件AutoAnnotated可視化網絡，把相似度>0.3的功能術語用邊連接，共產生14個功能集合術語，富集網絡圖顯示了功能集合內和集合間功能富集結果的相似性，網絡中的節點代表一個功能術語，邊代表功能之間的連接關系，同一顏色節點代表一個功能集合富集的過程，發現它們高度連接并聚成一個完整的網絡(圖5)，其中配子產生(gamete generation)包括：生殖細胞發育(germ cell development)、精子發生(spermatogenesis)、精子細胞分化(spermatid differentiation)、精子發育(spermatid development)、雄性配子產生(male gamete generation)和多細胞生物體內參與生殖的細胞過程(cellular process involved in reproduction in multicellular organism)。表明Blue模塊內基因是參與精子發生和精子細胞發育過程的重要的基因集合。

圖5 睪丸差異基因富集結果全局網絡圖Fig.5 Global network diagram of testicular DEGs enrichment result

2.2.4 關鍵基因的篩選 對2 058個DEGs通過PPI網絡的MCODE和WGCNA聯合分析，考慮到精子發生和精子細胞的發育過程相關基因主要分布在Blue模塊，因此利用Blue模塊內基因和PPI網絡關鍵的4個集合，確定了25個基因(表1、圖6)，發現僅來自不同營養水平的基因有19個：基質重塑相關8(matrix remodeling associated 8,MXRA8)、富含亮氨酸重復序列和免疫球蛋白樣結構域3(leucine rich repeats and immunoglobulin like domains 3,LRIG3)、神經細胞黏附分子2(neural cell adhesion molecule 2,NCAM2)、絲氨酸/蘇氨酸激酶17b(serine/threonine kinase 17b,STK17B)、V-集和跨膜結構域(V-set and transmembrane domain containing,VSTM4)、恩貝酸(embigin,EMB)蛋白、免疫球蛋白超家族成員3(immunoglobulin superfamily member 3,IGSF3)、CDC樣激酶(CDC like kinase 4,CLK4)、VRK絲氨酸/蘇氨酸激酶2(VRK serine/threonine kinase 2,VRK2)、淀粉樣前體蛋白(amyloid beta precursor protein,APP)、磷脂酰肌醇-4-磷酸3-激酶催化亞單位2β(phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 beta,PIK3C2B)、尼沙林(nischarin,NISCH)、分揀微管連接蛋白2(sorting nexin 2,SNX2)、紅細胞膜帶4.2蛋白(erythrocyte membrane protein band 4.1 like 2,EPB41L2)、SHC適配器蛋白(SHC adaptor protein 4,SHC4)、小核糖核蛋白多肽B2(small nuclear ribonucleoprotein polypeptide B2,SNRPB2)、SRA莖環相互作用的RNA結合蛋白(SRA stem-loop interacting RNA binding protein,SLIRP)、RNA結合基序單鏈相互作用蛋白2(RNA binding motif single stranded interacting protein 2,RBMS2)和RNA結合基序蛋白46(RNA binding motif protein 46,RBM46)；僅來自不同品種的基因有4個：生長因子受體結合蛋白14(growth factor receptor bound protein 14,GRB14)、酪氨酸激酶2(tyrosine kinase 2,TYK2)、異種核糖核蛋白A1樣 (heterogeneous nuclear ribonucleoprotein A1-like,LOC101109899)和ELAV類RNA結合蛋白1(ELAV like RNA binding protein 1,ELAVL1)；來自不同營養水平和品種的基因有2個：受體酪氨酸激酶3(receptor tyrosine kinases3,TYRO3)和肝細胞黏附因子(hepatic and glial cell adhesion molecule,HEPACAM)。25個基因分布于13條染色體上，且各基因間存在互作關系，如SLIRP與LOC101109899基因、RBM46、RBMS2和ELAVL1基因存在互作關系。

表1 Blue模塊25個關鍵基因

最外圈表示綿羊全基因組27條染色體的編號與長度，格越長表示染色體越長，同時外圈展示了25個關鍵基因在染色體上的位置；圈里的熱圖代表關鍵基因的表達量，表達量越高，顏色越深；最內圈的連線表示關鍵基因間互作關系The outermost ring indicates the number and length of 27 chromosomes in the whole genome of sheep，the longer the grid is,the longer the chromosome length is,meanwhile,the outermost ring shows the location of 25 key genes on the chromosome；The heat map in the circle represents the expression of key genes，the higher the expression,the darker the color；The innermost ring lines represent interactions between key genes圖6 關鍵基因Circos圖Fig.6 The key genes Circos map

SLIRP基因表達量最高為100.3，RBM46基因表達量為42.47，ELAVL1、RBM46和SLIRP基因顯著富集到AMPK信號通路和芳香族化合物分解代謝過程；TYK2基因表達量為39.18，主要參與白細胞介素-12、-23和-27介導的信號通路；TYRO3基因表達量為21.54，主要參與配子產生、精子發生、蛋白磷酸化。

3 討 論

PPI在任何生命的生長、繁殖和新陳代謝中都起著極其重要的作用，睪丸發育和精子發生是影響雄性綿羊生育力的重要因素。雖然已經有大量關于睪丸發育和精子發生的研究，但對睪丸PPI網絡的研究較少，本研究對睪丸表達的11 884個基因和2 058個DEGs分別構建綿羊PPI網絡，發現睪丸的PPI網絡可較為緊密地形成一個大的網絡，只有少部分關系稀疏，對2 058個DEGs提取互作關系，僅有1 644個DGEs對應2 020個蛋白質，產生46 169對PPI，通過MCODE、WGCNA和Metascape分析，發現Blue模塊是與精子發生相關的基因集。 Silva等[19]構建了人類睪丸PPT網絡，共鑒定出1 778個蛋白質及它們之間的32 187對相互作用，并有少量孤立的簇；富集分析顯示，最重要的生物學過程類別都與生殖相關，如有性生殖(P=3.71E-156)、精子發生(P=5.36E-147)和受精(P=4.33E-44)，與本研究結果相近。

睪丸是轉錄最復雜的器官，其表達的基因數量是所有哺乳動物器官中最多的，據報道，睪丸在人和其他物種中廣泛的轉錄超過所有蛋白質編碼基因的80%，長鏈非編碼RNA(lncRNA)基因和基因組的基因間序列，包括假基因和轉座因子等，向睪丸的傾斜表達更為明顯，對人類和小鼠精子發生的單細胞轉錄組研究發現，廣泛的轉錄是通過轉錄掃描機制來糾正DNA損傷，從而維持雄性生殖系中的DNA序列完整性。在精子發生過程中表達的基因在轉錄鏈上表現出較低的突變率，且在群體中具有較低的多樣性[20]。Soumillon等[21]對人、恒河猴、小鼠、負鼠和雞的不同組織器官進行轉錄組測序，發現睪丸分別表達20 322、20 078、21 819、18 198、14 574個蛋白質編碼基因，且常染色體蛋白質編碼基因在睪丸中的轉錄頻率高于其他器官。Clark等[22]對3頭成年公羊和3頭成年母羊的主要組織和5類細胞進行轉錄組測序，發現有19 921個基因至少在一個組織中表達。本研究通過嚴格的表達量篩選最終獲得11 884個睪丸表達基因，由于與參考基因組比對、過濾或表達量篩選時參數和條件與前人不同而產生差異，由此可見睪丸表達基因較為豐富。

哺乳動物睪丸是先天免疫的特殊器官，分泌多種免疫調節因子，對外來和自身抗原耐受[23]，Toll樣受體(Toll-like receptors，TLR)在支持細胞和生殖細胞中表達，在被配體激活后引發睪丸的先天免疫，Tyro3基因負調控支持細胞中的TLR3信號轉導[24]，降低炎癥因子表達水平，識別與結合凋亡細胞，控制睪丸對病原體的先天免疫應答[25]。本研究通過MCODE分析發現，Tyro3基因位于分值最高的集合1，表明其與周邊基因的密集程度很高，處于關鍵地位，且其在不同營養水平和不同品種組都是差異基因，表明其在精子發生中可能起關鍵作用。通過功能富集分析發現,Tyro3基因主要參與配子產生、精子發生、蛋白磷酸化等生物學過程。唐紅梅[26]研究發現，Tyro3基因在哺乳動物的神經、免疫、生殖和血管等組織中廣泛表達，并在神經發育、免疫調節、精子發生、NK細胞分化、血小板功能及血管發生中發揮著重要的作用。

ELAVL1基因位于1號染色體上，又名HuR或Hua，與mRNA的3′-非翻譯區結合，調節轉錄的穩定性。Ghosh等[27]研究發現，小鼠出生后ELAVL1基因的缺失可導致造血器官萎縮、腸絨毛丟失、阻塞性小腸炎，并在10 d內死亡。Chi等[28]對小鼠敲除和過表達HuR基因發現，HuR基因在原始生殖細胞中的失活與正常的減數分裂后細胞形成和精子細胞成熟是矛盾的，此外，在圓形精子細胞中過表達HuR基因延緩了精子細胞的分化，敲除ELAVL1基因使生殖細胞分化缺陷、精子完全失活，導致雄性不育[29]。本研究中，ELAVL1基因處于集合4，主要參與AMPK信號通路。研究發現，在支持細胞中AMPK信號通路調節能量代謝，連接復合體的穩定性和增殖，一旦平衡被打破，睪丸的微環境和精子的質量都會受到影響，如在α1 AMPK整體敲除小鼠中，精子表現為頭部異常、鞘彎曲和活動能力受損[30]。以上均表明ELAVL1基因與精子發生有關。

RBM46基因位于7號染色體上，是一種新的生殖細胞特異性因子，可調節減數分裂和細胞凋亡，在性腺組織的生殖細胞中特異表達。研究發現，RBM46基因突變的斑馬魚中，RBM46基因的特異性破壞導致突變體的雄性化和不育，雖然突變體的精原細胞基本正常，但精子發生受損，減數分裂未發生，表明RBM46基因在斑馬魚生殖細胞第一次減數分裂過程中起著至關重要的作用[31]。Zhai等[32]研究發現，RBM46基因分布在小鼠胚胎干細胞的細胞質中，通過靶向β-Catenin mRNA降解調節胚胎干細胞分化，β-Catenin是Wnt信號通路中的一個關鍵效應子，RBM46基因轉錄后調控在胚胎干細胞分化中起重要作用，RBM46基因過表達導致胚胎干細胞向滋養外胚層分化，而敲除RBM46基因導致胚胎干細胞向內胚層分化。本研究中RBM46基因表達量為42.47，處于集合4，顯著富集到AMPK 信號通路和芳香族化合物分解代謝過程，并與SLIRP、LOC101109899、RBMS2和ELAVL1基因存在互作關系，但其調控分子機制尚不清楚。SLIRP基因是核受體(NRs)抑制因子，其存在于核內與NR靶基因相關的復合物中，在骨骼肌、心臟、肝臟和睪丸中表達量最高[33]，本研究中SLIRP基因表達量為100.3，與此結論相一致。 Colley等[34]研究發現，當敲除SLIRP基因的純合雄性小鼠與野生型雌性雜交時，平均產仔數比野生型雄性和雌性小鼠雜交產仔數減少了1/3，且敲除小鼠的精子運動能力明顯低于野生小鼠，電鏡觀察發現敲除小鼠精子中段/環連接處斷裂，線粒體形態發生改變，表明SLIRP基因的缺失導致與精子結構和線粒體形態受損相關的弱精子癥。Shan等[35]研究表明，SLIRP基因調節雄性繁殖能力且與精子運動呈正相關，與正常精子相比，弱精子癥患者精子中SLIRP基因的mRNA和蛋白質表達量降低，氧化損傷增加，能量代謝降低。

TYK2基因位于5號染色體上，是JAK激酶蛋白家族中第一個鑒定的成員，JAK家族還有3個成員：JAK1、JAK2和JAK3[36]。TYK2基因在細胞因子轉導過程中起著關鍵作用，廣泛參與免疫、細胞增殖、生長、分化、遷移和細胞凋亡等過程[37]，本研究中TYK2基因顯著富集到參與白細胞介素-12、-23和-27介導的信號通路，與前人研究結果相似。D’Cruz等[38]通過免疫印跡和激光掃描共聚焦顯微鏡檢測到JAK家族的4個成員中僅TYK2基因在人精子尾部中表達，表明其可能調節精子活力。 本試驗中TYK2基因表達量為39.18，處于集合2，與RBM46、SLIRP等基因存在互作關系，但未對JAK1、JAK2和JAK3基因在綿羊精子尾部的表達情況進行驗證，仍有待進一步研究。

4 結 論

本研究通過對綿羊睪丸表達基因PPI網絡構建，以及不同營養水平、不同品種、不同發育階段差異基因的整合與構建PPI網絡和加權共表達網絡，最終找到與睪丸精子發生有互作關系的25個基因，進一步了解了睪丸的網絡調控機制，明確了影響綿羊精子生成中的關鍵因素，為綿羊繁殖研究提供了理論依據。