Deficiency of myeloid PHD proteins aggravates atherogenesis via macrophage apoptosis and paracrine fibrotic signalling

Abstract Aims Atherosclerotic plaque hypoxia is detrimental for macrophage function. Prolyl hydroxylases (PHDs) initiate cellular hypoxic responses, possibly influencing macrophage function in plaque hypoxia. Thus, we aimed to elucidate the role of myeloid PHDs in atherosclerosis. Methods and results Myeloid-specific PHD knockout (PHDko) mice were obtained via bone marrow transplantation (PHD1ko, PHD3ko) or conditional knockdown through lysozyme M-driven Cre recombinase (PHD2cko). Mice were fed high cholesterol diet for 6–12 weeks to induce atherosclerosis. Aortic root plaque size was significantly augmented 2.6-fold in PHD2cko, and 1.4-fold in PHD3ko compared to controls but was unchanged in PHD1ko mice. Macrophage apoptosis was promoted in PHD2cko and PHD3ko mice in vitro and in vivo, via the hypoxia-inducible factor (HIF) 1α/BNIP3 axis. Bulk and single-cell RNA data of PHD2cko bone marrow-derived macrophages (BMDMs) and plaque macrophages, respectively, showed enhanced HIF1α/BNIP3 signalling, which was validated in vitro by siRNA silencing. Human plaque BNIP3 mRNA was positively associated with plaque necrotic core size, suggesting similar pro-apoptotic effects in human. Furthermore, PHD2cko plaques displayed enhanced fibrosis, while macrophage collagen breakdown by matrix metalloproteinases, collagen production, and proliferation were unaltered. Instead, PHD2cko BMDMs enhanced fibroblast collagen secretion in a paracrine manner. In silico analysis of macrophage-fibroblast communication predicted SPP1 (osteopontin) signalling as regulator, which was corroborated by enhanced plaque SPP1 protein in vivo. Increased SPP1 mRNA expression upon PHD2cko was preferentially observed in foamy plaque macrophages expressing ‘triggering receptor expressed on myeloid cells-2’ (TREM2hi) evidenced by single-cell RNA, but not in neutrophils. This confirmed enhanced fibrotic signalling by PHD2cko macrophages to fibroblasts, in vitro as well as in vivo. Conclusion Myeloid PHD2cko and PHD3ko enhanced atherosclerotic plaque growth and macrophage apoptosis, while PHD2cko macrophages further activated collagen secretion by fibroblasts in vitro, likely via paracrine SPP1 signalling through TREM2hi macrophages.

Statistical analyses were performed using a student t-test (C) or two-way ANOVA with Bonferroni posthoc test (D). All results show mean ± SEM. *P<0.05 ***P<0.001.  56 and LysMCre transgenics, 57 were previously described. All PHD lines were crossed to low density lipoprotein receptor knockout (LDLRko) mice, obtained from an in-house breeding colony, originating from Charles River

ATHEROSCLEROSIS QUANTIFICATION AND IMMUNOHISTOCHEMISTRY
Mice were euthanized with a pentobarbital overdose (100 mg/kg i.p.) and blood was withdrawn via the right ventricle for flow cytometry, absolute white and red blood cell counts (Coulter Ac.T diff, Beckman Coulter) and total cholesterol analysis. Mice were perfused via the left cardiac ventricle with PBS containing sodium nitroprusside (0.1 mg/ml; Sigma-Aldrich, Seelze, Germany). Aortic arch, root and organs were subsequently excised and fixed in 1% paraformaldehyde overnight and paraffinembedded.
Aortic roots and arches were serially sectioned (4µm) and stained with hematoxylin and eosin (H&E, Sigma) for plaque area and lipid core content quantification. Five consecutive H&E sections at 20 µm intervals were analyzed blindly using computerized morphometry (Leica QWin V3, Cambridge, UK) and the sum of the three valves averaged per mouse. Necrotic core was defined as a-cellular and anuclear plaque area containing cholesterol clefts, and shown as the percentage necrotic content of the total plaque area. Sections within this 100 µm interval were used for remaining immunohistochemical stainings. Antigen retrieval was performed at pH 6 (Dako REAL target retrieval, Dako) (for MAC3, αSMA, and collagen type I), pH 9 (tris-EDTA, made in-house) (for platelet derived growth factor receptor beta (PDGFRβ)) or trypsin digestion (for CD31). Mouse atherosclerotic plaques were
Following antigen retrieval, PHD2 (Novus Biologicals NB100-2219) and CD68 (DAKO, M0814) were analyzed by non-fluorescent immunohistochemistry, followed by multispectral imaging to convert into pseudo-fluorescent images. Multispectral imaging (MSI) was performed to analyze human PHD1 and 2 expression, and PHD-CD68 co-localization. Spectral images were taken between 420-720 nm (10 nm interval) at a 10x (human) magnification using a Nuance spectral imaging system (Perkin Pseudo-colors were assigned to unmixed images and composite images showing co-localization were generated with the Nuance 3.0 software. Stainings were performed on serial sections. Images were converted to pseudo-fluorescence in Fiji.

PROLIFERATION ASSAY
Proliferation of SMC and 3T3 fibroblasts in response to macrophage conditioned medium (see above) was measured using CellTiter-glo luminescent cell viability assay (Promega, G7570) to determine ATP content of cells according to manufacturer's protocol. Primary SMC (2x104 cells) were isolated and seeded in a 96 wells plate and allowed to attach for 24 hours. After starvation for 24h in DMEM containing 0.5% FCS, BMDM-conditioned medium was added to the cells and incubated for 72 hours.
Luminescence was measured on a luminometer (Victor3, PerkinElmer) and proliferation calculated as difference on ATP between T0 and T72hrs. Proliferation of PHD2 WT and cko BMDMs was also measured on an ACEA xCELLigence (Roche). BMDMs (8x10 3 cells) were seeded on a gold electrode implemented in a 96 wells plate. Impedance was measured hourly and used to quantify proliferation (slope of impedance increment over time) using RCTA software (version 1.2, Roche).

MIGRATION ASSAY
Migration stimulation in primary murine SMC by conditioned macrophage medium of WT and PHD2cko was measured on an ACEA xCELLigence (Roche). SMC were starved in DMEM containing 0.5% FCS for 24 hours. Upper chambers of ACEA CIM 16 plates (ACEA, 20131122) were coated with 10ug/ml collagen G (Biochrom, L7213) for one hour per side prior to start of the experiment.
Subsequently, lower chambers were equilibrated for the respective conditioned mediums and controls.
BMDM-conditioned medium contained a final concentration of 1% FCS and 15% LCM. SMC (4x10 4 cells) were then added to the upper chamber and migration was monitored for 24h (hourly measurements), using the slope of the impedance increment over time.

INTRACELLULAR COLLAGEN CONTENT
Intracellular collagen content was measured using CNA35-FITC (Kindly provided by prof. Reutelingsperger, Biochemistry department Maastricht) shown to bind to collagen type I, III and IV. 60 SMCs and 3T3s were starved in DMEM containing 0.1% FCS for 24 or 48 hours, respectively. SMCs and 3T3s were subsequently treated with conditioned medium of either PHD2 WT or conditional knockout macrophages with or without a collagen producing stimulus (TGF-β1) (5 ng/ml, Biolegend, 763102) for 72 hours. Cells were fixed in 2% PFA for 15 min and permeabilized using 0.1% Triton X-100 in PBS for 15 min. Subsequently, cells were stained for internal collagen content with CNA35-FITC (1µM) and nuclei were stained with Hoechst (15µg/ml). Samples were analyzed using the BD Pathway 855 High Content Bioimager. Data was processed with Attovision and BD Diva software.

COLLAGEN SECRETION
After serum starvation, SMCs and 3T3s were treated with conditioned medium of either WT or PHD2cko macrophages for 72 hours. Culture medium of SMCs and 3T3s was collected after 72 hours and analyzed using Sircol soluble collagen assay as described by the manufacturer (Biocolor, S1000).
In comparable subsequent experiments, TGF-β1 was added to the conditioned medium or proteins were heat-inactivated (30 min, 85°C) prior to addition to 3T3 fibroblasts.

APOPTOSIS
BMDMs were stimulated with 50µM 7-ketocholesterol (Sigma, C2394) or 50µg/ml oxLDL (Isolated as described elsewhere) 61 for 24 hours to induce apoptosis. After stimulation, nuclei were stained with Hoechst (15µg/ml) and apoptotic cells with fluorescently labeled AnnexinxA5-FP488 (produced by Prof. Reutelingsperger, Maastricht) for 15 min. Samples were analyzed using a high-throughput, fluorescent reporter system, coupled to automated microscopy (BD Pathway 855 High Content Bioimager). Data was processed with Attovision and BD Diva software.

RNA SEQUENCING OF CULTURED CELLS
For RNA sequencing cells in vitro, RNA was isolated from triplicates of WT and PHD2cko BMDMs 24 hour after seeding, and from triplicate fibroblasts after 72 hours exposure to WT or PHD2cko conditioned medium. Bioanalyzer confirmed intact RNA (RNA Integrity number 10) for sequencing of 10µg RNA, on the NextSeq 500 system using v 2.5 chemistry, at ~15M single reads per sample by the c(Core Facility Genomics of the Medical Faculty Münster.
The gene-level expression of a total of 32,544 genes were quantified, 14,285 genes were retrieved for downstream analysis, 18,259 genes were discarded as lowly expressed genes.

BIOINFORMATICS ANALYSIS OF BULK RNA SEQUENCING DATA
Gene-level expression was quantified using Kallisto with the mouse genome (Mus Musculus GRCm38 assembly). 62 Principal component analysis was used for exploratory data analysis using the prcomp function in R (stats package, R 3.6.1 version). The limma R package (v3.40.6) was used to test for differential expression between conditions using the empirical Bayes method after voom transformation. 63 Lowly and non-expressed genes in the experiment were discarded from the analysis using the filterByExpr function (limma) to reduce potential false positives from the multiple testing.
Single-sample transcription factor activities were estimated using DoRothEA mouse regulons with A, B and C confidence classes. 64,21 Similar to differential gene expression analysis, the empirical Bayes method (limma) was also used to test for differential transcription factor activities using the TF activities (normalized enrichment scores) estimated by VIPER method (v.1.18.1). Genes and Transcription factor activities differentially dysregulated with FDR-adjusted p-values < 0.05 were considered significant. Pre-ranked Gene-Set Enrichment Analysis (GSEA) was performed using fgsea R package (v1.10.1) on the transcriptome-wide ranking of differential expression by the moderated tstatistics with the hallmark gene set collection from MSigDB and mouse gene sets from MatrisomeDB. 65,66 Human genes from the hallmark gene set collection were transformed to their orthologs in mouse using the biomaRt service from EnsEMBL. Gene sets with FDR < 0.05 were considered significantly enriched in the condition. Pathway analysis was performed using PROGENy 67,68 with the mouse model of pathway footprints of 100 genes, and 10,000 gene permutations of the ranking to build a null distribution for statistical estimations of significance. NicheNet 69 was used for the ligand-receptor analysis of stimulated pro-fibrotic fibroblasts by BMDM PHD2cko using the differentially over-expressed ligands from BMDMs PHD2cko (p-value < 0.05), and those target genes that were detected by the leading edge analysis in the significantly enriched MatrisomeDB gene sets from the pre-ranked GSEA in pro-fibrotic fibroblasts.

SINGLE CELL RNA SEQUENCING OF MURINE PLAQUES
PHD2 WT and PHD2 cko mice (n=11 and 9, respectively) were euthanized with a pentobarbital overdose (100 mg/kg i.p.) after 20 weeks of HCD and blood was withdrawn via the vena cava, followed by PBS perfusion via the left ventricle. The aortic root was subsequently excised and fixed in 4% paraformaldehyde overnight and paraffin-embedded.  Table S7).

BIOINFORMATICS ANALYSIS OF SINGLE-CELL RNA-SEQUENCING DATA
Raw sequencing data (FastQ files) were processed (alignment and gene-level expression quantification) using the CellRanger pipeline (10x Genomics, version 3.1) with the mouse genome (mm10 assembly). Seurat R package (v.3.1.0) was used to perform a standard analysis. 70 Quality control diagnostics were applied on library sizes, percentage of mitochondrial genes and gene detection coverage of the single cells. Cells with a gene coverage between 500 and 4,000 of genes expressed, and less than 7.5% of mitochondrial gene expression were retrieved to avoid bad quality cells, such as doublets and dead cells in downstream analysis. After data normalization using log transformation and applying a 10,000 scaling factor, the 2,000 most highly variable genes were selected using the variance stabilizing transformation method implemented in Seurat for each sample.
The first twenty principal components from Principal Component Analysis (PCA) applied on this selection of genes were used to find anchors for sample integration of the two conditions PHD2cko and WT pooled mice, and integrated using Canonical Correlation Analysis implemented in Seurat. 70 Myeloid leukocytes were identified following two rounds of unsupervised clustering. First, major clusters of cells expressing canonical myeloid leukocytes markers were selected. For this, PCA was applied on the scaled batch-corrected data to extract twenty-five principal components for the unsupervised clustering. Graph-based unsupervised clustering was performed using Shared-Nearest Neighbour algorithm. Louvain method was used to find clusters at resolution 0.1. Cluster of cells with positive expression of PTPRC, Lyz2, and CD68 expression, and absence of other vascular cell markers (CD3, CD19, MYH11, PECAM1) were detected as potential myeloid leukocytes. Second, this large and heterogeneous population of cells was selected for another iteration of unsupervised clustering to find more clusters with a higher resolution (at resolution 0.5), following the same workflow as described before. Distinct clusters were annotated based on the markers reported from Zernecke et al 2020, the top 10 marker genes for macrophages, monocytes, dendritic cells and top 10 overexpressed genes for Neutrophils as compared to the rest of cells. 53 . The main classes of myeloid leukocytes were identified among clusters, including cavity, IFNIC, inflammatory, Trem2-foamy and resident macrophages, mature-DC, moDC, monocytes and neutrophils. In addition, a cluster of proliferating cells was detected using the CellCycleScoring function from Seurat with the ortholog genes in mouse of the cell cycling genes, and another small cluster remained as not assigned (n.a.) due to the non-specific expression of cell-type markers. Uniform Manifold Approximation and Projection (UMAP) was used to reduce the dimensional space of the twenty five principal components to an embedding of two dimensions for visualization purposes with standard parameters in the Seurat package. Pathway and transcription factor activities were estimated using PROGENy (https://saezlab.github.io/progeny/) and DoRothEA (https://saezlab.github.io/dorothea/) in single-cell data as previously described (pathway footprints of 100 genes and TF regulons of A, B and C confidence classes, respectively). 68 The enrichment of the PHD2cko signature expression in singlecell basis was calculated using AUCell R package (v1.6.1) using the top 50 most up-regulated genes in BMDM upon PHD2cko for the PHD2 functional stratification of plaque-resident macrophages, prior discarding the main molecular players Bnip3 and Spp1 as potentially being defined as part of the signature. Correlations analysis between Hif1a activity, hypoxia response and BMDM PHD2cko expression signature was performed using Pearson correlation. For differential gene expression analysis, the two clusters of proliferating and n.a. cells were leave out, and the two clusters of resident macrophages found in the unsupervised clustering were merged as a single population of cells.
Differential gene expression was performed using Wilcoxon Rank Sum test. P-values from the Wilcoxon Rank Sum tests were adjusted for multiple testing using Bonferroni method for genome-wide multiple testing between two groups, and using FDR method for 2-group individual gene testing among cell types. Differentially expressed genes with adjusted p-values below 0.05 were considered statistically significant. R effect sizes from Wilcoxon Rank Sum test were calculated as Z divided by square root of the total observations. R effect size ranges from -1 to +1. Positive sign from r effect size relates to up-regulation in PHD2cko (group 1) as compared to WT condition (group 2). The greater the absolute r value is, the larger the effect size is. R effect size could be interpreted as small (r between 0.10 and 0.30), medium (r between 0.30 and 0.50) and large (r>0.50).