基于RNA测序技术鉴定和分析绵羊毛囊发育相关的可变剪接事件

袁晓春,王翊帆,王雅艳,孙昊然,孟 科,李新海

(宁夏大学 农学院,宁夏 银川 750021)

羊毛具有物理保护、分泌皮脂和汗液、调节体温,以及为皮肤干细胞提供稳态、再生和修复等重要功能。羊的毛囊(hair follicle, HF)是一个皮肤上的“微型器官”,是促进羊毛生长发育的重要组织来源[1]。绵羊毛囊发育、羊毛生长发育和脱落主要受环境、自身因素、季节和光照周期的影响[2],且经历从生长期、退行期、休止期的周期性变化。绵羊毛囊发育和羊毛生长的周期性发育受到WNT、FGF、MAPK、BMP等信号通路的调控[3-6],并有microRNA(miRNA)、长链非编码RNA(long noncoding RNA, LncRNA)、转录调控因子和可变剪接等的参与[7-9]。其中,可变剪接作为一种重要的调控方式,在绵羊毛囊发育、羊毛生长的周期性变化中具有非常重要的作用。

前体mRNA(pre-mRNA)是DNA的最初转录产物,pre-mRNA经历可变剪接,导致mRNA或非编码RNA(non-coding RNA, ncRNA)的成熟。由于选择外显子和剪接位点的不同,从而形成多种不同的mRNA,并可能被翻译成不同的蛋白质异构体。一个pre-mRNA通过不同的剪接方式(选择不同的剪接位点)产生不同的mRNA剪接异构体的过程称为可变剪接[10]。可变剪接是调控真核生物蛋白质多样性和遗传多样性的重要分子机制,在不同组织、发育阶段和环境等条件下控制真核生物的发育进程[10-11]。Zhang等[12]在奶山羊中鉴定出真核翻译延伸因子1 delta基因(eukaryotic translation elongation factor 1 delta,EEF1D)的新亚型EEF1Da和EEF1Dc,EEF1D和EEF1Da在脂肪中高表达,在调节脂肪发育和产奶中起着潜在的关键作用。Zhou等[13]在牛胚胎和成年牛中预测并验证了NFIX基因的5种亚型转录产物,在胚胎组织中,Nfix2仅在大脑、脾、肾和心脏中表达,Nfix3未在大脑中表达;在成年牛中,Nfix2未表达,Nfix1和Nfix4在所有组织中均被检测到,Nfix3仅在肌肉和肺中表达,Nfix5仅在肝、脾、肾和肌肉中检测到。这一结果揭示了牛NFIX基因在胚胎发育为成牛中所起的关键作用。在绵羊中,也发现可变剪接具有重要的调控作用。张嘉楠[14]利用克隆、噻唑蓝[3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide, MTT]等方法对IGF-Ⅰ基因Ⅰ类Ea型和Ⅰ类Eb型剪接异构体进行研究,结果表明绵羊IGF-Ⅰ基因的Ⅰ类Ea型和Ⅰ类Eb型对绵羊皮肤和胎盘成纤维细胞增殖、活性与凋亡等均存在影响。血管内皮生长因子(vascular endothelial growth factor, VEGF)基因及其剪接异构体在绵羊毛囊发育及其周期循环中具有重要的调控作用。张立春等[15]研究发现,VEGF基因具有多种剪接体,且不同剪接体之间对毛囊生长发育的调控存在差异。

基因经过可变剪接形成不同的剪接异构体影响羊毛囊的生长发育,但前人的研究集中于绒山羊毛囊发育和羊毛生长,而对绵羊毛囊发育和羊毛的整个循环周期过程,以及基因可变剪接体对绵羊毛囊发育影响的相关研究较少。为研究可变剪接在绵羊毛囊和羊毛生长周期性变化中的调控机制,本研究选择杜泊成年母羊,采集羊毛生长发育不同阶段的体侧中部皮肤进行转录组测序,对不同发育阶段的可变剪接事件进行比较分析,旨在挖掘绵羊毛囊发育和羊毛不同发育阶段相关基因可变剪接事件存在的差异,为解析其在绵羊毛囊发育和绵羊羊毛生长过程中的调控机制提供借鉴。

1 材料与方法

1.1 试验动物与样品处理

本研究所用杜泊成年母羊均来自宁夏回族自治区银川市中牧亿林畜产股份有限公司。试验羊只为饲养管理条件一致、体况良好、年龄相近(约2岁)的杜泊成年母羊5只,分别在2019年9月(生长旺期,S1)、2020年1月(休止期,S2)、2020年3月(生长初期,S3)对羊的体侧中部皮肤进行采样。采样依据:经过在试验羊场对羊只的长期观察,发现杜泊羊的毛囊发育和羊毛生长呈周期性变化,在体侧中部表现明显;且采样时羊只侧卧,可使体侧中部皮肤暴露,便于皮肤采样。采样标准:最后一根肋骨后缘与体中线交界处,取2 cm2皮肤组织,用液氮保存。采样流程:1)高温高压消毒弯形剪、手术剪刀、镊子等金属器具;2)用弯形剪根据采样标准把需采样部位羊毛剪掉;3)使用手术剪刀按采样标准剪取2 cm2羊皮组织;4)使用磷酸缓冲盐溶液(phosphate buffered saline, PBS)冲洗剪下的羊皮组织并对羊只采样皮肤处进行乙醇消毒;5)将采下的皮肤组织用PBS冲洗后快速装入冻存管,投入液氮罐保存。

1.2 总RNA的提取与转录组测序

将皮肤组织样在液氮中磨碎后,采用TRizol试剂提取3个时期杜泊成年母羊体侧中部皮肤组织的总RNA,利用NanoDrop 2000分光光度计、Agilent 2100 RNA 6000 Nano kit和琼脂糖凝胶电泳检测皮肤总RNA的纯度、浓度和完整性。经检验合格后委托基迪奥生物科技有限公司进行RNA-seq文库构建和测序。经FastQC软件质控,去除含接头(adapter)的reads、含N比例大于10%的reads、全部都是A碱基的reads和低质量的reads(质量值Q≤20的碱基数占整条read的50%以上),得到干净序列(clean reads)用于后续分析。

1.3 可变剪接事件鉴定

使用rMATS[16]软件对每个样本存在的可变剪接位点和相应表达量进行检测,rMATS软件是一款可以对转录组数据进行检测和分析差异可变剪接事件的有效计算工具。rMATS软件的统计模型是根据P值和错误发生率(false discovery rate,FDR)来鉴定差异的可变剪接。本试验以P<0.05和FDR<0.05作为标准鉴定差异可变剪接事件。可变剪接事件主要包括外显子跳跃(skipping exon, SE)、内含子滞留(retention intro, RI)、外显子互斥(mutually exclusive exons, MXE)、5’端可变剪接(alternative 5’splicing site, A5SS)、3’端可变剪接(alternative 3’splicing site, A3SS)[17]。

1.4 差异剪接基因的GO(gene ontology)注释和KEGG(Kyoto Encyclopedia of Genes and Genomes)通路分析

使用基迪奥生物云平台在线数据分析平台Omicshare工具(www.omicshare.com/tools)对差异剪接基因进行GO和KEGG富集分析。采用P<0.05作为显著富集标准。

1.5 实时荧光定量(qRT-PCR)检测

随机挑选6个差异表达AS相关基因(MT1C、KRT23、RPL13、TPRXL、S100A14、RPL6),利用qRT-PCR对AS事件进行验证。在NCBI数据库下载相关转录本外显子序列,使用Primer Premier 6.0软件设计引物,用NCBI的Primer BLAST进行引物特异性验证。以GAPDH作为内参基因,由生工生物工程(上海)股份有限公司合成所需引物(引物序列见表1)。利用TB Green® Premix ExTaqTM荧光定量试剂盒检测基因的表达量,qRT-PCR反应体系(20 μL):TB Green Premix ExTaqⅡ(Tli RNaseH Plus) 10 μL,正向引物(10 μmol·L-1)0.8 μL,反向引物(10 μmol·L-1)0.8 μL,ROX参比染料(ROX reference dye,50×)0.4 μL,cDNA 2 μL,RNase-free水(RNase-free water)6 μL。qRT-PCR反应程序:95 ℃ 30 s;95 ℃ 10 s,57.8 ℃ 30 s,65 ℃ 5 s,40个循环。

表1 实时荧光定量引物列表

1.6 统计分析

使用Excel 2019软件整理数据,采用2-ΔΔCT法计算基因相对表达量,使用SPSS软件(v20.0)对数据进行统计学分析,以P<0.05作为差异显著性的判断标准,用GraphPad软件绘图。

2 结果与分析

2.1 转录组测序数据质量分析

利用Illumina HiSeqTM 4000测序平台对杜泊成年母羊体侧中部皮肤样品总RNA进行转录组测序,每个体侧中部皮肤组织样本平均产生89 200 988个raw reads,去除含adapter的reads、含N比例大于10%的reads、全部都是A碱基的reads和低质量的reads后,平均每个组织样品筛选到88 932 430个clean reads(表2)。将得到的clean reads与绵羊的参考基因组(Oar_rambouillet_v1.0)进行比对,比对率均在93.65%以上,比对到基因组唯一位置的序列比例为77.34%～87.06%。高质量数据在测序数据中所占比例均高于99%,过滤后样品的GC含量均在48.51%～52.49%,且Q30值均大于92.4%。说明测序数据真实可靠,可用于进一步分析。

表2 测序数据质量评估及与参考基因组的对比

2.2 可变剪接事件的鉴定

使用rMATS对S1 vs S2、S1 vs S3和S2 vs S3进行AS事件差异分析,共鉴定出128 033个可变剪接事件,3个比较组共有3 973个差异表达的AS事件,3个比较组分别鉴定出1 952、862和1 559个差异可变剪接事件(表3)。3个比较组中,鉴定得到SE类型最多,共96 906个,占75.69%;其次是MXE类型,共14 467个,占11.30%;A5SS、A3SS和RI类型相对较少,分别有4 837、9 220和2 603个,所占比例分别为3.78%、7.20%和2.03%。对AS事件进行显著性分析发现,3个比较组中差异表达SE类型有2 588个,其中,1 077个上调,1 511个下调;差异表达MXE类型有350个,其中182个上调,168个下调;差异表达A5SS类型有299个,其中153个上调,146个下调;差异表达A3SS类型有524个,其中248个上调,176个下调;差异表达RI事件有212个,其中77个上调,135个下调。从上述数据中可以看出,不同发育阶段的绵羊毛囊中可变剪接事件的数量存在差异,且在S1 vs S2和S2 vs S3中数量较多,这与绵羊毛囊不同时期的生长发育相一致。

表3 绵羊毛囊不同时期AS 差异分析统计

2.3 毛囊发育的差异剪接基因鉴定

以FDR<0.05为标准对rMATS软件分析结果进行筛选,共筛选到3 448个显著的差异剪接基因(图1),S1 vs S2在SE、A5SS、A3SS、MXE、RI事件中鉴定到的差异剪接基因数分别为1 010、136、236、158、106;S1 vs S3在SE、A5SS、A3SS、MXE、RI事件中鉴定到的差异剪接基因数分别为518、74、91、58、35;S2 vs S3在SE、A5SS、A3SS、MXE、RI事件中鉴定到的差异剪接基因数分别为641、77、157、87、64。差异剪接基因中SE事件的占比最高(约62.90%),RI的占比最低(5.95%)。

S1,2019年9月(生长旺期)采样;S2,2020年1月(休止期)采样;S3,2020年3月(生长初期)采样。SE,外显子跳跃;RI,内含子滞留;MXE,外显子互斥;A5SS,5’端可变剪接;A3SS,3’端可变剪接。下同。S1, Sampling in September 2019 (growth peak period); S2, Sampling in January 2020 (telogen); S3, Sampling in March 2020 (early stage of growth). SE, Skipping exon; RI, Retention intro; MXE, Mutually exclusive exons; A5SS, Alternative 5’splicing site; A3SS, Alternative 3’splicing site. The same as below.图1 杜泊羊毛囊不同发育时期差异可变剪接基因数量Fig.1 Number of differentially spliced genes in Dorper wool follicles at different developmental stages

2.4 差异剪接基因的功能富集分析

为进一步了解杜泊羊中可变剪接基因的功能,以P<0.05为标准,对上述得到的差异可变剪接基因进行GO富集分析,生长旺期与休止期绵羊毛囊差异剪接基因富集结果见图2-A,生长旺期与生长初期绵羊毛囊差异剪接基因富集结果见图2-B,休止期与生长初期绵羊毛囊差异剪接基因富集结果见图2-C。3个比较组的差异剪接基因主要富集在细胞组分(cellular component)、生物学过程(biological process)和分子功能(molecular function)等,在细胞组分层面上主要富集在细胞(GO:0005623)、细胞部分(GO:0044464)、细胞器(GO:0043226)、微管组织中心(GO:0005815)、线粒体外膜(GO:0005741)、微管细胞骨架(GO:0015630);在分子功能层面上主要富集在结合(binding,GO:0005488)、细胞骨架蛋白结合(GO:0008092);在生物学过程层面上主要富集在细胞器组织(GO:0006996)、细胞成分组织或生物发生(GO:0071840)、RNA加工(GO:0006396)、mRNA剪接体(GO:0000398)。上述结果表明,各个发育阶段的可变剪接差异基因主要富集在细胞内物质的合成、分子代谢和基因的表达调控等过程。

对3个比较组的差异剪接基因进行KEGG通路富集分析发现,差异剪接基因主要显著富集在内吞作用、MAPK信号通路、剪接体、VEGF信号通路、Hippo信号通路、ECM-受体相互作用、mTOR信号通路、鞘脂代谢、PI3K-Akt信号通路、Rap1信号通路等信号通路上(图3、表4)。内吞作用、MAPK信号通路、剪接体、神经营养因子信号通路是3个比较组共同显著富集到的信号通路。S1 vs S2的差异剪接基因主要富集在VECF、Hippo、ECM-受体相互作用、mTOR等信号通路;S1 vs S3的差异剪接基因主要富集在VECF、Rap1等信号通路;S2 vs S3的差异剪接基因主要富集在PI3K-Akt、Rap1、NF-kappaB等信号通路。通过GO、KEGG功能富集结果,以及相关文献报道,筛选出毛囊的启动和发育(SMAD2)、毛囊的生长发育(FGFR2)、毛囊的发育和循环(BMP4)、细胞凋亡抑制剂(BCL2)等与毛囊发育和羊毛生长呈周期性变化的相关基因(表4)。

A, S1vs S2; B, S1vs S3; C, S2vs S3.图3 杜泊羊不同时期毛囊发育差异剪接基因KEGG富集结果前30条通路Fig.3 Top 30 pathways of splicing gene KEGG enrichment results in different stages of hair follicle development in Dorper sheep

表4 与毛囊发育相关的差异剪接基因KEGG通路富集

2.5 qRT-PCR结果

为验证测序结果的可靠性,随机选择6个基因进行qRT-PCR验证。结果(图4)显示,差异剪接基因表达量在各组中的表达趋势与转录组测序集中的结果基本一致,说明测序数据可靠。

图4 差异可变剪接基因的qRT-PCR和RNA-seq结果比对Fig.4 Comparison of qRT-PCR and RNA-seq results of differentially alternatively spliced genes

3 结论与讨论

近年来,通过转录组学对不同组织环境中剪接异构体的深入研究,证明了可变剪接是真核生物发育过程中普遍存在的基因表达重要调控机制[10]。在mRNA前体中单个基因通过不同的剪接方式产生不同的mRNA剪接异构体,增加蛋白质编码能力并以不同的组合组装,最终产生功能和结构多样性甚至功能相反的蛋白质[18]。例如,VEGF-A165b亚型是血管生成的抑制剂,而VEGF165可以诱导毛囊干细胞分化为血管内皮细胞[19-20]。本研究共鉴定出128 033个绵羊毛囊发育相关可变剪接事件,其中外显子跳跃的比例最高,约占75.69%。章会斌等[21]和张玲[22]分别对白藜芦醇致猪卵巢颗粒细胞、云南高峰牛肝组织和脾组织的可变剪接事件进行统计,均显示出SE类型的可变剪接事件所占数量最多,本试验结果与此一致。在S1 vs S2、S1 vs S3、S2 vs S3中分别鉴定到1 952、862、1 559个差异剪接基因,说明剪接基因的表达量伴随绵羊毛囊周期性发育的特点发生显著变化。曾东等[23]在研究黑毛白皮小香猪皮肤和毛囊中TGF异构体表达时发现,TGF-β1蛋白颗粒在生长后期和退行期的毛囊外根鞘和内根鞘细胞均存在阳性表达,TGF-β2蛋白颗粒在毛囊的退行期、毛囊外根鞘细胞、毛隆突细胞等存在阳性表达。

MAPK、PI3K-Akt、Hippo等信号通路与毛囊生长发育及其循环过程有密切的关系[24-25]。Li等[26]在研究绒山羊毛囊生长周期阶段特异性时,发现Hippo信号通路促进毛囊生长,Rap1、PI3K-Akt、NF-kappaB等信号通路则在退行期和休止期起作用。PI3K-Akt信号通路、ECM-受体相互作用在休止期至生长期阶段发挥重要作用,可作为休止期-生长期再生的候选生物标志物。以上研究主要是基于转录组学mRNA的差异表达分析得出的结论,并未进行可变剪接事件的深入分析。本研究中的差异剪接基因主要显著富集在与毛囊生长发育相关的VEGF、MAPK、PI3K-Akt、Rap1、Hippo、ECM-受体相互作用、NF-kappaB、神经营养因子等信号通路上,与毛囊发育相关的MAPK和神经营养因子信号通路在3个比较组中均被显著富集,表明MAPK和神经营养因子信号通路中的这些基因可能以不同的可变剪接方式参与整个毛囊的循环过程。Hippo、mTOR信号通路在S1 vs S2被显著富集,表明涉及基因的可变剪接可能参与调控使毛囊由生长期进入休止期。在S2 vs S3,Rap1、PI3K-Akt、NF-kappaB信号通路被显著富集,表明这些信号通路涉及基因可能以不同的可变剪接方式在毛囊休止期至毛囊进入下一个周期循环中起激活传导的作用,使毛囊重新进入生长初期,刺激绵羊毛囊中新毛干的生成,从而导致旧羊毛在毛囊中及时脱落。3个比较组中各通路富集的变化情况表明,可变剪接在绵羊毛囊生长发育和羊毛的生长周期中起着重要作用。Liu等[27]利用microRNA和mRNA关联分析对绒山羊毛囊周期的分子调控机制进行研究,结果表明,SMAD2基因参与毛囊的启动和发育;Ren等[28]和Kasuya等[29]分别通过生化指标、行为测定等建立模型的方法和全贴片测定法对老鼠毛发生长进行研究,结果均显示FGF2促进毛发的生长发育;Bai等[30]研究了辽宁绒山羊皮肤组织中BMP4基因在不同时期的转录模式和甲基化状态,发现BMP4参与毛囊的发育及其周期性循环;Geueke等[31]研究表明,小鼠表皮中BCL2的条件性敲除导致小鼠毛囊凸起干细胞中的细胞凋亡升高,表明异位BCL2表达会阻断毛囊消退期间的细胞凋亡,导致毛囊静止干细胞的积累并延迟小鼠毛囊生长。同样,上述基因在本研究中也被筛选为差异可变剪接基因,并且通过功能富集分析发现,上述基因富集到与绵羊的毛囊发育相关的内吞作用、Hippo信号通路、NF-kappaB信号通路上。进一步从富集通路中筛选出SMAD2、FGFR2、BMP4、BCL2等可能与绵羊毛囊发育、羊毛生长相关的差异剪接基因,预示差异可变剪接基因在绵羊毛囊发育和羊毛生长中同样具有重要意义,为后续开展绵羊毛囊发育和羊毛呈周期性生长过程中相关基因的筛选提供理论基础。

本研究发现,在绵羊毛囊发育和羊毛生长循环过程中,不同发育阶段间存在较多的差异AS事件,它们富集在与毛囊生长发育相关的VEGF、MAPK、PI3K-Akt、Rap1等通路上,可能在绵羊脱毛过程中起到重要的调控作用。通过富集功能分析和前人研究初步筛选出SMAD2、FGFR2、BMP4、BCL2等可能参与绵羊脱毛性状的差异可变剪接基因。