High-Resolution Transcriptome of Human Macrophages

Macrophages are dynamic cells integrating signals from their microenvironment to develop specific functional responses. Although, microarray-based transcriptional profiling has established transcriptional reprogramming as an important mechanism for signal integration and cell function of macrophages, current knowledge on transcriptional regulation of human macrophages is far from complete. To discover novel marker genes, an area of great need particularly in human macrophage biology but also to generate a much more thorough transcriptome of human M1- and M1-like macrophages, we performed RNA sequencing (RNA-seq) of human macrophages. Using this approach we can now provide a high-resolution transcriptome profile of human macrophages under classical (M1-like) and alternative (M2-like) polarization conditions and demonstrate a dynamic range exceeding observations obtained by previous technologies, resulting in a more comprehensive understanding of the transcriptome of human macrophages. Using this approach, we identify important gene clusters so far not appreciated by standard microarray techniques. In addition, we were able to detect differential promoter usage, alternative transcription start sites, and different coding sequences for 57 gene loci in human macrophages. Moreover, this approach led to the identification of novel M1-associated (CD120b, TLR2, SLAMF7) as well as M2-associated (CD1a, CD1b, CD93, CD226) cell surface markers. Taken together, these data support that high-resolution transcriptome profiling of human macrophages by RNA-seq leads to a better understanding of macrophage function and will form the basis for a better characterization of macrophages in human health and disease.


Introduction
Macrophages represent resident phagocytic cells in the tissue and are involved in tissue homeostasis and induction of inflammatory reaction towards pathogens by use of their broad range of pattern-recognition receptors [1]. In context of the respective immune response, macrophages are polarized to specific functional properties, often referred to as M1-like and M2-like phenotype. Human monocytes can be differentiated towards macrophages using either M-CSF or GM-CSF already resulting in differences in gene expression [2]. Classically polarized M1-like macrophages can be induced by IFN-c alone or together with LPS or TNF-a using M-CSF or GM-CSF [3]. M1-like macrophages are effector cells of classical inflammatory immune responses exerting an IL-12 high , IL-23 high and IL-10 low phenotype with secretion of inflammatory cytokines IL-1b, IL-6 and TNF-a. They display a phenotype characterized by the expression of CD86, CD64, and CD16 [4,5]. In contrast, macrophages that are activated by other mechanisms than IFN-c/LPS/TNF-a are grouped in the alternatively activated M2-like macrophage subset. Non-classically activated macrophages can be induced by cytokines including IL-4 and IL-13, but other stimuli have been described as well [4,5]. These cells share an IL-12 low and IL-23 low phenotype and express CD23. Over the last decade, phenotypic adaptations of macrophages to environmental stimuli have been linked to radical changes in transcriptional regulation mainly by applying microarray-based gene expression profiling [3,[6][7][8][9]. In fact, a large amount of data covering transcriptional reprogramming of macrophages has been accumulated, albeit not always systematic [3,[6][7][8][9][10][11]. However, molecular mechanisms controlling transcriptional reprogramming in macrophages are far from understood and it has been suggested that integrative analyses of epigenomic and transcriptomic data will be required to better understand how macrophages integrate the information they receive from their respective microenvironment [12], enabling the identification of specific transcription factor combinations being responsible for cellular macrophage programs.
The introduction of RNA sequencing (RNA-seq) to interrogate whole transcriptomes has challenged previously established gene expression profiling studies [13][14][15]. Advantages assigned to RNAseq over microarray analysis include increases in transcript quantity and quality, improved detection of alternative splicing events and gene fusion transcripts, and a larger dynamic range of detection [13][14][15].
To better understand polarization and integration of environmental signals by macrophages and to identify more specific markers for different functional states, high-resolution transcriptome data have been asked for [16]. Using M1 and M2 polarization as models we applied RNA-seq and compared the information content with data derived by microarray analysis. We provide new insights into human macrophage biology and determine several new markers associated with classical and alternative macrophage polarization in humans.

Ethics Statement
Buffy coats from healthy donors were obtained following protocols accepted by the institutional review board at the University of Bonn (local ethics vote no. 045/09). Informed written consent was provided for each specimen according to the Declaration of Helsinki.

Cell Isolation from Healthy Blood Donors
Peripheral blood mononuclear cells (PBMC) were obtained by Pancoll (PAN-Biotech, Aidenbach, Germany) density centrifugation from buffy coats from healthy donors. CD14 + monocytes were isolated from PBMC using CD14-specific MACS beads (Miltenyi Biotec) according to the manufacturers protocol (routinely .95% purity, Figure S1A).
RNA Isolation 5610 6 -2610 7 macrophages were harvested, subsequently lysed in TRIZOL (Invitrogen) and total RNA was extracted according to the manufactures' protocol. The precipitated RNA was solved in RNAse free water. The quality of the RNA was assessed by measuring the ratio of absorbance at 260 nm and 280 nm using a Nanodrop 2000 Spectrometer (Thermo Scientific) as well as by visualization the integrity of the 28S and 18S band on an agarose gel.
Quantitative PCR Conditions and Primer Sequences 500 ng RNA was reverse transcribed using the Transcriptor First Strand cDNA Synthesis Kit (Roche Diagnostics). qPCR was performed using the LightCyclerTaqman master kit with GAPDH as reference on a LightCycler 480 II (Roche). GAPDH was chosen as reference as both microarray analysis as well as RNA-seq did not show statistically significant differences in gene expression for GAPDH in the conditions assessed. Primer sequences and assay conditions were determined using the Universal Probe Library Assay Design Center (Roche). qPCR primer sequences are summarized in Table S1.
Isoform specific PCR to identify alternative splicing events were performed using the Maxima SYBR Green/Fluorescein qPCR Master Mix (Fermentas). The relative enrichment of each isoform relative to GAPDH was calculated using the 2 2DDCT method. Primer sequences and assay conditions were determined using Beacon Designer 7. qPCR primer sequences are listed in Table  S2.

Microarray-based Transcriptional Profiling and Bioinformatic Analysis of Microarray Data
Prior to array based gene expression profiling total RNA was further purified using the MinElute Reaction Cleanup Kit (Qiagen). Biotin labeled cRNA was generated using the Targe-tAmp Nano-g Biotin-cRNA Labeling Kit for the Illumina System (Epicentre). Biotin labeled cRNA (1.5 mg) was hybridized to Human HT-12V3 Beadchips (Illumina) and scanned on an Illumina HiScanSQ system. Raw intensity data were processed and exported with BeadStudio 3.1.1.0 (Illumina). Subsequent analyses were performed using the R programming language [17] with the Bioconductor software packages [18]. Quality of the array data was assessed using pairwise scatterplots whereby the correlation coefficient should account to r 2 $0.95. Moreover the present call rate was calculated and only samples with at least 30% present calls were included in further studies. Having passed the quality control check points, data was normalized using quantile normalization implemented in the ''limma'' package of Bioconductor. [18]. To remove non-informative genes, the data was filtered using the coefficient of variation. Hence, only genes with a coefficient of variation of at least 0.5 were kept for further analysis. From the resulting data sets we extracted a list of genes with a significant different expression in macrophage subtypes (fold-change $2.0, p-value ,0.05 Student's t-test with Benjamini & Hochberg false-discovery rate correction). Variable genes were plotted as heatmaps with hierarchical clustering using the correlation coefficient as a distance measure for the samples and  the average of each cluster for cluster formation of the genes using the ''amap'' package. Expression values are visualized with colors ranging from red (high expression) over white (indermediate expression) to blue (low expression). Principal component analysis (PCA) was performed using the ''pcurve'' package in R [19]. When visualizing PCA results, the first 2 principal components (coordinates) are shown z-transformed. Microarray data can be accessed under GSE35449.

RNA-seq and Data Analysis
Sequencing and analysis were performed individually on M1like and M2-like macrophages from 3 independent donors. Total RNA was converted into libraries of double stranded cDNA molecules as a template for high throughput sequencing using the Illumina CBot station and HiScanSQ following the manufacturer's recommendations using the Illumina TruSeq RNA Sample Preparation Kit. Shortly, mRNA was purified from 5-10 mg of total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in Illumina proprietary fragmentation buffer. First strand cDNA was synthesized using random oligonucleotides and SuperScript II. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/ polymerase activities and enzymes were removed. After adenylation of 39 ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 200 bp in length the library fragments were separated on a 2% (w/v) agarose gel. The corresponding gel-fraction for each library was excised and purified using the QIAquick gel extraction kit (Qiagen). DNA fragments with ligated adapter molecules were selectively enriched using Illumina PCR primer PE1.0 and PE2.0 in a 15 cycle PCR reaction. Products were purified (QIAquick PCR purification kit) and quantified using the Agilent high sensitivity DNA assay on a Bioanalyzer 2100 system (Agilent). After cluster generation, 100 bp paired-end reads were generated and analyzed using CASAVA 1.8. Alignment to the human reference genome hg19 from UCSC was performed stepwise. First, all reads passing the chastity filter were aligned to the reference genome. Next, reads were aligned to the RNA reference transcriptome. Based on these alignments the numbers of reads aligning to intragenic regions, or intergenic regions, respectively, were calculated. In addition the numbers of reads mapping to exonic and intronic regions as well as to splice sites were calculated based on the UCSC annotation file [20]. Reads per kilobase of exon model per million mapped reads (RPKM) values for Refseq genes were established using CASAVA 1.8. In order to identify reads spanning altered splicing events or gene fusion breakpoints we also analyzed reads using TopHat and Bowtie. Results were further processed using Cufflinks and Cuffdiff [21][22][23][24]. RNA-seq data can be accessed under GSE36952.

A Priori Information-based Network Analysis Using EGAN Software
To visualize connectivity between genes in high-throughput datasets contextual network graphs were generated based on a priori knowledge stored in literature, pathway, interaction, or annotation term databases by EGAN (exploratory gene association network) [25]. To visualize the transcriptional regulation of genes enriched in M1 respectively M2, array data were used and fold change differences calculated using unpolarized macrophages as comparison. Genes with a FC .2 for M1 and FC .1.65 for M2 were visualized; represented is the major network. Using the network topology established for M1-like macrophages the expression values for M2-like macrophages were plotted and vice versa. For comparison of network components and density between RNA-seq and array data, the network was first visualized for the RNA-seq data (FC .4 for M1 and FC .2.5 for M2). Keeping the network topology, genes were marked according to their fold change when visualizing the array-based network. Graphs for genes enriched in M1 respectively in M2 were generated independently.

Immunoblot Analysis
Cell lysates from human macrophages were prepared as previously described [26] followed by immunoblotting with APOL1 or APOL3 antibodies (Abcam) as well as human b-actin (C4; Millipore) as loading control.
(FC $2, p-value ,0.05 with FDR, diff .100) are depicted in red. (E) Left: network of genes highly expressed in M1-like macrophages (fold-change .2.0) in comparison to M0 macrophages identified by microarray analysis. Right: data for the comparison of M2-like versus M0 macrophages were loaded into the M1-network. (F) Right: network of genes highly expressed in M2-like macrophages (fold-change .1.65) in comparison to M0 macrophages identified by microarray analysis. Left: data for the comparison of M1-like versus M0 macrophages were loaded into the M2-network. All networks were generated using EGAN. doi:10.1371/journal.pone.0045466.g002

Statistical Analysis
For flow cytometry or qPCR data comparisons between groups were performed using the appropriate paired or unpaired Student's t-test after testing for equal variance or normal distribution of the data, respectively. P,0.05 was defined as statistically significant. The relationship between fold-changes in gene expression in RNA-seq and microarray analysis was investigated using linear regression. All statistical analyses were performed using the SPSS 19.0 statistical software package.

Generation of Human M1-and M2-like Macrophages as a Model System
To establish a high-resolution transcriptome of human macrophages as a result of specific polarization signals, we used classical (M1-like) and alternative (M2-like) polarization of human macrophages as a model system. Since both M-CSF and GM-CSF have been described to differentiate macrophages from blood-derived CD14 + monocytes with distinct gene expression profiles [2], we first compared the two different stimuli in respect to macrophage polarization and used expression of well-known macrophage markers as the initial readout [3,27]. For classical polarization we primarily used IFN-c as the model stimulus and IL-4 for alternative polarization. When assessing the macrophage surface marker CD11b, the total percentage of CD11b + cells under M1 and M2 polarization conditions was similar while the MFI was slightly higher in M2-like macrophages independent of the usage of GM-CSF or M-CSF (Fig. 1A). Further, we observed high expression of CD14 on all cells under M1 and M2 polarizing conditions irrespective of GM-CSF or M-CSF differentiation with a higher CD14 expression in M1-like macrophages (Fig. 1B). For both classical macrophage markers CD68 and MHC class II molecules ( Fig. 1C and D) we observed no differences in all four conditions tested. Similarly, on whole genome level GM-CSF and M-CSF induced M0 macrophages showed a very similar gene expression profile ( Figure S1A-B). Of note, when the IL-4 signal was provided to monocytes from the beginning of the differentiation period, immature dendritic cells were generated showing a distinct transcriptome ( Figure S2A-B) with a typical loss of macrophage markers such as CD14 or CD68 ( Figure S2C).
When assessing markers previously associated with M1 or M2 polarization [28], a selective induction of the M1 marker CD64 in M1-like macrophages was observed in cultures differentiated with both GM-CSF and M-CSF (Fig. 1E) while CD86 only showed an M1-specific expression in GM-CSF differentiated cells (Fig. 1F). Assessment of these markers following other M1-associated polarization signals, e.g. LPS, TNF-a or combinations thereof resulted in comparable results ( Figure S3). Inversely, strong induction of the M2-marker CD23 was observed in IL-4 polarized macrophages with significantly higher induction in GM-CSF polarized M2-like macrophages (Fig. 1G). As we were mainly interested in macrophage polarization under inflammatory conditions, we chose to differentiate monocytes with GM-CSF for 3 days for further experiments prior to polarization with either IFNc or IL-4 as model signals.

Characterization of M1-and M2-like Macrophages by Microarray-based Gene Expression Profiling
Most recently it has been suggested that assessment of macrophage polarization in humans cannot solely rely on few cell surface markers but should be accommodated by gene expression profiling [16]. Using one of the most recent array generations, gene expression profiling was performed on unpolarized and polarized macrophages derived from seven healthy donors. To determine sample relationships, PCA ( Fig. 2A) and hierarchical clustering (Fig. 2B) based on variable genes were performed and showed segregation of the samples by polarization state. Comparing our data with publically available datasets from M1-and M2-like macrophages generated with earlier array versions we observed concordant gene expression patterns ( Figure  S4) [3]. Heatmap visualization of known M1-and M2-like macrophage markers (Fig. 2C) further demonstrated that the genes encoding for the surface molecules FCGR1A and FCGR1B (both representing CD64) and CD86, the cytokine/chemokine genes CXCL10, CXCL9, IL-1B, IL-6, CXCL11, TNF, IL-23A, and the genes encoding for the intracellular protein GBP1 were increased in M1-like macrophages, results similar to what has been previously reported for M1 polarization [3,29,30]. Inversely, M2-associated genes encoding for the surface molecules FCER2 (CD23), IL27RA, and CLEC4A, the chemokine genes CCL22, CCL18, and CCL17, and the intracellular protein F13A1 were increased in the M2-like macrophages [29][30][31]. Using a foldchange difference of $2 and a p-value ,0.05 to determine differently expressed genes between M1 and M0 and M2 and M0 macrophages respectively, we observed 757 M1 and 436 M2 specific genes ( Figure 2D).
To further illustrate differential macrophage polarization, we generated a priori knowledge based M1-associated (Fig. 2E) and M2-associated (Fig. 2F) networks applying EGAN [25] using genes significantly upregulated in M1-(FC .2) respectively M2polarized cells (FC .1.65). When applying expression values from the comparison of M2-like with M0 macrophages on the M1associated network, the vast majority of genes showed either no change or even a reduction in expression, likely to represent an active repression of M1-associated genes in M2-like macrophages (Fig. 2E). Only few genes showed a simultaneous increase in expression, and these genes represented common cell cycle associated genes. Similarly, members of the M2-associated network were mostly not changed or even reduced in M1-like macrophages (Fig. 2F).

Increase in Overall Transcriptome Information by RNAseq
Gene expression profiling using microarrays has recently been suggested to be inferior to newer sequencing based technologies in providing genome-wide transcriptome information [14]. To address the potential information increase for macrophages, RNA-seq was performed on in vitro generated M1-and M2-like macrophages. After quality filtering, we obtained 15.062.8 million and 20.469.2 million reads for the M1-and M2-like macrophage cDNA libraries (Table 1). Consistent with RNA-seq data obtained from other eukaryotic cells [32] the majority of sequencing reads for M1-and M2-like macrophages mapped to annotated exons (Refseq transcripts). The remaining reads mapped to exon-intron boundaries, introns, or other uncharacterized intergenic regions (Table 1). RPKMs are measures of individual transcript abundance in RNA-seq datasets and have been shown to be highly accurate across multiple cell types [32]. We used CASAVA to assign RPKMs. To compare RNA-seq and microarray data we cross-annotated RNA-seq and microarray data using HGNC symbols. In human M1-and M2-like macrophages, 11317 and 11466 Refseq genes were expressed applying a previously defined optimal threshold (0.3 RPKM) for gene expression (Fig. 3A) [33]. The present call rate for Refseq genes for M1-(n = 10155) and M2-like macrophages (n = 10418) was only slightly lower when using microarrays (Fig. 3A). Furthermore, when assessing the levels of mRNA expression on a global scale a high correlation between RNA-seq and microarray data -similar to other cells [32] -was observed for M1-and M2-like macrophages (Fig. 3B).

RNA-seq Reveals Differential Expression at Significantly Higher Resolution
RNA-seq showed a significantly increased dynamic range over background mainly due to significantly lower background levels. This suggested that the assessment of differential expression using RNA-seq might lead to an improved resolution in comparison to array-based data. Applying standard filter criteria (FC $1.5, p,0.05, RPKM $0.3) revealed a total of 1736 genes elevated in M1-versus M2-like macrophages by RNA-seq, while 834 genes were recognized by array analysis (Fig. 3C). Similarly, 822 genes were identified as being elevated in M2-versus M1-like macrophages by RNA-seq, while 786 genes were detected by array analysis (Fig. 3D). More importantly, when categorizing differentially expressed genes according to their level of differential expression, RNA-seq data clearly revealed up to 4-fold more genes with FC .4 for M1-and M2-associated genes ( Fig. 3C and D), which was similarly true for M1-associated genes at lower levels (1.5, FC ,4). To reveal potential reasons for this difference on the single-gene level we utilized FC-FC plotting, correlating RNAseq and array-based data for individual genes (Fig. 3E). The majority of genes showed similar behavior in both RNA-seq and microarray experiments, albeit the relative differences were higher in RNA-seq data (Fig. 3E). Altogether, we observed a considerable increase in the dynamic range of fold-differences in RNA-seq data with differences spanning six orders of magnitude in contrast to only four orders of magnitude in the microarray data ( Fig. 3E and Figure S5). In addition, there was a significant number of genes showing differential expression in RNA-seq data (e.g. DUOX1, DUOX2, TBX21, GBP7) but not in the array data suggesting that the array data are not informative for this class of genes. As anticipated, when using Venn diagrams with a defined cutoff (-2$ FC $2, p,0.05, RPKM $0.3) for data presentation (Fig. 3F), both RNA-seq and array data revealed 595 genes to be differentially expressed, but RNA-seq revealed 900 additional genes of which 3 we validated on protein level ( Figure S6). Surprisingly, 308 genes were classified as being differentially expressed by array analysis alone (Fig. 3F). When further assessing these genes, it became apparent that these genes show a similar trend in the RNA-seq data but these differences did not yet reach statistical significance (Fig. 3G). In contrast, genes only identified by RNA-seq, clearly showed no differences when assessed by array analysis (Fig. 3H). Visualization of typical marker genes as depicted for array data in Fig. 2C demonstrated a comparable differentiation of M1-and M2-like macrophages when assessed by RNA-seq (Fig. 3I).

Exon Resolution Transcriptome Analysis of Known Macrophage Markers
Another advantage of RNA-seq is to resolve gene expression on the exon level (Fig. 4). For the macrophage related genes CD68 (Fig. 4A-D), CD64 (Fig. 4E-H) and CD23 (Fig. 4I-L), RNA-seq data were visualized for the genomic loci of the respective genes and compared with array, qPCR, and FACS data. For CD68, RNA-seq data revealed similarly high expression for M1 and M2 macrophages for all exons with little variance in expression levels between donors ( Fig. 4A and Figure  S7A). Slightly higher variance was observed for both microarray (Fig. 4B) and qPCR data (Fig. 4C) while protein levels showed equal expression in all samples analyzed (Fig. 4D). For CD64, RNA-seq revealed complete absence of expression for all exons in M2-like macrophages with high expression in M1-like macrophages ( Fig. 4E and Figure S7B), which was similarly observed by other technologies (Fig. 4F-H). For CD23, protein data suggest significantly elevated expression on M2-like macrophages with low level expression on M1-like macrophages (Fig. 4L), a result which was also observed for RNA-seq data ( Fig. 4I and Figure S7C) as well as array (Fig. 4J) and qPCR (Fig. 4K). Similar results were obtained for other marker genes (data not shown). Additionally, we were able to detect classical M1/M2-markers not accessible using microarrays ( Figure S8), suggesting that high-resolution RNA-seq data are predestined for exploration of genes not detectable using microarrays, novel marker genes, as well as biological principles of macrophage polarization.

RNA-seq Ameliorates Network-based Analysis in M1-and M2-like Macrophages
To understand if RNA-seq would also enhance the understanding of biological principles of macrophage polarization we applied network analysis based on a priori information assessing the information content of RNA-seq data in comparison to array data. Genes expressed at elevated levels in M1 RNA-seq data (FC .4) were used for network generation (Fig. 5A). This primary RNA-seq based M1 network was subsequently used to visualize array-based gene expression (Fig. 5B). When genes at a lower level of differential expression (FC .2) were included 73% of the network was revealed in the array data and central hubs of the network were also categorized as being highly (FC .4) differentially expressed. However, only RNA-seq data revealed two gene clusters of immunomodulating proteins highly enriched in the M1 network, namely apolipoproteins L (APOL) (Fig. 5A and Figure S9) and the leukocyte immunoglobulin-like receptor (LILR) family ( Fig. 5A and Figure  S10) [34][35][36]. As exemplified for LILRB1 and APOL3 both genes were clearly identified by RNA-seq, qRT-PCR, and flow cytometry respective western blotting ( Fig. 5E and F) but not by microarray analysis (data not shown). Applying the RNA-seq data-based M2 network (Fig. 5C) to the array data (Fig. 5D) revealed only 54% elevated genes and major network hubs were

Identification of Splice Variants and RNA Chimaera in Differentially Stimulated Human Macrophages
It has recently been suggested that cell differentiation can result in usage of alternative gene transcripts or isoform switching [21]. We applied Cufflinks and Cuffdiff to illuminate switches in dominant promoter usage, transcription start sites (TSS), and coding sequences (CDS) [21]. This analysis revealed 9 genes with alternative promoters (Table S3), 28 genes using alternative TSS (Table S4), and 20 genes with different CDS in M1-and M2-like macrophages (Table S5). We analyzed one of these genes in greater detail. For the gene encoding PDZ and LIM domain 7 (enigma) (PDLIM7) we observed a slight but significant increase in M1-like macrophages for the complete locus in RNA-seq data (Fig. 6A) while the probes on the microarray revealed no significant changes (Fig. 6B). Previous screening projects suggested three different CDS for PDLIM7. Applying Cufflinks and Cuffdiff to M1 and M2 RNA-seq data clearly revealed differential expression of individual CDS (Fig. 6C). While PDLIM7 v1 was mainly expressed by M1-like macrophages, M2-like macrophages mainly expressed PDLIM7 v2 while no difference was observed for PDLIM7 v4. We validated the usage of these variants by version-specific qPCR (Fig. 6D). Taken together, these new findings might open new avenues towards the role of alternative splicing in macrophages potentially linking alternative transcript usage with macrophage polarization.

New Markers for M1-and M2-like Macrophages Identified by Combined Transcriptome Analysis
In light of the still limited number of cell surface markers clearly distinguishing human M1-from M2-like macrophages, we interrogated the genes of the human surfaceome [37] for differential expression between M1-and M2-like macrophages. By this approach 475 surface molecules were found to be differentially expressed (Fig. 7A). As visualized in Fig. 7B, the cell surface molecules CD120b, TLR2, and SLAMF7 showed preferential expression in M1-like macrophages, which was true irrespective of differentiation of macrophages by GM-CSF or M-CSF ( Figure  S11A). Several surface molecules, including CD1a, CD1b, CD93 and CD226 were significantly increased in expression in M2-like macrophages ( Fig. 7C and Figure S11B). Taken together, screening higher-resolution transcriptome data established by RNA-seq allows for the identification of novel genes related to specific polarization programs in macrophages.

Discussion
Because of the enormous plasticity of human macrophages, the classification of polarization states on the basis of few cell surface markers will remain a substantial challenge [16]. Here, we addressed how RNA-seq based high-resolution transcriptome data can be utilized to better understand the biology of macrophage polarization. We observed a significant increase in dynamic range in RNA-seq data resulting in a significantly higher number of genes determined to be significantly differentially expressed. This was true despite the fact that we used seven biological replicates for array analysis but only three samples for RNA-seq. A priori information based network analysis further supported that the increased information content of RNA-seq data uncovered novel aspects of macrophage biology, which was illustrated by the recognition of differential expression of numerous family members of two gene families, namely the apolipoprotein L family and leukocyte immunoglobulin-like receptors. APOLs constitute a new class of apolipoproteins expressed by macrophages as they serve as lytic factors against invading pathogens, e.g. African trypanosomes inducing programmed cell death as well as inhibiting intracellular infection by Leishmania [34,35]. LILRs have been associated with balancing the effects of Toll-like receptor signaling, suggesting an important role of LILRs both in the initiation but also cessation of inflammatory responses mediated by macrophages [36]. Another aspect enhancing our knowledge about the polarization biology of macrophages was the identification of several genes with differential usage of alternative promoters and transcription start sites as well as differential splicing variants between M1-and M2like macrophages. As visualized for PDLIM7, an intracellular scaffold protein that contains a PDZ domain and three LIM domains linked to mitogenic signaling through actin cytoskeleton organization [38], regulating Tbx5 transcriptional activity [39], and suppressing p53 activity [40], RNA-seq revealed significant differences in splice variant usage for M1-and M2-like macrophages potentially linking p53 regulation with macrophage polarization [41]. Usage of splice variant-specific qPCR reactions supported these findings while this differential regulation was not revealed by microarray analysis. Altogether we detected differential promoter usage, transcription start site usage and splice variant usage in over 50 gene loci, a number that was surprisingly low taking into account that such mechanisms of transcriptional regulation have been suggested for the majority of gene loci in mammalian genomes [42].
While studies in other cell systems suggested that RNA-seq data will further improve cell characterization [13][14][15], the direct assessment of the new technology in macrophage polarization was necessary to estimate its potential information gain. Both, increased dynamic range and the identification of transcripts that were missed by microarray analysis were major reasons for the discovery of novel genes associated with either M1-or M2polarization. Nevertheless, despite a lower number of informative transcripts in the microarray data, 73% of the major M1-network was still revealed -at least when using transcripts defined to be enriched in M1-like macrophages. However, this rate dropped to only 54% in the M2-network and major hubs like MYC and TP53 where only revealed by RNA-seq data in M2-like macrophages. Overall these findings point towards an advantage of RNA-seq data, when the endpoint of the analysis is the identification of novel biological mechanisms.
An important aspect of genomic characterization is the identification of novel marker genes in macrophage polarization [16]. When focusing on genes being part of the human surfaceome in most cases RNA-seq data revealed larger differences between M1-like and M2-like cells when compared to microarray data. Nevertheless, some genes only reached significant differential expression in the array data clearly pointing toward the necessity to include a large enough number of biological replicates also like macrophages generated in the presence of GM-CSF with quantification shown in the graph at the right. Expression of (B) CD120b, TLR2, and SLAM7 as well as (C) CD1a, CD1b, CD93, and CD226. Isotype controls are depicted as dotted lines. *P,0.05 (Student's t-test). Numbers in plots indicate mean fluorescence intensity. Data are representative of nine independent experiments (B,C; mean and s.e.m.) each with cells derived from a different donor. doi:10.1371/journal.pone.0045466.g007 when applying RNA-seq. On the other hand, a subset of genes showed the well-known background noise effect in the microarray data resulting in non-significant differences between the two cell types. Irrespective of these different shortcomings of the two technologies, the overall differences between the two techniques in this defined gene space were less obvious suggesting that both technologies are well suited for cell surface marker identification. Taken together, we introduced several new marker genes for which we established FACS assays that can be used to distinguish between M1 and M2 polarization of macrophages and that can be combined with the analysis of common macrophage markers.
The identification of novel marker genes distinguishing human M1-like and M2-like macrophages opens new avenues towards understanding the biology of differentially polarized macrophages. One of the M1-marker identified in this study, namely CD120b (TNFR2) has been linked to cell survival, activation and even proliferation in other cell types such as T cells [43]. In contrast to TNFR1, TNFR2 preferentially leads to NFkB activation. Whether this is true in myeloid cells as well requires further investigation. However, earlier studies already suggested that production of TNF-a in macrophages might be interpreted as a phenotypestabilizing feed-forward loop [44] and TNFR2 might actually play an important role in such a process.
SLAMF7 was originally identified as a NK cell-associated surface molecule [45]. Subsequently, it was shown to be expressed on lymphocytes and monocytes [46]. More recently, a reduced expression on monocytes and NK cells with a simultaneous increase of SLAMF7 on B cells was observed in patients with lupus erythematosus [47]. The strongest link to SLAMF7 as an M1 marker gene comes from observations in intestine allograft rejection, demonstrating that tissue macrophages derived from patients rejecting the graft showed elevated levels of SLAMF7 [48]. It would be interesting to see if macrophages in other settings of transplant rejection are also enriched for this novel M1 marker gene. Considering the identification of single specific marker genes for macrophage polarization our findings clearly point to the necessity for multi-parameter analysis instead. This can be exemplified by the differential expression of CD1a and CD1b, two cell surface molecules that are mainly studied in context of antigen presentation by dendritic cells [49]. Previous reports suggested upregulation of CD1 proteins on human monocytes by GM-CSF [50]. However, we clearly present evidence that expression is induced in both M-CSF and GM-CSF driven macrophages and polarization towards M2-like macrophages is significantly increasing expression of CD1a and CD1b suggesting that they might be up-regulated on tissue macrophages in an M2driving environment. This is similarly true for CD93, which was originally identified to be expressed on early hematopoietic stem cells and B cells [51]. CD93 is involved in biological processes such as adhesion, migration, and phagocytosis [52,53]. CD93 expressed on myeloid cells can be shed from the cell surface and the soluble form seems to be involved in differentiation of monocytes towards a macrophage phenotype [54]. Since soluble CD93 has been implicated in inflammatory responses, it will be important to further elucidate how polarization-induced differential expression of CD93 contributes to specific inflammatory responses. Another surprising finding is the differential expression of CD226 between human M1-and M2-like macrophages, a molecule initially shown to be involved in cytolytic function of T cells [55]. Subsequently, it could be shown that CD226 has additional functions including the regulation of monocyte migration through endothelial junctions [56]. Similar to the other M2-associated markers, so far little is known about CD226 on polarized macrophages. Since CD226 expression levels on lymphocytes have been implicated in autoimmune diseases [57] further research is necessary to understand its role in the myeloid compartment during such processes.
Overall, by using RNA-seq we introduce a high-resolution transcriptome analysis of human macrophages unraveling novel insights into macrophage polarization. While previously established transcriptome datasets addressing macrophage biology are still very suitable to assess important biological and medical questions, a deeper understanding of transcriptional regulation during macrophage polarization will require higher resolution that is provided by current and future RNA-seq technologies. Moreover, the novel cell surface markers will help to better understand macrophage programs and functions in human disease.