不同绵羊品种发情期卵巢组织lncRNA和mRNA的差异表达分析

杨帆,姚晓磊,李康,李晓丹,唐轮,尤佩华,蔡庆贤,王锋*

(1.南京农业大学动物科技学院,江苏 南京 210095;2.广西大化瑶族自治县畜牧管理站,广西 大化 530800;3.江苏波杜农牧股份有限公司,江苏 南京 211800)

绵羊作为我国重要的经济动物,能够提供肉、奶、纤维等优质的农畜产品［1］。虽然我国绵羊品种资源丰富,但绝大部分表现为产羔率低、季节性发情,规模化舍饲的快速发展对种羊繁殖性能提出了更高的要求。因此,选育优良品种,提高绵羊繁殖力,一直是绵羊育种的关键。湖羊和小尾寒羊是我国优良的多胎绵羊品种,湖羊更是以性早熟、全年发情以及高繁殖力而闻名［2-3］。蒙古羊和陶赛特羊作为地方性品种,具有适应性强、产肉性能高等特点,但其繁殖力较低［4］。因此,开展绵羊不同品种间繁殖力差异的研究显得尤为重要。

绵羊的多胎性或高繁殖力现象早就引起科学界的高度重视,并已经成功鉴定出多个绵羊繁殖力相关基因BMPR1B［5］(骨形态发生蛋白受体1B)、GDF9［6］(生长分化因子9)、BMP15［7］(骨形态发生蛋白15)和B4GALNT2［8］(β-1,4-N-乙酰氨基半乳糖转移酶2)等,但数量还较为有限。湖羊和小尾寒羊是我国少有携带多胎主效基因FecB(fecundity booroola)的绵羊品种。但在实际生产中,携带纯合FecB基因的湖羊也存在产单羔现象,这种差异表明FecB与多胎性之间可能存在其他调控机制。

长链非编码RNA(long non-coding RNA,lncRNA)作为一种新的表观遗传调控分子,在家畜繁殖中直接或间接参与基因调控过程。为揭示lncRNA对哺乳动物繁殖力的调控机制,研究者围绕不同繁殖力母羊的下丘脑-垂体-卵巢-子宫轴开展了相关研究［9-13］。郑健等［14］对单多羔湖羊下丘脑进行转录组分析,发现下丘脑中神经信号传递通路和γ-氨基丁酸信号通路可能间接调控卵泡发育及排卵过程。从高、低繁绵羊垂体的lncRNA-mRNA调控网络中筛选到SMAD2(SMAD家庭成员2)、NMB(神经介素B)和EFNB3(肝配蛋白B3)等与繁殖相关的基因［15］。此外,在不同繁殖力湖羊卵巢的研究中,卵泡发育相关lncRNA FDNCR,可吸附颗粒细胞中的miR-543-3p,并阻止miR-543-3p与核心蛋白聚糖(DCN)基因的3′UTR区域结合,导致湖羊DCN基因转录激活和TGF-β通路抑制［16］。同时,研究者以发情期不同繁殖力的寒泊羊［17］、小尾寒羊与陶赛特［18-19］、安徽白山羊［20］卵巢为对象,通过转录组测序(RNA-seq)后均获得一定数量与繁殖相关的差异lncRNA。也有研究者发现lncRNA还参与绵羊卵巢［21］、睾丸［22］发育及初情期启动［11］等过程。以上研究表明,lncRNA在调控肉羊繁殖性能中发挥着关键作用。然而,上述研究大多数侧重于筛选同一品种不同繁殖力绵羊之间的差异lncRNA,而关于不同绵羊品种之间的研究相对较少。

本研究选取地方特色绵羊品种湖羊(HBB)和蒙古羊(Mon)为研究对象,在发情期采集卵巢组织进行RNA-seq,再结合NCBI数据库中已有的小尾寒羊(HanBB)和陶赛特羊(Dorset)发情期卵巢数据,比较分析不同绵羊品种发情期卵巢中差异表达的lncRNA和mRNA,构建mRNA-mRNA和lncRNA-mRNA调控网络,筛选关键基因并探究其表达规律,为绵羊繁殖调控及育种分子标记提供理论依据。

1 材料与方法

1.1 试验动物

本试验于江苏省泰州市姜堰某羊场进行,根据生产记录选择2～3岁体况相近、生殖生理状态相似的空怀湖羊［(49.99±1.40)kg］4只。同时,选择相同条件的空怀蒙古羊［(45.30±4.83)kg］4只。试验羊只每天08:00、18:00饲喂全混合(TMR)日粮,自由饮水。

1.2 试验材料

绵羊促卵泡激素(FSH)和促黄体生成素(LH)ELISA试剂盒(DRE-S0622c,DRE-S5741c)均购自上海卡迈舒生物科技有限公司;Trizol试剂盒购自Invitrogen;反转录试剂盒(R323-01)和荧光定量试剂盒(Q711-02)均购自南京诺唯赞生物科技有限公司;即用型免疫组化试剂盒(abs957)购自上海爱必信生物科技有限公司;基质金属蛋白酶2(MMP2)抗体(WL03224)购自沈阳万类生物有限公司。

1.3 试验设计

蒙古羊和湖羊各4 只,于11月份采用单栏方式饲养,于第1天18:00放拴处理,11 d后撤拴并注射氯前列烯醇,每只0.2 mg。次日用公羊试情,接受公羊爬跨视为发情。到下一次自然发情时,颈静脉密集连续采血,每12 min 1次,连续采集3 h,分离血清后置于-20 ℃保存备用。采血后随机挑选蒙古羊和湖羊各3 只进行屠宰,分别采集下丘脑、垂体、卵巢、子宫、心脏、肝脏、脾脏、肺脏、肾脏等组织样,单侧发育较好的卵巢放入液氮用于后续RNA-seq,另一侧卵巢以及垂体和子宫组织用布氏固定液固定,用于后续免疫组化检测。其余组织样经液氮速冻后于-80 ℃冰箱保存。

1.4 试验方法

1.4.1 生殖激素的测定对上述置于-20 ℃的血清进行FSH和LH含量测定,试验步骤严格按照FSH和LH ELISA检测试剂盒说明书进行。

1.4.2 总RNA提取及建库利用Trizol法对蒙古羊和湖羊全身组织提取总RNA。卵巢RNA的质量测定、文库构建和测序交由北京百迈客生物科技有限公司采用Illumina HiSeqTM2500 系统进行。同时,小尾寒羊(GSM2883548、GSM2883549、GSM2883550)和陶赛特羊(GSM2883554、GSM2883555、GSM2883556)的测序原始数据下载于NCBI(https://www.ncbi.nlm.nih.gov/)的GEO数据库。

1.4.3 转录组数据分析对测序的原始数据(raw reads)进行过滤,使用HISAT2软件将质控后序列(clean reads)与绵羊基因组(oar_v4.0)进行序列比对,获得clean reads在绵羊基因组上的位置信息(mapped reads)以及湖羊卵巢组织特有序列信息;使用StringTie软件对转录本和基因表达量进行分析;同时利用StringTie软件将mapped reads进行转录本拼接,以此对lncRNA进行鉴定。鉴定lncRNA包括筛选和潜在编码能力分析,选择长度 ≥ 200 bp,exon数 ≥ 2,FPKM(fragments per kilobase of transcript per million fragments mapped)≥ 0.1的转录本和利用CPC、CNCI、CPAT和PFAM等4种编码潜能的分析方法评估lncRNA;使用顺式作用(lncRNA邻近100 kb)和反式作用(LncTar软件)2种方法预测lncRNA靶基因;使用DESeq2方法对基因和转录本的表达量进行差异分析,获得差异lncRNA和mRNA。差异分组:DorsetvsHBB、HanBBvsHBB、MonvsDorset、MonvsHBB、MonvsHanBB,前者为对照组,限定条件:差异倍数(fold change,FC)≥ 2,错误发现率(false discovery rate,FDR)<0.01。使用GO seq R语言包中Wallenius对差异lncRNA的靶基因和mRNA在GO功能数据库中进行富集分析;使用KOBAS软件进行KEGG通路显著富集分析;使用STRING数据库(http://string-db.org/)预测并绘制mRNA-mRNA的蛋白互作(protein-protein interactions,PPI)网络,利用Cytoscape软件绘制lncRNA-mRNA互作图。

1.4.4 实时荧光定量PCR(RT-qPCR)为验证测序结果的准确性,从测序结果中随机挑选5个lncRNA和mRNA进行验证。根据GenBank公布的绵羊基因序列以及测序拼接后获得的lncRNA序列,采用Primer 5.0软件设计引物,以GAPDH为内参基因,引物由南京擎科生物有限公司合成(表1)。RT-qPCR按照诺唯赞公司试剂盒操作说明进行,反应体系为20 μL,采用2-ΔΔCT法计算基因的相对表达量。

表1 引物序列信息Table 1 Primer sequences information used for this study

1.4.5 免疫组化检测湖羊垂体、卵巢和子宫组织采用布氏固定液固定48 h后进行石蜡包埋和组织切片,切片厚度7 μm,37 ℃烤片24 h,免疫组化检测步骤参考本实验室的前期研究［23］。一抗为MMP2兔多克隆抗体,以1∶200进行稀释;二抗HRP酶标抗小鼠和兔二聚物为试剂盒即用试剂;对照组用兔血清代替一抗,其余步骤保持一致。

1.5 数据分析

所有数据使用Excel 2016软件初步处理后,使用Graphpad Prism 9软件进行差异显著性检验和单因素方差分析并绘图。数据结果以平均值±标准误表示。

2 结果与分析

2.1 发情期湖羊和蒙古羊生殖激素水平比较

对发情期湖羊和蒙古羊血清中FSH和LH激素含量进行检测,结果如图1所示,湖羊血清中FSH和LH含量均显著高于蒙古羊(P<0.05)。

图1 湖羊和蒙古羊血清中促卵泡激素(FSH)和促黄体生成素(LH)含量差异Fig.1 Differences of follicle stimulating hormone(FSH)and luteinizing hormone(LH)levels in serum between Hu sheep and Mongolian sheep*P<0.05. The same below.

2.2 不同绵羊品种卵巢中lncRNA与mRNA的比较分析

RNA-seq后对raw reads进行过滤,共得到150.89 Gb clean reads。如表2所示,数据与绵羊(oar_v 4.0)参考基因组比对,匹配效率68.60%～91.59%。Q30值均在87.05%以上,总碱基中G-C含量为50.30%,表明测序结果良好,clean reads质量合格,可进行下一步数据分析。本研究共鉴定到87 874个lncRNA(图2-A),其中基因间区lncRNA 44 457 个,占比50.60%,反义lncRNA 4 778 个,占比5.40%,内含子lncRNA 36 623个,占比41.70%,正义lncRNA 2 016 个,占比2.30%(图2-B);lncRNA在不同染色体中均有不同程度的分布,其中在1～3号常染色体上最多(图2-C),这与在绵羊下丘脑［23］上的研究结果一致;lncRNA表达量整体低于mRNA(图2-D);lncRNA的ORF长度在0～100 bp所占比例最高,而 mRNA的ORF长度主要分布在0～200 bp(图2-E);lncRNA和mRNA的序列长度主要分布在200～500 bp(图2-F);lncRNA主要包含2个外显子,其数量明显少于mRNA(图2-G)。

图2 不同绵羊品种卵巢组织中长链非编码RNA(lncRNA)和信使RNA(mRNA)的鉴定Fig.2 Identification of long non-coding RNA(lncRNA)and messenger RNA(mRNA) in the ovaries tissue of different sheep breeds A. lncRNA数统计的韦恩图;B. lncRNA的分类及数量;C. lncRNA的染色体分布［最外面的环代表染色体,中间从外向内分别是正义lncRNA(绿色)、基因间lncRNA(红色)、内含子lncRNA(蓝色)和反义lncRNA(灰色)］;D. lncRNA和mRNA相对表达量;E. lncRNA和mRNA ORF长度;F. lncRNA和mRNA长度;G. lncRNA和mRNA外显子数量。A. Venn diagram of lncRNA number;B. Number and classification of lncRNA;C. The distribution of lncRNA in different chromosomes［The outermost ring represents different chromosomes. From the outside toward the inside:sense-lncRNA(green),lincRNA(red),intronic-lncRNA(blue),and antisense-lncRNA(gray),respectively］;D. Relative expression level of lncRNA and mRNA;E. ORF length of lncRNA and mRNA;F. Length of lncRNA and mRNA;G. Exon number of lncRNA and mRNA.

表2 样品测序数据统计表Table 2 Sample sequencing data statistics

2.3 不同差异分组间lncRNA和mRNA差异表达分析

如图3所示,经过差异比对,在DorsetvsHBB、HanBBvsHBB、MonvsDorset、MonvsHBB、MonvsHanBB 5个差异分组中,分别获得6 920［↑(上调)2 735,↓(下调)4 185］、3 929(↑1 356,↓2 573)、5 505(↑3 966,↓1 539)、925(↑415,↓510)、3 661(↑2 712,↓949)个差异 lncRNA;12 251(↑5 126,↓7 125)、5 424(↑2 009,↓3 415)、11 548(↑6 982,↓4 566)、2 037(↑769,↓1 268)、4 791(↑3 299,↓1 492)个差异mRNA。为验证RNA-seq结果的准确性,分别随机选取5个lncRNA和mRNA进行RT-qPCR验证。结果表明,MSTRG.597282.1、MSTRG.645655.1和MSTRG.127771.15的lncRNA相对表达量在蒙古羊卵巢组织中显著高于湖羊(图4-A,P<0.05),这与测序中FPKM值的趋势一致(图4-B),而MMP2、GDF9、SFRP4(分泌型卷曲相关蛋白4)、LRP5(低密度脂受体蛋白5)和FZD7(卷曲同源物7)在内的5个mRNA表达差异也与测序相同(图4-C和D),证明测序结果可靠,可用于后续功能分析。

图4 湖羊和蒙古羊卵巢中lncRNA和mRNA的验证Fig.4 Verification of lncRNA and mRNA expression in the ovaries of Hu sheep and Mongolian sheep A. lncRNA的相对表达量Relative expression level of lncRNA;B. lncRNA的FPKM值Fragments per kilobase of exon model per million mapped fragments(FPKM)value of lncRNA;C. mRNA的相对表达量Relative expression level of mRNA;D. mRNA的FPKM值FPKM value of mRNA.**P<0.01. The same below.

2.4 不同差异分组中共有和特有差异基因分析及蛋白互作网络构建

如图5-A所示:在5个差异分组DorsetvsHBB、HanBBvsHBB、MonvsDorset、MonvsHBB、MonvsHanBB中存在129个共有差异表达基因(differentially expressed gene,DEG),去除新基因后经String网站预测,其中32个DEG存在互作关系,构建蛋白质互作(PPI)网络图(图5-B);功能富集分析进一步筛选到MMP2、CHAD(软骨黏附素)、HOXA10(子宫内膜同源盒基因A10)、GNB2(G蛋白β2)、COL4A1(Ⅳ型胶原α1)、ACTB(β-肌动蛋白)、COL6A3(Ⅵ型胶原蛋白α3)、CREB3L2(cAMP反应元件结合蛋白3样2)和LOC101115805等9个与繁殖相关的共有DEG;同时,5个组分别获得61、8、37、68、4个与繁殖相关的特有DEG(图5-A),MonvsHBB组中鉴定到GDF9、FZD7、FSTL3(卵泡抑素样蛋白3)等多个与繁殖相关的基因。

图5 共有差异基因(DEG)的蛋白互作网络图Fig.5 Protein interaction network of common differentially expressed gene(DEG)A. 各差异分组的DEG数Number of DEG in each differential group;B. 蛋白互作网络 Protein interaction network.

2.5 不同差异分组中共有和特有差异lncRNA分析及lncRNA-mRNA网络构建

进一步对各组共有和特有lncRNA进行筛选,共鉴定到60个共有差异lncRNA(图6-A),对其进行功能注释,从中选择与繁殖相关的差异lncRNA,主要富集到受精、卵母细胞成熟、配子发生等条目(图7-A)。在KEGG富集分析中,主要富集到嗅觉传导、核糖体、RNA转运等通路。繁殖相关差异lncRNA则富集到PI3K-Akt、Wnt、GnRH等信号通路(图7-B),其中,MSTRG.127771.15、MSTRG.29220.1和MSTRG.646338.1均靶向MMP2基因。进一步提高筛选条件后(|log2FC| ≥ 4,FDR<0.01),获得28个与雌性生殖相关的差异 lncRNA,构建其lncRNA-mRNA的互作网络(图6-B),其中存在多个靶向Wnt信号通路关键基因的lncRNA。

图7 与繁殖相关差异lncRNA靶基因的GO和KEGG分析Fig.7 Analysis of target genes of differentially lncRNA related to reproduction by GO and KEGG A. 繁殖相关差异lncRNA的GO分析GO analysis of lncRNA related to reproduction;B. 繁殖相关差异lncRNA的KEGG分析 KEGG analysis of lncRNA related to reproduction.

如图6-A所示:在5个差异分组中分别存在1 414、129、355、403、212个特有差异 lncRNA,经功能富集分析,获得1 139、79、313、325、133个特有差异lncRNA与繁殖相关。在5个差异分组中均存在靶向共有差异基因MMP2的特有差异 lncRNA。

2.6 基质金属蛋白酶2(MMP2)在湖羊和蒙古羊全身组织中的表达模式

如上所述,MMP2基因在各差异分组中均差异显著,因此,以MMP2为对象,检测其在蒙古羊和湖羊组织中的表达模式。结果(图8-A)表明,MMP2基因在湖羊和蒙古羊中表达趋势一致,均在卵巢组织中高表达,其次则是在肺脏中有较高表达。MMP2蛋白在湖羊腺垂体细胞(垂体)、颗粒细胞(卵巢)、子宫腔上皮细胞和腺上皮细胞(子宫)中均检测到阳性信号(图8-B)。

图8 基质金属蛋白酶2(MMP2)在湖羊和蒙古羊组织中的表达模式Fig.8 Expression pattern of matrix metallopeptidase 2(MMP2)in Hu sheep and Mongolian sheep tissue A. MMP2在蒙古羊和湖羊组织中的表达量;B. MMP2在湖羊垂体、卵巢和子宫中的表达定位(NC:阴性对照;AC:腺垂体细胞;GC:颗粒细胞;LE:腔上皮细胞;GE:腺上皮细胞)。A. Expression level of MMP2 in Mongolian sheep and Hu sheep;B. Expression location of MMP2 in pituitary,ovary and uterus of Hu sheep(NC:Negative control;AC:Adenohypophysis cell;GC:Granular cell;LE:Luminal epithelium;GE:Glandular epithelium).

3 讨论

多胎性是绵羊种畜最重要的评价指标,其不仅与遗传因素(品种)有关,还受到表观遗传调控的影响。研究表明,lncRNA参与卵巢发育［25］、卵子发生［26］、妊娠维持［27］等生殖过程,这为探究lncRNA调控绵羊繁殖性能提供契机。因此,本研究对不同绵羊品种发情期卵巢组织进行RNA-seq分析后发现,在5 个差异分组(DorsetvsHBB、HanBBvsHBB、MonvsDorset、MonvsHBB、MonvsHanBB)中,分别获得6 920、3 929、5 505、925、3 661个差异lncRNA和12 251、5 424、11 548、2 037和4 791个差异mRNA。有学者结合湖羊的体型以及蒙古羊的南迁历史,认为蒙古羊是湖羊的祖先［28］。而在MonvsHBB组中差异lncRNA和mRNA数最少,在一定程度上也证实了这个说法。因此,2个亲缘关系较近的绵羊品种,造成繁殖力差异的lncRNA和mRNA更值得我们关注。

经RT-qPCR验证后进一步对共有差异lncRNA和mRNA进行功能注释,发现与繁殖相关的共有差异lncRNA主要注释到受精、卵母细胞成熟、配子发生等条目,富集到PI3K-Akt、Wnt等已被证实与繁殖相关的通路中,并构建了lncRNA-mRNA互作网路图,这为后续深入探究绵羊繁殖力调控机制提供参考。其中,MSTRG.62667.21、MSTRG.64407.1等多个差异lncRNA均靶向Wnt通路相关基因,这与前人的研究结果一致［29］,提示lncRNA可能通过调控Wnt通路关键基因来影响绵羊繁殖力。同时,本研究筛选获得9个(HOXA10、MMP2、CHAD等)与繁殖相关的共有DEG,HOXA10基因被认为是子宫内膜容受性的关键调控因子［30］,CHAD蛋白也在小鼠胚胎和卵巢中较高表达［31］。在不同绵羊品种卵巢组织中均存在显著差异,表明这些DEG可能在绵羊卵巢的分子调控机制中发挥重要作用。此外,在MonvsHBB组中,筛选到GDF9、FSTL3、FZD7等与繁殖相关的特有DEG,GDF9已被证实与绵羊繁殖力相关［6］,FSTL3基因可以通过降低绵羊孕激素水平,调节绵羊的发情周期［32］。以上研究表明,lncRNA和mRNA可能参与绵羊繁殖力的调控过程,但具体机制有待揭示。

如上所述,在各差异分组中,Wnt信号通路关键基因MMP2及相关lncRNA的表达水平存在显著差异。MMP2蛋白能够降解Ⅳ型胶原蛋白促使卵泡形态发生改变,对哺乳动物卵泡发育发挥重要作用［33］。本研究结果表明,MMP2蛋白在湖羊垂体、卵巢和子宫组织中均检测到阳性信号,在蒙古羊和湖羊的卵巢组织中呈现高表达状态,提示MMP2基因可能参与绵羊繁殖调控过程。有学者认为MMP2的高表达可能是卵泡闭锁的标志,参与优势卵泡的选择［34］,这在一定程度上解释了本试验中低繁蒙古羊卵巢中MMP2表达量高于高繁湖羊的原因,但具体机制仍需进一步验证。

另外,FSH和LH是主导卵泡发育的2种主要激素。初情期后,FSH和LH对卵泡发育起主要调节作用,且二者相互协调,刺激卵泡成熟与排卵。由于多胎性与卵泡发育关系密切,且与排卵优势卵泡数量直接相关,因此,FSH与LH含量与绵羊的产羔性状存在一定的相关性。本课题组前期研究发现,高繁湖羊中LH水平显著高于低繁湖羊［35］。与此类似,不同绵羊品种中,发情期单羔美利奴羊中FSH和LH含量均显著低于多羔湖羊［36］,这在一定程度上验证了本研究中湖羊血清中FSH和LH含量显著高于蒙古羊的结论。但也有研究者认为这是高繁绵羊卵巢卵泡对FSH和LH的灵敏度更高造成的结果［37］。

综上,本试验绘制了不同绵羊品种卵巢组织中差异lncRNA和mRNA表达图谱,筛选到与繁殖相关的lncRNA和mRNA,探究候选基因MMP2在不同组织中的表达差异,并初步确定MMP2蛋白的表达定位,为后续揭示lncRNA参与绵羊繁殖力的调控机制奠定基础。