基于转录组数据挖掘外源褪黑激素影响水貂卵巢发育的分子机制

高娅薇,彭 弟,孙朝阳,晏子越,崔 凯,马泽芳

(青岛农业大学动物科技学院,青岛 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等基因促进卵巢的发育,研究结果为水貂的育种工作提供了理论基础。