PARP9 and PARP14 cross-regulate macrophage activation via STAT1 ADP-ribosylation

Despite the global impact of macrophage activation in vascular disease, the underlying mechanisms remain obscure. Here we show, with global proteomic analysis of macrophage cell lines treated with either IFNγ or IL-4, that PARP9 and PARP14 regulate macrophage activation. In primary macrophages, PARP9 and PARP14 have opposing roles in macrophage activation. PARP14 silencing induces pro-inflammatory genes and STAT1 phosphorylation in M(IFNγ) cells, whereas it suppresses anti-inflammatory gene expression and STAT6 phosphorylation in M(IL-4) cells. PARP9 silencing suppresses pro-inflammatory genes and STAT1 phosphorylation in M(IFNγ) cells. PARP14 induces ADP-ribosylation of STAT1, which is suppressed by PARP9. Mutations at these ADP-ribosylation sites lead to increased phosphorylation. Network analysis links PARP9–PARP14 with human coronary artery disease. PARP14 deficiency in haematopoietic cells accelerates the development and inflammatory burden of acute and chronic arterial lesions in mice. These findings suggest that PARP9 and PARP14 cross-regulate macrophage activation.

D espite medical advances, the global burden of ischaemic heart disease is increasing 1,2 . Pro-inflammatory macrophage activation plays key roles in the pathogenesis of many disorders, including arterial disease [3][4][5][6][7][8][9][10] . Some pathways associated with macrophage activation may contribute to the shared mechanisms of inflammatory diseases, as demonstrated previously 11,12 . Despite potent therapies such as cholesterol-lowering by statins, substantial residual cardiovascular risk remains 7,13,14 , which drives the active search for novel solutions against pro-inflammatory macrophage activation. Dissecting complex and intertwined mechanisms for macrophage activation requires well-defined mechanistic models. The evidence suggests that distinct types of macrophage activation are functionally different in disease pathogenesis, a classification that has helped to assess the heterogeneity of macrophages [15][16][17][18][19][20][21][22] . For instance, pro-inflammatory and anti-inflammatory phenotypes can oppose one another, develop in response to distinct cytokines, differ in the activating stimuli and produce different cytokines. A recently proposed nomenclature suggests that each macrophage subpopulation can be called based on a specific stimulator, for example, M(IFNg), M(LPS), M(IL-4), M(IL-10) 21 . This established paradigm demonstrates clear relationships between classical stimuli and their respective responses-interferon gamma (IFNg) for pro-inflammatory activation in settings such as atherosclerotic vascular disease and interleukin (IL)-4 for activation that can counter that of M(IFNg) or M(LPS) macrophages. Hence, we used this paradigm as a starting point to explore novel regulators through global proteomics.
Proteomics screening and bioinformatics in mouse and human data sets found that poly ADP-ribose polymerase 14 (PARP14), also known as ADP-ribosyltransferase diphtheria toxin-like 8 (ARTD8), and PARP9/ARTD9 both increased in M(IFNg) and decreased in M(IL-4) cells. The network analysis associated these PARP family members with human arterial disease. Sequence similarity to the PARP catalytic domain, which transfers ADP-ribose moieties from NAD to protein acceptors, characterizes the PARP family proteins 23 . The best-characterized member, PARP1/ARTD1, represents poly-ADP-ribosylation enzymes, which processively catalyse long and branching polymers of ADP-ribose additions starting from an initial post-translational modification, commonly of glutamate. Recent evidence also validates proteins that execute mono-ADP-ribosylation as having various functions 24 . PARP14/ARTD8 is an intracellular mono-ADP-ribosyltransferase. Previous reports indicated that PARP14 enhances IL-4-induced gene expression by interacting with the cytokine-induced signal transducer and activator of transcription 6 (STAT6) in B and T cells, thereby functioning as a transcriptional co-activator 25,26 that may mediate this effect. A recent study reported that PARP14 regulates the stability of tissue factor mRNA in M(LPS) in mouse 27 . Less information exists regarding the molecular function of PARP9/ARTD9. Although PARP9 appears to lack catalytic activity 28 , it increases IFNg-STAT1 signalling in B-cell lymphoma 29 .
This study employed a multidisciplinary approach, including proteomics, systems biology and cell and molecular biology to explore new mechanisms for modulating the functional profile elicited after macrophage activation. Mouse and human cell lines as well as primary macrophages were used for complementary analyses of PARP14-deficient mouse and human tissues. Ultimately, the analyses led to evidence that expression of PARP14 in haematopoietic cells restrains vascular inflammation in mouse models, which are not solely regulated by either IFNg or IL-4. Our findings suggest a novel mechanism for regulating the balance of macrophage phenotypes in vascular disease, and potentially other disorders in which macrophage activation has an impact on outcomes.

Results
Proteomics screening for regulators of macrophage activation. We used the tandem mass tagging (TMT) quantitative proteomics to identify regulators of pro-inflammatory and non/anti-inflammatory activation in mouse RAW264.7 and human THP-1 macrophage cell lines ( Supplementary Fig. 1a-c). In this paradigm of macrophage heterogeneity, IFNg and IL-4 promote distinctive subpopulations [15][16][17][18][19][20][21][22] . A pilot TMT proteomic study (Supplementary Figs 1d and 2) analysed the changes in the proteomes at 0, 12 and 24 h, and observed the expected increase and decrease in STAT1 in M(IFNg) and M(IL-4) cells, respectively, as determined by hierarchical cluster analysis ( Supplementary  Fig. 2). Within this pilot study, we first noted that PARP14 co-clustered with STAT1 in the M(IFNg) and M(IL-4) data ( Supplementary Fig. 2). To ascertain whether any changes in the M(IFNg) and M(IL-4) proteomes were not because of cell culture conditions, we performed a second, more in-depth study that included an unstimulated macrophage control for both RAW264.7 and THP-1 experiments, and extended the stimulation period to up to 72 h, sampling six time points for a more detailed timeresolved proteomic study ( Supplementary Fig. 1d).
In this latter proteomic study, we quantified 5,137 and 5,635 proteins in RAW264.7 and THP-1 cells, respectively, across the three conditions: unstimulated control, IFNg-stimulated macrophages and IL-4-stimulated macrophages-M(-), M(IFNg) and M(IL-4), respectively ( Fig. 1a and Supplementary Table 1). An overview of the protein intensities across the three conditions revealed that the magnitude of protein abundance levels for each of the IFNg and IL-4-stimulated RAW264.7 and THP-1 cells were generally higher than those of unstimulated cells (Fig. 1b), indicating that each stimulation promoted changes in protein abundance beyond those due to cell culture conditions alone.
To pursue one class of potential upstream regulators based on protein abundances, we used the following criteria: proteins exhibited (1) an early increase in M(IFNg) (within 24 h) followed by sustained levels until the later time points (up to 72 h), (2) a decrease in abundance in M(IL-4) and (3) no significant change in M(-). We employed two distinct informatics methods to explore proteins with such behaviours: data filtering (Method 1) and model-based clustering (Method 2; Supplementary Fig. 3). In Method 1, the M(-) data set permitted a subtraction of the background signal at all time points and thus facilitated data filtering ( Fig. 1b upper panels and Supplementary Fig. 3b). Three proteins in the data set from mouse RAW264.7 cells and 12 proteins in human THP-1 cells (Fig. 1b lower panels) met the three criteria above. PARP14 emerged from both data sets, and PARP9 appeared in the THP-1 data set ( Fig. 1b lower panels). In parallel, Method 2 ( Supplementary Fig. 3c) produced 15 and 20 clusters for the RAW264.7 and THP-1 data sets, respectively ( Supplementary Fig. 4a). Clusters mined for proteins whose abundances increased in M(IFNg) but decreased in M(IL-4) with respect to M(-) revealed 490 proteins in the RAW264.7 data set and 414 proteins in THP-1 fulfilling these criteria (Supplementary Fig. 4b and Supplementary Tables 2 and 3). The proteins were short-listed to those with similar time-resolved changes in both in human and mouse data sets, resulting in 38 candidate proteins (Fig. 1c). Both RAW264.7 and THP-1 data sets identified PARP9 and PARP14 ( Supplementary Fig. 4a). Collectively, while Method 2 identified PARP9 and PARP14 in RAW264.7 and THP-1 cells, PARP14 was the only common protein that both Methods 1 and 2 identified among over 5,000 proteins in mouse and human data sets.
The PARP9 and PARP14 network and coronary artery disease.
To understand the influence of PARP9 and PARP14 in a global ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12849 interaction network ('interactome') and predict their potential clinical impact in disease mechanisms, we applied a networkbased analysis ( Supplementary Fig. 5a). Increasing evidence suggests that disease genes are not distributed randomly on the interactome but work together in similar biological modules or pathways 30,31 . Moreover, gene products (for example, proteins) linked to the same phenotype likely interact with each other and cluster in the same network neighbourhood 31 . We thus postulated that, if PARP9 or PARP14 influences the network neighbourhood of a disease, its immediate neighbours should be close to a disease module compared with random expectation 30,31 . Using the random-walk method, we defined a set of genes as a human disease module for each of the cardiovascular and metabolic diseases, and IFNg-related diseases (see Methods). We then measured the average shortest distance of the immediate neighbours of the PARP9-PARP14 network to each disease module. An interactome mainly describes a set of physical or functional associations between proteins and does not provide  Hspd1  LRRC8D  DDX58  NFX1  MARK2  PNPT1  HSPD1  HIRA  NUDCD1  WDR75  NR3C1  HS2ST1   PELO  ANKLE2  CCDC50  NMI  MAP3K2  NUP188  CRNKL1  NUP155  TRRAP  SP110  OAS2  LCP2  ENSA  MRPL30  FCF1  MCM5  DNAJC9  FAM105A  PARP9  IFIH1  PARP14  RPL24  JUNB  IL1RN  HIST1H1C   Lower panels: extracted protein profiles of interest generated by data filtering. Red traces only graphs: extracted profiles of proteins whose abundances exceed the M(IFNg) threshold ( þ 0.13, maximum protein abundance in unstimulated control at 8 h, dashed line). From these M(IFNg)-filtered traces, those that decreased in IL-4 stimulation when compared with their unstimulated control are plotted to the right for RAW264.7 and THP-1 cells, respectively. (c) Hierarchical clustering of 38 proteins from that were identified in both RAW264.7 and THP-1 data sets. Each row corresponds to a protein gene ID. cell or tissue specificity. In our study, we thus examined the genes of immediate neighbours reported in macrophages in the public databases.
We included IFNg-related autoimmune diseases as positive controls, as we used this cytokine to promote pro-inflammatory macrophage activation. The PARP9-PARP14 network is significantly close to systemic lupus erythematosus, dermatomyositis and polymyositis, as expected (Fig. 2). The PARP9-PARP14 network had significantly greater proximity to the human coronary artery disease gene module compared with other cardiovascular and metabolic diseases ( Fig. 2 and Supplementary Fig. 5b). The analysis also linked PARP9-PARP14 with osteoporosis. Moreover, we quantified the closeness of PARP9-PARP14 to autoimmune diseases, coronary artery disease and other diseases in the form of the distribution of shortest distances ( Supplementary Fig. 5c) and demonstrated a clear separation between the distribution of diseases inside and outside the circle (Fig. 2), indicating an enrichment of shorter distances (d ¼ 1 and 2) for IFNg-related autoimmune diseases (yellow) and coronary artery disease and osteoporosis (blue), compared with other cardiovascular and metabolic diseases (red; Supplementary  Fig. 5c). These results may indicate the potential impact of PARP9 and/or PARP14 on the pathogenesis of arterial disease or the onset of its clinical complications.
In vitro validation in cultured macrophages. qPCR and western blot analysis validated the proteomic screening data on PARP9 and PARP14. Consistent with the proteomics data ( Fig. 3a and Supplementary Fig. 4a  PARP14 increased with IFNg and decreased by IL-4 ( Fig. 3b and Supplementary Fig. 6). At protein levels, PARP14 increased before PARP9 in response to IFNg (Fig. 3c).
PARP9 and PARP14 expression in plaque macrophages.
According to the network analysis that linked PARP9 and PARP14 with arterial disease, we performed immunohistochemistry in arterial lesions. Mouse (Fig. 3d, left) and human (Fig. 3d, right) atherosclerotic lesions exhibited PARP9 and PARP14 proteins, while they were less abundant in human carotid arteries with no apparent atherosclerotic changes ( Supplementary Fig. 7a). Immunohistochemistry localized PARP9 and PARP14 expression in the majority of macrophages (CD68) of human atherosclerotic plaques, while few if any smooth muscle cells (SMa-actin) and endothelial cells (CD31) stained positively for these PARPs ( Supplementary  Fig. 7b). These results suggest that macrophages are a major source of PARP9 and PARP14 in human atherosclerotic lesions. Pro-inflammatory PARP9 and anti-inflammatory PARP14. A series of subsequent experiments examined whether PARP9 or PARP14 plays a causal role in macrophage activation. Expression levels of gene products typical of the M(IFNg) phenotype (for example, IL-1b) and in the M(IL-4) phenotype (for example, MRC1) gauged the downstream effects of PARP silencing using small interfering RNA (siRNA) first in macrophage-like cell lines RAW264.7 and THP-1 since the original screening was performed in these cells, and then in human and mouse primary macrophages. In mouse RAW264.7 cells, PARP14 silencing enhanced the induction of tumour necrosis factor alpha (TNFa) and inducible nitric oxide synthase (iNOS) by IFNg, and suppressed the response of MRC1 to IL-4 ( Supplementary Fig. 8a). In human THP-1 cells, PARP14 silencing also increased TNFa and IL-1b mRNA induction in response to IFNg and decreased MRC1 induction by IL-4 ( Supplementary Fig. 8b). Increased levels of TNFa and IL-1b proteins in the supernatant of IFNg-treated THP-1 cells with PARP14 silencing (Supplementary Fig. 8c) supported these data. Silencing of PARP9 suppressed the induction of TNFa, IL-1b and CCL2/MCP-1 mRNA under IFNg stimulation in THP-1 cells, while MRC1 showed no significant change in IL-4 ( Supplementary Fig. 8d).
To obtain unambiguous evidence for the role of PARP9 and PARP14 in macrophage activation, we extended the in vitro validation studies to mouse and human primary macrophages. In human primary macrophages derived from CD14 þ peripheral blood mononuclear cells (PBMCs), PARP9 silencing suppressed the expression of TNFa, IL-1b and CCL2/MCP-1 in IFNg-treated macrophages, but exerted no significant effects on MRC1 induction by IL-4 ( Fig. 4a). In contrast, PARP14 silencing accelerated the induction of TNFa, IL-1b and CCL2/MCP-1 by IFNg, and suppressed MRC1 in IL-4-treated macrophages (Fig. 4a). In mouse bone marrow-derived macrophages (BMDMs), silencing PARP9 and PARP14 exerted the effects similar to those in human primary macrophages (Fig. 4b).
Overall, siRNA experiments provided consistent results in macrophage-like cell lines and primary macrophages. Neither PARP14 nor PARP9 showed significant effects on viability, proliferation or apoptosis of mouse primary macrophages ( Supplementary Fig. 8e). In addition, enforced expression of PARP14 in THP-1 cells suppressed the induction of TNFa, iNOS, TLR2 and TLR4 in M(IFNg) ( Supplementary Fig. 8f,g). Collectively, these findings indicate that PARP14 suppresses IFNginduced responses and augments IL-4-responses in macrophages. In contrast, PARP9 promotes responses to IFNg.
PARP9 and PARP14 regulate STAT1 and STAT6 activation. Although expression patterns of PARP14 and PARP9 in M(IFNg) and M(IL-4) were comparable, their respective siRNA experiments yielded opposing results, suggesting the involvement of distinct signalling mechanisms. IFNg signalling involves activation (phosphorylation) of pro-inflammatory STAT1, while the IL-4 pathway uses anti-inflammatory STAT6 phosphorylation 32,33 . Immunofluorescence staining demonstrated enhanced intracellular colocalization of PARP14 and STAT1 in IFNg-treated THP-1 cells compared with unstimulated cells (Supplementary Fig. 9). Moreover, in human primary macrophages derived from CD14 þ PBMCs, PARP14 silencing accelerated IFNg-induced STAT1 phosphorylation and suppressed IL-4-promoted STAT6 phosphorylation (Fig. 4c). PARP14 may suppress pro-inflammatory macrophage activation by modulating the IFNg-STAT1 axis, and promote the antiinflammatory IL-4-STAT6 pathway. In contrast, PARP9 silencing decreased STAT1 phosphorylation in IFNg-treated human macrophages (Fig. 4c), indicating that PARP9 may activate IFNg-STAT1 signalling and induce pro-inflammatory activation. In addition, siRNA experiments in THP-1 cells that we used for proteomic screening produced similar results on the role of PARP9 and PARP14 on the STAT1 and STAT6 pathways ( Supplementary Fig 10a-c). Evidence suggests the participation of STAT3 in immune responses in various contexts 34 . Neither PARP14 nor PARP9 silencing, however, exerted significant effects on STAT3 phosphorylation ( Supplementary Fig. 10d,e).
The introduction of mutations at Glu657 and Glu705 to prevent ADP-ribosylation helped to investigate the functional relevance of these two sites within STAT1a. Overexpression experiments in mouse BMDMs revealed that mutations at Glu657 and Glu705 enhanced Tyr701 phosphorylation of STAT1a in IFNg-treated cells as compared with wild-type STAT1a (Fig. 6c). Using the same mutant STAT1a, we have observed a similar response of STAT1a phosphorylation in HEK293 cells ( Supplementary Fig. 12a,b). As a functional consequence, mutant STAT1a in mouse BMDMs increased expression of pro-inflammatory iNOS, IL-1b and CCL2/MCP-1 (Fig. 6c). Collectively, these data indicate that PARP14-mediated ADP-ribosylation of Glu657 and Glu705 may control Try701 phosphorylation of STAT1a, a potential mechanism for the interplay between ADP-ribosylation and phosphorylation in macrophage activation.

PARP14 deletion enhances acute arterial lesion development.
To provide in vivo evidence that PARP14 participates in arterial lesion formation and macrophage activation, we further used PARP14 À / À mice. Peritoneal macrophages from PARP14 À / À or PARP14 þ / þ mice enabled the examination of the macrophage phenotype. PARP14 À / À macrophages expressed higher mRNA and protein levels of iNOS and TNFa under IFNg stimulation and lower levels of MRC1 and Arg1 mRNAs under IL-4 stimulation compared with PARP14 þ / þ cells (Fig. 7a,b). PARP14 deficiency also enhanced phosphorylation of STAT1 induced by IFNg and decreased STAT6 phosphorylation in IL-4-treated  The ratio of phosphorylated STAT1 and STAT6 protein levels to total STAT1 and STAT6 (pSTAT1/tSTAT1 ratio and pSTAT6/ tSTAT6 ratio) in human primary macrophages (n ¼ 6 and n ¼ 5, respectively) of the PARP9 and PARP14 silencing experiments. * Po0.05 and **Po0.01, respectively, by Student's t-test. Error bars indicate s.d.
Our network analysis closely linked the PARP9-PARP14 module with coronary artery disease (Fig. 2). To first examine whether PARP14 indeed plays a role in arterial diseases, we used two different models: (1) acute mechanical injury in femoral     arteries of PARP14 À / À mice and (2) acute injury in mice that underwent bone marrow transplantation from PARP14 À / À mice. In a model of acute arterial responses due to wire-mediated mechanical injury, PARP14 deficiency enhanced neointima formation (Fig. 7d, top panels) and increased macrophage accumulation (Fig. 7d, bottom panels). Laser capture microdissection (LCM) of the intima of injured femoral arteries followed by real-time PCR demonstrated that PARP14 deficiency increased the levels of TNFa and iNOS mRNA and decreased Arg1 mRNA (Fig. 7e). The spleen is a reservoir of monocytes/ macrophages that releases these cells in response to an event in remote organs 37 . Flow cytometry of splenic cells revealed that, in CD11b þ Ly6G À monocytes/macrophages ( Supplementary  Fig. 13b), acute arterial injury increased the Ly6C high population in PARP14 þ / þ mice, which was further expanded by PARP14 deficiency (Fig. 7f and Supplementary Fig. 13c).
To examine the relative contribution of PARP14 in the haematopoietic lineage to the lesion development after acute mechanical injury, we performed bone marrow transplantation from PARP14 þ / þ and PARP14 À / À mice. In lethally irradiated PARP14 þ / þ mice whose bone marrow was reconstituted by PARP14 À / À bone marrow (BMT PARP14 À / Àþ / þ mice), neointima formation after injury was accelerated, as compared with control PARP14 þ / þ mice whose bone marrow was reconstituted by PARP14 þ / þ cells (BMT PARP14 þ / þþ / þ mice; Fig. 7g). These results indicate that PARP14 derived from the haematopoietic cell lineage, the majority of which are macrophages in the neointima, plays an important role in the development of arterial disease.
Haematopoietic PARP14 deficiency enhances atherogenesis. To further examine the role of PARP14 in chronic arterial diseases, we used high-fat/high-cholesterol-fed low-density lipoprotein receptor-deficient (LDLR À / À ) mice, an established model of atherosclerosis. Lethally irradiated LDLR À / À mice whose bone marrow was reconstituted with PARP14 À / À cells (BMT PARP14 À / À -LDLR À / À ) or PARP14 þ / þ cells (BMT PARP14 þ / þ -LDLR À / À ) underwent high-fat/highcholesterol feeding to develop chronic atherosclerotic plaques. Sixteen weeks after the initiation of an atherogenic diet, BMT PARP14 À / À -LDLR À / À mice exhibited more plaque formation and macrophage accumulation in the aortic root compared with BMT PARP14 þ / þ -LDLR À / À mice (Fig. 8a). The aorta of BMT PARP14 À / À -LDLR À / À contained higher expression levels of the pro-inflammatory TNFa and CCL2/ MCP-1, while MRC1 expression tended to be lower (Fig. 8b). BMT from PARP14 À / À mice also increased PARP9 expression in the aorta (Fig. 8b). These results are consistent with other in vitro and in vivo data that we have reported in the present study. Taken together, these findings indicate that PARP14 derived from the haematopoietic lineage, the majority of which are macrophages in arterial lesions, plays a protective role against the development of arterial diseases, which verifies our prediction by network analysis.
Macrophage-rich human atheroma and PARP9. Seeking additional in vivo evidence for the potential role of PARP9 and PARP14 in arterial disease involved immunohistochemical analysis of human carotid atherosclerotic plaques surgically removed by endarterectomy. PARP9 and PARP14 signals were predominantly localized in plaque macrophages, as indicated by the overlapping CD68-positive signal ( Fig. 8c and Supplementary  Fig. 14). Quantitative analysis demonstrated that more macrophages were immunoreactive for PARP9 in macrophage-rich plaques than in macrophage-poor plaques, while there was no significant difference in PARP14-positive macrophages (Fig. 8c).
While some macrophages in the human plaques coexpressed PARP9 and PARP14, other cells stained positively for either PARP9 or PARP14 alone (Supplementary Fig. 14). These lines of in vivo evidence indicate that macrophages in arterial lesions are heterogeneous, which may reflect diverse levels of pro-inflammatory activation in individual cells.
Single-cell analysis of primary human macrophages. The diversity in expression patterns of PARP9 and PARP14 in human plaques necessitated the additional assessment of macrophage subpopulations using single-cell gene expression profiling 38 Supplementary Fig. 15a.
The expression patterns of 91 target genes (Supplementary Table 4) in M(-) (n ¼ 268) and M(IFNg) (n ¼ 252; n ¼ 520 in total) derived from the three different donors. All cells were evaluated against each other based on their gene expression dissimilarity (Manhattan distance, http://www.nist.gov/dads/). We present the distance matrix as a distance-based graph. By combining the two conditions-M(-) and M(IFNg)-across three donors (six conditions in total), M(-) cells (green) and M(IFNg) (red) are clearly segregated, forming distinct clusters (Fig. 9a). However, there are 'trails' of M(IFNg) that appear in the M(-) cluster, suggesting potential heterogeneity. Although this may be attributable to donor-to-donor variations, similar patterns were observed within individual donors as well ( Supplementary  Fig. 15b). As we are particularly interested in the expression profiles of PARP9 and PARP14, we examined variations/ correlations in their expression levels in M(-) and M(IFNg) human primary macrophages. In both phenotypes, PARP9 and PARP14 were correlated. However, M(-) cells showed lower levels of variation for PARP9 and PARP14 mRNA expression than did M(IFNg) (Supplementary Fig. 15c). These findings suggest that the 'IFNg-polarized' human primary macrophages are heterogeneous.

Discussion
The present study provides new evidence that PARP9 and PARP14 regulate macrophage activation. The specific novel findings demonstrated in this report include the following: (1) PARP9 promotes IFNg-induced responses in mouse and human macrophages; (2) PARP14 suppresses IFNg responses in mouse and human macrophages; (3) PARP14 induces IL-4triggered responses in mouse and human macrophages; (4) PARP9 and PARP14 appear to have physical and functional interactions; (5) PARP9 and PARP14 closely interact with components of IFNg signalling in macrophages; (6) PARP14induced mono-ADP-ribosylation of STAT1 inhibits its phosphorylation; (7) the PARP9 and PARP14 interactomes have significant proximity to the coronary artery disease module; (8) PARP14 deficiency indeed accelerates macrophage activation and lesion development in mouse models of acute and chronic arterial diseases; (9) PARP14 in the haematopoietic lineage exerts protective effects against arterial diseases; and (10) single-cell analysis revealed that IFNg-stimulated human primary macrophages derived from CD14 þ PBMC contain subpopulations, in which PARP9 and PARP14 are closely associated with genes of IFNg signalling. Our findings provide insight into new mechanisms for macrophage activation that play a critical role in the pathogenesis of inflammatory arterial diseases, a global health burden.
We aimed to establish the unambiguous evidence for the role of PARP9 and PARP14 in macrophage activation using in vitro and in vivo studies. In this study, M(IFN) and M(IL-4) as two multidimensional models of macrophage heterogeneity were used in an unbiased proteomics approach. We then validated our key results on the functionality of PARP9 and PARP14 in mouse and human primary macrophages. Bioinformatic analysis of singlecell gene profiling in human primary macrophages revealed close links among PARP9, PARP14 and IFNg pathway-related molecules (for example, JAK1, JAK2, STAT1 and STAT6), suggesting that these PARP family members contribute critically to the process of M(IFNg) activation.
We took advantage of the availability of a well-characterized paradigm of macrophage heterogeneity. The balance of diverse macrophage phenotypes (for example, pro-inflammatory versus anti-inflammatory subsets) may regulate normal homeostasis and disease mechanisms. Effective and safe anti-inflammatory therapies may require the fine-tuning of the imbalance of diverse macrophage phenotypes (for example, suppressing excessively activated pro-inflammatory macrophages without compromising protective functions or with enhancing anti-inflammatory activation). Our major goal was to identify molecules or pathways that may regulate the delicate balance of macrophage heterogeneity using global proteomics of M(IFNg) and M(IL-4). We thus used stringent criteria to choose proteins that increased during M(IFNg) activation and decreased in M(IL-4). Our in vitro experiments indicate that the interplay between PARP9 and PARP14 may indeed regulate the M(IFNg)/M(IL-4) balance. We further validated the in vitro results in human arterial lesions and mouse models of arterial disorders, neither of which is solely regulated by IFNg or IL-4, to provide clinically translatable evidence. Human and mouse atherosclerotic plaque macrophages express PARP9 and PARP14. We identified four macrophage subpopulations in human lesions: PARP9 þ /PARP14 þ ; PARP9 þ /PARP14 À ; PARP9-/PARP14 þ ; and PARP9 À / PARP14 À . Single-cell gene expression analysis further revealed that IFNg-stimulated human primary macrophages-M(IFNg)are heterogeneous, which may support an emerging concept of the multidimensional model of macrophage activation 21 . In addition, subpopulations within these IFNg-stimulated cells may have different functions, as shown by associations with MMP-9, CD36, TLR2 and TLR4. Future studies addressing the functional significance of the heterogeneity in human primary macrophages may lead to the development of personalized medical solutions.
To explore the evidence for the anti-atherogenic role of PARP14 beyond in vitro assays, we used PARP14 À / À mice. Genetic deletion of PARP14 indeed promoted macrophage activation and accumulation in the intima of mechanically injured arteries, offering in vivo proof of concept. Moreover, significant acceleration of acute arterial lesion formation and chronic atherosclerosis was observed in the mice reconstituted with PARP14 À / À bone marrow. These findings suggest that haematopoietic PARP14 plays a central role in arterial disease.
Our study further revealed the new biology for understudied members of the PARP family-PARP9 and PARP14. ADP-ribosylation assays demonstrated that PARP14 has the ability to mono-ADP-ribosylate STAT1 and mass spectrometric analysis identified two ADP-ribosylation sites proximal to its key phosphorylation site. Interestingly, PARP9 suppresses PARP14dependent mono-ADP-ribosylation of STAT1. Furthermore, mutations of STAT1a ribosylation sites for PARP14 enhanced its phosphorylation and pro-inflammatory gene expression in macrophages. The interplay of PARP9 and PARP14 in STAT1a ADP-ribosylation and phosphorylation may in part explain why PARP9 and PARP14 exert opposing effects on IFNg-and IL-4-induced macrophage activation. Supplementary Fig. 16 demonstrates an overview of regulatory mechanisms for the IFNg-STAT1 and IL-4-STAT6 pathways by PARP9 and PARP14.
In this study, unbiased global proteomics and bioinformatics of IFNg-and IL-4-treated macrophages implicated PARP9 and PARP14 in novel mechanisms of macrophage activation. The subsequent network-based prediction of the close relationship between the macrophage PARP14-PARP9 module and human coronary artery disease genes supported the premise that our proteomics screen would effectively identify regulators of macrophage activation in the context of cardiovascular diseases. The approach was then followed by a series of in vitro and in vivo analyses on mouse and human cells/tissues that demonstrated the novel concept that PARP9 promotes and PARP14 suppresses IFNg-induced activation of macrophages ( Supplementary  Fig. 16). The present study also represents our strategy of target discovery research assisted by global proteomics screening and subsequent validation studies (Supplementary Fig. 17). Collectively, our discoveries indicate that inhibition of PARP9 and/or activation of PARP14 may attenuate macrophage-mediated vascular diseases, and also provide new insight into the development of effective therapies for other inflammatory disorders.

Methods
Cell stimuli for cell culture and TMT sample preparation. In this study, we used mouse and human IFNg 10 ng ml À l and IL-4 10 ng ml À 1 (R&D systems) as stimuli for macrophage activation, respectively. The murine monocyte/macrophage cell line RAW264.7 was obtained from American Type Culture Collection (ATCC, TIB-71, Rockville, MD) and maintained in 10% fetal bovine serum (FBS, Life Technologies) containing DMEM (Sigma) supplemented with penicillin and streptomycin (Corning) at 37°C in humidified 5% CO 2 . Before stimuli (IFNg 10 ng ml À 1 and IL-4 10 ng ml À 1 ), cells were starved for 24 h with 0.1% FBS-containing media. THP-1 was also purchased from ATCC (TIB-202) and maintained in RPMI 1640 medium in 10% FBS with penicillin and streptomycin at 37°C in humidified 5% CO 2 . The macrophage-like state was obtained by treating the THP-1 monocytes for 48 h with PMA (200 ng ml À 1 , Sigma). Mycoplasma contamination test was routinely performed (once a month).
Each cell culture experiment, unstimulated, INFg and IL-4, was prepared for isobaric labelling using the 6-plex TMT strategy (Pierce). For sample preparation the cells were lysed and proteolysed (Lys-C, Wako Chemicals) using the in-solution urea þ RapiGest (Waters) strategy detailed previously (Supplementary Fig. 1a) 39 . Tryptic peptides were labelled with TMT 6-plex reagent (Pierce), combined and desalted using Oasis Hlb 1cc (10 mg) columns (Waters). The peptides were then fractionated into 24 fractions based on their isoelectric focusing point (pH range of 3-10) using the OFF-gel system (Agilent; Supplementary Fig. 1b). The fractions were dried using a tabletop speed vacuum (Fisher Scientific), cleaned with the Oasis columns and resuspended in 40 ml of 5% acetonitrile (Fisher Scientific) and 5% formic acid (Sigma-Aldrich) for subsequent analysis by liquid chromatography/ mass spectrometry (LC/MS).
Liquid chromatography tandem mass spectrometry. The high-resolution/ accuracy LTQ-Orbitrap Elite (Thermo Scientific) analysed TMT peptide samples, and the Q Exactive (Thermo Scientific) and Elite analysed in vitro ribosylated peptides. Both mass spectrometers are fronted with a Nanospray FLEX ion source, and coupled to an Easy-nLC1000 HPLC pump (Thermo Scientific). The peptides were subjected to a dual column set-up: an Acclaim PepMap RSLC C18 trap column, 75 mm Â 20 cm (50 mm Â 15 cm on the Q Exactive), and an Acclaim PepMap RSLC C18 analytical column 75 mm Â 250 mm (Thermo Scientific). For TMT analysis the analytical gradient was run at 250 nl min À 1 from 10 to 30% Solvent B (acetonitrile/0.1% formic acid) for 90 min, followed by 5 min of 95% Solvent B. Solvent A was 0.1% formic acid. For ribosylated peptides the gradient was run at 250 nl min À 1 from 5 to 28% Solvent B for 10 or 30 min, followed by 5 min of 95% Solvent B. All reagents were HPLC-grade. For TMT analysis, LTQ-Orbitrap was set to 120 K resolution and the top 20 precursor ions (within a scan range of 380-2,000 m/z) were subjected to higher-energy collision-induced dissociation (HCD, collision energy 40%, isolation width 3 m/z, dynamic exclusionenabled, starting m/z fixed at 120 m/z and resolution set to 30 K) for peptide sequencing (MS/MS). The Q Exactive was set to 140 K resolution with a top 10 precursor selection method (scan range of 380-1,500 m/z). HCD was set to a stepped normalized collision energy of 25±10%, isolation width of 1.6 m/z, dynamic exclusion-enabled and resolution set to 17.5 K for MS/MS. Ribosylated peptide candidates were screened in the MS/MS scan by the m6 peak of 348.1 (refs 40,41). Unmodified forms were calculated by subtracting the mass of the ADP-ribose (541.06 Da) from the observed precursor. Modified and unmodified m/z values and corresponding retention time windows were submitted to an inclusion list and analysed in using the data-independent acquisition module of the Q Exactive (R ¼ 35 K).
The MS/MS data were queried against the mouse or human UniProt database (downloaded on 27 March 2012) using the SEQUEST search algorithm via the Proteome Discoverer (PD) Package (version 1.3, Thermo Scientific) 42 , using a 10 p.p.m. tolerance window in the MS1 search space and a 0.02 Da fragment tolerance window for HCD. Methionine oxidation was set as a variable modification, and carbamidomethylation of cysteine residues and 6-plex TMT tags (Thermo Scientific) were set as fixed modifications. The peptide false discovery rate (FDR) was calculated using Percolator provided by PD: the FDR was determined based on the number of MS/MS spectral hits when searched against the reverse, decoy mouse or human database 43,44 . Peptides were filtered based on a 1% FDR. Peptides assigned to a given protein group, and not present in any other protein group, were considered as unique. Consequently, each protein group is represented by a single master protein (PD Grouping feature). Master proteins with two or more unique peptides were used for TMT reporter ratio quantification (Supplementary Table 2 contains a summary of peptides and peptide-spectrum matches (PSMs) for PARP9 and PARP14). Ribosylation spectra were manually annotated 40,41 .
Proteomics normalization and filtering strategy. For each PSM the TMT ion channel intensities were normalized to the time-zero channel (reference normalization, for example, R i ¼ xi x1 ; where i ¼ 1, 2, y, 6 and x 1 is the time-zero abundance). The protein's abundance was then calculated by taking the median of its corresponding PSM ratios 45 . To extract the proteins that increase in IFNgstimulated (M(IFNg)) but decrease in IL-4-stimulated condition (M(IL-4)), we applied a simple filtering logic that exploited the available unstimulated control data set. In this study, unstimulated control, M(-), serves as a control for the biological signal owing to the cell culture condition; that is, protein traces M(IFNg) and M(IL-4) that exceed the edges of the baseline are more likely to be bona fide IFNg-and IL-4-specific responses (Supplementary Fig. 3b). Protein profiles whose abundances surpassed the baseline, which was defined as maximum relative abundance at 8 h of M(-) ( þ 0.13, log10 of relative abundance), after supplement of INFg established the general threshold for the entire time M(IFNg) data set. This cutoff value was the same that was used for both RAW264.7 and THP-1 M(IFNg) data sets. Moreover, the proteins extracted from the IFNg-stimulated filtering step were cross-referenced not only to M(-) but also to M(IL-4) data sets where proteins were expected to possess opposite profiles with respect to M(IFNg) ( Supplementary  Fig. 3b). Therefore, the final list of candidate proteins (Fig. 1b) had profiles whose abundances increased in M(IFNg) and decreased in M(IL-4), respectively, beyond their baseline controls.
Proteomics data clustering. Clustering was performed using the model-based algorithm 46 in R, which is based on finite mixture models; as such, it can successfully be applied to time series data analyses 44 such as those acquired for M(IFNg) and M(IL-4) 47 (Supplementary Fig. 3c). In this approach, each time series (protein profile) yi, i ¼ 1, 2, y, N (where N is the number of TMT channels) is considered to be a single entity connected by a line. Clustering is achieved for a traditional finite mixture model by assigning each time series, yi, to a cluster. The non-supervised model-based clustering uses the Expectation Maximization algorithm to assign the profile to a specific cluster. MCLUST has the advantage of using the Bayesian Information Criteria to determine the number of clusters that best partition the data set by maximizing the intradata set variability. Finally, stronger covariance (and thus also dependencies) between sets of two time points is enabled by sum normalization 48 . After clustering was performed ( Supplementary Fig. 4a) we focused on proteins shared in all three data sets (unstimulated control, IFNg-stimulated and IL-4-stimulated) for both RAW264.7 and THP-1 cells, with the purpose of identifying proteins whose profiles increased in IFNg-stimulated condition and decreased in IL-4-stimulated, but maintained at basal levels in unstimulated control condition. We inspected all clusters looking for an increase in the relative abundance of proteins with respect to time zero at any time point in IFNg-stimulated, and subsequently cross-referenced the IL-4-stimulated and unstimulated control clusters for proteins whose profiles decreased and remained within the baseline, respectively. The proteins that fulfilled these three criteria are shown in the heat maps for RAW264.7 and THP-1 cell data sets ( Supplementary Fig. 4b). Finally, to further narrow the list of proteins of interest, we prioritized those that were common to both species (Fig. 1c), as determined by identical Uniprot protein IDs. Heat maps were used only with the aim of ordering the protein expression levels according to the row obtained Z-score (that is, z ¼ (x-mean)/s.d.). The proteins were clustered in the horizontal direction hierarchically using Euclidean distance and average linkage methods ( Fig. 1c and Supplementary Fig. 4b). In these data matrices, each column represents the TMT time point analysis 0, 8, 12, 24, 48 and 72 h, and each row corresponds to a protein gene ID.
PARP14-PARP9 network analysis. As a parallel approach to further investigate the candidacy of PARP9 and PARP14 as cardiovascular and metabolic diseases, we turned to in silico or network-based prediction methods under the premise that a potential PARP14-PARP9 interactome would be close in vicinity to its pertinent vascular disease network/module.
To evaluate the impact of PARP14-PARP9 neighbours, we used the HumanNet interaction database 49 . We then used the random-walk methodology 50 to construct the disease modules from the functional database using gene-disease associations extracted from genome-wide association studies and from the OMIM (http://www.ncbi.nlm.nih.gov/omim) and the MalaCards (http://www.malacards.org) databases for: Cardiomyopathy, Coronary Heart disease, Heart Failure, Hypercholestermia, Hypertension, Metabolic Traits, Osteoporosis, Cardiovascular Risk Factors, Sudden Cardiac Arrest, Systemic lupus erythematosus, Sjögren's syndrome, polymyositis and dermatomyositis. The PARP9-PARP14 module was determined by first neighbours in the functional database. Further, we restricted the first neighbours to those that are expressed in macrophages based on gene expression. The motivation of this network analysis is based on two recent studies of ours that investigated the localization of disease in the network 51 and measure the separation of disease modules in the network 52 , both of which have proved instrumental in identifying disease modules. These network topology-based methods have been rigorously validated in both of these peer-reviewed publications by means of biological enrichment as well as compared with other network-based methodologies. On the basis of this selection, we considered 55 first neighbours for PARP9 and 149 for PARP14 that we annotated collectively as the PARP9-PARP14 module. Further, we measured the closeness of the first neighbours using the shortest-path topology measure with the disease modules above. For each disease module, we calculated the shortest-path distances to the PARP9-PARP14 module and compared these distances to the random distance distribution with the same module size (Wilcoxon test and Benjamini-Hochberg correction for multiple testing). Finally, we verified that the module size does not affect the significance of the P values for the Coronary Heart Disease and Osteoporosis modules by recalculating the random distances by reducing the module size to the top 35 genes from random walk.
Mouse peritoneal macrophages. Four days after intraperitoneal injection of 2.5-3 ml of 4% thioglycollate (Fisher Scientific), peritoneal macrophages were obtained by injecting 5 ml of ice-cold PBS into the peritoneal cavity using 25 G needle, followed by collection of the fluid using 23 G needle. In all, 5 Â 10 5 cells were cultured on 24-well plates (Corning) with RPMI containing 10% FBS. Nonadherent cells were discarded after 16 h. After washing with PBS, cells were incubated with IFNg 10 ng ml À 1 and IL-4 10 ng ml À 1 for 24 h until harvesting.
LCM and RNA amplification. LCM was performed on the Leica LMD6500 Microdissection System. Neointima were cut using the following LCM parameters: power, 50 mW; pulse duration, 2 ms; and spot size, 20 Hm. RNA was isolated using the PicoPure RNA Isoaltion Kit, followed by RNA amplification using the RiboAmp HS Plus RNA Amplification Kit (both Arcturus, Mountain View, CA, USA), according to the manufacturer's protocol. PCR array was performed using the Fluidigm PCR system.
ELISA. The amounts of TNFa and IL-1b released into the culture media after stimulation were measured using an ELISA kit following the manufacturer's instructions (Duoset Kit, R&D). The culture medium of unstimulated macrophage was used as the negative control. Standard, control or sample solution was added to the ELISA-well plate, which had been pre-coated with specific monoclonal capture antibody. After being shaken gently for 3 h at room temperature, the polyclonal anti-TNFa antibody, conjugated with horseradish peroxidase, was added to the solution and incubated for 1 h at room temperature. A substrate solution containing hydrogen peroxidase and chromogen was added and allowed to react for 20 min. The levels of cytokines were assessed by a plate reader at 450 nm and normalized with the abundance of standard solution.
Nitrate quantification. To quantify nitric oxide in cell culture media of macrophages, mimicking iNOS concentration, the Griess Reagent Kit (G-7921, Life Technologies) was used.
Co-IP. Cells were lysed in IP lysis buffer (Thermo Scientific). A volume of 100 mg of protein was incubated with PARP14 antibody (5 mg, Invitrogen) and Dynabeads streptavidin (Life Technologies) by rotation overnight at 4°C, followed by washing three times with PBS/Tween 20 (0.02%), using a magnet to collect the beads after each wash. Five per cent of the precipitated protein sample was subjected to SDS-PAGE. Protein expression was detected using Pierce ECL western blotting substrate Reagent and ImageQuant LAS 4000.
ADP-ribosylation assays. Recombinant human PARP14 and PARP9 (BPS Bioscience Inc.) proteins and BSA (Sigma-Aldrich) were incubated with recombinant human STAT1a (OriGene Technologies Inc.) or STAT6 protein (Sino Biological Inc.) at a final concentration of 5 ng ml À 1 in the presence of 100 mM b-nicotinamide adenine dinucleotide hydrate (NAD; Sigma-Aldrich) or 6-biotin-17-NAD (Trevigen Inc.) in 50 mM Tris-HCl buffer (pH 7.4) for 1 h at room temperature. Ribosylation of STAT1a and STAT6 was detected by liquid chromatography tandem mass spectrometry (LC-MS/MS) after trypsin digestion (using the biotin-free NAD reaction) or by western blotting using Streptavidin-HRP (Abcam) after SDS-PAGE. Quantification of the relative abundances of ribosylated STAT1a peptides was completed by calculating the area under the curve (AUC) of the extracted ion chromatograms of the monoisotopic peaks of the modified versus unmodified peptides. The ratios were reported as AUC mod./(AUC mod. þ AUC unmod).
The raw real-time PCR reads for each array were transformed into n Â m matrices using Python's Pandas libraries (http://pandas.pydata.org/; where n ¼ cell index and m ¼ gene). Each data matrix was then processed and analysed using an in-house developed platform.
We first run a check for each column (genes) to see whether there are undetected values (CT ¼ 999) interspersed among positive reads (CTo25). If less than 10% of the reads are positive, we substitute those values with 999, and consider the entire gene undetectable for this array. On the other hand, if more than 10% of the genes are positive, their corresponding reads with undetected values are substituted with the average value of the positive reads. Granted, more sophisticated missing-value imputation (MVI) techniques exist 59 ; however, we may not have enough features (96 features) to fully benefit from MVI nor is it clear whether this necessarily leads to improvement in signal. Moreover, genes requiring missing-value estimation tend to fall near the limit of detection, and are unlikely to benefit fully from MVI. For the purpose of having a fully populated matrix with no missing values, the averaged value therefore should suffice. The final missing-value adjusted reads are converted into log2exp via the following equation (equation 1), where LOD stands for the limit of detection and set at recommended default value of 24.
Although we do not normalize using housekeeping genes (see below for normalization method), they can be good indicators of the overall read quality for a given cell. Cells without the housekeeping gene expression (that is, GAPDH) were removed from analysis. Next, we calculate the mean expression value of the GAPDH gene (average log2exp of all cells for GAPDH). Cells with outlier GAPDH expression (more than 3 s.d.'s from the mean) were also excluded from the analysis. Individual cells may exhibit extreme reads because of transcriptional burst, or because of high-instrument sensitivity [5]. Such extreme values need to be contained. Trimming methods (removal of extreme outliers based on the 5% quantile on either side of the read distribution) are unsuitable, as they would alter the dimensions of the matrix as well as lead to the loss of data. For example, given 100 genes with 100 observations each. If there is a 5% chance for every gene that at least one observation would contain an extreme outlier, then via the process of trimming, only 95 of the observations could be used for analysis since the outliers, and hence the entire column (containing useful information as well) would be discarded. To circumvent this, and maintain the overall integrity of the matrix, we performed Winsorization 60 , where we set boundaries of the extreme values to the values of the 5 and 95%th quartiles, respectively.
To compare the chips, individual genes are converted into z-scores by subtracting the mean from the log2exp, and division by the s.d.
Arrays are compared on a per-donor basis (Donor_i unstimulated control, IFNg-stimulated, where i ¼ 1,2,3). For each pair, the arrays are merged into a combined matrix. We calculate the Manhattan distance 61 between all cells based on their gene parameters, and represent the clustering using the minimum spanning tree 62 .
Statistical analysis. Data are given as mean ±s.d. Moreover, 'n' indicates the number of independent experiments or number of animals/samples. Tests with a P value less than 0.05 were considered statistically significant. Pairwise group NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12849 ARTICLE comparisons were performed using a Student's t-test (GraphPad prism 5, Prism Software Inc. (La Jolla, CA)). If F test showed the variance was significantly different, unpaired t-test with Welch's correction was performed. Exclusion criteria were set by Grubbs' test. The experiments were not randomized. No statistical method was used to predetermine sample size. The experiments were not randomized.
Data availability. The data that support the findings of this study are available within the article and its Supplementary Information Files.