刘 宇,段晓燕
(1.河北北方学院 实验动物中心,河北 张家口 075000;2河北北方学院 教务处,河北 张家口 075000)
随着测序技术的不断发展,大量哺乳动物基因组已经完成了测序并进行了组装。在哺乳动物基因组中,mRNA被认为是最重要的元件,但是其仅占基因组或转录组的一小部分,其余为非蛋白编码基因包括siRNAs、miRNAs、piRNAs以及lncRNA等等,这些非编码RNA也参与生命活动各个方面的调控过程[1-2]。
不同于小RNA分子,lncRNA是一种类似于mRNA长度大于200个碱基,在转录水平、转录后水平和染色体重塑等多种层面调控基因表达的的保守性较低的非蛋白编码RNA分子[3-4]。非编码RNA的作用包括且不仅限于X染色体失活、P53介导的细胞凋亡、癌细胞转移、诱导干细胞的重编程以及其他的一些过程[5-6]。随着第二代测序技术的发展以及lncRNA数据库的建立,为大规模筛选与性成熟调节过程相关的lncRNA提供了技术依据。
家禽的性成熟过程受到不同转录本构建的分子调控网络共同作用,而卵巢是一个独特又重要的生殖器官和内分泌腺,蛋鸡的卵巢发育程度直接关系着其产蛋性能的优劣,卵巢组织有关基因或非编码RNA表达的差异类型及其表达量的变化在分子水平上影响蛋鸡繁殖性能[7]。但相关研究,特别是关于lncRNA仍报道甚少。
本研究基于二代测序技术对不同产蛋阶段蛋鸡卵巢组织差异表达lncRNA进行筛选与鉴定,进一步挖掘影响蛋鸡繁殖性能的关键调节基因以及功能元件,研究成果有助于研究人员更好的理解蛋鸡性成熟过程中的分子调控机制,从而为蛋鸡育种工作服务。
以张家口地区饲养的京红蛋鸡为研究对象。整个试验周期内,所有蛋鸡均自由取食、饮水,执行常规光照和免疫等程序。于20周龄、30周龄分别随机选取产蛋鸡,称重后取最接近平均值的3羽产蛋鸡,屠宰并取出卵巢组织,放入液氮中保存备用(20周龄:LH_5M,30周龄:LH_6M)。获取的卵巢组织采用Trizol一步法提取总RNA,采用Bioanalyzer 2100及RNA 6000 Nano Lab Chip Kit RNA对RNA质量与纯度进行检测。采用Ribo-Zero Gold Kit对10μg总RNA进行核糖体RNA去除,离子打断法将总RNA裂解为小片段,逆转录法合成cDNA文库,以PCR法检测文库大小及浓度,采用Illumina Hiseq 4000(杭州联川)双末端测序法上机测序。
对测序原始数据去除低质量数据,得到的有效数据(Valid Data)采用Bowtie2[8]和物种的参考基因组(Gallus_gallus_5.0)进行比对,同时根据基因组注释文件(gtf和gff)指定的基因位置信息分别进行统计。对获得的转录本去除已知的mRNA和小于200 bp的转录本,再对剩下的转录本进行lncRNA预测。预测软件为CPC(Coding Potential Calculator)和CNCI(Coding-Non-Coding Index)。
采用StringTie和Ballgown对测序获得的转录本进行表达水平计算,基因的表达水平主要采用每百万测序碱基中每千个转录子测序碱基中所包含的测序片断数 FPKM度量基因表达的丰度值。差异表达lncRNA使用R语言的Ballgown包进行差异分析,选择阈值矫正P值Padj<0.05。
对筛查到的差异表达lncRNA进行的功能富集分析依据其靶基因进行。根据lncRNA位置确定lncRNA靶基因。筛选获得的靶基因用于功能富集分析,包括Gene ontology和KEGG通路分析。采用David软件进行功能富集分析,ggplot2作散点图展示。
对高通量测序结果筛查到的差异表达lncRNA,利用实时荧光定量验证表达量差异性。qRT-PCR引物由Primer 5.0软件设计(表1)。以3-磷酸甘油醛脱氢酶基因(GAPDH)作为内参。反应体系包括95 ℃ 3 min、95 ℃15 s以及60 ℃和72 ℃各30 s,共计40循环,熔解阶段由95 ℃ 30 s和60 ℃ 1 min组成,运用2-ΔΔCt法计算不同组别lncRNA的倍数变化,使用SPSS 26.0单因素方差分析进行分析,设定P<0.05为差异有统计学意义。
表1 引物信息Table 1 Primer information
各样品测序数据产出如表2所示,共得到57.95 Gb Clean data,单个样品Clean data大于8.29 Gb。GC 含量在48%~51%范围内,Q30碱基百分比大于94.71%;说明测序数据获取完整,可用于后续分析。参考基因组比对结果显示检测样品的Reads与参考基因组的比对效率介于76.90%~79.61%之间,大部分转录本在参考基因组中均是唯一比对位置。
表2 转录组测序及序列比对结果Table 2 Result of RNA-seq and sequence align
对测序获得的数据,按照筛选依据,在20周龄和30周龄京红蛋鸡卵巢组织筛选到1 283个lncRNA。对lncRNA的分类(图1A)表明,筛选到的lncRNA多为未知或基因间转录本,其次为反义链外显子重叠区域。对筛选获得的lncRNA进行基因组定位,结果见图1B。对lncRNA与mRNA的结构特征和表达模式分析发现,lncRNA长度与mRNA类似,多数分布在>1 000 bp范围(图1C)。而表达量上,lncRNA要明显低于蛋白质编码基因(图1D),log10(FPKM)分布在0~1之间,同时数量也少于蛋白质编码基因。
图1 lncRNA转录组特征A.lncRNA分类统计;B.lncRNA染色体分布统计;C.lncRNA长度统计;D.lncRNA表达量与数量统计Fig.1 Transcriptom characteristics of lncRNAA. Classified statistic of lncRNA;B. Chormsome distribution of lncRNA;C. Statistic of lncRNA length;D. Expression level and amount of lncRNA
在本研究中,采用矫正P值≤0.05作为阈值,共筛选到10个差异表达lncRNA,其中7个在20周龄蛋鸡卵巢组织中高表达,3个在30周龄蛋鸡卵巢组织中高表达(图2A)。火山图展示了两组间最具显著性的差异基因(|log2FoldChange|≥1且Padj≤0.05)(图2B),可以看到筛选到的10个差异表达lncRNA表达量均较低,符合lncRNA表达模式,但大量基因处于差异不显著状态。对筛选到的10个差异表达lncRNA的qRT-PCR结果显示与转录组测序结果相同(图2C)。
图2 差异表达lncRNA分析A.差异表达lncRNA数量;B.火山图;C.qRT-PCR结果Fig. 2 Analysis of differential expressed lncRNAsA. amount of differential expressed lncRNA;B. volcano map;C. qRT-PCR result
基于位置关系,1 283个lncRNA潜在的靶基因有2 960个,结合表达模式分析,10个差异表达lncRNA中有2个lncRNA具有顺式调控靶基因,分别为MSTRG.14606对应的肌动蛋白γ2基因(actin gamma 2,ACTG2)和MSTRG.15311对应的肌间线蛋白基因(desmin,DES),同时lncRNA与mRNA表达量存在负相关。对lncRNA的功能富集分析结果如图3所示。由图3可见,显著富集的前5条GO通路包括:心率调节(GO:0086091)、结构分子激活(GO:0005198)、蛋白激酶结合(GO:0019901)、离子通道结合(GO:0044325)以及基底外侧质膜(GO:0016323)。显著富集的前5条KEGG通路包括:焦点粘连通路(4510)、mTOR 信号通路(4150)、孕酮介导的卵母细胞成熟(4914)、卵母细胞成熟(4114)以及p53 信号通路(4115)。
图3 功能富集分析A.GO功能富集分析;rich factor表示位于该GO的差异基因个数/位于该GO的总基因数;B.KEGG通路分析; rich factor表示位于该KEGG的差异基因个数/位于该KEGG的总基因数Fig.3 Functional enrichment analysisA. GO enrichment analysis; rich factor represents the DE gene number/total gene number in this GO term;B. KEGG pathway analysis; rich factor represents the DE gene number/total gene number in this KEGG pathway
高通量测序和组学技术研究方式的发展使得研究人员能够对不同的生物过程中的分子调控网络进行筛查,这已经成为解析生物学过程的金标准[9]。本研究利用RNA-seq技术对不同生长发育阶段的蛋鸡卵巢组织的lncRNA进行筛选与鉴定,在20周龄和30周龄京红蛋鸡卵巢组织中,共鉴定到1 283个lncRNA表达,并对lncRNA进行了转录组学特征分析,lncRNA在不同染色体分布不同,且表达量与mRNA存在较大差异,这与其他研究结果一致[10]。同时qRT-PCR技术也验证了RNA-seq的结果,提示RNA-seq在lncRNA研究中可以发挥重要作用。
20周龄为京红蛋鸡的产蛋初期,30周龄京红蛋鸡达到产蛋高峰期,这期间伴随着蛋鸡生长与繁殖过程,存在着复杂的分子调控网络,特别是基因表达的变化,往往会受到非编码RNA的影响[11]。本研究在20周龄和30周龄京红蛋鸡卵巢组织中筛查到10个差异表达lncRNA,提示这两个生长发育阶段lncRNA表达谱差距较小,但在本研究中发现差异表达lncRNA MSTRG.14606和MSTRG.15311对ACTG2和DES可能存在顺式调节作用机制,且lncRNA表达模式与mRNA存在负相关。Felicioni等[12]报道了ACTG2基因与猪繁殖性能有关,MSTRG.14606是否能够通过影响ACTG2表达量来影响蛋鸡繁殖性能,MSTRG.14606和MSTRG.15311能否发挥分子海绵作用机制或是存在其他调控机制,仍有待于深入研究。
GO将采用lncRNA靶基因进而对lncRNA功能进行预测[13-14]。本研究发现,lncRNA在心率调节、结构分子激活、蛋白激酶结合、离子通道结合以及基底外侧质膜等通路中富集。蛋鸡繁殖性能与多种因素有关。蛋白激酶激活通路中AMP激活的蛋白激酶与蛋鸡生产性能、血液生化指标等相关[15],提示被研究发现的差异表达lnRNA可能通过下丘脑轴调节性激素分泌从而影响蛋鸡繁殖性能。细胞组分方面,基底外侧质膜位于前列,该成分与细胞代谢有关,参与代谢物的转运,对维持细胞稳态发挥重要作用[16-17]。蛋鸡繁殖过程中代谢活动旺盛,后续研究中,可围绕lncRNA与蛋鸡代谢过程的关系开展。KEGG通路侧重分析基因产物及其在代谢途径作用中的关系[18-19],本研究发现差异表达lncRNA的靶基因富集在与卵母细胞成熟有关的通路中。蛋鸡性成熟过程中,卵黄前体物:极低密度脂蛋白VLDL和卵黄生成素VTG能够经由卵泡膜受体介导卵母细胞成熟[20]。这一过程受到多因素调节,本研究发现的差异表达lncRNA可能通过调节相关基因表达参与到上述生物学过程中。两种功能富集结果表明,lncRNA可能参与到蛋鸡不同生长发育阶段卵巢组织各种生物学途径中,但是具体的作用机制仍有待于功能研究、体外试验等进一步验证。
本研究表明,从20周龄和30周龄蛋鸡卵巢组织筛选到了lncRNA,其转录组特征与mRNA存在差异,筛选获得的差异表达lncRNA显著富集在蛋鸡繁殖性能相关的生物学过程,细胞组分,分子功能以及多种信号通路中。