Transcriptional profiling reveals developmental relationship and distinct biological functions of CD16+ and CD16- monocyte subsets

Background Human peripheral blood monocytes (Mo) consist of subsets distinguished by expression of CD16 (FCγRIII) and chemokine receptors. Classical CD16- Mo express CCR2 and migrate in response to CCL2, while a minor CD16+ Mo subset expresses CD16 and CX3CR1 and migrates into tissues expressing CX3CL1. CD16+ Mo produce pro-inflammatory cytokines and are expanded in certain inflammatory conditions including sepsis and HIV infection. Results To gain insight into the developmental relationship and functions of CD16+ and CD16- Mo, we examined transcriptional profiles of these Mo subsets in peripheral blood from healthy individuals. Of 16,328 expressed genes, 2,759 genes were differentially expressed and 228 and 250 were >2-fold upregulated and downregulated, respectively, in CD16+ compared to CD16- Mo. CD16+ Mo were distinguished by upregulation of transcripts for dendritic cell (DC) (SIGLEC10, CD43, RARA) and macrophage (MΦ) (CSF1R/CD115, MafB, CD97, C3aR) markers together with transcripts relevant for DC-T cell interaction (CXCL16, ICAM-2, LFA-1), cell activation (LTB, TNFRSF8, LST1, IFITM1-3, HMOX1, SOD-1, WARS, MGLL), and negative regulation of the cell cycle (CDKN1C, MTSS1), whereas CD16- Mo were distinguished by upregulation of transcripts for myeloid (CD14, MNDA, TREM1, CD1d, C1qR/CD93) and granulocyte markers (FPR1, GCSFR/CD114, S100A8-9/12). Differential expression of CSF1R, CSF3R, C1QR1, C3AR1, CD1d, CD43, CXCL16, and CX3CR1 was confirmed by flow cytometry. Furthermore, increased expression of RARA and KLF2 transcripts in CD16+ Mo coincided with absence of cell surface cutaneous lymphocyte associated antigen (CLA) expression, indicating potential imprinting for non-skin homing. Conclusion These results suggest that CD16+ and CD16- Mo originate from a common myeloid precursor, with CD16+ Mo having a more MΦ – and DC-like transcription program suggesting a more advanced stage of differentiation. Distinct transcriptional programs, together with their recruitment into tissues via different mechanisms, also suggest that CD16+ and CD16- Mo give rise to functionally distinct DC and MΦ in vivo.

A dramatic increase in circulating CD16 -Mo has been reported in inflammatory pathologies such as sepsis, HIV infection, tuberculosis, and asthma [11,24,25]. Studies of patients infected with Mycobacterium leprae demonstrated that CD16 + and CD16 -Mo differentiate into DC-SIGN + MΦ and CD1b + DC-SIGN -DC, respectively, and the presence of CD1b + DC-SIGN -DC in M. leprae lesions was associated with healing [26]. An increased frequency of CD16 + Mo was associated with non-healing Leishmania chagasi lesions [27]. Thus, CD16 + and CD16 -Mo differentiation into MΦ or DC subpopulations with distinct phenotypes influences host defenses in infectious disease. Consequently, there is interest in developing therapeutic strategies that target specific Mo subpopulations [6,8,28,29].
Mo heterogeneity is conserved across mammalian species [7,8,19,30]. In mice, Gr1 + CX3CR1 low Mo (homolog of human CD16 -Mo) are recruited into the peritoneal cavity or draining lymph nodes under inflammatory conditions by mechanisms dependent on CCR2 and CD62L, and subsequently differentiate into DC [19,[31][32][33]]. In contrast, Gr1 -CX3CR1 high (homolog of human CD16 + Mo) are constitutively recruited into peripheral tissues including spleen, gut, lungs, and brain [19]. Gr1 -CX3CR1 high patrol vascular endothelium by mechanisms involving LFA-1 and CX3CR1, and are rapidly recruited into inflamed tissues where they differentiate into MΦ expressing the transcription factors cMaf and MafB and transiently producing TNF-α [33]. Studies in CX3CR1/ApoE double knockout mice suggest that CX3CR1 high Mo play a critical role in development of atherosclerotic lesions [34,35]. CX3CR1 + Mo may be precursors for lamina propria DC, which depend on CX3CR1 to form transepithelial dendrites, enabling direct sampling of luminal antigens [36]. Furthermore, adoptive transfer studies in rats demonstrated that CCR2 low CX3CR1 high Mo are constitutively recruited into the gut where they give rise to intestinal lymph DC [37]. Studies [17,43], or IL- 10 [44], and upregulate CX3CR1 expression upon CCL2 stimulation via CCR2 [45]. Mo differentiation is associated with decreased CCR2 expression and increased CCL2 production [46]. Engrafted Gr1 high CX3CR1 low Mo in peripheral blood traffic to the bone marrow, differentiate into Gr1 low CX3CR1 high Mo, and contribute to mucosal, but not splenic, generation of DC [5]. These findings suggest that human CD16 + Mo and mouse Ly-6C low CCR2 low CX3CR1 high Mo differentiate from CD16 -Mo and Gr1 -Ly-6C high CCR2 high CX3CR1 low Mo, respectively.
Here, we investigate the developmental and functional relationship between CD16 + and CD16 -Mo subsets. Whole genome transcriptome analysis suggests that these Mo subsets originate from a common myeloid precursor, with CD16 + Mo being at a more advanced stage of differentiation and having a more MΦ -and DC-like transcription program. Upregulation of the transcription factors RARA and KLF2 in CD16 + Mo coincided with the absence of cutaneous lymphocyte associated antigen (CLA) expression, indicating potential imprinting for non-skin homing in CD16 + Mo. These results define distinct transcriptional profiles of CD16and CD16 + Mo subsets suggesting different stages of myeloid differentiation, new markers to distinguish these Mo subpopulations, and unique roles in immune responses and inflammatory diseases.

Flow cytometry analysis
Blood from healthy individuals was collected with informed consent and IRB approval from Dana-Farber Cancer Institute. PBMC isolated from peripheral blood by Ficoll-Paque gradient density centrifugation were stained with fluorochrome-conjugated Abs and analyzed by multi-color flow cytometry (BD FACSCalibur or LSRII).

RNA isolation and microarray analysis
Total RNA from Mo pellets was isolated by Trizol extraction and purified using RNeasy columns (Qiagen). The quality of RNA was assessed by visualization of intact bands corresponding to 18S and 28S rRNA on formaldehyde agarose gels. Total RNA (10 μg) from matched CD16 + and CD16 -Mo samples isolated from 4 different healthy donors was quality tested using an Agilent 2100 Bioanalyzer chip, reverse transcribed, and hybridized on the GeneChip ® Human Genome U133 Plus 2.0 Array (Affymetrix), which includes 54,000 probe sets on a single array (i.e., 47,000 transcripts and variants, including 38,500 well-characterized human genes). Primary data analysis performed using GeneSpring software (Biopolymer core facility, Harvard Medical School) generated Excel spreadsheets with relative gene expression values for the 4 matched CD16 + and CD16 -Mo subsets.

Microarray data analysis
A total of 16,328 probe sets were detected in these 8 samples (present calls, defined as probe sets detected in at least 3 samples). Normalization was performed as described [48] to account for variation between microarrays. Missing value estimation was performed using a modified KNN algorithm [49,50]. T-test was used to identify probe sets differentially expressed in CD16 + and CD16 -Mo (p < 0.05). These genes were sorted according to their t-statistics and fold change ratios, which were calculated by computing the mean expression in CD16 + and CD16 -Mo. Clustering analysis using fuzzy-c-means [50,51] was performed based on the genes selected by Ftest (n = 2,759 probe sets). False discovery rates (FDR) [52] were estimated using dChip software (build date: Jan 27 2009) [53] by performing 100 random permutations using all 8 samples with p-values < 0.05, expression ratio cut-off = 2.0-fold, and present call cut-off of 20%, yielding a median FDR of 0.07. Expression ratios for differentially expressed probe sets were calculated in CD16 + versus CD16 -Mo (cut-off 2-fold; p < 0.05). Heat maps for biological function categories were generated by dChip software using signal values from each of the 8 samples for genes that were > 2-fold upregulated or downregulated in CD16 + Mo compared to CD16 -Mo. The entire microarray dataset and technical information requested by Minimum Information about a Microarray Experiment (MIAME) are available at the Gene Expression Omnibus (GEO) database under accession number GSE16836 (Transcriptional profiling of CD16 + and CD16peripheral blood monocytes from healthy individuals) http://www.ncbi.nlm.nih.gov/geo.

Gene set enrichment analysis (GSEA)
Gene set enrichment analysis (GSEA) and Molecular Signature DataBase (MSigDB) http://www.broad.mit.edu were used to identify differentially expressed gene sets [54]. GSEA is a computational method that determines whether an a priori defined set of genes shows statistically significant concordant differences between two biological states (e.g. phenotypes). MSigDB contains more than 3000 gene sets for use with GSEA. An enrichment score (ES) that represents the difference between the observed and expected rankings from phenotype correlation was calculated for every gene set. A nominal p-value for the specific ES was then estimated from an empirical permutation-based null distribution that preserves the complex correlation structure of the gene expression data. Multiple testing was corrected via the FDR, with FDR less than 10% considered statistically significant.

Quantitative Real time RT-PCR
One step SYBR Green real time RT-PCR (Qiagen) was carried out in an iCycler BioRad EN270 PCR machine according to manufacturer's recommendations. Absolute quantification of target gene expression was performed using a 10-fold serial dilution of purified PCR products as described [55]. Briefly, 25-50 ng total RNA was reverse transcribed in 25 μl 1× SYBR Green mix (Qiagen) containing 0.5 μM primers, and 10 nM fluorescein calibration dye (Bio-Rad). Agarose gel electrophoresis was used to determine the size of amplification products (100-200 bp) and allowed cDNA purification (QIAquick Gel Extraction Kit; Qiagen) for standard curve preparation (i.e., 200, 20, 2, 0.2, and 0.02 fg cDNA). Primers spanning one or multiple exons were purchased from Qiagen (i.e., SIGLEC10, MafB, C1QR, C3AR1, CDKN1C, CSF1R, CSF3R, FcγRIII, TNFRSF8, ICAM-2 QuantiTect primer sets). Samples without template and reverse transcriptase were used as negative controls. The concentration of each gene was normalized to the 28S ribosomal RNA (RRN28S) internal control [55]. Each RT-PCR reaction was performed in triplicate.

Distinct gene expression profiles in CD16 + and CD16monocytes
To define transcriptional profiles of monocyte subsets in vivo, we performed genome wide transcriptome analysis of matched CD16 + and CD16 -Mo subsets in peripheral blood of four healthy individuals. We identified 2,759 probe sets that were differentially expressed and 13,569 genes that were similarly expressed in these Mo subsets ( Figure 1A) (GEO database accession number GSE16836, http://www.ncbi.nlm.nih.gov/geo). Clustering analysis separated the 8 samples into 2 groups that perfectly matched CD16 + and CD16 -Mo, with 1,402 genes downregulated and 1,357 genes upregulated in CD16 + compared to CD16 -Mo ( Figure 1B). Calculation of expression ratios for 2,759 differentially expressed probe sets showed that 250 probe sets were downregulated (corresponding to 166 genes and 23 unknown transcribed sequences) and 228 probe sets were upregulated (corresponding to 153 genes and 19 unknown transcribed sequences) in CD16 + compared to CD16 -Mo (cut-off 2-fold; p < 0.05) ( Figure  1C-D, Additional files 1, 2). These 2-fold lists of differentially expressed genes included known markers for CD16 + (i.e., FCGR3A/CD16, CX3CR1, ITGAL/LFA-1, and CD31/ PECAM1) and CD16 -Mo (i.e., CD14, CCR2, SELL/ CD62L, FCGR1/CD64) [18][19][20] (Additional files 1, 2), providing initial validation of microarray results. Signature transcripts for other blood cell lineages were absent in both CD16 + and CD16 -Mo (i.e., CD3 and CD8 for T cells, CD56 for NK cells, CD19 for B cells, and DC-SIGN and CD1c for DC), consistent with results obtained by flow cytometry demonstrating the purity of sorted Mo (>98%) and absence of DC (i.e., CD1c), NK cell (i.e., CD56), B cell (i.e., CD19), T cell (i.e., CD3), neutrophil (i.e., CD16b and CD66b) markers on CD16 + and CD16 -Mo. The difference in relative expression of some probe sets for donor #1 probably reflects normal donor-to-donor variability, since post-sort cell viability, RNA quality, and MicroArray Quality Controls were similar for the four donors. A more stringent analysis was performed where in addition to a cut-off >2-fold and p-value < 0.05, probe sets with expression levels >3-fold higher than background were selected; by this approach, we identified 132 downregulated and 183 upregulated probe sets in CD16 + compared to CD16 -Mo (data not shown). These genes were further selected for those with the highest levels of expression (>10,000 AU (arbitrary units), cut-off >2-fold; p-value < 0.05) and two lists of top genes were generated, with 30 and 31 transcripts upregulated in CD16 + and CD16 -Mo, respectively (Tables 1 and 2). Other genes were differentially expressed with a difference <2-fold. In CD16 + compared to CD16 -Mo, downregulated markers included the early myeloid markers CD13 (2,971 ± 1,753 versus 5,656 ± 2,392; ratio 0.53, p-value < 0.05) and CD33 (1,860 ± 703 versus 3,478 ± 686; ratio 0.52, p-value < 0.05). Thus, despite a high level of transcriptional similarity (approximately 83%), a subset of probe sets were significantly downregulated (n = 250) or upregulated (n = 228) in CD16 + compared to CD16 -Mo, suggesting that these Mo subsets represent different stages of myeloid differentiation and have distinct biological functions in vivo.

Biological functions of differentially expressed genes
Differentially expressed genes, corresponding to 250 downregulated and 228 upregulated probe sets in CD16 + compared to CD16 -Mo, were classified into eight functional categories using Gene Ontology. Heat maps for biological function categories ( Figure 5A-H) showed clear distinctions in patterns of gene expression between the Mo subpopulations for these categories.
Cell cycle, proliferation and differentiation CD16 + Mo were distinguished from CD16 -Mo by upregulation of the cell cycle related genes cyclin-dependent kinase inhibitor 1C (CDKN1C, p27, or KIP2, which is a negative regulator of cell proliferation [66] induced by TGF-β [67]), and metastasis suppressor 1 (MTSS1, a transcript involved in cytoskeleton organization missing in metastasis [68]). CD16 -Mo preferentially expressed mRNA for genes encoding the CD1d antigen (member of the MHC family that mediates presentation of primarily lipid/glycolipid antigens to T cells), and myeloid cell nuclear differentiation antigen (MNDA, which is expressed in human monocytes and granulocytes and earlier stage cells in the myeloid lineage [69]) ( Figure 5F). These results provide evidence that CD16 + Mo represent a more advanced stage of differentiation compared to CD16 -Mo.
Transcription factors CD16 + Mo expressed significantly higher levels of mRNA for several transcriptional factor genes including the macrophage transcription factor v-maf musculoaponeurotic fibrosarcoma oncogene homolog B (MafB, an essential determinant of the monocytic program in hematopoietic cells [7,70,71]), the pleckstrin homology, Sec7, and coiled-coil protein-binding protein (PSCDBP or CYBR, a cytohesin-1-binding protein expressed in NK cells stimulated with IL-2 and IL-12 that plays a role in integrinmediated cell adhesion [72]), the Kruppel-like factor 2 (KLF2, reported to license mature T-cells for trafficking from the thymus and recirculation through secondary lymphoid tissues [73]), and retinoic acid receptor, alpha (RARA, expressed in dendritic cells [74] and involved in myeloid differentiation [75] and imprinting for gut homing [76,77]). In contrast, CD16 -Mo preferentially expressed the aryl hydrocarbon receptor (AHR, a ligand dependent E3 ubiquitin ligase [78] and modulator of anti-viral immunity [79]) transcript ( Figure 5H). Both CD16 + and CD16 -Mo lacked PU.1 expression, a transcription factor upregulated in DC [71], as demonstrated by microarray analysis and RT-PCR (data not shown). Thus, CD16 + and CD16 -Mo express distinct transcription factors that may differentially regulate biological functions in vivo.

Gene set enrichment analysis (GSEA)
To extract further meaning from differentially expressed genes in CD16 + and CD16 -Mo, GSEA, a knowledge based approach for interpreting genome-wide expression profiles [54], was applied to test for sets of genes that share common biological functions. Enrichment scores (ES), nominal p-values, false discovery rate (FDR), and family wise-error rate (FWER) values were generated for a large number of gene sets for GSEA available on the Molecular Signatures Database (MSigDB) of the Broad Institute. One gene set was significantly enriched in CD16compared to CD16 + Mo: HADDAD_HPCLYMPHO_ENRICHED (p < 0.001; both FDR q-value and FWER p-value-< 0.2). According to MSigDB, this set includes genes enriched in CD45RA hi Lin -CD10 + versus CD45RA int CD7and CD45RA hi CD7 hi hematopoietic progenitor cells [80]. Four uncharacterized open reading frames upregulated in hematopoietic progenitor cells (i.e., C18ORF1, CYORF15B, C6ORF62, and C6ORF111) [80] were significantly enriched in CD16 -Mo.

Potential imprinting for non-skin homing in CD16 + monocytes
We demonstrated increased expression of RARA mRNA in CD16 + compared to CD16 -Mo ( Figure 5H and 6A). SLP-76 (Src-homology 2 domain-containing leukocyte specific phosphoprotein of 76 kDa), a RA-induced target [75], was significantly upregulated in CD16 + compared to CD16 -Mo (6363 ± 611 versus 3882 ± 565; CD16 + /CD16ratio 1.64; p = 0.005), indicative of RARA pathway activation in CD16 + Mo. Activation of the RARA transcription factor pathway leads to loss of skin homing potential in lymphocytes via downregulation of the cutaneous lymphocyte-associated antigen (CLA, an epitope on PSGL-1) [83] and imprinting for mucosal homing [77]. CLA expression was quantified on CD16 + and CD16 -Mo by FACS on PBMC from healthy individuals. CLA expression was undetectable on all CD16 + Mo and a fraction of CD16 -Mo ( Figure 6B). The frequency of CLA + CD16 -Mo was negatively correlated with the frequency of CD16 + Mo ( Figure 6C). CD16 + Mo express M-DC8 (an epitope on PSGL-1) [13], which was previously reported to be expressed by a subset of mucosal DC [84] (Figure 6D). CX3CR1 expression on CD16 + Mo ( Figure 4D) [18] may contribute to recruitment of these cells into CX3CL1 expressing tissues including the gut. These results indicate that CD16 + Mo, similar to RA-stimulated T-cells [77], lack expression of the skin-homing addressin CLA and therefore are potentially imprinted for non-skin homing. Because retinoic acid (RA) is an important factor driving myeloid differentiation [75], this reprogramming of CD16 + Mo homing potential may be in part a consequence of RARA pathway activation.

Discussion
In this study, we define transcriptional profiles of human CD16 + and CD16monocytes (Mo) and provide new insights into their developmental relationship and biological functions. Despite remarkable transcriptional similarity (approximately 83%), a significant number of transcripts were differentially expressed (n = 2,759), with 228 and 250 >2-fold upregulated and downregulated, respectively, in CD16 + compared to CD16 -Mo. Differentially expressed genes related to cell-to-cell adhesion and trafficking, immune responses and inflammation, metabolism and stress response, signaling and signal transduction, cell cycle, proliferation, and differentiation, cytoskeleton, and regulation of transcription. Gene set enrichment analysis (GSEA) demonstrated that CD16 + Mo are enriched in genes related to NK-mediated cytotoxicity, inositol phosphate metabolism, actin binding, and oxidative stress, while CD16 -Mo are enriched in genes related to hematopoietic cell lineage, receptor-mediated endocytosis, arginine and proline metabolism, NTHi pathway, and lipid binding. The transcriptional profiles suggest that CD16 + and CD16 -Mo subsets originate from a common myeloid precursor, with CD16 + Mo being at a more advanced stage of myeloid differentiation and having distinct biological functions in vivo.
CD16 + Mo expressed high levels of transcripts for RARA, which controls transcription of genes involved in cell trafficking and mucosal homing. RA imprints lymphocytes with non-skin mucosal homing properties by decreasing cutaneous lymphocyte-associated antigen (CLA, an epitope on PSGL-1) expression [83] and increasing expression of CCR9 and integrin beta 7, two mucosal addressins [77]. RA also controls reciprocal differentiation of Th17 and regulatory T cells [101], modulates myeloid gene expression and differentiation [75], and regulates survival and antigen presentation by DC [74]. Consistent with our hypothesis that the RARA pathway is activated in CD16 + Mo, we demonstrated CLA downregulation on these cells, together with upregulation of two RA-induced targets: SLP-76 [75] and CXCL16 [102]. RA induces mucosal-type DC, which produce TGF-β and thereby imprints T-cells for gut homing by inducing CCR9 and integrin beta 7 [76]. CD16 + Mo-derived MΦ and DC constitutively produce TGF-β [23,100], but whether they also instruct T-cells for gut homing remains to be determined.
KLF2 mRNA is expressed at very high levels and significantly upregulated in CD16 + compared to CD16 -Mo. KLF2 belongs to a family of zinc-finger transcription factors that is induced by PI3K signaling [103] and controls expression of several genes including those coding for CD62L, CCR7, integrin beta7, sphingosine-1-phosphate receptor (S1PR1) [73], and lymphotoxin beta [104]. CCR7 and CD62L are essential for migration into lymph nodes, S1P1 regulates T-cell thymic egress and recirculation [105], and integrin beta 7 mediates cell recruitment into Peyer's patches and mesenteric lymph nodes [106]. Together, these findings raise the possibility that preferential expression of KLF2 in CD16 + Mo may confer an increased potential for trafficking.
CD16 + compared to CD16 -Mo express very high levels of transcripts for cyclin-dependent kinase inhibitor 1C (CDKN1C or p57/KIP2) (18.4-fold increase) and metastasis suppressor 1 MTSS1 (5.7-fold increase). CDKN1C is a potent inhibitor of several G1 cyclin-dependent kinase (cdk) complexes, and negative regulator of G1/S cell cycle transition and cell proliferation [66]. CDKN1C [66] and MTSS1 [68] are candidate tumor suppressor genes, and their high expression is consistent with the inability of Mo to proliferate [7]. CDKN1C is induced by TGF-β [67], a cytokine known to induce CD16 + Mo differentiation [17,43]. Thus, our results are consistent with a potential link between TGF-β pathway activation and CD16 + Mo differentiation in vivo.
The CD16 + Mo subset includes two subsets with distinct levels of CD14 expression: CD14 high CD16 + and CD14 low CD16 + [11,112]. CD14 high CD16 + Mo exhibit a phenotype intermediate between that of CD14 high CD16 neg and CD14 low CD16 + Mo in terms of adhesion molecule (e.g., CL62L) and chemokine receptor expression (e.g., CCR2, CXCR2, and CX3CR1) [18]. Both CD14 high CD16 + and CD14 low CD16 + Mo contributed to the transcriptional profile of CD16 + Mo in this study. The expression of some genes we identified as markers for CD16 + Mo may be distinct on CD14 high CD16 + and CD14 low CD16 + Mo. Consistent with this prediction, we demonstrated intermediate expression of CD115 and CD114 on CD14 high CD16 + Mo compared to CD14 high CD16 neg and CD14 low CD16 + Mo, and high expression of CD93 and C3aR1, similar to that on CD14 high CD16and CD14 low CD16 + Mo, respectively (Additional file 3). These findings suggest a developmental relationship between these Mo subsets in which CD14 high CD16 neg Mo, CD14 high CD16 + Mo, and CD14 low CD16 + Mo represent sequential stages of monocyte differentiation [41].

Conclusion
Comparative transcriptome analysis of CD16 + and CD16 -Mo indicates that CD16 + Mo represent a more advanced stage of myeloid differentiation with a more MΦ -and DC-like transcription program, whereas CD16 -Mo are more closely related to a common myeloid precursor. Given the ability of CD16 + and CD16 -Mo to be recruited into specific tissues via distinct mechanisms, these Mo subsets are likely to give rise to DC and MΦ subpopulations with distinct phenotypes and roles in immunity and disease pathogenesis. Further studies to characterize phenotypic differences between CD16 + and CD16 -Moderived DC and MΦ are relevant for development of DCbased vaccines, and will also provide a better understanding of their functional roles in immune responses, inflammation, and disease pathogenesis.