近日，香港城市大学Patrick Lee团队在环境领域著名学术期刊The ISME Journal上发表了题为“Chronic exposure to polycyclic aromatic hydrocarbons alters skin virome composition and virus–host interactions”的研究论文。该研究分析了124名中国女性脸颊皮肤拭子的宏基因病毒组，发现长期暴露于不同水平的多环芳烃会显著改变皮肤病毒组的构成和功能，以及病毒与宿主之间的相互作用。这一发现对于理解空气污染如何影响皮肤健康具有重要意义，并为开发新的皮肤保护策略提供了科学依据。
研究意义
环境健康关联：皮肤作为人体与环境的物理屏障，其微生物群落受环境因素影响。多环芳烃PAHs 作为普遍的空气污染物，研究其对皮肤病毒的影响，能深入了解环境污染物与皮肤健康的关系。
微生物群落复杂性探索：揭示皮肤微生物群落中病毒的多样性、组成和功能，以及与细菌宿主的相互作用，有助于理解皮肤微生物群落的生态动态和功能机制。
疾病预防与治疗的潜在应用：研究结果为开发针对性的皮肤健康维护策略提供依据，例如利用病毒 - 宿主相互作用来防治因环境因素导致的皮肤疾病，或开发基于噬菌体的治疗方法。
研究结果
不同水平的多环芳烃暴露影响了皮肤病毒组成
从124个皮肤样本中获得267个高质量的病毒（完整度 ≥ 50%），聚类分析得到238个物种水平的病毒操作分类单元，多数为有尾噬菌体目，这表明噬菌体是皮肤病毒的重要组成。此外，以裂解性生活方式的噬菌体占主导。先前的宏基因组研究显示，皮肤微生物组可以分为1型和2型两种皮肤分型，它们在多样性、分类学和功能组成方面存在显著差异。2型与高水平的多环芳烃暴露相关，具有更高的细菌多样性和较低的痤疮丙酸杆菌相对丰度，发现于干燥和暗淡的皮肤表型中，而1型与低水平的多环芳烃暴露相关，具有较低的细菌多样性和较高的痤疮丙酸杆菌相对丰度，发现于正常皮肤表型中。2型受试者头发中10种高分子量多环芳烃的平均浓度显著高于1型受试者，而两种类型之间低分子量多环芳烃的平均浓度没有差异。偏最小二乘路径建模显示多环芳烃与细菌和病毒均显著正相关，表明细菌在多环芳烃和病毒之间的关系中扮演了关键角色。主坐标分析证实，皮肤病毒群落也可以被分为两个与皮肤细菌群落相对应的不同分型。
病毒功能与皮肤分型有关
研究发现，与皮肤细菌类似，1型和2型皮肤病毒的分类学和功能是一致的，表明具有相似分类学组成的样本倾向于具有相似的功能。通过分析病毒操作分类单元，发现皮肤病毒群落包含高度多样化的功能，尤其是编码DNA复制、修复和重组功能的蛋白质。在1型和2型皮肤分型中，共有755个KEGG KOs表现出差异富集，其中1型中有675个，2型中有80个。1型中富集的KOs涉及多种与生物合成相关的代谢功能，而2型则以一组更受限的代谢功能为特征。与精氨酸代谢途径相关的12个KOs在1型和2型皮肤病毒群落中表现出差异富集，其中11个在1型中显著富集。这些KOs与精氨酸代谢相关的两种代谢物显示出显著的正相关性。在1型中发现的一个有尾噬菌体病毒含有五个KOs，这些KOs编码参与精氨酸生物合成途径的蛋白质，表明1型中可能存在一个功能性的精氨酸生物合成操纵子。在2型皮肤病毒群落中，高水平的多环芳烃暴露富集了外源物质降解功能。研究还发现，与外源物质降解相关的九个KOs在病毒群落中也表现出富集，且这些KOs在2型中更为普遍，其中五个在1型中完全缺失。这些KOs与苯甲酸降解和多种化学物质的生物降解有关，且与皮肤样本中检测到的与苯甲酸相关的代谢物显示出显著的正相关性。
病毒-宿主相互作用与不同水平的多环芳烃暴露有关
病毒-宿主相互作用数量在两种皮肤分型中不同，皮肤分型1病毒-宿主相互作用数量多于皮肤分型2。在两种皮肤分型中，大多数细菌属的平均病毒-宿主比率大于1，表明病毒倾向于利用裂解感染来靶向和裂解大多数细菌。然而，包括丙酸杆菌和微球菌在内的少数细菌属的平均病毒-宿主比率接近1。在2型中，大多数细菌属水平（包括一些具有生物降解能力的细菌属）的病毒-宿主比率显著低于1型，且病毒与其宿主的相对丰度之间的相关性在两种类型之间有所不同，其中皮肤分型1符合 “Piggyback - the - Winner” 生态模型，皮肤分型2符合 “Piggyback - the - Persistent” 生态模型。在高水平多环芳烃暴露下，细菌中CRISPR-Cas系统的数量增加。这表明，2型中的细菌宿主可能进化出了更多策略来保护自己免受病毒感染。尽管宿主拥有CRISPR/Cas系统来抵御病毒感染，但它们似乎仍然被噬菌体感染，这可能是由于噬菌体中抗CRISPR蛋白的存在。然而，在1型和2型病毒群落中发现的抗CRISPR蛋白数量相似，表明多环芳烃暴露水平并不影响抗CRISPR蛋白的系统发育分布。
噬菌体编码的辅助代谢基因与高水平多环芳烃暴露相关
研究发现了153个辅助代谢基因，编码61种代谢功能，其中最常见的辅助代谢基因与氨基酸代谢相关，并且与裂解或溶原生活方式的噬菌体有关。部分辅助代谢基因在2型皮肤病毒群落中的平均覆盖度高于1型，例如编码磷酸盐/磷酸转运蛋白和queuosine生物合成的基因，它们位于相同的基因组区域并共享相同的启动子和终止子。研究发现与转运蛋白相关的五个辅助代谢基因富集在2型中。此外，在2型中富集的噬菌体中发现一个参与外源物质生物降解的辅助代谢基因，其编码芳基醇脫氫酶，且被整合到了宿主基因组中。同样，参与精氨酸生物合成（编码氨基酸N-乙酰转移酶）和脯氨酸代谢（编码吡咯啉-5-羧酸还原酶）的辅助代谢基因也存在于感染2型中富集的细菌的噬菌体中。
文中部分图表：
Figure 1. Differentiation of skin microbiomes into two cutotypes after exposure to low and high levels of polycyclic aromatic hydrocarbons (PAHs). (A) Workf low of sample processing. The illustration was made by biorender. (B) The skin phenotypes associated with cutotypes 1 and 2 after exposure to low and high levels of PAHs. The illustration was drawn based on the results from the reference (Leung et al., 2023).
Figure 2. Overview of cheek skin virome. (A) Taxonomy and lifestyle of viruses. (B) Relationships between PAHs, bacteriome, and virome based on partial least squares path modeling. The numbers next to the lines and in the brackets denote statistically significant direct and indirect path coefficients, respectively (∗P < .05 and ∗∗∗P < .001). (C) Principal coordinate analysis of the Bray–Curtis dissimilarity matrix of vOTUs in all cheek skin samples. Cutotypes 1 and 2 were defined based on the microbes in the samples. Ellipses represent 95% confidence intervals. (D) Pearson’s correlation between forehead sebum concentration and the first dimension of the viral community compositional principal coordinate axis (i.e. PCoA1). The linear regression line and coefficient of determination (R2) are shown.
Figure 3. Functional profile of skin viromes. (A) The prevalence of enriched KEGG Orthologies (KOs) involved in arginine biosynthesis (left panel) and xenobiotic degradation (right panel) pathways between the two cutotypes, as analyzed using MaAsLin2. The vertical bars indicate the enrichment of these KOs in cutotype 1 and cutotype 2. (B) KOs and their related metabolites involved in the arginine biosynthesis pathway are shown. KOs highlighted in color are those that are enriched in their respective cutotypes. Similarly, metabolites shown in color indicate significant positive correlations with the enriched KOs in their respective cutotypes. (C) Schematic illustration of genes in the arginine biosynthesis pathway carried by a Caudoviricetes virus in cutotype 1. Different colors represent various functional modules, and arrows indicate predicted terminators.
Figure 4. Numbers of virus–host interactions in the two cutotypes. (A) Left panel: Average relative abundances of viruses and their linked bacterial hosts. Right panel: Heatmap of the regression coefficient from the general linear model of MaAsLin2 showing differentially abundant hosts and viruses between the two cutotypes (Wilcoxon rank sum test (WRST), two-sided, *padj < .05, **padj < .01, ***padj < .001). Only genera enriched or depleted in cutotypes 1 and 2 are shown. A positive coefficient indicates the degree of enrichment, whereas a negative coefficient indicates the degree of depletion in cutotype 2. (B) Pearson’s correlations of relative abundances (log10scale) between viruses and their linked hosts of various Cutibacteriumspecies (including eight C. acnes strains) and Micrococcusspecies (including nine M. luteus strains) in cutotypes 1 and 2. The inset boxplots show the slopes of linear regression for strains of each species. (C) Virus-to-host ratio (VHR) at the bacterial genus level for each virus–host pair. The statistical significance of the difference in VHRs between cutotypes 1 and 2 is indicated on the x-axis (WRST, two-sided, *padj < .05, **padj < .01, ***padj < .001). A VHR > 1 indicates that the viruses tended to have a lytic lifestyle, whereas a VHR < 1 indicates that the viruses tended to have a lysogenic lifestyle. (D) Pearson’s correlations of the relative abundances (log10 scale) between viruses and their linked host (i.e. a Marmoricola sp.) in the two cutotypes. (E) Pearson’s correlations between the average VHR of selected Marmoricolaand Janibacter species and the measured concentrations of polyaromatic hydrocarbons (PAHs; only the VHR–PAH pairs with the two largest Pearson’s correlation coefficients (r) are shown). Linear regression lines and coefficients of determination (R2) are shown in panels d and e. In each box-and-whisker plot, the box indicates the median, first quartile, and third quartile; the whiskers span 1.5 times the interquartile range; and the diamond indicates the mean. Points that lie beyond the whiskers indicate outliers.
Figure 7. Schematic illustration of the potential mechanistic roles of viruses in cutotypes 1 and 2. (A) Differences in bacterial and viral diversity and functions. (B) Mechanisms of virus–host coevolution.
（本文由论文第一作者杜世聪博士撰写）
