摘要
为探索内分泌干扰物DEHP在母鼠怀孕期间暴露后对子鼠生殖系统的影响，为DEHP诱导雄性生殖毒性机制提供生物信息学参考。采用美国国立生物技术信息中心的GEO(Gene expression omnibus)数据库筛选DEHP组与对照组差异表达的mRNA并进行GO和 KEGG通路富集分析。使用STRING数据库、Cytoscape软件对差异表达基因构建蛋白质互作网络(PPI),并根据CytoHubba插件的MCC算法筛选出前十个关键基因。利用GSEA软件对基因集进行基因集富集分析(Gene set enrichment analysis，GSEA)，筛选在DEHP暴露后有统计学意义的基因集。本研究首先基于差异表达筛选出Cyp51，Sqle，Sc5d，Fdps，Dhcr7，Acss2，Insig1，Idi1，Scd1等关键基因。富集分析显示这些mRNA可能通过影响胆固醇合成代谢、类固醇合成代谢等相关途径发挥作用。此外，GSEA分析结果也清楚地表明DEHP暴露后会诱导小鼠胆固醇合成、mTORC1通路、雄性激素响应相关途径发生改变。研究发现，在母鼠怀孕期间暴露DEHP可能会影响雄性子鼠mTORC1通路，胆固醇合成相关途径引起小鼠雄性激素紊乱导致精子发生障碍，干扰小鼠生殖系统。
关键词
Abstract
To explore the effects of endocrine disruptor DEHP on the reproductive system of offspring rats after exposure during pregnancy, and provide bioinformatics references for the mechanism of DEHP-induced male reproductive toxicity.The GEO (Gene expression omnibus) database of the National Center for Biotechnology Information is used to screen differentially expressed mRNAs between the DEHP group and the control group, and GO and KEGG pathway enrichment analysis are performed. The differentially expressed genes are mapped into a PPI protein interaction network using the STRING database and Cytoscape software, and the top ten key genes are selected based on the MCC algorithm of the CytoHubba plugin. The gene set enrichment analysis (GSEA) is performed using GSEA software to screen gene sets with statistically significant differences after DEHP exposure.This study first screens out key genes such as Cyp51,Sqle,Sc5d,Fdps,Dhcr7,Acss2,Insig1,Idi1,Scd1 based on differential expression. Enrichment analysis shows that these mRNAs may play a role by affecting cholesterol synthesis and metabolism, steroid synthesis and metabolism, and other related pathways. In addition, the results of GSEA analysis also clearly shows that DEHP exposure induces changes in cholesterol synthesis, mTORC1 pathway, and male hormone response-related pathways in mice.This study finds through bioinformatics that exposure to DEHP during pregnancy in female mice may affects the mTORC1 pathway and cholesterol synthesis-related pathways in male offspring, causing male hormone disorders in mice, resulting in impaired spermatogenesis and interfering with the reproductive system of mice.
Keywords
邻苯二甲酸二（2-乙基己基）酯（DEHP）是邻苯二甲酸盐的主要代表，是一种环境内分泌干扰物，常作为聚氯乙烯商品中的增塑剂使用，如建筑材料、涂料、儿童玩具和医疗器械等[1]。由于邻苯二价酸酯在添加使用时是与原料混合而非化学结合的，所以在使用过程中会向环境，空气，食物中扩散[2]。人们可通过食品、皮肤、药械等途径接触到这种物质[3]。邻苯二甲酸二乙基己酯因其内分泌活性和致相关代谢性疾病而引起人们的广泛关注，目前研究表明邻苯二甲酸酯暴露可对雄性生殖系统造成损伤、具有雄性生殖毒性，而孕妇与发育期儿童对邻苯二甲酸酯的毒性机制更为敏感[4]。有大鼠实验研究表明雌性大鼠在怀孕期间暴露邻苯二甲酸盐可致雄性后代会出现生殖系统异常，如输精管、前列腺畸形和尿道下裂等[5]。与这些影响密切相关的是，在雄激素正常控制下生殖道发育的关键窗口期，邻苯二甲酸酯暴露后胎儿睾丸睾酮的产生显著减少[6]。邻苯二甲酸酯暴露后睾酮产生减少的原因可能是邻苯二甲酸酯导致胆固醇运输障碍，其相关机制不清楚。因此为了进一步研究DEHP暴露后致胎儿生殖系统毒性损伤的作用机制，本次采用生物信息学的方法，从基因表达综合数据库（Gene expression omnibus，GEO）寻找DEHP染毒相关的表达谱，鉴定差异表达基因并进行功能和通路富集分析，探寻DEHP暴露后影响胎儿生殖系统的关键基因及其相关机制。
1 材料和方法
1.1 基因芯片数据下载
在GEO数据库（https://www.ncbi.nlm.nih.gov/geo）中以“DEHP”和“Testis”为关键词，检索数据集，从中筛选获得了GSE41901数据集。该研究是基于GPL7294平台的芯片，其在原始实验设计阶段皆以500 mg/（kg·d）的剂量口服灌胃于妊娠母鼠，小鼠胎儿出生后立即麻醉处死并分离胎儿睾丸，提取RNA制备基因组芯片样本。经筛选，GSE41901数据集包含4例正常对照样本与4例DEHP处理样本，本文选取此样本数据进行研究分析。
1.2 差异表达基因（Differentially expressed genes，DEGs）的筛选
使用R软件（ Version 4.3.0）对GSE41901矩阵数据进行归一化处理，首先使用“corrplot”R包和OmicShare Tools（https://www.omicshare.com/tools）分别对基因矩阵进行相关性分析和PCA主成分分析，以评估样本之间的重复性与相关性。接下来使用“Limma” R包，对基因矩阵进行差异表达基因筛选，|log2FC|＞0.5和P＜ 0.05作为DEGs的筛选条件[7]。
1.3 DEGs的功能富集分析
使用R 软件（Version 4.3.0）的“clusterProfiler”R包对DEGs进行基因本体（Gene Ontology，GO）功能富集分析；同时使用“org.Mm.eg.db”R包进行京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）信号通路分析，以P＜0.05为分界点。
1.4 蛋白质互作用网络分析
将获得的DEGs上传至STRING 11.5版本数据库（https:// STRING-db.org）构建蛋白质互作用网络（PPI），物种限定为“Mus musculus”，以互作评分＞0.4作为阈值条件。使用Cytoscape（Version 3.9.1）可视化网络，运用Cytohubba插件的MCC算法对网络图中的相互作用基因进行筛选，选取得分最高的前10个基因视为核心基因（Hub 基因）。
1.5 基因集富集分析GSEA
使用GSEA软件（Version 4.3.2）和MSigDB数据库（www.gsea-msigdb.orgs）进行基因集富集分析[8]，输入归一化处理后的基因表达矩阵后对基因进行排序，以确定基因在各基因集中的状态，判断该基因集在实验处理后是否被抑制或激活。使用MSigDB数据库中的hallmark基因集用作富集输入，并将富集结果中标准化富集分数|NES|＞1，Pvalue＜0.05，FDR＜0.25的基因集作为有统计学意义的基因集。
1.6 Hub基因的表达验证分析
使用DISCO数据库（https://www.immunesinglecell.org/）[9]验证筛选出的Hub基因在睾丸组织中的细胞成分分布，从而判断雄性生殖系统遭到破坏的具体机制。
2 结果
2.1 DEGs的分析与筛选
使用R软件（Version 4.3.0）中的“Limma”包对GSE41901选择的数据集进行分组处理及DEGs分析。首先为了评估两组样本之间的差异性，进行了pearson相关系数分析，如图1（a）所示，样本间最小相关系数为0.94，最大相关系数为0.99，表明我们得到的数据具有很好的相关性和保真性，其次可以看到在正常组与正常组间的相关性系数，要显著高于正常组与DEHP组间的相关性，说明对照组与DEHP组在实验后基因表达可能已经产生了差异。除此之外，主成分分析（PCA）也能更好的体现出两组样本间的差异程度，PCA作为一种简化数据集的技术，可通过线性变换将原始数据变换为一组各维度线性无关的表示，利用降维思想将原始数据从一堆多项指标转化为少数几个综合指标，即为主成分。主成分分析经常用于减少数据集的维数，同时保持数据集对方差贡献最大的特征，本次研究利用PCA对基因组所有样本的表达数据进行降维分析，在二维坐标平面内得出样本降维后的分布和离散情况，验证本数据集的分组合理性，如图1（b）所示，X轴是第一主成分贡献值PC1，Y轴是第二主成分贡献值PC2，两组样本中有39.9%的差异来自PC1，两组样本均被PC1很好的隔离开，表明在DEHP暴露后，正常组与DEHP组小鼠的转录组表达可能发生了差异。差异分析处理后，共有86个mRNA被鉴定为上调或下调，其中上调基因14个，下调基因72个（|log2FC| ＞ 0.5，P＜0.05）。将差异基因在正常组与DEHP组中以热图形式展示（图1（c）），这些差异基因在正常组与DEHP组之间呈现了很好的聚类。如图1（d）所示，将差异基因展示出来，其中左上侧区域表示下调基因，右上侧区域表示上调基因，中间区域表示其它基因。
2.2 DEGs的功能富集分析
为探究在DEHP暴露后是否有对小鼠的生理功能产生影响，对得到的DEGs进行了GO及KEGG富集分析（表1）。图2（a）为丰富度靠前的KEGG通路富集结果与各基因之间的流动，包括了皮质醇的合成和分泌、皮质醇增多症、类固醇生物合成、卵巢类固醇生成、醛缩酮的合成和分泌、萜类骨架生物合成。差异基因富集到丰富度靠前的GO术语如图2（b）所示，揭示DEGs主要参与了类固醇代谢过程、类固醇合成过程、胆固醇代谢过程、皮质醇合成与代谢等相关生物调节过程，这些都通常被认为与血清睾酮水平紊乱密切相关[10]。
图1差异分析结果
Fig.1Difference analysis results
注：（a）样本相关性分析;（b）PCA主成分分析;（c）差异表达热图;（d）差异表达火山图.
表1GO与KEGG富集分析
Table1GO and KEGG enrichment analysis
图2KEGG通路与GO富集分析通路图
Fig.2KEGG pathway and GO enrichment analysis pathway diagram
注：（a）KEGG富集分析（P＜ 0.05）;（b）GO富集分析（P＜ 0.05）.
2.3 蛋白质相互作用（PPI）网络构建与核心基因的筛选鉴定
使用STRING数据库对筛选的DEGs构建PPI网络，将网络节点信息导入到Cytoscape软件进行网络图的可视化展示，并运用CytoHubba插件中的MCC算法对网络图中的相互作用基因进行评分筛选，选取了前10个得分最高的基因作为Hub基因。如图3（a）所示，得到的Hub基因依次为Hmgcs1，Cyp51，Sqle，Sc5dl，Fdps，Dhcr7，Acss2，Insig1，Idi1，Scd1。采用Pearson相关性分析，对10个Hub 基因的表达水平进行相关性分析。由于筛选出的10个Hub基因均为表达组中的下调基因，所以10个Hub基因的表达均呈正向相关趋势，且相关性高，表明10个基因之间联系密切（图3（b））。如图3（c）所示，与正常组相比，10个Hub基因在DEHP组表达水平均呈显著下调趋势。在10个Hub基因中，根据其相关性程度与之前功能富集通路和基因之间的联络关系，可以选择Hmgcs1，Cyp51，Dhcr7，Insig1，Idi1作为之后的重点研究对象。
图3关键mRNA筛选分析图
Fig.3Key mRNA screening analysis chart
注：（a）10个Hub基因的蛋白质网络;（b）10个Hub基因之间的相关性分析；（c）10个Hub基因在正常组与DEHP组的表达水平分析.
2.4 基因集富集分析
对正常组与DEHP组差异分析后进行了基因富集分析，以研究一些关键编码通路基因集在本研究的基因表达数据中处于什么样的位置，以便于观察通路在实验前后的激活或抑制情况。本次GSEA所采用的hallmark基因集由定义生物状态和进程的标志基因组成，可传达特定的生物状态或过程并显示连贯的表达，由此来反映本次数据样本的核心基因参与了哪些生物学功能。本次根据筛选阈值选出了3个基因富集显著的基因集，图4（a）、图4（b）、图4（c）分别为胆固醇恒定（NES=2.66，FDR=0.00）、MTORC1通路（NES=1.57，FDR=0.05）和雄性激素响应（NES=1.45，FDR=0.085），此3类基因集皆与雄性生物的精子合成密切相关。图中折线图的峰值若出现在前端，则峰值前的基因为核心基因；若峰值出现在后端，则峰值后的基因为核心基因。热图红色部分对应的基因在正常组中高表达，蓝色部分对应的基因在实验组中高表达。图中可以看到3组基因集的峰值皆在前端，故峰值前的基因为核心基因，核心基因皆在正常组中呈高表达，表明三组基因集在实验后表达有所下调。图4（d）为在3组基因集中选取共同基因的维恩网络图，在3组基因集中选出了Idi1，Hmgcs1，Hmgcr三个共同基因。其中Idi1和Hmgcs1也在我们前面筛选出的Hub基因中。
2.5 关键mRNA的细胞成分表达
DISCO数据库包含来自158种疾病的107种组织的11 265个样本中超过4 800万个细胞测序数据，本研究使用DISCO数据库分析了筛选出的关键枢纽基因如Cyp51，Dhcr7，Insig1，Idi1在睾丸组织内不同细胞类型中表达情况。如图5所示，揭示了每个基因在睾丸内各组织细胞的表达情况，而其中每个基因在睾丸支持（Sertoli）细胞和睾丸间质（Leydig）细胞内的表达情况都要显著高于其他细胞，表明这几个基因在Sertoli细胞和Leydig细胞内起着关键作用。
图4基因集富集分析
Fig.4Gene Set Enrichment Analysis
注：（a）胆固醇恒定（|NES|＞1，P＜0.05，FDR＜0.25）;（b）MTOCR1通路（|NES|＞1，P＜0.05，FDR＜0.25）;（c）雄性激素响应（|NES|＞1，P＜0.05，FDR＜0.25）;（d）GSEA共基因网络.
图5Hub基因在睾丸组织内不同细胞类型中表达
Fig.5Expression of hub genes in different cell types in testicular tissue
注：（a）Dhcr7;（b）Cyp51;（c）Insig1;（d）Idi1.
3 讨论
邻苯二甲酸二（2-乙基己）酯（DEHP）是一种全球性的环境污染物，广泛分布环境介质，具有生殖毒性。邻苯二甲酸二（2-乙基己）酯已被确定为可以引起机体生殖激素紊乱，导致生殖系统凋亡等，母体在怀孕期间暴露于DEHP可能导致雄性后代睾丸损伤和生育能力障碍。然而，这种损伤背后的详细机制尚不清楚。因此，基于生物信息学RNA转录组分析进行了本研究，深入分析DEHP暴露后致胎儿生殖系统毒性损伤的作用机制，为DEHP的生殖毒性研究提供进一步理论依据。
有研究发现青春期前小鼠在DEHP暴露量为500 mg/（kg·d）时，会引发其精子发生失调、睾丸间质细胞受损[11]。而本研究在探究母鼠妊娠期间暴露邻苯二甲酸二（2-乙基己）酯的影响时也参考了500 mg/（kg·d）的剂量使用。于是首先通过基因表达综合数据库获取了对妊娠母鼠每日灌胃500 mg/（kg·d）的DEHP剂量后诞生的雄性小鼠睾丸RNA转录组测序数据。分析后鉴定出了14个上调基因和72个下调基因。对差异基因进一步富集分析表明，胆固醇代谢和类固醇生物合成、代谢及皮质醇生物合成代谢相关过程在GO和KEGG分析中显著富集，这些途径都是睾酮合成所不可或缺的，表明雄性小鼠的睾丸发育可能已经受到干扰。由PPI网络分析得到了10个Hub基因，而在KEGG与GO术语的关联分析中显示，几个关键基因如Cyp51，Insig1，Dhcr7在富集的功能通路都涉及到了多个通路术语，在类固醇代谢、胆固醇合成等雄性激素分泌功能中呈多样性分布，表明这几个基因可能在DEHP暴露后引起的雄性小鼠睾丸生殖系统破坏起着关键性作用。胆固醇是高等动物中主要的甾醇类物质，在大脑功能发育中起着至关重要的作用，胆固醇的前身是7-脱氢胆固醇（7-dehydrocholesterol，7-DHC），其通过DHCR7酶转化为胆固醇，因此当DHCR7遭到破坏时，会导致7-DHC无法转化为胆固醇而累积[12]。而7-DHC作为一种易氧化的脂质会因其不稳定性而使自身氧化，形成了7-DHC衍生的氧化甾醇，其具有毒性且会影响细胞生长和分化。已有研究发现某些药物都可以抑制DHCR7的功能从而导致7-DHC的累积和毒性氧化甾醇的产生[13]。甾醇14 α-去甲基化酶（Sterols 14 α-Demethylase，Cyp51）是一种关键的类固醇、胆固醇生成酶，Cyp51在哺乳动物的所有组织中表达，在睾丸中表达水平最高[14]。
为了避免单一的按照差异倍数筛选差异基因会遗漏具有微小变化的重要基因，通过GSEA对小鼠睾丸转录组数据进行了全面分析以探查DEHP引起胎鼠生殖功能损伤的潜在机制。研究发现，在GSEA分析中，不止睾酮合成方面途径受损，小鼠的雄激素响应和mTORC1通路也呈显著性富集，其中大多数核心基因都已在DEHP处理后呈下调。哺乳动物雷帕霉素靶蛋白（Mammalian target of rapamycin，mTOR）是一种重要的细胞生长调节因子，其参与了各种营养物质及蛋白质、脂质等生理功能相关的调控，并在调节细胞生长、分化、凋亡过程中起着核心调控作用。MTORC1是MTOR的一种复合体，主要负责细胞生长、增殖和蛋白质合成。雄激素由睾丸分泌产生，睾酮就是一种天然的雄性激素，而且睾酮在精子的发生中发挥着关键的作用，影响着生精细胞的成熟和分化，胆固醇是合成性激素的重要物质，加快睾酮发生，促进精子生成。而mTORC1信号通路在精子发生过程中会调节精原干细胞的增殖分化与自我更新，在mTORC1持续活化后会引起精原干细胞耗竭而导致精子发生障碍，因此mTORC1通路被认定为在精子发生和生殖障碍方面有着很大影响[15]。Hmgcs1，Hmgcr与Idi1是本研究在GSEA筛选出的三个显著性基因集中都共同存在的三个关键基因。3-羟基-3-甲基戊二酰辅酶A合酶1（human 3-hydroxy-3-methylglutaryl-Coenzyme A synthase1，HMGCS1）是胆固醇合成甲羟戊酸途径中的关键酶，可以通过影响胆固醇的含量来调节睾酮的合成[16]。3-羟基3-甲基戊二酰辅酶A还原酶（3-hydroxy-3-methylglutaryl-CoA reductase，Hmgcr）也是胆固醇生物合成途径的主要酶，Gao等[17]发现内分泌干扰物双酚A（bisphenol A，BPA）及其替代品双酚S（bisphenol S，BPS），双酚F（bisphenol F，BPF）在雄性小鼠暴露后会引起睾丸内胆固醇相关基因Hmgcs1，Hmgcr等下降，这表明DEHP可能与双酚类环境干扰物有相似的睾丸毒性诱导作用。人异戊烯二磷酸δ异构酶1（Isopentenyl-diphosphate delta isomerase1，Idi1）和人异戊烯二磷酸δ异构酶2（Isopentenyl-diphosphate delta isomerase2，Idi2）是异戊烯基二磷酸异构酶（Isopentenyl-diphosphate delta isomerase，IDI）的一种亚型，Nakamura等[18]发现Idi1可能与神经元中胆固醇代谢物的产生有关。通过单细胞数据库分析了差异分析与基因集富集分析选取的关键基因在睾丸组织细胞层次方面的表达分布情况，发现几个基因大部分都广泛分布在睾丸间质细胞与睾丸支持细胞中，另外之前的分析表明了睾丸激素水平可能受到了影响，因此推测DEHP对雄性小鼠诱导的生殖毒性很可能是在睾丸间质细胞受损发生的。综上所述，本研究发现在孕鼠DEHP灌胃暴露后可能会引起睾丸间质细胞受损影响Dhcr7，Cyp51，Insig1，Idi1等基因的表达引起mTORC1通路失调，导致了雄激素的下降致使精子发生紊乱及生殖系统破坏。
4 结论
通过转录组表达分析探究了DEHP破坏子鼠生殖系统的机制过程，基于差异表达分析和基因集富集分析筛选出了潜在关键基因及其破坏机制，推测DEHP暴露可能会破坏睾丸间质细胞引起胆固醇合成相关基因的表达降低从而影响mTORC1通路失调导致生殖系统紊乱及精子发生失调，缩小了未来的研究范围，为研究DEHP的生殖毒性暴露机制提供了一定数据参考。但本研究作为生物信息学研究，仍存在不完善之处，核心基因与DEHP诱导的生殖毒性机制仍需要进一步的体内外实验验证。
