雌雄文昌鱼差异表达基因的转录组分析❋

李浩毅， 徐高倩， 刘振辉

(中国海洋大学海洋生命学院， 海洋生物多样性与进化研究所， 山东 青岛 266003)

文昌鱼(Amphioxus 或Lancelet)属于脊索动物门头索动物亚门(Chordate，Cephalocordate)，其在终生成长过程中均具有背神经管、脊索和咽鳃裂[1]。在进化上，文昌鱼是无脊椎动物进化到脊椎动物过程中的一个重要过渡类群，是研究脊椎动物起源与演化的理想的模式生物[2]。

文昌鱼为雌雄异体动物，文昌鱼的性腺成熟后卵巢呈淡黄色，精巢呈乳白色。除此之外，雌雄个体无外形上的差异[3]。陈大元等[4]观察到青岛文昌鱼(B.belcheri)精子中部有多个线粒体存在，且出现了终环和隐窝等高等脊椎动物精子才存在的结构。而在文昌鱼卵子皮层以内的细胞质中含有许多核糖体和线粒体，宋裕昌等曾发现文昌鱼卵黄颗粒内存在一种由9根直径为500～700 μm的小管呈环形排列而成的类似于“微管”的结构[5-8]。与文昌鱼生殖有关的生理活动如性腺的发育、成熟，以及其排卵、释精等过程受温度、光照强度、海水盐度以及性激素等多种因素的影响[9-10]。张致一等在文昌鱼发育成熟的性腺中，检测到雌激素和雄激素，这一结果证明性腺发育至成熟的过程中有性类固醇激素的参与。方永强等发现文昌鱼生长发育从小生长期至大生长中期的过程中，17β-雌二醇、睾酮和孕酮水平呈上升趋势，大生长中期之后开始下降，表明在文昌鱼卵黄生成和精子发生的过程中，性类固醇激素参与了调控[11-14]。另外，有通过对文昌鱼染色体的核型与带型分析发现，青岛文昌鱼的染色体有36条，其中第2对染色体可能是性染色体[15]。Shi等[16]则通过性染色体连锁基因nodal突变体繁育后代的性别比例分析显示，佛罗里达文昌鱼的性染色体属于ZW型。

虽然研究文昌鱼领域的学者们从不同方面报道了文昌鱼精子和卵子的发生、胚胎发育以及器官的形成[17-27]，但关于文昌鱼的性腺发育和性别分化的分子调控机制等科学问题还知之甚少，甚至连不同性别的文昌鱼在基因表达上有哪些差异也不清楚。本文通过转录组测序的手段，探讨了雌雄文昌鱼基因的差异表达情况，发现雌雄差异表达的基因主要富集在神经活性配体-受体相互作用(Neuroactive ligand-receptor interaction)和细胞粘附分子(Cell adhesion molecules)通路中，为从基因水平上研究文昌鱼雌雄性别的差异乃至性别分化的分子机制提供了基础信息。

1 材料与方法

1.1 实验材料样品采集

青岛文昌鱼(Branchiostomajaponicum，曾用名B.belcheritsingtauense)取材于山东省青岛市崂山区沙子口附近海区。6月份采集的文昌鱼处于性腺成熟期，可肉眼观察到雌性文昌鱼体内黄色颗粒状卵巢和雄性文昌鱼体内乳白色颗粒状精巢。根据文昌鱼性腺特征将雌雄文昌鱼分开喂养，并做好标记，将细沙用灭菌海水淘洗多次以除去细沙中残存的杂质，将细沙铺至收纳箱中至一定厚度，再倒入过滤的海水，加入充氧泵，分别放入雌雄文昌鱼，并标记如下：排卵前雌鱼(F1组)、释精前雄鱼(M1组)、排卵后雌鱼(F2组)、释精后雄鱼(M2组)。

1.2 方法

1.2.1 RNA提取 每个处理组分别取3只青岛文昌鱼全鱼样本，将其置于DEPC处理后的1×PBS 缓冲液中清洗3遍，在液氮中研磨全鱼，研磨至粉末状后，用Invitrogen公司的TRIzol®Reagent试剂盒提取样品总RNA。每组将3个 RNA 样品等浓度混合后用于文库构建。提取的 RNA 纯度均满足OD值(260/280)在1.8～2.2之间，且含量大于3.5 μg，浓度大于200 ng/μL。

1.2.2 RNA测序和质量控制 RNA样品送到安诺优达公司(中国北京)完成测序。待RNA样品质检合格后进行制备测序文库。测序得到的某些原始下机序列，会含有测序接头序列以及低质量序列。为了保证信息分析数据的质量，作者对原始序列进行过滤，去除接头污染的Reads(Reads中接头污染的碱基数大于5 bp；去除低质量的Reads(Reads中质量值Q ≤ 19的碱基占总碱基的50%)；去除含N比例大于5%的Reads；对于双端测序，一端测序不符合上述要求，则去掉两端Reads。得到质量较高的Clean Reads，再进行后续分析，后续分析都基于Clean Reads。

1.2.3 差异表达基因GO和KEGG富集分析 由于雌雄文昌鱼差异表达基因众多，本文设置了较高的基因差异表达的标准，将基因显著差异表达的参数设置为|log2FoldChange|≥5和q<0.05。转录组序列通过计算RPKM(Reads Per Kilobase Millon Mapped Reads)值来确定基因的表达量[28]。RPKM计算公式为：

式中：RPKM(A)为基因A的表达量；R为唯一比对到基因A的Reads数；N为唯一比对到基因的总Reads数；L为基因A的长度。

此外，利用Gene Ontology(GO)[29]、Kyoto Encyclopedia of Genes and Genomes(KEGG)[30]等数据库对差异表达基因进行功能注释。利用DAVID在GO数据库中, 通过计算差异表达基因参与生物学过程、细胞组分、分子功能的超几何分布，进行功能分类注释和富集分析。同时，用KEGG数据库获得Unigenes的代谢通路注释信息。以上分析均P<0.05作为统计学显著性的阈值。

2 结果

2.1 测序结果的组装

对于4个组别F1、F2、M1、M2的测序结果进行分析。测序组装后的所有Reads 数分别为64 379 874、61 564 104、63 850 758、60 502 630。去掉低质量序列(当总碱基的50%的Reads是质量值Q≤19的碱基时)和接头污染后(Reads中接头污染的碱基数大于5 bp)的Reads数为63 179 798、60 446 482、62 606 300和59 280 314。所有组别的组装错误率均小于等于0.04%，Q30均大于94%(见表 1)，说明测序与组装质量较好。

表1 测序组装结果统计

2.2 组间差异表达基因分析

对各组间差异表达的基因进行聚类分析，文昌鱼排卵与释精前组(F1组和M1组)聚为一类，文昌鱼排卵与释精后组(F2组和M2组)聚为一类，文昌鱼排卵与释精前后的基因表达差异明显(见图 1)。

(根据差异表达基因在每个样品里的表达量，取以2为底的对数后，计算欧氏距离，再利用系统聚类法(Hierarchical Cluster)，最终得到样品的整体聚类结果。在图中，表达量的变化用颜色的变化表示，蓝色表示表达较低，黄色表示表达量较高。According to the expression level of differentially expressed genes in each sample, the logarithm base 2 is taken to calculate the Euclidean distance, and then the Hierarchical Cluster method is used to obtain the overall clustering result of the sample. In the figure, the change of expression level is represented by the change of color, with blue indicating the low expression level and yellow indicating the high expression level.)

分别对排卵和释精前的F1组与M1组、排卵和释精后的F2组与M2组的转录组测序结果中的显著差异表达基因进行分析, 发现F1组与M1组有13 434个显著差异表达基因，其中上调基因5 737个，下调基因7 697个；F2组与M2组有5 957个显著差异表达基因，其中上调基因2 512个，下调基因3 445个。通过火山图将差异表达基因可视化(见图 2)。

(横坐标代表同样品组中的基因表达倍数变化，纵坐标代表基因表达量的变化显著程度，不同颜色表示不同的分类(上调、下调和不变)。The abscissa is the fold change of gene expression in different sample groups, and the ordinate is the statistically significant degree of gene expression change. Different colors indicate different classifications (up-regulated, down-regulated and unchanged).)

2.3 GO功能分析

对差异表达的基因进行GO(Gene Ontology) 富集分析，分别对F1、M1 组和F2、M2中差异表达基因的前66个GO term进行列举，并对富集到GO中的term进行显著性分析，发现F1与M1 组有295个GO term 差异富集，其中9个GO term极显著差异富集，而F2与M2组有376个GO term 差异富集，其中3个GO term极显著差异富集。

在F1与M1 组中的差异表达基因，在生物过程方面，大部分富集在生物粘附、生物调节等过程，如DNA的复制，细胞粘附；在分子成分方面，集中富集在细胞器、细胞膜及分子复合物等组分；在分子功能方面，主要富集在细胞过程、代谢过程，如跨膜转运、钙离子结合、白细胞的跨内皮迁移等。F2与M2 组的差异表达基因，生物过程部分主要以单一生物过程、对刺激的反应、细胞表面受体信号通路等功能为主；细胞组分中主要以细胞膜、细胞器等组成成分以及细胞外区为主；分子功能主以转移酶活性、RNA结合、金属离子结合为主(见图 3)。

(①细胞过程；②单一生物体过程；③代谢过程；④生物调节；⑤发育过程；⑥刺激应答；⑦多细胞生物过程；⑧细胞组分组成或生物合成；⑨定位；⑩生物粘附；生殖过程；运动；多生物过程；免疫系统过程；生长；信号；行为；生物阶段；生物间歇过程；生殖；参与化学突触传递的突触前过程;细胞聚集；解毒作用；细胞杀伤；细胞成分；膜成分；细胞器；细胞膜；细胞器成分；大分子复合体；细胞外区域成分；细胞外区域；细胞连接；突触成分；细胞膜内腔；超分子复合体；突触；其它生物成分；病毒体成分；类核；病毒体；线粒体相关粘附复合体；病毒包埋体；细胞；合胞体；其它生物；结合；催化；分子传递；运输；信号传递；结构分子；细胞核结合转录因子；分子功能调节因子；电子传递；转录因子活性，蛋白结合；抗氧化剂；翻译因子；化学排斥物；化学引诱物；形态发生；金属伴侣蛋白；丙氨酸酰基载体蛋白；营养储存；蛋白标签。横坐标为GO的三个大类及其具体的子类，左侧纵坐标为差异基因在该子类中所占的比例，右纵坐标为该子类富集到的基因数。不同的颜色代表不同的组别。The abscissa is the three major classes of GO and their specific subclasses. The left ordinate is the proportion of different genes in this subclass, and the right ordinate is the number of genes enriched in this subclass. Different colors represent different groups. )

2.4 KEGG聚类分析

对于差异表达的基因进行KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析，得到众多差异富集的通路。F1与M1 组的差异表达基因在33个KEGG通路中较为富集，其中有7个极显著差异富集通路，分别为神经活性配体-受体相互作用、细胞粘附分子、造血细胞系、肥大型心肌病、细胞外基质相互作用、PI3K-Akt信号通路、金黄色葡萄球菌感染；而F2与M2 组的差异表达基因在32个KEGG通路中显著富集，其中有3个极显著差异富集通路，包括神经活性配体-受体相互作用、细胞粘附分子、病毒性心肌炎(见图 4)。

(①病毒性心肌炎；②弓形虫病；③金黄色葡萄球菌感染；④小细胞肺癌；⑤沙门氏菌感染；⑥核糖体；⑦类风湿性关节炎；⑧肌动蛋白骨架调节；⑨癌症的蛋白聚糖；⑩PI3k-Akt信号通路；吞噬体；百日咳；氮代谢；尼古丁成瘾；神经活性配体-受体相互作用；疟疾；白细胞跨内皮迁移；利什曼病；军团病；肥大心肌病；低氧诱导因子通路;造血细胞系；谷氨酸突触；氨基丁酸突触；细胞外基质的相互作用；DNA复制；扩张型心肌病；补体系统；细胞黏附分子；上皮细胞细菌入侵；轴突导向；右室心肌病；抗原加工与递呈；b: F2、M2组的差异表达基因 KEGG 通路。①病毒性心肌炎；②弓形虫病；③金黄色葡萄球菌感染；④小细胞肺癌；⑤沙门氏菌感染；⑥唾液分泌；⑦类风湿性关节炎；⑧肌动蛋白骨架调节；⑨癌症的蛋白聚糖；⑩PI3k-Akt信号通路；吞噬体；百日咳；癌症通路；破骨细胞分化；神经活性配体-受体相互作用；疟疾；亚油酸代谢；白细胞跨内皮迁移；利什曼病；肥大心肌病；低氧诱导因子通路;造血细胞系；糖鞘脂生物合成神经节；细胞外基质相互作用；扩张型心肌病；细胞粘附分子；钙信号通路；上皮细胞细菌感染；轴突导向；右室心肌病；花生四烯酸代谢；黏附连接。纵坐标为通路名称，横坐标为富集因子，圆圈颜色代表 q value 值，颜色越红代表显著富集性越可靠，圆圈越大代表富集的基因数目越多。The ordinate is the name of the pathways， the abscissa is the enrichment factors.The circle color represents q value.The color is redder， the significant enrichment is more reliable.The circle is larger means the number of genes are greater.)

此外，作者也分析了排卵前后(F1组与F2组)和释精前后(M1组与M2组)的显著差异表达基因，结果表明排卵前后的文昌鱼(F1组与F2组)有11 288个显著差异表达的基因，释精前后的文昌鱼(M1组与M2组)有8 476个显著差异表达的基因。对这些差异表达基因进行KEGG通路聚类分析，发现文昌鱼排卵前后(F1组与F2组)，差异表达的基因极显著地富集到4个KEGG通路中，分别为神经活性配体-受体相互作用、DNA复制(DNA replication)、核糖体合成和细胞粘附分子。而文昌鱼释精前后(M1组与M2组)，差异表达的基因极显著地富集到3个KEGG通路中，包括细胞粘附分子、金黄色葡萄球菌感染和神经活性配体-受体相互作用。

3 讨论

性别决定和分化是由许多不同基因控制的复杂过程，尽管人们对果蝇、小鼠等模式动物的性别分化基因有较深入的了解, 但由于不同的物种在性别分化机制上存在很大的不同，如今还很缺乏对从低等动物到人类性别分化和发育机制的进化规律的认识。文昌鱼作为从无脊椎动物到脊椎动物进化的一个过渡类群, 在研究性别分化和发育的分子机制方面有其独特的优势。然而目前关于文昌鱼的性腺发育和性别分化的分子调控机制等科学问题还知之甚少，甚至连不同性别的文昌鱼在基因表达上有哪些差异也不清楚。为此, 作者对雌雄文昌鱼进行了转录组测序，鉴定雌雄文昌鱼差异表达基因，为进一步研究文昌鱼的性别分化和性腺发育的分子机制及其进化奠定基础。

通过分析转录组测序数据，作者获知不同性别的文昌鱼至少有5 957个显著差异表达的基因。通过对几个在其它物种中已有报道的基因进行查验，验证了转录组测序结果的正确性。比如，促性腺激素释放激素基因(Gnrh)可刺激促黄体生成素的产生，进而促进排卵和精子的释放[31]，其在雌性文昌鱼中要比在雄性文昌鱼中的表达高，这与斑马鱼中Gnrh3的表达模式相一致[32]。再如，锌指蛋白Glp-1(GATA-like protein 1)是卵子发育命运的重要细胞决定因子，其在雌性小鼠中的表达显著高于在雄性中的表达[33]。本研究的转录组数据也显示该基因在雌性文昌鱼中的表达显著高于其在雄性文昌鱼中的表达。

本研究利用 KEGG 数据库分别对成熟期雌雄文昌鱼的差异表达基因进行聚类后分析发现，这些在文昌鱼不同性别中差异表达的基因，还与鱼的生长发育、免疫应答和新陈代谢等生命活动密切相关。排卵和释精前，雌雄文昌鱼的差异表达基因主要聚类在细胞和代谢过程，其中神经活性配体-受体相互作用富集的差异基因数最多，细胞粘附分子、细胞外基质相互作用、PI3K-Akt信号通路富集的差异基因数较多；而排卵和释精后，雌雄文昌鱼的差异表达基因主要集中在代谢和组织系统过程，包括白细胞迁移、轴突导向、吞噬体、神经活性配体-受体相互作用、细胞粘附分子等通路上，其中神经活性配体-受体相互作用和细胞粘附分子富集的差异基因数最多。

KEGG通路分析显示，无论排卵和释精前还是排卵和释精后，雌雄文昌鱼的差异表达基因极富集的共同信号途径是神经活性配体-受体相互作用和细胞粘附分子信号通路。在鱼类、小鼠等脊椎动物中，有一类性分化关键基因AMH等属于TGF-β细胞因子超家族，其也是通过配体-受体相互作用发挥作用的，说明性别分化的分子机制在进化上存在保守性[34-35]。再如，在斑马鱼中，Gnrh3可通过MAPK途径调控PGCs(原始生殖细胞)的增殖，在早期性别分化中承担重要功能[36]，这整个过程离不开一系列配体-受体的相互作用。另一在雌性中高表达的绒毛膜促性腺激素受体(LHCGR)基因的纯合子错义突变会导致空卵泡综合征和XY性发育障碍。研究表明，LHCGR启动子区域的去甲基化是调节卵泡发育过程中细胞类型特异性分化的关键机制[37-38]。

对文昌鱼雌雄差异表达的相关基因进一步分析发现，与血管内皮生成相关的基因如F11R、JAM2、CDH5等，无论排卵和释精前还是排卵和释精后，雌性都比雄性文昌鱼呈显著的高表达(见图 5a)；而与白细胞免疫相关的基因如ITGB1、SELP、SELE等，以及与神经系统的神经突触相关的基因如CDH2、NCAM、L1CAM、PTPRF、SDC1等，排卵和释精前，其在雌性要比雄性文昌鱼中表达高，而排卵和释精后，其在雌性要比雄性文昌鱼中表达低(见图 5b)。其中，F11R、CDH2等基因主要参与细胞粘附过程，并未见它们在性别分化过程中发挥作用的报道，而本研究结果暗示其在文昌鱼性别分化及决定过程中可能有一定的作用。

(①免疫系统；②抗原递呈细胞(直流巨噬细胞)；③T细胞；④肽；⑤T细胞受体信号通路；⑥毒性T细胞；⑦靶细胞；⑧辅助T细胞；⑨B细胞；⑩上皮细胞；紧密连接；白细胞跨内皮迁移；白血球；血小板；补体系统；神经系统；神经元(突触前)；神经元(突触后)；神经元(生长锥)；神经元(轴突)；施万细胞；节点；副节点；近节点；施万细胞(髓磷脂)；周期线；少突细胞(髓磷脂)；其它系统；黏附连接；精细胞；足细胞；纤毛体；色素上皮；无色素上皮；成肌细胞；凯尼海撒实验室；细胞粘附分子。a: F1与M1组间差异表达基因富集到的通路图注释示例。b : F2与M2组间差异表达基因富集到的通路图注释示例。红色表示注释到该KO的基因为上调基因，绿色表示注释到该KO的基因为下调基因。a :Examples of annotated pathway diagrams in which differentially expressed genes between F1 and M1 groups are enriched. b : Examples of annotated pathway diagrams in which differentially expressed genes between F2 and M2 groups are enriched. The red indicates that the gene annotated to the KO is up-regulated, and the green indicates that the gene annotated to the KO is down-regulated.)

本研究通过转录组分析获得的雌雄文昌鱼差异表达的基因，为从基因水平上研究文昌鱼雌雄性别的差异乃至性别分化的分子机制提供了基础信息。