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

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

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

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

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

1 材料與方法

1.1 供試材料

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

1.2 RNA提取

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

1.3 miRNA芯片檢測

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

1.4 轉錄組測序

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

1.5 差異表達基因的篩選

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

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

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

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

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

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

2 結果

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

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

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

2.2 家蠶5齡幼蟲精巢和卵巢轉錄組數據質量分析與基因比對

測序結果表明，精巢樣本(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軟件對樣品的有效序列數進行比對，結果表明：精巢樣本中能比對到家蠶基因組上的序列數占94.45%，其中在參考序列上有唯一比對位置的序列數占91.38%；卵巢樣本中能比對到家蠶基因組上的序列數占93.60%，其中在參考序列上有唯一比對位置的序列數占89.70%。將比對到基因組上的序列分布情況進行統計，定位區域分為外顯子、內含子和基因間隔區。精巢樣本中，88.77%的序列定位在外顯子區，3.93%的序列定位在內含子區，7.30%的序列定位在基因間隔區。卵巢樣本中，85.94%的序列定位在外顯子區，4.82%的序列定位在內含子區，9.24%的序列定位在基因間隔區。

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

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

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

圖3 家蠶5齡幼蟲精巢和卵巢轉錄組中差異表達基因的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芯片數據的qRT-PCR驗證

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

圖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個上調miRNA Eight up-regulated miRNAs;B:12個下調miRNA Twelve down-regulated miRNAs.圖中數據為平均值±標準誤。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 [細胞質多聚腺苷酸化元件結合蛋白(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個靶基因與轉錄組測序數據一致(圖5:B)。

圖5 家蠶5齡幼蟲精巢和卵巢間預測的差異表達的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 討論

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

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