家蚕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)，可以作为进一步功能验证的重点候选基因。