基于高通量测序技术筛选与山羊精子发生有关的lncRNAs

代子诺,何其杰,乔 蕾,何莉娜,黄思艺,赵中权 (西南大学 动物科技学院,重庆 北碚400715)

大量的转录本通常超过200个核苷酸,这些转录本通常是多聚腺苷酸化的,并且没有明显的开放阅读框(ORF),它们被定义为长链非编码RNA(lncRNAs)[1]。lncRNAs参与了许多重要的生理过程,如X 染色体失活、遗传印记、染色质修饰和细胞周期调节[2-3]。由于lncRNAs的功能众多且重要,对于lnc RNAs的进一步探索已成为生物信息学的一个重要方向。

之前的研究已经证实,lnc RNAs在睾丸发育和精子发生中起着重要的调节作用。ANGUERA等[4]的研究结果表明了lncRNAtsx在粗线期精母细胞中有特异性表达,但在精原细胞和圆形精母细胞中不表达。小鼠tsx基因的敲除促进了更多精子的凋亡,证实了其对精子发育中的重要作用。

山羊是人类重要的经济家畜,主要提供肉、奶、羊绒、皮等多种产品。与其他山羊相比,大足黑山羊具有较高的繁殖率。大足黑山羊母羊的平均窝产羔数为2.72只,公山羊具有性成熟早、生长发育快、抗病性强、遗传稳定等特点。

本研究采用高通量测序方法分析了成熟和不成熟山羊睾丸中差异表达的lncRNAs,并对靶基因进行了预测。此外,还对这些靶基因进行了分类注释、GO功能分析和KEGG 通路分析,筛选出了与精子发生有关的候选靶基因及lncRNA。这些数据有助于研究人员深入了解睾丸发育过程以及生殖过程,为进一步了解lncRNAs的功能提供了许多有价值的信息。

1 材料与方法

1.1 试验样品及初情期鉴定方法本研究选取9只大足黑山羊按照年龄平均分为3组。ET 组包括3只出生时的山羊(ET1、ET2和ET3),PT 组包括3只初情期山羊(PT1、PT2和PT3),C 组包括3只成年期山羊(AT1、AT2和AT3)。山羊麻醉后,立即收集睾丸,一部分送入深圳华大基因公司进行测序,一部分置入液氮中进行保存。按初情期的定义来确定公畜的初情期在实践中较为困难。目前常采用体质量、睾丸大小的评定来做估测。初情期的鉴定方法:初情期与体质量的关系通常比年龄更为密切,当公羊初次表现出求偶、勃起、爬跨、交配等行为时,对这些公羊进行称重发现其体质量为成年体质量的40%。此外,对这些公羊的睾丸大小进行测量发现,睾丸直径约为6 cm。综合以上特征,判断其为处于初情期的公羊。

1.2 苏木精伊红染色采用4%甲醛溶液和常规组织学方法对3个不同发育阶段的睾丸组织进行固定处理。之后移入二甲苯与酒精混合液处理5 min,依次浸入100%,95%,85%,70%酒精,分别处理2～5 min,蒸馏水转入苏木精染液进行染色5～15 min,蒸馏水洗涤后,转入0.1%～0.5%伊红染液染色1～5 min后经脱水、脱色、封片后加盖玻片封固,在显微镜下进行观察。

1.3 RNA提取以及文库的构建通过TRIzol Plus RNA 纯化试剂盒(Invitrogen,美国)提取总RNA。用NanoDrop-1000分光光度计测定总RNA 的浓度和纯度。利用生物素标记的特异性探针去除核糖体r RNA。纯化之后,将RNA 在一定的温度和离子环境下片段化。随后使用TruSeq®标准试剂盒中随机引物和反转录酶合成cDNA 第1链,之后使用DNA聚合酶I和核糖核酸酶H 来合成双链的cDNA。双链的cDNA 产物随后进行了加“A”碱基和接头连接。连接产物将被扩增,纯化后即得到了最终的cDNA 文库。最后将构建好的测序文库进行测序。

1.4 数据的过滤组装以及编码能力预测过滤掉含有接头,大量低质量碱基,以及含N 过多的读数,得到清洁读数。之后使用比对软件HISAT[5]将清洁读数比对到山羊参考基因组上并用String Tie[6]进行组装。在获得每个样品所有转录本的序列后用Cuffcompare[7]将这些转录本与已知的m RNA 及lncRNA 进行比较,获得它们相互位置关系的信息。并使用蛋白质数据库pfam[8]以及预测软件CPC[9],txCdsPredict,CNCI[10]对转录本的编码进行打分以区分m RNA 和lncRNA。匹配上pfam 的转录本则认为是mRNA,否则为lncRNA。4种方法至少有3种方法结果一致,才会确定该转录本的类型。3款预测软件的打分阈值分别为CPC threshold＝0,CNCI threshold＝0,txCdsPredict threshold＝500。

1.5 基因定量分析及组间差异分析使用RSEM[11]计算转录本的表达量并对基因的表达量进行FPKM 标准化处理。利用差异分析软件DEGseq[12]进行组件差异分析。显著性差异基因的过滤条件为:Fold Change≥2.000和FDR≤0.001。

1.6 IncRNA的靶基因以及GO 功能和KEGG 富集分析计算lnc RNA 与m RNA 的2个相关系数,斯皮尔曼相关性系数与皮尔森相关性系数,要求斯皮尔曼相关性系数≥0.6且皮尔森相关性系数≥0.6。lncRNA 在mRNA 上游10 000内或在mRNA 下游20 000内则判定为顺式调控作用。超出这范围的,将用RNAplex[13]分析lncRNA 与m RNA 的结合能,若结合能＜－30,则判定为反式调控。利用BLAST2GO[14]和Diamond[15]对lnc RNA 显著性差异的靶基因进行GO 和KEGG 功能注释,显著富集的筛选条件:FDR≤0.01。

1.7 逆转录与定量聚合酶链反应根据lnc RNAs序列设计引物避开与靶基因有重叠的区域,引物由上海生工生物公司合成。GAPDH 被用作内参基因。我们使用TaKaRa PrimescriptTMRT 试剂盒从基因组中去除gDNA,并进行RNA 的逆转录。在Bio-Rad CFX96 实时PCR 检测系统中使用TBGrenTMPremix ExTaqTMⅡ(Ta KaRa)进行反应,如下所示:95℃持续30 s,随后40 个95℃循环持续5 s,60℃持续30 s。用2－ΔΔCt法测定lncRNAs的相对表达水平。每个选择鉴定的lncRNAs均进行3次生物学重复,试验所用引物如表1所示。

1.8 统计分析所有qPCR 数据均表示为±s。使用单因素方差(one-way ANOVA)对所有数据进行分析,通过Levene检验数据方差的均匀性,然后进行最小显著性差异(s x) 检验。

2 结果

2.1 不同发育阶段睾丸组织形态我们将收集到的一部分睾丸组织用横切和纵切的方式制作切片,然后用苏木精染色,放置显微镜下进行形态学观察(图1)。我们观察到这3组睾丸的形态学差异很大。在放大200倍的情况下,出生期睾丸的曲细精管直径明显小于初情期和成年期;曲细精管排列较近,空间内结缔组织较小;精原细胞数量稀少,未有精子细胞的产生。初情期内精原细胞数量明显增多,已有精母细胞和精子细胞的产生。在成年期可见大量精原细胞已经发育为初级精母细胞和次级精母细胞,次级精母细胞又发育为精子细胞,部分精子细胞已逐渐分化为了精子。结果表明,3 组处在不同发育阶段的睾丸在形态学上有显著性差异,收集的3组睾丸与我们预期的发育时期相吻合。

图1 不同发育阶段睾丸的显微形态观察 A.出生期山羊睾丸组织横切切片形态(200×);B.出生期山羊睾丸组织纵切切片形态(200×);C.初情期山羊睾丸组织横切切片形态(200×);D.初情期山羊睾丸组织纵切切片在形态(200×);E.成年山羊睾丸组织横切切片形态(200×);F.成年山羊睾丸组织纵切切片形态(200×)

2.2 测序数据概述ET、PT 和AT 组共计分别产生了654 200 616,658 526 698和720 412 344个原始读数,过滤掉含接头和大量低质量碱基以及含N过多的读数后,分别产生了619 375 206(94.68%),620 322 606(94.20%)和688 471 996(95.57%)的清洁读数,这些清洁读数将用于后续分析。总共将近90%的清洁读数被定位到山羊的参考基因组中,表明文库质量良好。

2.3 山羊睾丸中lncRNAs的识别及结构特征我们从CPC,txCdsPredict,CNCI和pfam 分析结果的交叉点中鉴定产生了8 183个新发现的lncRNA(图2)。同时,我们计算了这些lncRNAs转录本的平均长度为8 676.50 bp。然而,计算m RNA 的平均长度为10 376.40 bp。此外,lncRNAs外显子的平均数为2.76,其中62.65%的转录本仅含有1个或2个外显子。然而,47 193个蛋白质编码转录本中的外显子平均数为11.99,远远大于lncRNAs的转录本。表明lnc RNA 编码能力较弱,主要是从表观修饰的层面上对机体进行调控(图3)。

图2 8 183个新的lncRNAs从CPC、txCdsPredict、CNCI和pfam 分析结果的交集中获得

图3 外显子分布图

2.4 各组间lncRNAs的差异分析通过定量分析,我们发现在预测的8 183个新的lncRNAs中,7 808个lncRNAs在样本中表达。3组间lncRNAs的统计分析表明,AT 和ET 组相比共有7 508个lnc RNAs表达差异,其中5 465 个表达差异显著。AT和PT 组相比共有7 513个不同表达的lncRNAs,其中5 392个表达差异显著。PT 和ET 组相比共有5 593个lnc RNAs差异表达,其中只有1 167个表达差异显著(图4)。

2.5 qPCR验证差异表达的lncRNAs我们从差异表达的数据中随机选择9个lnc RNAs,用qPCR 验证其相对表达量。如图5所示,大多数lncRNAs在体内3个阶段表达,表达趋势以及差异显著与我们的测序分析数据高度一致。

2.6 IncRNAs的靶基因预测IncRNA 的功能主要通过顺式或反式作用于靶基因来实现。顺式作用靶基因预测的基本原理认为lncRNAs的功能与其临近的编码蛋白基因有关[16],对lnc RNAs 以及mRNA 的位置关系进行了详细的分类并进行统计(图6),计算lncRNAs与m RNA 的相关系数并判定在m RNA 上游10 000 内或在mRNA 下游20 000 内的lncRNAs。共有4 111个lnc RNAs参与到了顺式调控作用,5 477个m RNA 被lnc RNAs顺式调控。在此范围之外,用RNAplex分析lncRNAs和mRNA 的结合能,最终发现只有445个lnc RNAs参与了反式作用。

2.7 差异表达lncRNAs靶基因的富集分析及与精子发生相关lncRNA的筛选与所有可能被调控的靶基因相比,我们更关注的是组件显著性差异的靶基因功能。为了探索这些靶基因的功能,我们进行了GO 功能和KEGG 通路富集分析。在GO 功能分析中,精子发生、精细胞发育、精细胞分化等功能模块显著富集(图7),提示这些差异表达的靶基因在雄性生殖中可能扮演着重要角色。此外对于KEGG 通路的富集分析发现靶基因主要富集于MAPK 信号通路,HTLV-I信号通路,醛固酮的合成分泌、酮体信号通路,其中MAPK 信号通路,醛固酮的合成分泌、酮体信号通路均与精子发生密切相关。57个候选m RNA 富集在精子发生功能中共靶向到了50个新发现的lnc RNAs(表2),表明这50个lncRNAs可能是调控精子发生的关键基因。

2.8 与精子发生相关lncRNAs的初步鉴定从筛选出的50个与精子发生相关的lncRNAs选择出3个进行qPCR定量分析(图8),结合mRNA 测序与已知的研究结果,对这些lncRNAs及其靶基因进行关联分析。选出的3 个lncRNAs LTCONS_00002707、LTCONS_00052201、LTCONS_00072416在成年期的表达量比其他时期极显著增加,其表达趋势与其靶基因和测序数据基本一致,表明在出生期和初情期精子发生相关的基因表达量比较低,在成年期性成熟时则显著表达。推测这3个lncRNAs有可能是通过促进精子发生相关的mRNA的表达而发挥作用。

图4 不同组中差异表达的lncRNAs和m RNA 数量 A.过滤前差异表达基因数量;B.过滤后差异表达基因数量;1.AT-VSET;2.AT-VS-PT;3.PT-VS-ET

表1 试验所用引物

图5 qPCR 验证9个差异表达的lncRNAs ∗.P＜0.05;∗∗.P＜0.01。下同

图6 lncRNAs和靶基因的位置分布

图7 靶基因的GO 功能富集图谱

3 讨论

在山羊睾丸的研究中,更多的研究集中在蛋白质编码RNA 和miRNA 上,而lncRNAs作为一种重要的调节因子鲜有研究进行测序发掘。在本研究中,我们利用深度测序方法分析了不同发育阶段山羊睾丸中lncRNAs的表达,通过靶基因预测以及组件差异分析寻找出组件差异显著的lncRNAs及其靶基因,将这些显著性差异的靶基因通过GO 功能与KEGG 富集分析以筛选出与精子发生相关的lncRNAs。

表2 筛选出的部分与精子发生有关的lncRNAs及靶基因

图8 与精子发生相关lncRNAs的qPCR 验证

本研究共获得了12.68 Gb的数据,90%的基因映射到了山羊参考基因组上。发掘出新的lncRNAs 8 183个,并筛选获得了50个新发现与精子发生有关的lncRNAs及对应的57个候选靶基因,大多数候选靶基因都位于lnc RNAs附近,通过顺式作用被lncRNA 所调控,提示顺式调控可能是lncRNAs对靶基因的主要调控方式,这与之前报道的一致[16]。这些候选靶基因在精子的形成过程中扮演着重要角色。例如tdtp(睾丸发育相关蛋白)之前发现在睾丸精原细胞中表达量很高,与生殖细胞的减数分裂有关。当tdtp缺失时将会导致精子发生功能停止,此外临床研究也表明,在患有精子生成障碍症的男性中,tdtp基因几乎不表达[17]。zpbp2(透明带结合蛋白)是在精子顶体中特异表达并高度保守的蛋白。之前的研究认为zpbp2 是透明带在精子上的重要受体,主要在精子与卵子的相互作用中行使功能,最近的研究表明zpbp2在精子发生过程中对精子顶体的形成起重要作用[18]。敲除zpbp2会导致生育力下降和精子头部细微畸形,这与zpbp2沿精子吻脊的离散定位相对应。odf2(外周致密纤维2)是精子尾部的促成成分,是维持精子结构的骨架蛋白,在精子运动过程中可以保护鞭毛不受剪切力的作用从而维持精子的运动[19]。此外,tssk家族也是我们筛选出的候选靶基因,它可以动态磷酸化odf2来维持精子的完整性与odf2协同调节精子结构和运动[20]。

从筛选出的lnc RNAs中经qPCR 验证发现它们的表达趋势与靶基因一致,即在出生期和初情期表达量较低,在成年性成熟期表达量显著增高,推测这些lnc RNAs可能促进它们靶基因的表达,即通过促进精子发生相关的m RNA 的表达发挥作用。本研究的睾丸组织染色结果也表示在出生期和初情期几乎未发现明显的精子形态,而在成年期观察到了大量精子的产生,表明这些lncRNAs可能在精子发生中扮演着重要角色。

总之,本研究的结果极大地缩小了精子发生相关lnc RNAs的研究范围,为进一步了解lnc RNAs在睾丸发育和精子发生中的作用提供了重要的数据支持。