家蠶5齡幼蟲精巢和卵巢microRNA芯片及轉(zhuǎn)錄組比較分析

張 冰,李 娜,闞云超,2,*

(1.南陽師范學(xué)院,河南省伏牛山昆蟲生物學(xué)重點實驗室,昆蟲生物反應(yīng)器河南省工程實驗室,河南南陽 473061;2.河南大學(xué)生命科學(xué)學(xué)院,河南開封 475004)

微小RNA (microRNAs,miRNAs)是一類高度保守的長約為22 nt的非編碼RNA，通過其種子序列以完全或不完全匹配的方式結(jié)合靶基因的3′或5′端非翻譯區(qū)(untranslated region,UTR)，在轉(zhuǎn)錄或轉(zhuǎn)錄后水平上實現(xiàn)對靶基因表達的調(diào)控(Saliminejadetal.,2019)，參與調(diào)節(jié)信號傳導(dǎo)以及細胞增殖、凋亡和分化等生命活動(Bartel,2004)。Lee等(1993)首次在秀麗隱桿線蟲Caenorhabditiselegans中發(fā)現(xiàn)小分子RNAlin-4可能通過一個反義的RNA-RNA互作來抑制lin-14基因的翻譯。2000年，研究人員在線蟲中揭示了攜帶lin-41基因3′UTR的報告基因以依賴于小分子RNAlet-7的方式進行調(diào)節(jié)表達(Reinhartetal.,2000)。近幾年的研究表明，miRNA在昆蟲中廣泛存在，在昆蟲變態(tài)發(fā)育、卵子發(fā)生與胚胎發(fā)生過程中發(fā)揮重要的生物學(xué)功能(Lucasetal.,2015)。miR-100/let-7/miR-125簇與半變態(tài)和全變態(tài)發(fā)育的物種翅形態(tài)建成密切相關(guān)(Belles,2017)。Dicer 1 (Dcr1)和Argonaute 1 (Ago1)兩個關(guān)鍵酶通過調(diào)節(jié)miRNA的產(chǎn)生加工過程調(diào)控其功能(F?rstemannetal.,2007;Czech and Hannon,2011)。miR-184,miR-14以及miR-2/miR-13a/miR71簇等依賴于Ago1產(chǎn)生的miRNA是東亞飛蝗Locustamigratoria中保幼激素介導(dǎo)的卵黃發(fā)生、卵母細胞成熟和卵巢發(fā)育所必需的(Songetal.,2013)。敲減miRNAlet-7導(dǎo)致雌性果蠅Drosophila生殖力和產(chǎn)卵能力出現(xiàn)中度缺陷(Sokoletal.,2008)。Dcr1缺失突變體通過miR-280調(diào)節(jié)卵母細胞的成熟過程(Nakaharaetal.,2005)。

家蠶Bombyxmori是鱗翅目的模式物種，miRNA參與其發(fā)育、新陳代謝以及疾病發(fā)生和生殖過程的調(diào)節(jié)(Wangetal.,2014)。miR-let7通過與FTZ-F1和Eip74EF (E74)靶基因互作調(diào)控家蠶蛻皮以及變態(tài)發(fā)育(Lingetal.,2014)。miR-14通過維持蛻皮激素的動態(tài)平衡以促進家蠶的蛻變(Liuetal.,2018)。miR-1a-3p通過下調(diào)bmVMP23蛋白增加卵子的結(jié)構(gòu)完整性(Chenetal.,2013)。在家蠶質(zhì)型多角體病毒(Bombyxmoricytoplasmic polyhedrosis virus,BmCPV)感染的細胞中microRNA Novel-31*與溶血素基因5′UTR結(jié)合，上調(diào)溶血素基因的表達(施莉莉等,2016)。bmo-miR-2739和miR-167協(xié)調(diào)調(diào)控家蠶卵子發(fā)生過程中卵黃蛋白原受體的表達(Chenetal.,2020)。然而家蠶精巢組織發(fā)育中是否有miRNA參與卻尚未見報道。本研究通過收集家蠶5齡幼蟲的精巢和卵巢樣本進行miRNA芯片和轉(zhuǎn)錄組分析，比較精巢和卵巢組織中差異表達的miRNAs及預(yù)測的相關(guān)差異表達靶基因。為深入研究miRNA調(diào)控家蠶性腺發(fā)育提供基礎(chǔ)。

1 材料與方法

1.1 供試材料

供試家蠶為大造p50品系，幼蟲用新鮮桑葉飼養(yǎng)，溫度為25±3℃，相對濕度為60%±10%，光周期為12L∶12D。

1.2 RNA提取

對收集的1.1節(jié)家蠶5齡幼蟲進行解剖獲得精巢和卵巢樣品(分別命名為Test和Control)，每種樣品3組，每組來自30頭個體，分別用液氮研磨，取約50 mg粉末加入1 mL Trizol(Invitrogen，美國)試劑，按照說明書進行RNA抽提，利用NanoDrop ND-2000對提取RNA的濃度和純度進行檢測，瓊脂糖凝膠電泳檢測RNA的完整性。

1.3 miRNA芯片檢測

為了找到家蠶精巢和卵巢組織差異表達miRNAs，我們對1.2節(jié)提取得到的RNA進行加尾孵育后使用FlashTag標(biāo)記試劑盒(Genisphere，美國)進行生物素標(biāo)記miRNA，與Affymetrix miRNA 4.0芯片進行雜交(GeneChip Hybridization Oven 645)，然后進行清洗染色和掃描(GeneChip Scanner 3000 7G)，最后使用AGCC軟件(Affymetrix GeneChip Operating Software)將芯片的熒光掃描圖像保存成DAT文件。Agilent Feature Extraction軟件處理DAT原始數(shù)據(jù)，GeneSpring軟件將原始數(shù)據(jù)進行歸一化分析并以Excel格式輸出數(shù)據(jù)。利用SAM(significance analysis of microarray)R程序包分析差異表達miRNA，篩選標(biāo)準(zhǔn)是：P<0.05且變化倍數(shù)log2(fold change,FC)≥2；GraphPad Prism 6軟件用于圖形展示。將差異表達miRNA 做無監(jiān)督層次聚類分析(hierarchical clustering)，以熱圖(heatmap)形式顯示。此節(jié)實驗由北京中康博生物有限公司輔助完成。

1.4 轉(zhuǎn)錄組測序

對1.2節(jié)提取得到的RNA進行樣品檢測，利用去核糖體試劑盒去除rRNA 并將RNA片段化，逆轉(zhuǎn)錄合成單鏈cDNA，再合成雙鏈cDNA，之后進行末端修復(fù)并連接測序接頭。最后進行PCR富集和文庫質(zhì)檢，使用Illumina HiSeq-pe150上機測序。測序得到的原始圖像文件經(jīng)堿基識別(base calling)分析轉(zhuǎn)化為測序序列(sequenced reads)，即原始序列數(shù)據(jù)(raw reads)。對原始序列數(shù)據(jù)進行過濾、去接頭序列以及去除rRNA(28S rRNA,18S rRNA,5.8S rRNA和5S rRNA等)序列后，得到高質(zhì)量的序列數(shù)據(jù)(clean data)。使用Hisat2軟件將測序數(shù)據(jù)與參考基因組進行比對，對測序數(shù)據(jù)覆蓋區(qū)域及覆蓋深度做出綜合評價。之后將比對到基因組上的reads分布情況進行統(tǒng)計，定位區(qū)域分為外顯子內(nèi)含子和基因間隔區(qū)域。

1.5 差異表達基因的篩選

采用RPKM計算和歸一化基因表達量。利用DEGseq軟件篩選TestvsControl的差異表達基因(differentially expressed genes,DEGs)，標(biāo)準(zhǔn)為q≤0.05且|log2(fold change)|≥1。q值為使用DEGseq對P值進行校正所得，可以有效去除總體假陽性。q值越低，基因表達差異越顯著。當(dāng)q值小于默認閾值0.05時，表示基因表達差異非常顯著。

1.6 差異表達miRNA及其潛在靶基因的qRT-PCR驗證

從TestvsControl比較組中隨機選取20個差異表達miRNA，包含8個上調(diào)miRNA和12個下調(diào)miRNA，以家蠶U6基因(GenBank登錄號：AY649381)為內(nèi)參基因。miRNA qRT-PCR驗證的上游引物為成熟miRNA序列(表1)，下游引物為反轉(zhuǎn)錄試劑盒中的通用引物。1.2節(jié)中提取所得的部分RNA采用All-in-oneTMmiRNA qRT-PCR Detection Kit(GeneCopoeia)反轉(zhuǎn)錄試劑盒取2 μg合成cDNA第1鏈，于-20℃保存?zhèn)溆?。為了找到差異表達miRNA的靶基因，我們用miRanda，pita和RNAhybrid 3種預(yù)測軟件在轉(zhuǎn)錄組數(shù)據(jù)中分析潛在的靶基因，預(yù)測出來的靶基因進行KEGG通路分析，依據(jù)KEGG通路校準(zhǔn)后的P<0.05篩選存在差異的靶基因，結(jié)合轉(zhuǎn)錄組數(shù)據(jù)篩選與miRNA表達趨勢相反或表達趨勢一致的靶基因。使用GeneTool軟件設(shè)計不同基因的qRT-PCR引物，選擇家蠶肌動蛋白ActinA3基因(actinA3)為內(nèi)參基因(表1)。qRT-PCR采用KAPA SYBR?FAST qPCR Kit Master Mix (2×)Universal熒光定量PCR試劑盒進行。PCR反應(yīng)條件：95℃預(yù)變性5 min；95℃變性20 s，58℃退火20 s，72℃延伸30 s，共40個循環(huán)。技術(shù)重復(fù)3次。qRT-PCR結(jié)果采用2-ΔΔCt法(Livak and Schmittgen,2001)計算相對表達量。采用SPSS19.0軟件中的ANOVA法進行單因素方差分析。所得結(jié)果均以平均值±標(biāo)準(zhǔn)誤表示，并利用Excel進行圖形分析。

表1 qRT-PCR引物信息Table 1 Primer information for qRT-PCR

1.7 DEGs及差異表達miRNA的靶基因的功能注釋

對1.5節(jié)篩選到的DEGs和1.6節(jié)篩選到的差異表達miRNA的靶基因進行KEGG信號通路分析，探索靶基因參與的生物學(xué)過程和涉及的信號通路。Cytoscape V3.2.1軟件用來解讀差異性基因相關(guān)的信號通路及功能，顯著性的檢驗水準(zhǔn)為P<0.05。通過Fisher Exact Test計算P值，以P<0.05為顯著性閾值分別得到相對于背景具有統(tǒng)計意義的信號轉(zhuǎn)導(dǎo)及疾病通路，從而得到基因集合在KEGG類別上的分布信息和顯著性情況。

2 結(jié)果

2.1 家蠶5齡幼蟲精巢和卵巢中miRNA表達譜

基因芯片檢測結(jié)果顯示有68個差異表達miRNAs，與卵巢組織相比，在精巢組織中高表達的miRNA有36個；相反，在精巢組織中低表達的miRNA有32個。為了全面直觀地展示卵巢和精巢基因芯片數(shù)據(jù)之間的關(guān)系及差異情況，將差異表達miRNA做無監(jiān)督層次聚類分析，以熱圖形式顯示(圖1)，在精巢組織中高表達的miR-2817,miR-2796-3p和miR-2774c聚在一枝。在精巢組織中低表達的miR-998,miR-263a-5p和miR-3340聚在一起。

圖1 家蠶5齡幼蟲精巢vs卵巢轉(zhuǎn)錄組中差異表達miRNAFig.1 Differentially expressed miRNAs in the testes vs ovary transcriptomes of the 5th instar larvae of Bombyx moriTest:精巢Testis;Con:卵巢Ovary.紅色代表上調(diào)，綠色代表下調(diào)。Red represents up-regulated and green represents down-regulated.

2.2 家蠶5齡幼蟲精巢和卵巢轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量分析與基因比對

測序結(jié)果表明，精巢樣本(Test)共獲得43 551 550條序列，其中GC含量為43.27%，Q20為98.78%，Q30為96.37%。卵巢樣本(Control)共獲得39 399 194條序列，其中GC含量為45.07%，Q20為98.73%，Q30為95.58%。

使用Hisat2軟件對樣品的有效序列數(shù)進行比對，結(jié)果表明：精巢樣本中能比對到家蠶基因組上的序列數(shù)占94.45%，其中在參考序列上有唯一比對位置的序列數(shù)占91.38%；卵巢樣本中能比對到家蠶基因組上的序列數(shù)占93.60%，其中在參考序列上有唯一比對位置的序列數(shù)占89.70%。將比對到基因組上的序列分布情況進行統(tǒng)計，定位區(qū)域分為外顯子、內(nèi)含子和基因間隔區(qū)。精巢樣本中，88.77%的序列定位在外顯子區(qū)，3.93%的序列定位在內(nèi)含子區(qū)，7.30%的序列定位在基因間隔區(qū)。卵巢樣本中，85.94%的序列定位在外顯子區(qū)，4.82%的序列定位在內(nèi)含子區(qū)，9.24%的序列定位在基因間隔區(qū)。

2.3 家蠶5齡幼蟲精巢和卵巢轉(zhuǎn)錄組的DEGs及KEGG通路富集

家蠶5齡幼蟲精巢和卵巢轉(zhuǎn)錄組得到3 991個DEGs，其中上調(diào)基因和下調(diào)基因分別為2 033和1 958個(圖2)。KEGG通路富集分析結(jié)果顯示，精巢和卵巢樣本中有972個DEGs富集在21條通路，涉及代謝(metabolism)(15)、遺傳信息處理(genetic information processing)(4)、細胞進程(cellular processes)(1)和有機系統(tǒng)(organismal systems)(1)5個大類，其中富集基因數(shù)最多的是新陳代謝(metabolism)(339)、核糖體(ribosome)(87)、氧化磷酸化(oxidative phosphorylation)(76)、碳代謝(carbon metabolism)(59)和剪接體(spliceosome)(56)(圖3)。

圖2 家蠶5齡幼蟲精巢vs卵巢轉(zhuǎn)錄組中差異表達基因Fig.2 Differentially expressed genes (DEGs)in the testes vs ovary transcriptomes of the 5th instar larvae of Bombyx mori紅色代表上調(diào)，綠色代表下調(diào)。Red represents up-regulated and green represents down-regulated.

圖3 家蠶5齡幼蟲精巢和卵巢轉(zhuǎn)錄組中差異表達基因的KEGG通路富集Fig.3 KEGG pathway enrichment of differentially expressed genes (DEGs)of the testis and ovary transcriptomes of the 5th instar larvae of Bombyx mori

2.4 家蠶5齡幼蟲精巢和卵巢差異表達miRNA芯片數(shù)據(jù)的qRT-PCR驗證

為了驗證miRNA芯片的結(jié)果，在精巢和卵巢比較組中隨機選取20個差異表達miRNAs，包含8個上調(diào)miRNA和12個下調(diào)miRNA(表1)，進行qRT-PCR驗證，結(jié)果如圖所示，8個上調(diào)miRNA在精巢組織中表達均高于卵巢組織中的表達(圖4:A)，12個下調(diào)miRNA在精巢組織中表達均低于卵巢組織中的表達(圖4:B)。這些實驗驗證均與miRNA芯片測定的結(jié)果一致。

圖4 qRT-PCR驗證20個差異表達miRNA在家蠶5齡幼蟲精巢和卵巢組織中的表達Fig.4 Expression of 20 differentially expressed miRNAs in the testis and ovary of the 5th instar larvae of Bombyx mori verified by qRT-PCRA:8個上調(diào)miRNA Eight up-regulated miRNAs;B:12個下調(diào)miRNA Twelve down-regulated miRNAs.圖中數(shù)據(jù)為平均值±標(biāo)準(zhǔn)誤。Data in the figure are mean±SE.下同The same below.

2.5 家蠶5齡幼蟲精巢和卵巢差異表達miRNA的靶基因和qRT-PCR驗證

最后找到了4組表達趨勢相反和1組表達趨勢一致的miRNA與靶基因(圖5:A)。分別是bmo-m iR-2774a與LOC101745556[磷脂酰肌醇蛋白聚糖(glypican-4)基因]，bmo-miR-92b與LOC101735954 [細胞質(zhì)多聚腺苷酸化元件結(jié)合蛋白(cytoplasmic polyadenylation element-binding protein 2)基因]；bmo-miR-3266與LOC733130[線粒體磷酸烯醇丙酮酸羧激酶(mitochondrial phosphoenolpyruvate carboxykinase)基因]和LOC778467[果糖1,6-二磷酸醛縮酶(fructose 1,6-bisphosphate aldolase)基因]以及bmo-miR-3321和LOC101744895(功能未知)。其中LOC101745556參與Wnt信號通路，LOC101735954參與背腹軸形成通路，LOC733130和LOC778467參與糖酵解/糖異生通路，LOC101744895參與脂代謝。通過qRT-PCR驗證找到的這5個靶基因與轉(zhuǎn)錄組測序數(shù)據(jù)一致(圖5:B)。

圖5 家蠶5齡幼蟲精巢和卵巢間預(yù)測的差異表達的miRNA的靶基因及qRT-PCR驗證Fig.5 Predicted target genes of differentially expressed miRNAs between the testis and ovary of the 5th instar larvae of Bombyx mori and qRT-PCR validationA:4組表達趨勢相反和1組表達趨勢一致的差異表達microRNA與潛在靶基因Four groups of differentially expressed miRNAs and their potential target genes showing the opposite expression trend and one group of differentially expressed miRNA and its target gene showing the consistent expression trend;B:5個潛在靶基因的qRT-PCR驗證Verification of the five potential target genes by qRT-PCR.

3 討論

精巢中的精母細胞通過前期的同源染色體配對、聯(lián)會、交叉、重組實現(xiàn)遺傳物質(zhì)的交換，產(chǎn)生帶有不同遺傳信息的精子細胞；卵巢中的卵母細胞停留在前期Ⅰ，直至產(chǎn)卵或受精之后再繼續(xù)減數(shù)第2次分裂。細胞學(xué)的觀察表明家蠶精巢和卵巢組織發(fā)生過程中存在著顯著差異。最近的研究發(fā)現(xiàn)雄蠶最早在2齡幼蟲時即可觀察到部分精母細胞進入細線期(Yamamotoetal.,2017)，至5齡幼蟲時，可以觀察到第1次減數(shù)分裂的諸多時期，不同時期所占比例不同。

精巢和卵巢組織中除了性母細胞減數(shù)分裂的細胞學(xué)差異之外，同樣存在諸多蛋白和表觀調(diào)控修飾。線蟲卵母細胞發(fā)育過程中，至少有23個miRNAs以依賴或不依賴于Drosha蛋白的方式調(diào)控卵子發(fā)生的關(guān)鍵過程(Minogueetal.,2018)。毛立明等(2007)通過對家蠶蛹期第2天的雌雄生殖腺蛋白質(zhì)雙向電泳比較分析，找到了55個雌蠶特異性和73個雄蠶特異性蛋白斑點。我們通過對家蠶5齡幼蟲精巢和卵巢組織進行轉(zhuǎn)錄組分析，找到3 991個DEGs(圖2)，主要涉及代謝以及遺傳信息處理相關(guān)通路(圖3)。結(jié)合miRNA芯片分析，找到4組對應(yīng)的miRNA與靶基因(圖5)。bmo-miR-2774a和bmo-miR-3266尚未見相關(guān)功能報道。bmo-miR-92b可以通過負調(diào)控Mef2調(diào)節(jié)果蠅肌肉發(fā)育(Chenetal.,2012)，還可以通過下調(diào)jigr1基因調(diào)節(jié)神經(jīng)母細胞的自我更新(Yuva-Aydemiretal.,2015)，最近的研究還發(fā)現(xiàn)miR-92b直接與crebA的3′UTR結(jié)合提高沃爾巴克氏菌感染果蠅的學(xué)習(xí)記憶能力(Bietal.,2019)。本研究中bmo-miR-92b的預(yù)測靶基因為細胞質(zhì)多聚腺苷酸化元件結(jié)合蛋白2基因，其在精巢中的表達遠遠低于卵巢(圖5:A)，可以作為進一步功能驗證的重點候選基因。