近年来，环境内分泌干扰物(environmental endocrine disruptors，EEDs)污染产生的生物毒性日益引起关注，EEDs能干扰DNA甲基化水平、影响细胞增殖，具有生殖毒性，并通过阻碍β-淀粉样蛋白合成、神经递质代谢及改变离子通道等方式造成神经系统损伤[1-2]。EEDs主要包含双酚A(bisphenol A，BPA)、壬基酚和多氯联苯等。BPA作为溴代阻燃剂广泛应用于塑料、纺织和电子等产品中，现已发现其可增加人类患癌症、心脏病和糖尿病的风险[3]。目前，在工业生产中用BPA类似物取代一部分BPA，试图减少BPA的毒性风险，但据Ao等[4]报道，BPA类似物双酚S(bisphenol S，BPS)同样具有内分泌干扰效应和生殖毒性。BPA广泛存在于水环境和沉积物中，已有学者在人类和水生生物样本中检测到一定水平的BPA[5-6]。相对于陆生动物和淡水生物，有关BPA对海洋生物的毒性报道较少。
单环刺螠(Urechis unicinctus)属于螠虫动物门(Echiura)刺螠属(Urechis)，俗称海肠，是中国、俄罗斯、日本和韩国沿海重要的底栖经济动物。单环刺螠肉味鲜美、营养丰富且富含活性物质，已成为开发价值较大的海洋生物资源，但其自然资源日趋枯竭，现已成为新兴的海水养殖对象。关于BPA对单环刺螠的毒性研究仅见丁子媛等[7]检测了急性毒性效应，并发现随着胁迫时间的延长，低浓度BPA暴露对单环刺螠蛋白酶和脂肪酶活力由诱导转为抑制作用。
转录组学是研究基因表达调控的有效方法，可筛选特定环境暴露下的机体差异表达基因(differentially expressed genes，DEGs)，探究生物体或细胞对外界应激的响应规律，实现对污染物的早期预警[8]。目前，水产动物转录组研究多集中于生长发育和病害相关基因的筛选，以及温度、盐度和溶氧等环境因子对基因表达谱的影响[9]，而关于EEDs的影响报道有限。本研究中，采用Illumina Hiseq 2500高通量测序平台，建立了单环刺螠对照组和BPA暴露组的转录组文库，分析了BPA胁迫对单环刺螠基因表达和信号通路的影响，以期为单环刺螠抗逆机制及环境内分泌干扰物对海洋生物的毒性效应研究提供有益参考。
1 材料与方法
1.1 材料
试验用单环刺螠购自连云港海宁路水产品市场，挑选体表完整、活力较强的个体作为试验对象，体长为(13.6±2.7)cm，体质量为(17.3±1.5)g。
RNA提取试剂、反转录和荧光定量PCR试剂盒均购自大连TaKaRa公司。
1.2 方法
1.2.1 试验设计 单环刺螠于室温(15±1)℃、pH为7.9±0.5、盐度为25±1条件下暂养1周，每日更换2/3海水、连续充气，投喂小球藻(Chlorella pacifica)。根据预试验结果所得96 h无死亡发生的最高浓度，设置对照组(海水)和BPA胁迫组，每组设置3个平行，每个平行随机放10条刺螠，试验期间条件与暂养时相同；BPA胁迫组在水体中放入终浓度为3 mg/L的 BPA(助溶剂用量为体积分数为0.15%的无水乙醇，对机体无影响[7])，每天换水后维持BPA浓度。处理15 d后，分别从对照组和BPA胁迫组随机取3条单环刺螠，用针管先吸取抗凝剂再抽取体腔液(抗凝剂与体腔液体积比为1∶1)，4 ℃下以1 000 r/min离心2 min，取体腔液细胞沉淀置于液氮中速冻备用。
1.2.2 总RNA提取与文库构建 采用Trizol法分别提取对照组和BPA胁迫组单环刺螠体腔细胞总RNA。用DNaseⅠ去除基因组DNA，用核酸检测仪Agilent 2100(Agilent Technologie，美国)分析RNA的质量和浓度。RNA-Seq文库由百迈克生物科技有限公司构建，测序平台为Illumina Hiseq 2500。
1.2.3 生物信息学分析 采用Trinity软件对过滤后的测序数据组装得到单环刺螠的单基因序列(Unigene)库。采用BLAST软件将Unigene分别在NR、GO、KOG、Swiss-Prot、COG、eggNOG4.5、KEGG和Pfam等数据库进行序列比对，得到Unigene的功能注释和通路富集信息。采用DESeq2筛选BPA胁迫组和对照组的DEGs，筛选标准为错误发现率(false discovery rate，FDR)<0.01且差异倍数|log2(fold change，FC)|≥1，对筛选出的DEGs进行GO功能分类、COG直系同源分类、KEGG Pathway注释及通路富集分析。
1.2.4 差异表达基因验证 为验证转录组测序(transcriptome sequencing，RNA-Seq)结果的准确性，挑选8个DEGs(4个上调和4个下调基因)进行qRT-PCR检测。qRT-PCR所用特异性引物采用Primer Premier 6.0软件设计(表1)，以β-actin作为内参基因。反应体系(20 μL)：2×One Step TB Green RT-PCR Buffer Ⅳ 10 μL，PrimeScript 1 Step Enzyme Mix Ⅱ 0.8 μL，正、反向引物(20 μmol/L)各0.8 μL，Rox Reference Dye(50×)0.4 μL，cDNA模板2 μL，ddH2O 5.2 μL。反应程序：42 ℃下合成cDNA 5 min，95 ℃下变性10 s；95 ℃下循环变性5 s，60 ℃下延伸30 s，共进行40个循环；95 ℃下变性15 s，60 ℃下退火1 min。采用2-ΔΔCt法计算基因的相对表达量。试验设置3个平行，采用SPSS 25.0软件对数据进行统计分析。
表1 qRT-PCR特异性引物序列
Tab.1 qRT-PCR specific primer sequences
基因 gene引物序列 primer sequence(5′-3′)内参基因 β-actinF:CTCAACCCCAAGGCTAACR:GAAACACCATCTCCGGAA碱性磷酸酶基因phoAF:CTCACAAAGAAGGCAACGR:CTCAAACAAGCCAAGCAA谷胱甘肽过氧化物酶基因GPxF:GGTTCTTGCGTTCCCTTGCR:TTCCGTTGACGATGACTTTACT巨噬细胞迁移抑制因子MIFF:TCYCAGTGCCTTGTTGCTR:ATCCCACATCTTGTTTCG谷胱甘肽巯基转移酶基因GSTF:TATGTCTCCTTGGCTTGTR:CATTTGGGTTTCTTATCC动力蛋白重链8基因DHC8F:CTGCACGCTACCTTCTAATAR:ATCTTGAGCATACCCTCCC锌指CCHC结构域含蛋白7基因ZCCHC7F:AAAACAGCCGTCTACAAGTATCACR:GGAGCCGACACCATCAATA细胞色素P450基因CYP15C1F:GCCGACCTTCATTTAACACCR:ACAACTCAGCCGCCATTTT核仁GTP结合蛋白2基因NOG2F:GGAAATGGGCAAAGTTATAR:TTCTATGTGCTTGGACCTG
2 结果与分析
2.1 转录组测序与序列组装
采用转录组测序技术，分别构建对照和BPA胁迫单环刺螠体腔细胞的cDNA文库，获得13.74 Gb高质量数据，clean data质量值≥30的碱基所占的百分比(Q30)≥93.02%(表2)。
表2 单环刺螠对照组与BPA胁迫组测序数据
Tab.2 Sequencing data of Urechis unicinctus in control group and BPA stress group
组别groupread总数read number总碱基数base numberGC含量/%GC contentQ30/%对照 control22 214 9736 606 590 62247.0193.20BPA胁迫BPA stress23 831 0717 128 491 28847.4293.02
使用Trinity软件对测序数据进行组装，共得到108 262个转录本和30 882个Unigene，转录本和Unigene的N50长度分别为3 187、2 549 bp。Unigene的平均长度为1 380.30 bp，其中，长度为300～500 bp的Unigene有11 288个，占总数的36.88%；长度为500～1 000 bp的Unigene有7 321个，占总数的23.71%；长度为1 000～2 000 bp的Unigene为5 293个，占17.14%；长度>2 000 bp的Unigene有6 880个，占22.28%。
2.2 Unigene功能注释
将Unigene分别与NR、GO、KOG、Swiss-Prot、COG、eggNOG4.5、KEGG和Pfam等8个数据库进行序列比对，共有14 581个Unigene得到注释(表3)。NR数据库注释到的Unigene最多(14 185个)。有7 242个Unigene注释到KEGG数据库，归类到284条代谢通路，基因数量排名前20的通路及Unigene数量依次为核糖体(139)、内吞作用(139)、嘌呤代谢(135)、氧化磷酸化(122)、RNA转运(118)、溶酶体(117)、内质网中的蛋白质加工(115)、碳代谢(113)、泛素介导的蛋白质水解(113)、剪接体(103)、过氧化物酶体(93)、吞噬体(84)、Mtor信号通路(83)、氨基酸生物合成(77)、嘧啶代谢(74)、磷脂酰肌醇信号转导系统(69)、Wnt信号通路(66)、甘油磷脂代谢(61)、真核生物核糖体生物合成(60)和FoxO信号通路(56)。
表3 Unigene注释统计
Tab.3 Unigene annotation statistics
数据库database注释Unigene数量number of annotated Unigene长度≥300 nt的Unigene数量number of Unigene with length≥300 nt长度≥1 000 nt的Unigene数量number of Unigene with length≥1 000 ntCOG4 6478113 836GO6 2291 5624 667KEGG7 2421 6065 636KOG9 1421 9757 167Pfam11 1692 4968 673Swiss-Prot7 9181 6576 261eggNOG4.512 5313 2869 245NR14 1854 15810 027合计 total14 5814 42710 154
通过与NR数据库进行比对发现，与单环刺螠注释匹配的物种中，相似序列比例最高的是海蠕虫(Capitella teleta)7 383个，占序列总数的52.10%，随后依次是鸭嘴海豆芽(Lingula anatina)1 656个(11.68%)、虾夷扇贝(Mizuhopecten yessoensis)582个(4.11%)、太平洋牡蛎(Crassostrea gigas)445个(3.14%)、美洲牡蛎(Crassostrea virginica)392个(2.77%)、厦门文昌鱼(Branchiostoma belcheri)276个(1.94%)、囊舌虫(Saccoglossus kowalevskii)227个(1.60%)、霸王莲花青螺(Lottia gigantea)220个(1.55%)、鲁氏水蛭(Helobdella robusta)217个(1.53%)、双斑蛸(Octopus bimaculoides)212个(1.50%)和其他物种2 562个(18.08%)。
2.3 差异表达基因
BPA胁迫下单环刺螠产生1 073个DEGs，包括649个上调表达基因和424个下调表达基因(图1)。共有838个DEGs得到注释，其中NR为819个，eggNOG为725个，Pfam为691个，KOG为528个，Swiss-Prot为471个，KEGG为399个，GO为350个，COG为321个。
图1 BPA胁迫组和对照组差异表达基因火山图
Fig.1 Volcano map of differentially expressed genes in BPA stress group and control group
2.3.1 差异表达基因的GO功能富集 GO富集分析显示，共有350个DEGs在GO二级功能中被注释分类(表4)。根据其参与的生物学过程(biological process，BP)主要富集在代谢过程(141)、细胞过程(138)和单个有机体过程(110)；细胞组分(cellular component，CC)主要富集在细胞(109)、细胞组分(107)和生物膜(106)；分子功能(molecular function，MF)富集最显著的是催化活性(167)和结合(154)，其他还涉及抗氧化活性、转运体活性、信号转导和电子载体活性等功能。
表4 BPA胁迫组和对照组差异表达基因GO富集分析
Tab.4 GO enrichment analysis of differentially expressed genes in BPA stress group and control group
1级条目entry hierarchy 12级条目entry hierarchy 2注释基因总数total annotated gene number差异基因数DEGs number登录号accession IDcellular process3 310138GO:0009987metabolic process2 859141GO:0008152single-organism cellular process1 84175GO:0044763biological regulation1 00435GO:0065007 localization88837GO:0051179response to stimulus67030GO:0050896cellular component organization or biogenesis56532GO:0071840signaling45715GO:0023052developmental process1186GO:0032502生物学过程biologicalprocessmulticellular organismal process1086GO:0032501 locomotion381GO:0040011 reproduction353GO:0000003biological adhesion304GO:0022610multi-organism process263GO:0051704 reproductive process223GO:0022414 immune system process162GO:0002376 growth161GO:0040007rhythmic process51GO:0048511cell2 358109GO:0005623cell part2 331107GO:0044464membrane1 916106GO:0016020 organelle1 54471GO:0043226membrane part1 12563GO:0044425 macromolecular complex94430GO:0032991organelle part88341GO:0044422endomembrane system25010GO:0012505mitochondrion22513GO:0005739membrane-enclosed lumen21011GO:0031974细胞组分cellularcomponent microtubule associated complex1104GO:0005875chromosome783GO:0005694cell junction331GO:0030054 binding2 820154GO:0005488catalytic activity2 814167GO:0003824 organic cyclic compound binding1 75993GO:0097159 transporter activity39115GO:0005215 structural molecule activity1814GO:0005198 signal transducer activity1416GO:0004871 molecular transducer activity1416GO:0060089nucleic acid binding transcrip-tion factor activity912GO:0001071分子功能molecularfunction electron carrier activity622GO:0009055 metallopeptidase activity593GO:0008237 transcription factor activity, pro-tein binding331GO:0000988antioxidant activity238GO:0016209
2.3.2 差异表达基因COG分类 利用COG数据库对单环刺螠差异表达基因序列进行直系同源分类，富集基因数量较多的功能主要有信号转导机制、核糖体结构与生物合成、翻译后修饰、蛋白质转换与分子伴侣、碳水化合物运输与代谢、脂质转运与代谢及防御机制等(图2)。
图2 BPA胁迫组和对照组差异表达基因COG注释分类
Fig.2 Classification of COG annotations for differentially expressed genes in BPA stress group and control group
2.3.3 差异表达基因KEGG分类及富集 KEGG分析显示，163个DEGs被注释到151条KEGG信号通路，分为细胞过程、环境信息处理、遗传信息处理、代谢和有机系统5大类通路，以代谢和遗传信息处理2大类通路富集的差异基因较多(图3)。在这些代谢通路中，富集差异基因较多的通路为真核生物核糖体生物合成(11)、碳代谢(9)、谷胱甘肽代谢(8)、内质网蛋白质加工(8)、溶酶体(7)、氧化磷酸化(7)、细胞色素P450对外源物质代谢(6)、范可尼贫血通路(6)、内吞作用(6)、药物代谢-细胞色素P450(5)和花生四烯酸代谢(5)等。选取前20位显著富集通路进行KEGG富集分析，结果见图4。富集程度最显著的3条通路为真核生物核糖体生物合成、谷胱甘肽代谢和细胞色素P450对外源物质代谢。与解毒相关的通路主要有谷胱甘肽代谢、细胞色素P450对外源物质代谢、药物代谢-细胞色素P450通路等；与免疫和疾病相关的通路主要有RIG-I样受体信号通路、癌症通路、单纯疱疹感染、Ras信号通路、范可尼贫血通路、丙型肝炎、Rap-1信号通路和HIF-1信号通路等；与代谢相关的通路有真核生物核糖体生物合成、类固醇生物合成、PI3K-Akt信号通路、甾类激素生物合成、初级胆汁酸生物合成和花生四烯酸代谢等；与生殖有关的通路有卵母细胞减数分裂等。
图3 BPA胁迫组和对照组差异表达基因KEGG分类
Fig.3 KEGG classification map of differentially expressed genes in BPA stress group and control group
图4 BPA胁迫组和对照组差异表达基因KEGG通路富集分析(前20)
Fig.4 KEGG pathway enrichment analysis of differentially expressed genes in BPA stress group and control group(top 20)
从表5可见，上调的DEGs主要有碱性磷酸酶、谷胱甘肽过氧化物酶、巨噬细胞迁移抑制因子、细胞色素P450 CYP3A4、谷胱甘肽巯基转移酶和羧酸酯酶等；下调的DEGs主要有细胞色素P450 CYP46A1、多药耐药相关蛋白1、动力蛋白重链8、核仁GTP结合蛋白2、蛋白激酶结构域和脂肪酰基转移酶等。
表5 BPA胁迫下单环刺螠主要差异表达基因
Tab.5 Major differentially expressed genes of Urechis unicinctus under BPA stress
注：—，基因在该数据库无注释。
Note：—，there is no annotation of the gene in this database.
基因编号gene ID基因geneCOG分类注释COG classification annotationKOG分类注释KOG classification annotationlog2(FC)c36888.graph_c1碱性磷酸酶离子转运和代谢离子转运和代谢7.02c33306.graph_c0谷胱甘肽过氧化物酶防御机制;脂质转运和代谢翻译后修饰;蛋白质周转;分子伴侣5.73c40915.graph_c0巨噬细胞迁移抑制因子次生代谢物生物合成、运输和代谢防御机制 5.37c33882.graph_c0细胞色素P450 CYP3A4—次生代谢物生物合成、运输和代谢4.33c38287.graph_c0谷胱甘肽巯基转移酶能量产生和转换翻译后修饰;蛋白质周转;分子伴侣3.96c42783.graph_c0细胞色素P450 10235次生代谢物生物合成、运输和代谢;防御机制 —3.61c28726.graph_c0镧硫氨酸合成酶c样蛋白防御机制 防御机制 3.58c37280.graph_c0羧酸酯酶1C脂质转运与代谢—3.26c29142.graph_c0谷胱甘肽过氧化物酶防御机制;脂质转运和代谢翻译后修饰;蛋白质周转;分子伴侣3.11c36334.graph_c0ADP-糖焦磷酸酶亚型x2 防御机制复制、重组和修复2.81c32162.graph_c0细胞色素P450 CYP46A1次生代谢物生物合成、运输和代谢;防御机制 —-6.40c42774.graph_c0多药耐药相关蛋白1防御机制 次生代谢物生物合成、运输和代谢-4.79c37980.graph_c0动力蛋白重链8 —细胞骨架 -4.62c28201.graph_c0谷胱甘肽巯基转移酶C末端结构域翻译后修饰;蛋白质周转;分子伴侣翻译后修饰;蛋白质周转;分子伴侣-4.60C40050.graph_c0锌指CCHC结构域含蛋白7—翻译后修饰;蛋白质周转;分子伴侣-4.51c38477.graph_c0细胞色素P450 CYP15C1次生代谢物生物合成、运输和代谢;防御机制次生代谢物生物合成、运输和代谢-4.02c41942.graph_c0核仁GTP结合蛋白2翻译、核糖体结构和生物合成—-3.81c34207.graph_c0蛋白激酶结构域信号转导机制信号转导机制-3.47c32808.graph_c0脂肪酰转移酶脂质转运和代谢脂质转运和代谢-3.22c39556.graph_c0类谷氨酰胺转氨酶超家族细胞周期调控、细胞分裂和染色体分区细胞周期调控、细胞分裂和染色体分区-3.05
2.4 差异表达基因qRT-PCR验证
运用qRT-PCR测定挑选的8个基因表达量，并与转录组测序结果进行对比。结果显示，单环刺螠暴露于BPA 15 d后，各基因荧光定量表达结果与转录组测序所得趋势相一致(图5)，表明转录组测序数据可靠。
图5 差异表达基因的qRT-PCR验证
Fig.5 qRT-PCR validation of differentially expressed genes
3 讨论
3.1 单环刺螠体腔细胞转录组组装与注释
由于缺乏相应物种基因组数据作为参考，导致非模式生物转录组研究较为困难，而RNA-Seq技术可通过相似物种同源序列比对进行生物基因信息的解码，弥补了其基因信息少的不足，是分析非模式生物中差异表达基因的有效方法[10]。邢月婷等[11]通过RNA-seq技术从日本医蛭(Hirudo nipponia)唾液腺组织中获得145 981个Unigene，均能得到注释。江红霞等[9]从克氏原螯虾(Procambarus clarkii)肌肉、卵巢和肝胰腺组织中共得到53 006个Unigene，其中有15 956个(30.1%)Unigene得到注释，且与美洲钩虾(Hyalella azteca)(31.28%)序列一致性最高。聂鸿涛等[12]从大竹蛏(Solen grandis)鳃组织中得到190 856个Unigene，有63 337个(33.18%)Unigene得到注释，且与长牡蛎(Crassostrea gigas)序列一致性最高(53.3%)。高琼[13]从刺参(Apostichopus japonicus)体腔细胞得到89 891个Unigene，其中有20 060个(22.32%)Unigene得到注释，且与紫色球海胆(Strongylocentrotus purpuratus)的序列一致性最高(46.7%)。本研究表明，从单环刺螠体腔细胞共获得30 882个Unigene，其中有14 581个Unigene得到注释(47.2%)，且与海蠕虫序列一致性最高(52.1%)。单环刺螠未能获得注释的Unigene有16 301个，表明基因数据库中单环刺螠基因资源仍较缺乏，需要进一步研究予以补充。
3.2 单环刺螠响应BPA胁迫的合成、免疫及解毒代谢相关基因分析
3.2.1 差异表达基因 单环刺螠为开管式循环，与刺参和光裸星虫(Sipunculus nudus)类似，其体腔液相当于血淋巴液，对机体起着重要的防御保护功能[14]。史伟卜等[15]发现，吐脏后的刺参体腔细胞有922个DEGs显著富集到GO分类中的蛋白质水解、细胞黏附、肽链内切酶活性、细胞外区、转化因子-β和FoxO等通路，可能是通过体腔细胞的迁移和黏附，并分泌大量细胞外基质和生长因子，进行吐脏后内脏器官的再生。低氧-复氧胁迫下脊尾白虾(Exopalaemon carinicauda)鳃组织的DEGs显著富集到核糖体、碳代谢、谷胱甘肽代谢和RNA降解等通路[16]。壬基酚处理后的红鲫(Carassius auratus)胚胎DEGs大部分与生长、发育、细胞骨架和血液循环相关，主要富集在黏着斑和转化生长因子-β信号等通路，基因表达异常导致胚胎发育畸形[17]。本研究表明，BPA胁迫下单环刺螠体腔细胞产生1 073个DEGs，包括649个表达上调基因和424个表达下调基因，KEGG富集分析显示，富集程度最显著的3个通路为真核生物核糖体生物合成通路(DEGs 1个上调、10个下调)、谷胱甘肽代谢通路(DEGs 7个上调、1个下调)和细胞色素P450对外源物质代谢通路(DEGs 5个上调、1个下调)。由此可见，不同环境条件对不同物种及不同组织器官的转录组影响效应差异较大。
3.2.2 合成代谢 核糖体在细胞内的作用是合成蛋白质。真核生物核糖体生物合成发生在核仁和细胞质中，其过程涉及200多种在rRNA合成、加工和组装中起协调作用的蛋白质，与细胞生长、增殖密切相关[18]。有研究表明，核糖体生物合成贯穿肿瘤的发生和发展过程，与肿瘤细胞增殖相关的各种因子及癌细胞基因产物也会调节核糖体生物合成过程[19]。细胞应激源的存在会阻断蛋白质翻译过程，并诱导应激颗粒(stress granule，SG)形成，SG的组装有助于应激状态下细胞的存活，因此，SG的产生和动态变化被用来确定应激反应关键因素的起点，H/ACA sno核糖核蛋白复合物具有其他非核功能，可参与SG形成，有助于应激反应[20]。本研究表明，BPA胁迫下核糖体生物合成通路中富集的11个DEGs仅有H/ACA核糖核蛋白复合物亚基3表现为上调，其余均表现为下调，与上述研究结果相一致，表明BPA胁迫会干扰核糖体生物合成和阻断蛋白质翻译。
3.2.3 解毒和防御 谷胱甘肽代谢和细胞色素P450对外源物质代谢两条通路与机体的解毒和防御功能相关。谷胱甘肽能维持机体正常的免疫机能，具有一定的抗氧化作用[21]。GST家族广泛存在于有机体中，在细胞面临氧化胁迫、失衡及死亡时起到重要作用[22]。细胞色素P450酶系在对外源物质的转化和药物代谢中起着重要作用，其中CYP3A亚家族参与了人类50%以上的药物代谢，CYP3A4是肝脏和肠道中含量最丰富的CYP酶，具有广泛的底物特异性[23]。张文强等[24]通过光谱法与计算机模拟方法证明CYP3A4和BPA能发生相互作用，BPA可使CYP3A4发生荧光增强及蛋白质构象发生变化。Saetan等[21]报道沙尖鱼(Silver sillago)在缺氧时，DEGs主要富集于类固醇与氨基酸生物合成、谷胱甘肽代谢、细胞色素P450对外源物质代谢和药物代谢-细胞色素P450途径等通路，且细胞色素P450酶系和GST家族基因表达显著上调。本研究表明，BPA胁迫下单环刺螠DEGs主要富集的通路为真核生物核糖体生物合成、谷胱甘肽代谢和细胞色素P450对外源物质代谢，且后两个通路涉及的大部分基因表达呈现显著上调，表明谷胱甘肽代谢通路和细胞色素P450对外源物质代谢通路及相关基因在单环刺螠对抗BPA胁迫过程中发挥重要作用。
本研究中GO功能分类显示，DEGs主要涉及催化活性、结合、抗氧化活性、转运体活性、信号转导和电子载体活性等分子功能，COG数据库直系同源分类显示，富集DEGs较多的功能主要为信号转导机制、核糖体结构与生物合成、翻译后修饰、蛋白质转换与分子伴侣、碳水化合物运输与代谢、脂质转运与代谢及防御机制等，表明BPA暴露对单环刺螠机体代谢和免疫等产生显著影响。朱龙[25]发现，BPA可诱导稀有鉤鲫(Gobiocypris rarus)精巢产生氧化应激、炎症反应和组蛋白DNA甲基化水平异常，造成精子质量下降、生殖功能异常。Canesi等[26]报道，BPA可影响紫贻贝(Mytilus galloprovincialis)肝胰腺的基因表达，导致参与氧化还原平衡的酶(谷胱甘肽相关酶、过氧化氢酶等)活性和溶酶体膜稳定性降低，免疫功能显著下降。外源性BPA可通过调节人细胞周期蛋白D1(Cyclin D1)和细胞周期蛋白依赖激酶抑制物P21的表达水平，促进甲状腺乳头状癌细胞KTC-1的增殖或凋亡[27]。本研究表明，BPA胁迫下单环刺螠产生的1 073个DEGs中除了649个上调基因外，还有424个基因表达下调，包括细胞色素P450 CYP46A1和多药耐药相关蛋白1基因等。Pikuleva[28]认为，细胞色素P450 CYP46A1酶的功能是消除神经系统内的胆固醇，可将其转化为24-羟基胆固醇，与神经退行性疾病的发生发展密切相关。多药耐药相关蛋白属于ATP结合盒(ATP-binding cassette，ABC)转运蛋白超家族的成员，是质膜上氨基酸、糖、磷脂和肽的转运蛋白，能将代谢废物及药物、毒素等排出细胞外[29]。本研究表明，BPA胁迫下单环刺螠细胞色素P450 CYP46A1和多药耐药相关蛋白1基因表达显著下调，表明BPA降低了机体的脂类代谢及物质的跨膜转运作用，对神经调控及解毒功能都有一定影响。本研究中，仅分析了单一剂量的BPA暴露对单环刺螠转录组的影响，今后的研究中需开展不同剂量BPA和暴露时间的多组学分析，以明确这些DEGs在机体解毒代谢中的调控机制。
4 结论
1)本研究中获得了单环刺螠体腔细胞转录组数据，但未能得到注释的基因数量超过半数，说明基因数据库中单环刺螠基因信息仍较缺乏。
2)通过分析BPA胁迫下单环刺螠体腔细胞转录表达差异，发现差异表达基因主要涉及催化、结合、抗氧化、转运体、信号转导和电子载体等分子功能，说明BPA可影响单环刺螠体腔细胞的生化反应和生理功能。
3)差异基因主要富集在真核生物核糖体生物合成、谷胱甘肽代谢和细胞色素P450对外源物质代谢等通路，说明BPA胁迫对单环刺螠的合成、免疫及解毒代谢产生显著影响，本试验结果为揭示海洋动物BPA胁迫分子调控机制提供了有益参考，也为环境内分泌干扰物的海洋生物毒性研究提供了科学依据。
[1] ISMANTO A,HADIBARATA T,KRISTANTI R A,et al.Endocrine disrupting chemicals (EDCs) in environmental matrices:occurrence,fate,health impact,physio-chemical and bioremediation technology[J].Environmental Pollution,2022,302:119061.
[2] 徐彤.TBBPA-DHEE暴露对斑马鱼的神经毒性、易感性及分子机制研究[D].镇江:江苏大学,2021.
XU T.The neurotoxicity,susceptibility and mechanism of TBBPA-DHEE exposure in zebrafish[D].Zhenjiang:Jiangsu University,2021.(in Chinese)
[3] SHANKAR A,TEPPALA S,SABANAYAGAM C.Bisphenol A and peripheral arterial disease:results from the NHANES[J].Environmental Health Perspectives,2012,120(9):1297-1300.
[4] AO J J,LIU Y J,TANG W F,et al.Bisphenol S exposure induces intestinal inflammation:an integrated metabolomic and transcriptomic study[J].Chemosphere,2022,292:133510.
[5] AINA I G, LYDIA D T, MARY I A.Polycarbonate plastic monomer (bisphenol-A) as emerging contaminant in Nigeria:levels in selected rivers,sediments,well waters and dumpsites[J].Marine Pollution Bulletin,2022,176:113444.
[6] LI Z M, MAO W L, YAO L, et al.First report on occurrence of bisphenol A isomers in human serum and whole blood[J].Journal of Hazardous Materials,2022,424:127549.
[7] 丁子媛,刘顺,王思婕,等.双酚A(BPA)对单环刺螠急性毒性及消化酶活力的影响[J].江苏海洋大学学报(自然科学版),2022,31(1):9-15.
DING Z Y,LIU S,WANG S J,et al.Effects of bisphenol A (BPA) on acute toxicity and digestive enzyme activities of Urechis unicinctus[J].Journal of Jiangsu Ocean University (Natural Science Edtion),2022,31(1):9-15.(in Chinese)
[8] CHUNG J Y,KIM Y J,KIM J Y,et al.Benzo[a] pyrene reduces testosterone production in rat leydig cells via a direct disturbance of testicular steroidogenic machinery[J].Environmental Health Perspectives,2011,119(11):1569-1574.
[9] 江红霞,刘慧芬,马晓,等.转录组测序筛选克氏原螯虾卵巢发育、免疫和生长相关基因[J].水产学报,2021,45(3):396-414.
JIANG H X,LIU H F,MA X,et al.Transcriptome analysis of Procambarus clarkii to screen genes related to ovary development,immunity and growth[J].Journal of Fisheries of China,2021,45(3):396-414.(in Chinese)
[10] KIM C H,GO H J,OH H Y,et al.Transcriptomics reveals tissue/organ-specific differences in gene expression in the starfish Patiria pectinifera[J].Marine Genomics,2018,37:92-96.
[11] 邢月婷,管峰,李诗语,等.日本医蛭唾液腺对饥饿胁迫响应的转录组比较[J].水产学报,2020,44(5):754-766.
XING Y T,GUAN F,LI S Y,et al.Comparative transcriptome analysis of salivary glands of Hirudo nipponia in response to starvation[J].Journal of Fisheries of China,2020,44(5):754-766.(in Chinese)
[12] 聂鸿涛,姜力文,郑梦鸽,等.大竹蛏高通量转录组测序数据组装和分析[J].大连海洋大学学报,2017,32(6):658-663.
NIE H T,JIANG L W,ZHENG M G,et al.Transcriptome data assembly and analysis of grand jackknife clam Solen grandis through high-throughput sequencing[J].Journal of Dalian Ocean University,2017,32(6):658-663.(in Chinese)
[13] 高琼.刺参体腔细胞转录组测序与抗灿烂弧菌感染相关免疫基因的筛选及表达研究[D].青岛:中国海洋大学,2015.
GAO Q.Transcriptome sequencing of coelomocyte and identification and expression of immune-related genes after Vibrio splendidus infection[D].Qingdao:Ocean University of China,2015.(in Chinese)
[14] 郝瑞娟,王庆恒,焦钰,等.光裸星虫体腔液细胞的显微和超微结构[J].海洋湖沼通报,2015(4):40-48.
HAO R J,WANG Q H,JIAO Y,et al.Microstructure and ultrastructure of coelomocytes of Sipunculus nudus[J].Transactions of Oceanology and Limnology,2015(4):40-48.(in Chinese)
[15] 史伟卜,王轶南,王浩文,等.转录组测序解析刺参波里氏囊腔与体腔中体腔细胞对吐脏胁迫的响应差异[J].海洋学报,2021,43(2):116-125.
SHI W B,WANG Y N,WANG H W,et al.Transcriptome analysis provides insights into the response of coelomocytes in Polian vesicle and coelomic cavity of sea cucumber Apostichopus japonicus to evisceration[J].Haiyang Xuebao,2021,43(2):116-125.(in Chinese)
[16] 史文军,王盼,胡润豪,等.低氧-复氧胁迫下脊尾白虾鳃组织差异表达基因趋势分析[J].南方农业学报,2022,53(3):735-747.
SHI W J,WANG P,HU R H,et al.Trend analysis of differentially expressed genes in gill of Exopalaemon carinicauda under hypoxia-reoxygenation stress[J].Journal of Southern Agriculture,2022,53(3):735-747.(in Chinese)
[17] 田雨苏,孙远东,欧密,等.壬基酚致红鲫发育畸形的机制[J].水产学报,2020,44(10):1619-1636.
TIAN Y S,SUN Y D,OU M,et al.Preliminary studies on the mechanism of nonylphenol-induced malformation of Carassius auratus red var[J].Journal of Fisheries of China,2020,44(10):1619-1636.(in Chinese)
[18] LARIZZA L,CALZARI L,ALARI V,et al.Genes for RNA-binding proteins involved in neural-specific functions and diseases are downregulated in Rubinstein-Taybi iNeurons[J].Neural Regeneration Research,2022,17(1):5-14.
[19] DERENZINI M,MONTANARO L,TRER
D.Ribosome biogenesis and cancer[J].Acta Histochemica,2017,119(3):190-197.
[20] BELLI V,MATRONE N,SAGLIOCCHI S,et al.A dynamic link between H/ACA snoRNP components and cytoplasmic stress granules[J].Biochimica et Biophysica Acta Molecular Cell Research,2019,1866(12):118529.
[21] SAETAN W,TIAN C X,YU J W,et al.Comparative transcriptome analysis of gill tissue in response to hypoxia in silver sillago (Sillago sihama)[J].Animals:an Open Access Journal from MDPI,2020,10(4):628.
[22] AL TAMEEMI W,DALE T P,AL-JUMAILY R M K,et al.Hypoxia-modified cancer cell metabolism[J].Frontiers in Cell and Developmental Biology,2019,7:4.
[23] ALAK G,YELTEKIN A Ç,TAS I H,et al.Investigation of 8-OHdG,CYP1A,HSP70 and transcriptional analyses of antioxidant defence system in liver tissues of rainbow trout exposed to eprinomectin[J].Fish &Shellfish Immunology,2017,65:136-144.
[24] 张文强,刘红艳,唐琳,等.细胞色素P450酶CYP3A4与双酚A的相互作用研究[J].分析科学学报,2020,36(6):869-873.
ZHANG W Q,LIU H Y,TANG L,et al.Study on the interaction between cytochrome CYP3A4 and bisphenol A[J].Journal of Analytical Science,2020,36(6):869-873.(in Chinese)
[25] 朱龙.BPA对稀有鮈鲫雄性生殖和子代发育的影响及分子机制研究[D].杨凌:西北农林科技大学,2021.
ZHU L.Impact of BPA on the male reproduction and offspring development in rare minnow (Gobiocypris rasus) and the underlying molecular mechanisms[D].Yangling:Northwest A &F University,2021.(in Chinese)
[26] CANESI L,BORGHI C,CIACCI C.Bisphenol-a alters gene expression and functional parameters in molluscan hepatopancreas[J].Molecular and Cellular Endocrinology,2007,276(1/2):36-44.
[27] 汪洁英,宁博,陈钰玮,等.双酚A对甲状腺乳头状癌细胞KTC-1增殖的影响[J].环境与职业医学,2020,37(7):712-718.
WANG J Y,NING B,CHEN Y W,et al.Effects of bisphenol A on proliferation of thyroid papillary carcinoma KTC-1 cells[J].Journal of Environmental and Occupational Medicine,2020,37(7):712-718.(in Chinese)
[28] PIKULEVA I A.Targeting cytochrome P450 46A1 and brain cholesterol 24-hydroxylation to treat neurodegenerative diseases[J].Exploration of Neuroprotective Therapy,2021,1(3):159-172.
[29] LIU X D.SLC family transporters[J].Advances in Experimental Medicine and Biology,2019,1141:101-202.