家养绵羊Ovis aries属于鲸偶蹄目Cetartioda ctyla反刍亚目Ruminantia牛科Boridae绵羊属Ovis，作为主要家畜在人类的经济生活中占据着重要地位，经过人工选育出现了不同的品种(张立岭，2004)。绵羊的繁殖力极大地影响着绵羊养殖的经济性，对绵羊繁殖力进行研究是提升中国绵羊生产水平的必由之路(吕晓曼等，2012)。每胎产羔数是影响绵羊繁殖力的重要因素，提高绵羊的每胎产羔数可在短期内较大地增加绵羊养殖的经济效益，因此，绵羊的产羔性状研究一直受到人们的重视(杨广礼，罗玉柱，2004)。发现和利用绵羊多羔性状相关遗传资源，阐明绵羊多羔性状形成的分子调控机理，筛选和鉴定绵羊多羔性状相关基因，逐渐成为绵羊多羔研究中的热点问题，有利于人们从遗传角度理解和利用绵羊的多羔性状(付绍印等，2017)。畜牧业是青藏高原地区重要的经济支柱型产业，其中绵羊畜牧业占据着举足轻重的地位。彭波半细毛羊是西藏自治区选育的家养绵羊新品种，也是西藏第一个国家级家畜新品种，目前在西藏地区大量养殖，具有极为重要的经济价值(央金，2012)。自然条件下彭波半细毛羊每胎以产单羔为主，偶有产双羔现象。目前虽已发现MMPR-1B、GDF9、BMP15等基因的突变与特定品种绵羊的产羔数相关，但彭波半细毛羊之前的研究未发现此类基因与产羔数的相关性(Yan et al., 2001；Chong et al., 2019；穆方方等，2019；Zheng et al., 2019)。并且，研究也证实绵羊的产多羔性状是一个低遗传性状，每胎产羔数受品种、饲养、年龄、遗传等多因素影响(Tang et al., 2019)。因此，彭波半细毛羊产双羔的机理需要从更全面的角度进行研究。

绵羊产多羔现象与其卵巢卵泡发育及排卵数等因素直接相关(付绍印等，2017)，对产双羔彭波半细毛羊的卵巢特异表达基因进行研究，有利于更好地理解产双羔彭波半细毛羊的分子特异调控机制，也有可能找到与其产双羔性状相关的候选基因并加以利用，从而促进彭波半细毛羊的人工快速繁育工作。到目前为止，有关彭波半细毛羊产双羔相关基因的研究还极为缺乏。本研究基于转录组测序(RNA-Seq)对彭波半细毛羊发情期的卵巢特异表达基因进行了研究，重点关注产双羔彭波半细毛羊卵巢的基因特异表达模式，分析和筛选了可能对彭波半细毛羊卵泡发育和排卵数产生影响的特异表达基因。相关研究结果不仅有利于人们加深对彭波半细毛羊这一高原藏绵羊品种产双羔机理的理解，也有利于产双羔候选基因的筛选利用研究工作。

1 材料与方法 1.1 材料

彭波半细毛羊养殖于西藏自治区农牧科学院畜牧兽医研究所位于西藏自治区拉萨市林周县的良种养殖基地，选取年龄3~4岁且有2次连续产双羔记录的母羊3只，为双羔组；选取年龄3~4岁且有2次连续产单羔记录的母羊3只，为单羔组。各组绵羊体质量相近，健康状况良好。

1.2 方法 1.2.1 绵羊卵巢总RNA的提取制备及测序文库构建

双羔组及单羔组绵羊于养殖基地现场屠宰，快速分离出卵巢，经磷酸缓冲盐溶液(PBS)清洗后置于干冰中冷冻保存及运输。实验室超净工作台中将卵巢于液氮中研磨粉碎，之后使用Trizol进行卵巢组织总RNA的提取。总RNA质量检测合格后，每个样品取3 μg总RNA使用Illumina的NEBNextⓇ Ultra^TM RNA Library Prep Kit建库试剂盒进行文库构建。通过Oligo(dT)磁珠富集mRNA，将mRNA随机打断片段化后以随机寡核苷酸为引物合成cDNA，并使用RNaseH降解RNA链。纯化后的双链cDNA经过末端修复、加A尾并连接测序接头后进行PCR扩增，PCR扩增产物经纯化后得到文库。文库构建完成后，使用Agilent 2100 bioanalyzer对文库的insert size进行质量检测(杨永刚等，2019)。

1.2.2 文库测序和差异基因表达分析

不同文库分别采用Illumina Hiseq平台进行测序，测得的序列数据(reads)经过滤后得到clean reads。使用NCBI数据库中Ovis aries基因组数据(https://www.ncbi.nlm.nih.gov/genome/?term=Ovis+aries)作为参考基因组，利用HISAT2(2.0.5)将配对末端clean reads与参考基因组进行比对。采用subread中的featureCounts工具对基因进行表达水平的定量，根据基因长度计算每个基因的FPKM值。双羔组和单羔组之间的基因差异表达分析使用DESeq2R(Love et al., 2014)进行统计分析，以|log₂(FoldChange)|>1和多重假设检验校正后的P值(padj) < 0.05为差异表达基因的筛选标准。

1.2.3 转录组数据的荧光定量PCR(qPCR)验证

以双羔组及单羔组各样品总RNA为模板逆转录成cDNA，设计10个基因的定量PCR引物，使用荧光定量PCR仪进行检测，检测结果使用2^-△△Ct分析方法进行差异表达量计算。

1.2.4 差异表达基因的功能分类及代谢通路富集分析

差异基因GO功能富集和KEGG通路富集通过clusterProfiler R进行，分别基于GO数据库(http://www.geneontology.org)和KEGG数据库(http://www.genome.jp/kegg/pathway.html)进行富集分析，上调和下调差异表达基因的GO和KEGG分析分别进行。

2 结果与分析 2.1 差异基因筛选与验证

各组样品总RNA提取后，检测浓度和完整性，合格后进行文库构建和转录组测序分析。对2组样品所建文库进行测序，每个样品的测序量约15 G，测序所得clean reads进行基因组比对后，比对率均达到90%以上。

所测数据通过与参考基因组数据进行比对后，共发现21 509个基因。以|log₂(FoldChange)|>1和padj < 0.05为筛选阈值，与单羔组相比，双羔组共有2 489个基因出现差异表达，其中上调基因1 021个，下调基因1 468个。以-log₁₀(padj)为纵坐标，log₂(FoldChange)为横坐标，绘制双羔组相对于单羔组的卵巢差异表达基因分布情况(图 1)，当padj=0.05时，-log₁₀(padj)= 1.301。

为了验证转录组数据的可靠性，选取了10个差异表达基因进行qPCR验证(表 1)，以绵羊beta-actin基因作为内参(引物序列5'-3'分别为：GAGCATCCCCAAAGTTCTAC、AATCACCTCCCCTGTGT-G)。与单羔组相比，双羔组基因NFIC、MAST1、GLB1L3、MYCBPAP、RBMX2的表达下调，基因CA5A、SLC7A5、CMTM6、SCD、TMED5的表达上调，且基因的差异表达趋势与转录组数据基本一致(图 2)。上述数据表明转录组测序所得到的数据具有较高的可靠性。

2.2 差异表达基因的GO功能富集分析

为了对双羔组的卵巢差异表达基因进行功能分析，本研究首先使用GO数据库对差异基因的功能进行了富集分析。GO分析分别从分子功能、细胞组分、生物过程3个层面对差异表达基因进行了功能富集(图 3)。如以padj < 0.05为显著性富集的阈值，可发现部分上调基因发生了显著性的富集。其中，发生显著性富集的上调基因在生物过程层面主要富集于跨膜转运，在细胞组分层面主要富集于质膜蛋白复合体、质膜、质膜部分，在分子功能层面主要富集于跨膜受体蛋白激酶活性(表 2)。显著富集的上调基因均与细胞膜及膜相关功能有关，表明双羔组的卵巢细胞中有较活跃的跨膜运输或者信号转导等膜结构相关的生理活动，而这种生理活动有可能与双羔组卵巢的卵泡发育相关。

以padj < 0.05为阈值，双羔组下调基因在生物过程层面未发生显著性的GO富集，但仍可发现下调的基因主要富集于翻译、物质代谢以及细胞呼吸等生物过程。在细胞组分层面上，双羔组下调基因显著富集于细胞内核糖核蛋白复合体和核糖核蛋白复合体。这预示着双羔组卵巢细胞中核糖核蛋白复合体有所变化，同时细胞内的物质代谢及蛋白翻译等活动有所减缓。

2.3 差异表达基因的KEGG富集分析

KEGG富集分析有助于更深入了解双羔组差异表达基因对卵巢调控的影响，差异表达基因经KEGG富集(图 4)：上调基因多富集于细胞粘附以及膜相关信号转导通路，下调基因多富集于产能产热及物质代谢相关通路，这与GO功能分析的结果相互印证。

以padj < 0.05为阈值，-log₁₀(padj)＞1.301时表示发生了KEGG显著性富集，-log₁₀(padj)值越大，富集越显著。双羔组卵巢上调基因出现了较显著的富集(图 4；表 3)，下调基因未见显著性富集。

双羔组的上调基因显著富集于与卵巢内细胞粘附及运动、发育相关的通路。其中细胞粘附分子、粘附斑、细胞骨架调节、细胞外基质受体相互作用均与细胞粘附、细胞运动直接相关，而细胞粘附与运动是贯穿于卵泡发育及排卵整个过程的生命活动。

另一方面，双羔组卵巢内大量上调基因显著富集于TGF-beta、cGMP-PKG、Hippo、cAMP等信号通路中，这些信号通路与细胞的生长、分化、发育、迁移以及激素调节等功能有较紧密的关系(Sandoval et al., 2017；Gasior et al., 2019)，如TGF-beta信号通路可影响动物卵巢内卵泡的发育调节过程(Zheng et al., 2008)。

结合GO和KEGG的富集分析，发现双羔组彭波半细毛羊发情期内卵巢中有大量基因表达出现显著性的上调，上调基因的功能主要与细胞膜相关的生理活动相关，包括细胞粘附、运动以及质膜相关的信号转导通路。

3 讨论

绵羊产双羔性状受多方因素影响，但与绵羊卵巢中卵泡的发育及排卵数直接相关。因此，产双羔或多羔绵羊卵巢的调控机制研究对于理解其特殊的产双羔或多羔性状具有重要的意义。RNA-Seq作为一种高通量差异基因表达研究方法，可以为相关研究提供更快速和全面的研究手段(兰道亮等，2014)。对于绵羊等动物繁殖力的研究，转录组测序能够更有效和全面地帮助人们发掘繁殖性状相关的候选基因以及内在的分子机制(Miao et al., 2014)。彭波半细毛羊作为藏绵羊的代表新品种，其产双羔性状可能具有特殊的机理，但目前研究资料极为缺乏，严重阻碍青藏高原地区特有绵羊品种的快速繁育研究。

绵羊卵巢中，卵泡的发育实际上是卵母细胞和周围的颗粒细胞及卵泡膜细胞互相调控的过程，卵泡的形态和功能也随着发育时间而不断变化(Matikainen et al., 2001；肖国宏，2014)。卵泡的发育过程由原始卵泡的启动开始，原始卵泡启动机制可能与卵泡内外分泌的各种因子、神经递质或卵母细胞-颗粒细胞间的信号传递等有关(Kidder & Mhawi，2002)。碱性成纤维细胞生长因子、干细胞因子、白血病抑制因子、神经生长因子、骨形态发生蛋白、胰岛素和角化细胞生长因子等，均能诱导原始卵泡的启动性募集，从而发育为成熟卵泡(Nilsson & Skinner，2004；应诗家，2012)。卵泡中的卵母细胞周围包裹着各种体细胞，体细胞合成的相关激素物质会通过卵泡液运送到卵母细胞中，以调控和促进卵母细胞的发育成熟和启动减数分裂(Matzuk et al., 2002；Wigglesworth et al., 2015)。可见，卵泡的发育及卵子的成熟过程均需要卵泡组织及卵母细胞膜相关信号转导通路和物质跨膜运输功能的配合，而这些过程均与卵巢内各细胞的细胞膜结构相关，这预示着卵巢内的细胞膜结构在卵泡发育及排卵过程中担负着举足轻重的作用。

转录组测序结果表明，双羔组卵巢组织中的基因在转录层面出现明显的差异。GO分析发现双羔组的卵巢组织细胞中有较活跃的跨膜运输等膜结构相关的生理活动，KEGG富集分析则发现双羔组卵巢中的细胞粘附、细胞运动以及跨膜信号转导等相关生理活动显著增强，而细胞粘附、运动、信号转导等通路的实现则离不开细胞的跨膜物质运输以及细胞膜结构。考虑到绵羊等哺乳动物卵泡发育以及排卵的特点，双羔组卵巢组织中细胞膜相关生理活动的活跃，预示着双羔组卵巢内有更为活跃的卵泡细胞分化、发育及排卵活动，这可能是双羔组卵巢内独特的调控机制。

付绍印等(2017)报道了巴美肉羊下丘脑、垂体、卵巢组织的转录组研究数据，发现粘附斑以及细胞外基质受体相互作用等通路的基因在上述组织中均出现了差异表达；而细胞粘附因子通路基因仅在卵巢组织中出现差异表达。彭波半细毛羊卵巢中也出现了细胞粘附因子通路基因的上调表达及富集，因此细胞粘附因子可能是产双羔彭波半细毛羊潜在的筛选标志物。另外，本研究发现了双羔组卵巢中TGF-beta信号通路基因的上调和富集，而绵羊繁育力相关基因MMPR-1B、GDF9、BMP15均为TGF-beta超家族成员，在多羔小尾寒羊卵巢中也发现特异表达TGF-beta信号通路蛋白(贾建磊，2015)。这表明虽然彭波半细毛羊未发现产羔性状与上述基因的突变有关，但TGF-beta信号通路的活跃与产羔性状可能有关。此外，研究结果中发现大量溶质运载蛋白家族(solute carrier family，SLC)基因出现明显上调，SLC7A5基因的qPCR结果也验证了相关的转录组结果。溶质运载蛋白一般位于细胞膜等质膜上，承担着各种生物分子转运的作用。溶质运载蛋白基因在双羔组的卵巢中大量上调表达，表明双羔组卵巢细胞中的跨膜物质运输机能显著高于单羔组，从而促使卵巢内发生一系列最终有利于卵泡分化、发育及排卵的生理过程。这预示着卵巢细胞中溶质运载蛋白家族基因的上调表达亦可作为彭波半细毛羊产双羔的一个信号，与此信号相关的基因或者激素等小分子可以作为彭波半细毛羊产双羔性状的潜在标志物进行深入研究。

4 总结

本研究基于转录组测序研究了产双羔彭波半细毛羊卵巢的特异表达基因，发现上调基因1 021个，下调基因1 468个。对差异基因进行功能分析后发现，双羔组显著上调的基因多与细胞膜相关的细胞活动有关，如细胞的粘附、运动、跨膜信号转导及物质转运。双羔组彭波半细毛羊卵巢内细胞的物质运输和信号转导频率加快，从而促进卵泡的分化和发育；同时，双羔组卵巢内细胞的粘附及运动功能增强，更有利于卵母细胞的成熟和成熟卵子的排出，从而诱导绵羊产双羔性状的出现。上述结果预示着卵巢内细胞膜相关生命活动的增强是彭波半细毛羊产双羔性状的一个重要因素。此外，细胞粘附因子、TGF-beta信号通路以及溶质运载蛋白家族基因的上调表达可作为彭波半细毛羊产双羔性状的潜在信号进行更深入的研究。本研究结果不仅有助于对彭波半细毛羊产双羔性状内在机制的了解，更有助于藏绵羊快速繁育相关问题的深入研究。