三唑酮对大型溞代际影响的转录组学分析

侯琳，徐建,2，梁雪芳，金小伟

1. 中国环境科学研究院环境健康风险评估与研究中心，北京 100012 2. 中国环境科学研究院国家环境保护化学品生态效应与风险评估重点实验室，北京 100012 3. 中国环境监测总站，北京 100012 4. 内蒙古大学生态与环境学院，呼和浩特 010021

三唑酮是一种三唑类杀菌剂，占全球三唑类杀菌剂消费量[1]的30%以上。三唑酮在水环境中的普遍存在及其对水生生物的潜在不利影响已经引起广泛关注。由于农药的有效利用率极低[2]，大量残留的三唑酮进入水环境，其使用区域的临近水体检测出了较高浓度的三唑酮[1,3]。在中国的九龙江、太湖和白洋淀等地均有三唑酮的检出，最高检出浓度为12 μg·L-1[3-6]。三唑酮具有稳定性强，流动性和吸附性好等特点，可在水生生物体内累积并发生生物转化[7-9]，对水生生物产生毒性效应。例如，三唑酮可影响生物体的发育[10-12]、繁殖[13-14]，引起神经毒性[15]、肝毒性[16]、心血管毒性[17]和氧化应激等[18-19]。三唑酮主要影响两栖动物的变态发育，农药的使用可能是两栖动物种群数量下降的原因之一[20]。三唑酮通过改变下丘脑-垂体-甲状腺轴(hypothalamic-pituitary-thyroid, HPT)相关基因的表达来干扰甲状腺激素的分泌，影响胚胎的孵化率，降低鱼类[13]和两栖动物[10]的繁殖能力。三唑酮影响两栖动物和鱼类生殖能力的另一途径是抑制睾酮向雌二醇的转化[12,21]。目前，关于三唑酮的毒性机制研究主要集中在高营养级的脊椎动物，如两栖动物和鱼类，而其对低营养级水生无脊椎动物的毒性效应研究相对较少，刘娜等[22]关于三唑酮毒性效应的研究发现三唑酮对大型溞蜕皮产生毒性作用，且其敏感性高于生长和繁殖，并进一步整理出三唑酮的有害结局路径[23]，结果表明三唑酮可能通过抑制P450酶的活性从而降低甲壳类动物蜕皮激素的活性。目前仍缺少关于三唑酮对大型溞代际影响的研究以及分子机制上的探究。

由于基因组序列的可获得性、较短的生命周期和易于处理等优点，大型溞成为生态毒性基因组学研究的优秀模式生物[24]，可用于评估化学品对水生生态系统的毒性影响[25]。大型溞处于藻-溞-鱼水环境食物链的中间位置，占据相对重要的生态位[26]，因此，有关大型溞的毒性效应研究可能对整个水生环境食物链具有一定的参考作用。

本研究在整理归纳三唑酮的毒性效应和机制的基础上，开展大型溞的多代试验，探索三唑酮对大型溞的生长和繁殖毒性以及相应的毒性作用机制和途径。试验期间记录并分析大型溞的体长及产溞数量等指标，使用HiSeq平台进行转录组测序，分析处理组与对照组的差异表达基因及其功能富集通路，为探索三唑酮对水生无脊椎动物的毒性机制提供参考。

1 材料与方法(Materials and methods)

1.1 培养条件

试验所用大型溞为本实验室(中国环境监测总站)连续稳定培养5代以上的单克隆品系，个体之间差异较小。按照OECD 211 (Organization for Economic Co-Operation and Development)[27]的要求，大型溞培养使用Elendt M4培养基(配制所需化学药品均购自上海国药集团化学试剂有限公司)，置于恒温培养箱中，培养温度为(21±1) ℃，每天光照16 h，黑暗8 h，光照强度在1 110～1 480 lx。每日以羊角月牙藻喂养，投饵密度为3.0×105～4×105cells·mL-1。藻类按照OECD 201[28]要求培养，培养箱温度设为(23±1) ℃，每天光照10 h，黑暗8 h，光照强度在4 440～8 880 lx，培养4～5 d后在5 000 r·min-1转速下离心浓缩5 min，使用灭菌水稀释以降低营养盐浓度。

1.2 试验设计

三唑酮购于阿拉丁有限公司(Aladdin，中国上海)，CAS号为43121-43-3，纯度为99.8%。丙酮购自康林科技有限公司(中国北京)。根据ECOTOX数据库(https:// cfpub.epa.gov/ecotox/search.cfm)中三唑酮毒性效应浓度和以往研究基础[22]，设置6个浓度梯度，分别为5、12.5、25、50、100和200 μg·L-1(浓度从低到高依次标记为3～8)，每个浓度设置5个重复。以丙酮(<0.1 mL·L-1)为助溶剂，使用M4培养基配制三唑酮母液。试验空白对照分别为M4培养基和含有0.02 mL·L-1丙酮的M4培养基(分别标记为1和2)。对照组未检出三唑酮。

采用图1中描述的试验设计进行三唑酮暴露的多代试验(F0、F1、F2)。试验开始前，选取20只个体大、游泳能力强且怀卵多的母溞置于1 L装有M4培养基的烧杯中，驯化培养。试验开始时，随机选取上述母溞所产<24 h的第三胎幼溞32只分别置于含有80 mL暴露液的100 mL烧杯中，每个烧杯一只溞，记为F0代。F0代所产<24 h的第三胎幼溞32只被选作子代F1，F1代所产<24 h的第三胎幼溞32只被选作子代F2。选择第三胎是因为早期的幼溞(第一胎和第二胎)存活不稳定。每代按照OECD 211[27]慢性毒性试验标准进行21 d半静态暴露试验，每天更换培养液取出新生儿并记录F0代大型溞的每胎产溞数量和时间以及每代相隔时间，暴露21 d后在显微镜下测定大型溞体长(不包括尾刺)，并取1 mL暴露液过0.22 μm滤膜，随后保存在-20 ℃下待测。F1代和F2代出生7 d后进行转录组测试和基因表达分析。

1.3 三唑酮暴露液浓度检测

暴露21 d后取暴露液1 mL，经过前处理去除杂质和影响物质，使用LC-MS/MS-8040 (岛津，日本)检测暴露液中三唑酮实际浓度。液相色谱柱为ACQUITY UPLC®BEH C18 (1.7 μm, 50 mm×2.1 mm, Waters, MA, USA)，进样量为5 μL，流动相由0.1%甲酸水溶液和乙腈组成。流动相流速为0.3 mL·min-1，柱温为40 ℃。洗脱梯度如表1所示。

定量分析使用LC-MS/MS-8040 (岛津，日本)完成，采用多反应监测(multiple reaction monitoring, MRM)与电喷雾正离子源(ESI+)模式。氮气作为脱溶剂和雾化气。母离子294 (m/z)，子离子197 (m/z)。名义浓度为5.00、12.50、25.00、50.00、100.00和200.00 μg·L-1的三唑酮暴露液实测浓度分别为4.72、11.36、32.02、63.24、116.79和229.00 μg·L-1。

1.4 转录组测序

每个处理组随机挑选约40只大型溞作为一个生物学重复，使用TRIzol(Invitrogen，美国)提取大型溞总RNA，并使用DNaseⅠ(TaKara，大连，中国)去除基因组DNA。利用2100 Bioanalyser(Agilent，美国)、ND-2000微量分光光度计(Thermo Scientific，美国)对所提取RNA的浓度和纯度进行检测，保证所有RNA样品的完整性(OD260/OD280=1.8～2.2, OD260/OD230≥2.0, RIN≥8.5, 28S/18S≥1.0)并进行转录组测序。mRNA测序基于HiSeq平台，采用Illumina TruSeqTMRNA sample prep Kit方法进行文库构建，主要操作流程为：提取总RNA(>1 μg)，富集mRNA，mRNA片段化，反转成cDNA，adaptor连接，Illumina测序。

图1 多代试验示意图Fig. 1 Diagram of multi-generational experiment

1.5 差异表达基因

使用RSEM(https://deweylab.biostat.wisc.edu/resm/)对基因的表达水平进行定量分析，获得基因/转录本的Read Counts数后，使用DEGseq (http://bioconductor.org/packages/stats/bioc/DESeq/)软件进行样本间基因的表达差异分析，以FPKM(fragments per kilobases per millionreads)衡量差异表达的显著性，使用错误发现率(false discovery rate, FDR)和差异倍数(fold change, FC)作为评判标准，当一个基因同时满足FDR<0.05和|log2FC|>1时，视该基因为差异表达基因(differentially expressed gene, DEG)。

1.6 差异表达基因的KEGG富集

使用KOBAS(https:// kobas.cbi.pku.edu.cn/home.do)对不同组别的差异表达基因进行功能富集分析，在基因功能水平上阐述三唑酮暴露对大型溞可能的作用机制和途径。利用KEGG数据库，可将基因按照参与的通路或行使的功能分类，能够识别出与生物现象最相关的生物学过程。使用BH(Benjamini and Hochberg)方法进行多重检验矫正P值，P≤0.05，认为此GO富集功能或KEGG通路存在显著富集情况。

1.7 统计分析

本研究涉及的数据统计分析使用Origin 2021(OriginLab Corporation, Northampton, MA, USA)，The R Programming Language 4.0.3和IBM SPSS Statistics 26(BM Corporation, Armonk, NY, USA)。使用单因素方差分析进行显著性检验和方差同质性检验，重复组取算术平均值，表示为平均值±标准差(standard deviation, SD)。P<0.05表示两样本间存在显著性差异。

2 结果(Results)

2.1 三唑酮对大型溞F0代生长和繁殖的影响

F0代的体长是大型溞的生长指标，暴露结束后，在显微镜下测定不同三唑酮暴露浓度下大型溞的体长(不包括尾刺)，结果显示，仅在200 μg·L-1浓度时体长明显降低(P=0.000)，其他浓度下体长变化不显著(图2(a))。本研究大型溞为孤雌生殖，将F0代大型溞21 d平均每胎产溞量作为大型溞的繁殖指标，结果显示，在浓度为100 μg·L-1和200 μg·L-1时21 d平均每胎产溞量均明显降低(P=0.014,P=0.000)，而其他浓度下的繁殖能力差异不明显(图2(b))。除体长和繁殖外，本研究也记录了不同浓度下F0代产下一代F1以及F1代产下一代F2所需要的时间间隔。通过组间方差分析发现，浓度为5、12.5、25、50和100 μg·L-1时F0从出生到产下F1的间隔时间和F1从出生到产下F2的间隔时间增加不明显，而当浓度为200 μg·L-1时两代世代间隔时间均显著增加(P=0.008，P=0.000)(表2)。比较不同浓度下两代世代间隔时间的差异可知(表2)，与对照组相比，12.5 μg·L-1(P=0.033)和100 μg·L-1(P=0.008)浓度处理下世代间隔时间显著增加。

2.2 差异表达基因

F1和F2代大型溞在三唑酮中暴露7 d后，进行转录组测试分析，共获得162.94 Gb Clean Data，且各样品Clean Data均达到6.38 Gb以上，Q30碱基百分比在93.5%以上。转录组测试共检测到20 187个基因。分别将各样品的Clean Reads与指定的参考基因组(https://www.ncbi.nlm.nih.gov/gen-ome/10953?genome_assembly_id=434305)进行序列比对，比对率为79.02%～91.76%。基于各浓度下基因表达量，进行样本间相关性分析和PCA分析，分别如图3和图4所示，样本间相关性分析表明样本间相关性较高，PCA分析可以发现高浓度(200 μg·L-1)处理组与其他浓度组间有较大差异，同时，F1代和F2代间的转录组响应也存在显著差异。因此我们对F1代、F2代以及F1 vs F2的基因表达情况进行了进一步分析。

表1 目标物质的洗脱梯度Table 1 Elution gradient of target substance

表2 三唑酮暴露对大型溞世代间隔时间的影响Table 2 Influences of triadimefon exposure on generation interval time of Daphnia magna

基于上述表达量分析，F1代对照组(空白对照和溶剂对照)与处理组(试验设置的6个浓度组)的差异表达基因共376个，其中上调基因322个，下调基因54个(图4(a))。F2代对照组与处理组的差异表达基因共422个，其中上调基因87个，下调基因335个(图4(b))。而F1和F2代的差异表达基因共2 604个，其中上调基因1 101个，下调基因1 503个(图4(c))，由火山图可知，F1代上调基因明显高于下调基因，而F2代则相反。相比之下，F1代与F2代比较的差异表达基因上调与下调基因数基本相似。F1代对照组(空白对照和溶剂对照)与所有处理组(试验设置的6个浓度组)共有差异表达基因36个，F2代对照组与所有处理组共有差异表达基因82个，将F1代的基因与F2代的基因对比得到两者共有差异表达基因共336个(如图4(d)～(f)所示)。

图3 样本间相关性分析和PCA分析注：F1_1～F1_8分别代表F1代空白对照组、丙酮对照组以及5、12.5、25、50、100、200 μg·L-1三唑酮处理组；F2_1～F2_8分别代表F2代空白对照组、丙酮对照组以及5、12.5、25、50、100、200 μg·L-1三唑酮处理组；(d)和(e)中红色虚线圈出是200 μg·L-1处理组，(f)中红色虚线圈出是F1代，蓝色虚线圈出是F2代。Fig. 3 Correlation analysis and PCA analysis among samplesNote: F1_1～F1_8 represent control group, acetone control group, and 5, 12.5, 25, 50, 100, 200 μg·L-1 triadimefon treatment group, respectively; F2_1～F2_8 represent control group, acetone control group, and 5, 12.5, 25, 50, 100, 200 μg·L-1 μg·L-1 triadimefon treatment group; circled by a red dotted line in (d) and (e) is treatment group of 200 μg·L-1; circled by a red dotted line and a blue dotted line in (f) represent F1 and F2, respectively.

图4 差异表达基因火山图(a)、(b)、(c)和Venn图(d)、(e)、(f)注：F1_0和F2_0分别代表F1代和F2代对照组，F1_3～F1_8分别代表F1代5、12.5、25、50、100、200 μg·L-1三唑酮处理组；F2_3～F2_8分别代表F2代5、12.5、25、50、100、200 μg·L-1三唑酮处理组。Fig. 4 Volcano figures of differentially expressed genes (a), (b), (c) and Venn figures of differentially expressed genes (d), (e), (f)Note: F1_0 and F2_0 represent control group in F1 and F2, respectively; F1_3～F1_8 represent 5, 12.5, 25, 50, 100, 200 μg·L-1 triadimefon treatment group in F1, respectively; F2_3～F2_8 represent 5, 12.5, 25, 50, 100, 200 μg·L-1 triadimefon treatment group in F2.

2.3 关键通路和关键基因

通过上述基因表达差异的分析，对F1代和F2代的所有差异表达基因进行功能聚类和通路分析(表3和表4)。结果表明，F1代差异表达基因显著富集的主要通路包括甾类激素生物合成、视黄醇代谢、酮体的合成与降解等代谢通路、氧化应激、内质网应激以及与寿命调节和消化系统相关的通路。其中寿命调节和蛋白质消化吸收是富集最为显著的通路。F2代显著富集通路主要有酮体的合成与降解、类固醇生物合成和谷胱甘肽代谢等代谢通路、抗原的处理和呈递以及造血细胞谱系等免疫系统通路、内分泌系统和氮代谢。其中氮代谢和造血细胞谱系是富集最显著的通路。对富集到关键通路的显著差异表达基因在不同浓度水平下进行表达量分析(图5)。F1代中富集到关键通路的绝大部分关键基因在浓度100 μg·L-1和200 μg·L-1时显著上调(图5)，仅基因trypsin alpha-1-like和基因zinc metalloproteinase nas-12-like在浓度200 μg·L-1时显著下调，且其他与蛋白质消化吸收相关的基因在该浓度下与对照组之间无显著差异，而在较低浓度时上调。这说明在较低浓度下，机体通过上调蛋白质消化吸收相关基因应对外界胁迫，而200 μg·L-1可能已经超过其基因水平的毒性阈值，从而引起不利的基因表达。而与药物代谢相关的CYP3A (CYP3A是一种重要的CYP450酶系，可催化许多外源药物的代谢)和UGT(葡萄糖醛酸转移酶)相关基因的表达在不同浓度下呈现不同的表达模式，这些基因富集在甾类激素生物合成通路，这可能导致内分泌干扰作用[29-30]。3-羟基丁酸脱氢酶、全反式视黄醇13,14还原酶相关基因在各处理浓度下均上调。不同于F1代，在F2代中，大部分关键基因显著下调。LIPA(脂肪酶)、CD13(氨肽酶类)、CD63(四磷酸腺苷)相关基因在高浓度处理组(50、100和200 μg·L-1)显著下调，CA(碳酸酐酶)、CTSL(组织蛋白酶L)和BDH(3-羟基丁酸脱氢酶)相关基因在各处理组显著下调。F1代和F2代基因表达的差异说明F1代的暴露会影响其子代的基因表达，加剧对子代的毒性作用。

图5 关键差异表达基因表达量(log(TPM+1))分布图注：F1_0和F2_0分别代表F1代和F2代对照组，F1_3～F1_8分别代表F1代5、12.5、25、50、100、200 μg·L-1三唑酮处理组；F2_3～F2_8分别代表F2代5、12.5、25、50、100、200 μg·L-1三唑酮处理组。Fig. 5 Expression distribution (log(TPM+1)) of key differentially expressed genesNote: F1_0 and F2_0 represent control group in F1 and F2, respectively; F1_3～F1_8 represent 5, 12.5, 25, 50, 100, 200 μg·L-1 triadimefon treatment group in F1, respectively; F2_3～F2_8 represent 5, 12.5, 25, 50, 100, 200 μg·L-1 triadimefon treatment group in F2.

表3 F1代差异表达基因KEGG富集主要信号通路Table 3 Main pathways of KEGG enrichment of differentially expressed genes in F1

表4 F2代差异基因KEGG富集主要信号通路Table 4 Main pathways of KEGG enrichment of differentially expressed genes in F2

3 讨论(Discussion)

本研究依据OECD 211[27]测试方法开展大型溞多代实验，评估三唑酮对大型溞的毒性效应。结果表明：(1)三唑酮浓度为200 μg·L-1时，显著影响大型溞体长。刘娜等[22]探究了三唑酮对大型溞不同测试终点的毒性效应，发现当三唑酮浓度≤100 μg·L-1时，大型溞体长与对照组无明显差别，计算以体长为生长指标的最大无影响浓度(no observed effect concentration, NOEC)为100 μg·L-1；(2)当三唑酮浓度≥100 μg·L-1时，大型溞繁殖能力显著降低，说明该浓度可能已经达到造成繁殖毒性的阈值；(3)三唑酮的暴露可能导致世代间隔时间延长。胡方华等[31]关于三唑酮对大型溞的慢性毒性效应的研究表明，三唑酮对第2代大型溞染毒的影响比对第1代的影响更大。这说明三唑酮对大型溞存在长期毒性效应，有必要考虑三唑酮对大型溞的传代效应。

对比F1代和F2代的KEGG通路和关键基因，KEGG通路均与代谢和免疫相关，但F2代内分泌干扰作用和免疫反应更加明显，F1代更多表现在氧化应激和消化系统相关的通路。F1代上调基因占比高而F2代下调基因占比高，在F1代中显著上调的BDH相关基因在F2代中显著下调。这些差异可能是由于F1代的暴露降低了三唑酮对F2代的毒性阈值，从而加剧了F2代的毒性作用，长期暴露后，氧化损伤最终导致细胞凋亡和机体功能损失等。F1代在低浓度暴露时，生物体通过提高基因表达来适应污染环境，而高浓度时可能超过了部分基因表达毒性阈值从而引起基因表达量的变化。F2代在低浓度暴露时基因下调，可能是因为其母代F1暴露于三唑酮从而加剧了对子代的毒性。因此，随着时间和代际延长，三唑酮对大型溞的毒性作用可能会加强，大型溞对某一阈值浓度的三唑酮有一定的调控能力，但超过一定浓度，可能会存在传代效应。

综上所述，低浓度三唑酮长期暴露大型溞会引起慢性毒性，≥200 μg·L-1的三唑酮显著影响大型溞的体长和繁殖。三唑酮可能会引起大型溞的氧化应激反应、神经毒性、免疫毒性和内分泌干扰，存在传代效应，长期暴露会下调免疫反应等相关基因的表达，降低毒性阈值。

通讯作者简介：金小伟(1985—)，男，博士，正高级工程师，主要研究方向为生态毒理学和生态风险评价。