加权基因共表达网络分析全氟和多氟烷基化合物对人间充质干细胞的毒性靶点

潘一帆，秦会，刘薇

大连理工大学环境学院 工业生态与环境工程重点实验室，大连 116024

全氟辛烷磺酸(perfluorooctane sulfonate, PFOS)的生产和使用逐渐被禁止后，新型替代品用量显著升高，目前已有超过3 000种全氟多氟化合物(per- and polyfluoroalkyl substances, PFASs)被投入全球市场[1-2]。新型替代品的环境效应和健康风险未知，给PFASs的监管带来了极大的挑战。氯代多氟醚基磺酸(chlorinated polyfluoroalkyl ether sulfonates, Cl-PFESAs)作为铬雾抑制剂在中国使用长达30余年，因其与PFOS相似的化学结构和理化性质，被视作PFOS的替代品[3]。初步研究发现，Cl-PFESAs具有肝脏毒性和发育毒性等有害效应，有必要进行深入的毒理学研究，揭示其敏感作用靶点[4-6]。传统的动物实验无法针对大量的化学品开展高通量研究，且与人体健康效应存在种间差异。因此，为揭示大量PFASs的毒性效应与敏感靶点，亟需敏感可靠的高通量毒理学模型。

人源干细胞由于具有增殖分化潜能，且其分化过程良好地模拟了人体的生长发育，是一类研究化学品潜在发育毒性的适宜模型[7-8]。基于人体解剖和动物实验的研究表明，PFASs可于骨骼中蓄积，这提示骨骼可能是PFASs的作用靶组织[9]。流行病学研究指出，暴露于PFASs与骨密度降低和骨关节炎的发生有一定的关联[10-11]。因此，来源于人骨髓的间充质干细胞(human bone marrow mesenchymal stem cells, hBMSCs)可能是PFASs的作用靶细胞。笔者在前期研究中发现，PFOS暴露干扰了hBMSCs的成骨和成脂分化[12]。这些线索提示，hBMSCs是PFASs毒性靶点研究的一种灵敏、准确的毒理学模型。

转录组学分析是分子毒理学研究的重要手段，通常按照表达规律对基因进行富集归类，然后对已分类基因开展功能或调控分析。加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)基于多基因间的表达相关性，不仅可鉴定高度协同变化的基因模块，同时也可根据基因模块与特定性状或表型之间的关联鉴定生物标记基因和敏感靶点，广泛应用于疾病发生与基因关联分析，是筛选疾病治疗靶点的出色工具[13]。本研究选择hBMSCs体外模型，基于其基因表达谱，结合WGCNA方法探究Cl-PFESAs、PFOS、全氟己烷磺酸(perfluorohexane sulfonate, PFHxS)和全氟辛酸(perfluorooctanoic acid, PFOA)这4种PFASs暴露与特定基因模块间的关联，分析PFASs干扰的生物学过程和敏感靶基因。对PFOS和PFOA这2种全氟化合物已有较为深入的毒理学研究，而PFHxS和Cl-PFESAs则被选作2种代表性替代品。本研究将其与PFOS/PFOA暴露下细胞基因组学特征相比较，有利于预测和评估替代品的潜在毒性效应及分子靶点，以期为该类化学品的毒理学研究、安全性评价与监管提供依据。

1 材料与方法(Materials and methods)

1.1 试剂

工业级Cl-PFESAs产品购自上海中科合臣股份有限公司，该产品为8碳和10碳Cl-PFESA的混合物，商品名为F-53B(Cl(CF2)6O(CF2)2SO3K, CAS 73606-19-6)，纯度大于98%。纯度≥98% PFOS (C8F17KO3S)、纯度为95% PFOA (C8F15KO2)和纯度≥98% PFHxS (C6F13SO3K)购自Sigma-Aldrich。二甲基亚砜(DMSO)，纯度大于99.5%，购自Panreac-Applichem(德国)。PFASs暴露液配制过程中均使用终浓度为0.1%的DMSO作为助溶剂。

1.2 细胞培养及处理

由广州赛业生物提供的人骨髓间充质干细胞，取自健康供体，细胞常规检测结果表明，细菌、真菌和支原体显阴性，内毒素小于10 EU。细胞鉴定检测结果表明，细胞复苏贴壁率大于95.1%，可分化为成骨、成脂和成软骨细胞。细胞表面标记分子检测结果显示，CD105、CD29、CD73和CD44呈阳性，CD34、CD45和CD11b呈阴性。

细胞置于37 ℃、5% CO2的培养箱中，利用hBMSCs基础培养基(广州赛业生物)进行培养及传代。取第6代细胞接种于细胞培养板中，1 d后将基础培养基换为含100 nmol·L-1的Cl-PFESAs、PFOS、PFHxS或PFOA的暴露培养基，以及含0.1% DMSO的溶剂对照培养基，持续暴露7 d，每2～3天更换培养基。

1.3 细胞RNA抽提、纯化和芯片检测

暴露7 d后，利用TRIzol(Life technologies，美国)抽提细胞总RNA。利用核酸定量仪测定RNA纯度和浓度，并用琼脂糖凝胶电泳判断质量，然后进行芯片检测。芯片杂交和扫描由上海伯豪生物公司完成。芯片杂交按照Affymetrix表达谱芯片(美国)的标准流程进行，使用GeneChip®Hybridization和Wash and Stain Kit试剂盒(Affymetrix)，在滚动杂交炉(Hybridization Oven 645, Affymetrix)中45 ℃滚动杂交16 h，杂交完成后在洗涤工作站(Fluidics Station 450, Affymetrix)进行芯片的洗涤。扫描芯片结果(GeneChip®Scanner 3000, Affymetrix)，用Command Console Software 4.0 (Affymetrix)读取原始数据。原始数据已上传于GEO数据库(GSE115836)。芯片测试共分为5组(DMSO溶剂对照组、Cl-PFESA暴露组、PFOS暴露组、PFHxS暴露组和PFOA暴露组)，每组3个平行。

1.4 加权基因网络共表达分析

利用R语言软件包WGCNA，按照默认参数构建加权基因共表达网络，筛选代表性PFASs共同影响的生物学功能和敏感靶基因[13]。首先对芯片中总计49 293个探针进行初步筛选，排除多个探针对应同一基因的情况，共获得20 014个基因。为筛选差异表达倍数差异较高的基因，计算了15个细胞样本中20 014个基因各自差异表达倍数的方差，并选取方差较大的前25%的基因，共得到了5 004个基因用于后续分析[14-15]。对样本进行皮尔逊相关性分析，筛查是否存在离群样本。确保不存在离群样本后，在基因间进行皮尔逊相关性分析，按照无尺度网络标准构建邻接矩阵。为降噪和更好地衡量基因网络连通性，将邻接矩阵转化为拓扑重叠矩阵。然后利用动态剪切枝法划分基因模块，并将相似分支进行合并。

1.5 基因模块与PFASs暴露相关性分析

为探究哪些基因模块与PFASs暴露存在关联，计算模块特征向量(module eigengene, ME)，并利用ME计算基因模块与溶剂对照组以及PFASs暴露组细胞基因表达谱间的相关系数(correlation coefficient,r)，将P≤0.05定义为显著相关。

1.6 基因模块功能注释与枢纽基因筛选

提取基因模块中的基因，利用Database for Annotation, Visualization, and Integrated Discovery(DAVID)进行功能注释[16]，并利用cytoscape对WGCNA构建的包含PFASs潜在靶点的基因共表达网络进行可视化，依据连通性筛选网络中的枢纽基因。

2 结果(Results)

2.1 PFASs作用下hBMSCs的基因共表达网络构建

利用15个细胞样本方差前25%的5 004个基因构建共表达网络。首先用flashclust函数执行层次聚类分析，无明显离群值。按无尺度网络标准建立临近矩阵，为使网络节点连接度的对数与该节点出现概率的对数相关系数≥0.9(图1(a))，计算确认临近函数加权参数，可得软阈值(β)为16，进行网络构建，得到14个模块(图1(b))，此外有938个基因未能分配到任何一个模块(灰色)。随后对14个模块进行层次聚类分析(图1(c))，判别不同基因模块间的关联。

图1 软阈值确定(a)、基因聚类树和模块颜色划分(b)和模块聚类分析(c)Fig. 1 Identification of soft threshold (a), clustering dendrograms of genes together with assigned module colors (b) and module cluster analysis (c)

2.2 基因模块与PFASs暴露的相关性分析

对基因模块和PFASs暴露组样本进行相关性分析，溶剂对照组与Brown(r=0.63,P=0.01)、Greenyellow(r=0.66,P=0.001)和Blue模块(r=0.66,P=0.007)显著正相关(图2)。Cl-PFESA暴露组与Greenyellow模块显著负相关(r=-0.61,P=0.02)，PFOS暴露组与Blue模块显著负相关(r=-0.62,P=0.01)。PFOA和PFHxS暴露组与Greenyellow和Blue模块无显著相关性。这表明，PFASs暴露组中Brown、Greenyellow和Blue模块基因表达模式与对照组有较大差异，提示相关模块中基因调控的生物学过程可能受到4种PFASs的影响，且模块中可能包含4种PFASs共同作用的靶基因。此外，聚类分析结果发现，black模块和red模块相关度较高(图1(b))，并与PFOS暴露显著正相关(图2)，提示其调控的生物学过程可能受到PFOS的干扰。

图2 基因模块与全氟多氟化合物(PFASs)暴露细胞样本间的相关系数注：Cl-PFESAs表示氯代多氟醚基磺酸，PFHxS表示全氟己烷磺酸，PFOS表示全氟辛烷磺酸，PFOA表示全氟辛酸；* P≤0.05，显著相关。Fig. 2 Correlation coefficients between modules and per- and polyfluoroalkyl substances (PFASs) exposed samplesNote: Cl-PFESAs, chlorinated polyfluoroalkyl ether sulfonates; PFHxS, perfluorohexane sulfonate; PFOS, perfluorooctane sulfonate; PFOA, perfluorooctanoic acid; *significance of correlation at P≤0.05 level.

2.3 PFASs影响基因模块的生物学功能分析

为了解每个模块所表征的生物学功能，笔者提取了每个模块中的基因，并利用DAVID进行生物学过程富集分析。每个模块中最显著的3个生物学过程如图3所示。Greenyellow模块中显著富集“骨髓白细胞分化负调控”以及“嘧啶和嘌呤核苷酸的代谢过程”，提示PFASs可能干扰核酸代谢以及免疫细胞分化。与PFOS暴露显著相关的Red和Black模块中分别富集出与糖酵解和线粒体氧化磷酸化相关的生物学过程，提示PFOS可能影响细胞呼吸过程。Blue模块中富集出蛋白质合成的tRNA氨酰化以及胆固醇和异戊二烯的生物合成，提示PFASs可能干扰蛋白质合成以及脂质代谢等。此外，在Blue模块中显著富集出“成骨分化”，提示PFASs对hBMSCs成骨分化的干扰作用。Brown模块中显著富集“骨化”，并且与PFOS显著相关的Tan模块中显著富集出“成骨分化正调控过程”，进一步提示PFASs对骨形成的干扰作用。

2.4 识别PFASs暴露干扰的模块枢纽基因

由于Blue模块中显著富集出与脂质代谢相关的生物学过程，而干扰脂质代谢是PFASs的主要毒性效应之一[17]，因此，进一步探究PFASs影响的Blue模块的基因共表达，以及筛选枢纽基因。利用WGCNA输出的模块网络文件结合Cytoscape的Blue模块基因共表达网络图(基因节点间权重>0.1)，并根据连通性筛选枢纽基因(图4)。在Blue模块中，连通性最高的4个枢纽基因分别为Delta(24)-sterol reductase (DHCR24)、Squalene monooxygenase (SQLE)、Laminin subunit alpha-1 (LAMA1)、3-beta-hydroxysteroid-delta(8), delta(7)-isomerase (EBP)。这4个枢纽基因均显著下调(P≤0.5)(表1)。

图4 Blue模块基因共表达网络图注：SQLE表示枢纽基因鲨烯环氧酶，LAMA1表示枢纽基因层粘连蛋白α1亚基，DHCR24表示枢纽基因24-脱氢胆固醇还原酶，EBP表示枢纽基因甾醇δ异构酶；基因节点大小、颜色与基因节点连通性呈正相关。Fig. 4 Gene co-expression network of Blue moduleNote: SQLE, squalene monooxygenase; LAMA1, hub gene laminin subunit alpha-1; DHCR24, hub gene delta(24)-sterol reductase; EBP, hub gene 3-beta-hydroxysteroid-delta(8), delta(7)-isomerase; size and color of gene node positively correlate with connectivity of gene node.

3 讨论(Discussion)

利用WGCNA分析了4种PFASs暴露所影响的hBMSCs的关键生物学过程，并通过筛选相应生物学过程中的枢纽基因识别靶标基因。

基于PFASs暴露后hBMSCs的基因表达谱信息构建的14个基因模块中，Brown、Blue和Greenyellow模块中的基因表达模式在PFASs暴露组细胞中与对照组间存在较大差异。Greenyellow模块中富集的最显著的生物学过程是骨髓白细胞分化的负调控。Brieger等[18]研究了PFOS和PFOA体外暴露对人外周白细胞的影响，发现PFOA显著增强了骨化三醇诱导的HL-60细胞株的单核细胞分化，并增加了巨噬细胞IL-6和TNF-α的释放，而在脂多糖刺激后PFOS暴露降低了促炎细胞因子TNF-α的释放。此外，Corsini等[19]发现全氟丁烷磺酸(perfluorobutane sulfonic acid, PFBS)、全氟癸酸(perfluorodecanoic acid, PFDA)、PFOS和氟调聚醇(fluorotelomer, 8:2 Telomer)均抑制了脂多糖诱导的TNF-α生成。动物实验和流行病学研究均已发现PFOS等全氟烷酸化合物暴露与免疫功能异常之间的关联[17]，但目前尚无关于Cl-PFESA等新型替代品的免疫毒理学研究。本研究通过PFASs对转录组的影响和WGCNA分析，证明免疫功能是PFASs毒性效应的敏感靶标。骨髓是人体重要的免疫器官，hBMSCs在调节骨髓微环境及其免疫功能中起到关键作用[20]。因此，以hBMSCs为毒理学研究模型，对于揭示免疫干扰作用和机制具有重要意义。

Blue模块中最显著富集的生物学过程为胆固醇生物合成。此外，其基因共表达网络中连通性最高的枢纽基因DHCR24、SQLE和EBP均与脂质代谢相关。DHCR24可催化甾醇中间体δ-24双键的还原，此外还可通过降低氧化应激诱导的细胞凋亡过程中caspase 3的活性来保护细胞免受氧化损伤[21]。只要上游胆固醇合成过程未被完全抑制，DHCR24下调很可能导致胆固醇累积[22]。SQLE可催化甾醇生物合成中的第一个氧化步骤，是该途径的限速酶之一，研究表明SQLE受到影响可干扰胆固醇稳态[23]。EBP可催化δ(8)-甾醇转化为其相应的δ(7)-异构体[24]。PFASs暴露导致以上几个基因表达下调，为进一步验证DHCR24、SQLE和EBP为PFASs作用的潜在敏感靶基因，本研究利用比较毒理组学数据库对PFOS和PFOA的互作基因进行了检索[25]。结果表明，PFOA的2 731个互作基因和PFOS的2 001个互作基因中包含本研究富集出的胆固醇的生物合成过程14个基因中的13个，并且在PFOS和PFOA互作基因中都包含EBP和SQLE，而DHCR24包含于PFOA互作基因中。此外，来源于大鼠[26]、小鼠[27]和人肝癌细胞[28]等体内体外模型的研究表明，PFOS或PFOA可干扰DHCR24、SQLE和EBP基因或蛋白的表达，变化趋势与暴露模式和受试生物的种类有关。比较毒理组学数据和实验结果均支持本研究利用WGCNA分析hBMSCs基因组数据识别出的枢纽基因。

表1 PFASs暴露对Blue模块中枢纽基因表达水平的影响Table 1 Effects of PFASs exposure on hub genes expression from Blue module

PFASs干扰啮齿类动物脂质代谢过程的主要机制是激活PPARα这一脂质代谢的关键调控因子[17,29]。在脂肪干/祖细胞中，PFASs主要通过激活PPARγ促进细胞成脂分化[5,30]。hBMSCs是脂肪细胞和成骨细胞共同的祖细胞，其成骨成脂分化平衡异常与代谢疾病和发育相关疾病密切相关[31]。与PFOS显著相关的Tan模块中“成骨分化正调控过程”显著富集(r=0.51,P=0.05)，并且Blue和Brown模块中分别显著富集“成骨分化”和“骨化”，提示PFOS可干扰成骨/成脂分化平衡。前期笔者在基于细胞分化的分子标志物和表型标志物的研究中，发现PFOS在人体暴露相关浓度下即可促进hBMSCs成脂分化并抑制其成骨分化，部分原因可能是由于PFOS对Wnt信号通路的干扰，但可能仍存在其他机制，需要更深入研究[12]。而本研究结果提示，若继续探究PFASs干扰成骨或成脂分化的机制，可围绕本研究识别的相关枢纽基因展开。基于hBMSCs基因表达谱的WGCNA分析，有助于识别PFASs毒性作用敏感靶标，揭示PFASs干扰脂质代谢的分子机理。

综上，本研究基于hBMSCs的基因表达谱，利用WGCNA方法对4种典型PFASs毒性作用干扰的敏感生物学过程和靶基因进行了筛选，发现免疫调控和脂质代谢是PFASs的敏感靶标，DHCR24、SQLE和EBP可能是PFASs干扰脂质代谢的敏感靶基因。可进一步利用hBMSCs体外模型研究PFASs的相关毒性作用和机制，有助于建立该类化学物毒性预测和筛选的生物标志物，为其安全性评价和监管提供科学依据和新方法。