|Home | About | Journals | Submit | Contact Us | Français|
Sustained and dysfunctional macrophage activation promotes inflammatory cardiometabolic disorders, but the role of long intergenic noncoding RNA (lincRNA) in human macrophage activation and cardiometabolic disorders is poorly defined. Through transcriptomics, bioinformatics, and selective functional studies, we sought to elucidate the lincRNA landscape of human macrophages.
We used deep RNA sequencing to assemble the lincRNA transcriptome of human monocyte‐derived macrophages at rest and following stimulation with lipopolysaccharide and IFN‐γ (interferon γ) for M1 activation and IL‐4 (interleukin 4) for M2 activation. Through de novo assembly, we identified 2766 macrophage lincRNAs, including 861 that were previously unannotated. The majority (≈85%) was nonsyntenic or was syntenic but not annotated as expressed in mouse. Many macrophage lincRNAs demonstrated tissue‐enriched transcription patterns (21.5%) and enhancer‐like chromatin signatures (60.9%). Macrophage activation, particularly to the M1 phenotype, markedly altered the lincRNA expression profiles, revealing 96 lincRNAs differentially expressed, suggesting potential roles in regulating macrophage inflammatory functions. A subset of lincRNAs overlapped genomewide association study loci for cardiometabolic disorders. MacORIS (macrophage‐enriched obesity‐associated lincRNA serving as a repressor of IFN‐γ signaling), a macrophage‐enriched lincRNA not expressed in mouse macrophages, harbors variants associated with central obesity. Knockdown of MacORIS , which is located in the cytoplasm, enhanced IFN‐γ–induced JAK2 (Janus kinase 2) and STAT1 (signal transducer and activator of transcription 1) phosphorylation in THP‐1 macrophages, suggesting a potential role as a repressor of IFN‐γ signaling. Induced pluripotent stem cell–derived macrophages recapitulated the lincRNA transcriptome of human monocyte‐derived macrophages and provided a high‐fidelity model with which to study lincRNAs in human macrophage biology, particularly those not conserved in mouse.
High‐resolution transcriptomics identified lincRNAs that form part of the coordinated response during macrophage activation, including specific macrophage lincRNAs associated with human cardiometabolic disorders that modulate macrophage inflammatory functions.
Long noncoding RNAs (lncRNAs), defined as >200 nucleotides in length and often 3′ polyadenylated,1 are increasingly implicated in cardiovascular diseases.2 Compared with mRNAs, lncRNAs are less abundant, mostly spliced but with fewer exons, and more species‐ and tissue‐specific,1, 3 emphasizing the importance of studies on human‐ and cell‐specific lncRNAs in human physiology and diseases.
As a critical component of the innate immune system, macrophages demonstrate remarkable plasticity and wide‐ranging states of activation. Sustained and dysfunctional macrophage activation promotes inflammatory cardiometabolic diseases (CMDs) such as atherosclerosis and metabolic dysregulation.4, 5 Macrophage activation to M1 (classic inflammatory activation by lipopolysaccharide and IFN‐γ [interferon γ]) and M2 (alternatively activated by IL‐4 [interleukin 4]) phenotypes are well characterized in vitro models for study of human and murine macrophage biology.6, 7 Although the protein‐coding transcriptome of human macrophages has been well characterized,7, 8, 9 the lncRNA landscape in human macrophage biology remains elusive. Long intergenic noncoding RNAs (lincRNAs) lincRNA‐Cox210 and lincRNA‐EPS11 have been shown to be critical regulators of inflammation in murine macrophages, but both lack human orthologs, limiting translational relevance to human. A handful of studies have mapped human macrophage lncRNAs, most using microarray12, 13, 14 and THP‐1 monocyte‐derived macrophages (THP‐1Φ),12, 13 yet THP‐1Φ is karyotypically abnormal and immature and thus may differ from primary human macrophages.15 Recent FANTOM5 (functional annotation of the mammalian genome) cap analysis of gene expression data sets have profiled transcription start site (TSS) and enhancer elements of lncRNAs in human monocyte‐derived macrophages,16, 17 yet deep RNA sequencing (RNA‐seq) of human macrophages has been lacking and is required to provide genomewide assembly of lncRNAs and to facilitate prioritization of promising lncRNAs for functional validation.
We have previously generated deep RNA‐seq data sets of human peripheral blood mononuclear cell–derived macrophages (HMDMs)18 with thorough characterization of their coding transcriptome18 and alternative splicing events during M1 and M2 activation.19 In this study, we focused on lincRNAs, a major subset of lncRNAs, using the same RNA‐seq data sets.18 LincRNAs do not overlap annotated protein‐coding regions, facilitating experimental validation.3 Because most genetic signals for complex traits are in intergenic regions, functional genetic variation in lincRNAs are likely to contribute to the intergenic genomewide association study (GWAS) signals for complex traits.20 Through de novo transcriptome assembly, we (1) report a comprehensive lincRNA catalog (31% are newly annotated); (2) identify specific lincRNA expression patterns that correspond to distinct M1‐ and M2‐activated phenotypes; (3) stratify macrophage lincRNAs based on synteny, conservation, tissue enrichment, and regulatory features defined by essential macrophage transcription factors (TFs) as well as histone H3 lysine 4 monomethylation (H3K4me1) and histone H3 lysine 4 trimethylation (H3K4me3) chromatin immunoprecipitation sequencing profiles21, 22; (4) use GWAS data to identify macrophage lincRNAs related to human complex CMDs; (5) perform initial functional validation of MacORIS, a lincRNA that harbors single nucleotide polymorphisms (SNPs) associated with central obesity; and (6) characterize human induced pluripotent stem cell–derived macrophages (IPSDMs) as a model for functional assessment of human lincRNAs in macrophage biology. Our findings constitute a unique translational proof of principle and resource for the comprehensive interrogation of human macrophage lincRNAs in macrophage differentiation, inflammatory and metabolic functions, and relationship to human CMDs.
All human protocols for this work were approved by the University of Pennsylvania and Columbia University Medical Center human subjects research institutional review boards, and all participants provided written informed consent.
Human macrophage preparation, RNA‐seq library preparation, sequencing, and data analysis were described previously18 and in Data S1. Briefly, human peripheral blood mononuclear cells were differentiated to macrophages using 100 ng/mL macrophage colony‐stimulating factor, and activation was induced by 18‐ to 20‐hour incubation with 100 ng/mL lipopolysaccharide and 20 ng/mL IFN‐γ for M1‐like activation or 20 ng/mL IL‐4 for M2‐like activation.18 Strand‐specific, poly(A)+ libraries underwent deep sequencing at 100‐bp paired‐end reads to obtain in macrophages ≈130 million filtered reads per sample with >95% mapping rate18 and in monocytes ≈280 million filtered reads per sample with >93% mapping rate.23 RNA‐seq reads were aligned with the hg19 reference genome, and transcript abundance was measured in FPKM (fragments per kilobase of transcript per million fragments mapped) using Cufflinks 126.96.36.199 De novo assembly was performed using Cufflinks 188.8.131.52 Differential expression, defined as false discovery rate–adjusted (FDR‐adjusted) P<0.01 and a fold change >2, was tested with Cuffdiff. RNA‐seq data are available under the accession number GSE55536.18 (Table S1 shows participant demographics, and Figure S1 shows correlation between biological replicates). The bioinformatics pipeline for the annotation of the human macrophage lincRNA catalog, including the long intergenic transcript filters, coding potential filters and reliable expression filters, is outlined in Figure 1 and described in detail in Data S1.
Synteny is defined as conserved gene order along the chromosomes of different species. We examined the synteny of macrophage lincRNAs in mouse using HomoloGene release 68, then further subdivided syntenic lincRNAs as annotated or not annotated in mouse, using GENCODE M4 annotation. To evaluate sequence conservation for syntenic lincRNAs, the human lincRNA sequence was queried against the mouse genome with an E‐value cutoff of 1×10−10 using BLASTN. Sequence hits in the mouse within the syntenic region were then searched in human samples with the same E‐value cutoff. Sequences that passed the reciprocal steps were considered conserved.
The fractional expression value of each lincRNA and mRNA was calculated by dividing the FPKM value in HMDMs by total FPKM values across HMDMs and 16 tissues from Human BodyMap RNA‐seq data sets3 (GSE30611).25
Histone H3 lysine 4 monomethylation (H3K4me1),21, 22 histone H3 lysine 4 trimethylation (H3K4me3),21 histone H3 lysine 27 acetylation,22 and PU.1 and C/EBPβ22 chromatin immunoprecipitation sequencing data sets for human HMDM were downloaded from GSE31621 22 and GSE58310.21 Data were quantified using computeMatrix from deepTools v184.108.40.206
We first explored the overlap of macrophage‐enriched lincRNAs with trait‐associated SNPs that reached a significance level of P<1×10−5 using data from the comprehensive NHGRI (National Human Genome Research Institute) GWAS Catalog.27 To further interrogate SNPs within macrophage lincRNAs for specific association with the 13 cardiometabolic traits (Table S2), 63 586 genotyped and imputed (HAPMAP28) SNPs were mapped to macrophage lincRNAs (±1 kb) and interrogated using either the minimum P value for the corresponding SNPs within each lincRNA (Bonferroni‐adjusted threshold of P<0.05)29 or a class‐based method, GenCAT (Genetic Class Association Testing),30 to test the overall impact of all SNPs within the region. Significantly associated lincRNAs were further prioritized to include only those that contained the strongest SNP‐level P value in the region (±500 kb of the lincRNA) or if it was in low‐linkage disequilibrium (r 2<0.3) with a stronger single SNP in the region, suggesting an independent signal at the lincRNA locus.
Quantitative reverse transcriptase–polymerase chain reaction (qRT‐PCR) was used to validate lincRNA expression, and primers are listed in Table S3. Knockdown of a top GWAS‐associated lincRNA, MacORIS, was performed in THP‐1Φ by transfection of single‐stranded antisense oligonucleotides (Exiqon) and small interfering RNA (Dharmacon) using Lipofectamine RNAiMAX Transfection Reagent (Thermo Fisher Scientific). Western blotting was used to determine expression and phosphorylation of JAK2 (Janus kinase 2; Try1008) and STAT1 (signal transducer and activator of transcription 1; Tyr701). Flow cytometry was used to determine the expression of IFNGR1 (IFN‐γ receptor 1).
Specific analyses of RNA‐seq and genomic data are described within each section. For analysis of gene ontology pathways in RNA‐seq data, significant enrichment was declared at FDR‐adjusted P<0.05 using the Benjamini and Hochberg method.29 Nonsequencing data were analyzed with GraphPad Prism 6 (GraphPad Software). Differences between 2 groups were assessed by Student t tests (2‐tailed). One‐way ANOVA followed by the Dunnett test was used to correct for multiple comparisons. Results were declared significant if P<0.05.
We interrogated a stringent set of known multiexon human lincRNAs (>200 bp, no overlap with a protein‐coding gene within ±1 kb of lincRNAs) collated from (1) the “Cabili” Human BodyMap “stringent” data set (4273 lincRNAs)3 and (2) the GENCODE V19 data set (7114 lincRNAs; Figure 1).31 By combining the 2 sets, a catalog of 8045 known multiexon lincRNAs was generated. Next, we performed de novo transcriptome assembly by Cufflinks v2.1.1 and excluded previously annotated multiexon lincRNAs, as described in Data S1. We filtered single‐exon lincRNAs because of greater probability of transcriptional noise.32 We applied a coding potential filter on newly annotated lincRNAs using iSeeRNA33 and HMMER‐3 on Pfam34 (Figure 1). To identify lincRNAs that were robustly and reliably expressed in human macrophages, we included only lincRNAs expressed at >1% FPKM based on FPKMs for all lincRNAs and mRNAs in at least 50% of subjects in all HMDM samples. Through this conservative, multilayered analysis, we identified 2766 distinct multiexon lincRNAs that are reliably expressed in human macrophages, of which 861 were previously unannotated (Figure 1). Coding potential for all 2766 lincRNAs was further validated by PhyloCSF,35 as described in Data S1, and the lincRNAs with scores higher than the threshold cutoff were listed in Table S4.
Among the 2766 lincRNAs, 1282 lincRNAs were found in all 6 M0, M1, or M2 HMDMs, and 562 lincRNAs were expressed in all 18 unique macrophage samples (Table S5). More than 50% (1407) of lincRNAs were found across all M0‐, M1‐, and M2‐HMDM activation states, whereas a smaller portion of lincRNAs were highly specific to M0, M1, or M2 HMDMs (Figure S2A). These latter “activation state”–specific lincRNAs were more likely to be previously unannotated lincRNAs (Figure S2B and S2C), underscoring the importance of interrogating lincRNAs within cell‐specific and functional contexts.
As in other cell types,3 compared with protein‐coding genes, macrophage lincRNAs were generally expressed with less abundance, were shorter, and had fewer exons than mRNAs (Figure 2A and Figure S3). Expression levels of newly annotated macrophage lincRNAs were lower than those of known lincRNAs (Kolmogorov–Smirnov test, P<2.2×10−16; Figure 2A).
Many functional lincRNAs are suggested to have synteny (genomic regions flanked by homologous protein‐coding genes)1 despite low sequence similarity across species.36 Of the 2766 macrophage‐expressed lincRNAs, 61% (1678) were syntenic, yet only 400 of these 1678 lincRNAs were annotated as expressed in mouse GENCODE M4. At this time, ≈85% (2366) of human macrophage lincRNAs are nonsyntenic or are syntenic but not expressed in mouse (Figure 2B). Among the 1678 syntenic lincRNAs, only 24% (395) demonstrated significant sequence conservation (Figure 2C) between human and mouse.37 Compared with previously known lincRNAs, a lower proportion of newly annotated macrophage lincRNAs was conserved (28% versus 15%, respectively; P=3.49×10−9), suggesting a higher level of species specificity (Figure 2C). Although lincRNAs that were syntenic and annotated as expressed in mouse showed higher expression in human macrophages than other lincRNAs (Figure 2B), sequence conservation per se was not associated with expression levels in human macrophage (Figure 2C).
Enrichment of lincRNAs in macrophages, relative to other tissues, may suggest their specific roles in macrophage biology. Consequently, we determined the tissue enrichment of macrophage lincRNAs by calculating their expression in macrophages relative to the sum of expression across 16 tissues in data from the Human BodyMap RNA‐seq.3 Applying a fractional expression of >0.2 to define “enriched” lincRNAs,3 595 lincRNAs within the 2766 macrophage lincRNAs (21.5%) demonstrated enriched expression in HMDMs (Table S5). Relative to protein‐coding genes, macrophage lincRNAs were proportionally more macrophage enriched (eg, 15.3% versus 9.8% in M0 HMDM; Figure 3A). Certain lincRNAs were specifically enriched only in M1 or M2 HMDM (Figure 3B, Table S5). Expression levels of macrophage‐enriched lincRNAs were higher than those of nonenriched lincRNAs (Figure 3C).
PU.1 and C/EBPβ are essential TF regulators of macrophage differentiation.22 Leveraging public data sets (GSE31621),22 we discovered that TF occupancy was significantly higher around (±2 kb) lincRNA TSS and gene bodies for macrophage‐enriched lincRNAs than for non–macrophage‐enriched lincRNAs; indeed, the majority of the enriched lincRNAs demonstrated PU.1 or C/EBPβ binding (Figure 3D). In comparing M0 HMDM to our human monocyte data (RNA‐seq of 6 age‐ and race‐matched subjects23; Table S1), numerous lincRNAs were differentially expressed (DE; fold change >2 and FDR <0.01) during monocyte transition to HMDM (Figure S4A and S4B and Table S6); many of the upregulated lincRNAs harbored PU.1 and C/EBPβ binding sites, and most (72 of 114) were also macrophage enriched (Figure S4C and S4D). This highlights the potential roles of a subset of highly macrophage‐specific lincRNAs in macrophage maturation and function.
Regulatory features at lincRNA loci increase the likelihood of biological and functional roles. Many macrophage lincRNAs overlap macrophage enhancer marks. Using public human macrophage chromatin immunoprecipitation sequencing data sets,21 the majority of protein‐coding genes (15 201 genes) displayed punctate binding of the H3K4me3 promoter mark around the TSS (±1.5 kb). In contrast, lincRNA TSS intervals showed relatively weaker signals, and only a small subset displayed high H3K4me3 density. In contrast, binding of H3K4me1, an enhancer mark, at lincRNAs was greater than at protein‐coding genes (Figure 4A). Transcription from putative enhancer regions characterized by high levels of H3K4me1 relative to the H3K4me3 is a major feature of enhancer RNAs,38 and lncRNAs that act as enhancer RNAs have been shown to modulate monocyte immune response.39 For macrophage lincRNAs, we used H3K4me1/H3K4me3 ratios of >1.2 and <0.8 to define enhancer and promoter states,39 respectively. In contrast to the predominance of promoter features at mRNAs (Figure 4A and Table S5), the majority of macrophage lincRNAs (60.9%, 994 of a total of 1632 M0‐HMDM lincRNAs with H3K4me1 and H3K4me3 signals) exhibited an H3K4me1/H3K4me3 ratio of >1.2, suggestive of “elincRNAs” with potential enhancer functions based on bioinformatics prediction (Table S5). Furthermore, elincRNAs were more macrophage enriched but had lower expression levels compared with promoter‐associated lincRNAs, or “plincRNAs.” Consistent with previous findings in mouse erythroblasts,40 the polyadenylated macrophage elincRNAs in our catalog were more likely to have unidirectional transcription (Figure S5A through S5C), unlike most nonpolyadenylated bidirectional enhancer RNAs.38
We hypothesized that elincRNAs that overlap enhancer signatures (H3K4me1/H3K4me3 >1.2), would have greater coregulation with nearby protein‐coding genes during M1 activation than lincRNAs with classical promoter signatures (H3K4me1/H3K4me3 <0.8). The correlation of fold change in expression during M1 activation (M0 versus M1) of elincRNAs and their closest protein coding genes was significantly stronger than that for plincRNAs (Z test, P=0.0039; Figure 4B). Gene ontology analysis of the 194 protein‐coding genes nearest the M1‐induced elincRNAs revealed enrichment for mRNAs involved in immune response and response to bacteria as well as transcriptional regulation (Figure S5D and Table S7). Ingenuity pathway analysis also suggested M1 activation related signaling and biological functions (Figure S5E and S5F and Tables S8 and S9). Indeed, a number of elincRNAs were paired with neighboring protein‐coding genes known to have roles in the macrophage inflammatory response, including Acyl‐CoA Synthetase Long‐Chain Family Member 1 (proximal to RP11‐701P16.5), Interleukin 6 (proximal to AC002480.2), C‐C Motif Chemokine Ligand 8 (proximal to CUFF.135177), Interferon Induced Transmembrane Protein 3 (proximal to RP11‐326C3.12), and Heparan Sulfate‐Glucosamine 3‐Sulfotransferase 3B1 (proximal to AC022816.2), or novel GWAS genes, such as Transmembrane Protein 171 (proximal to linc‐ZNF366‐6) associated with serum urate and gout,41 and Chromosome 6 Open Reading Frame 223 (proximal to CUFF.59743 and CUFF.59265), associated with circulating levels of vascular endothelial growth factor42 (see Table S10 for complete list). The genomic structure and the potential cis‐regulatory effects of these elincRNAs on their coding gene partners require further functional investigation.
Thus, macrophage enrichment and TF binding patterns facilitate prioritization of macrophage lincRNAs that are more likely to be functional, whereas enhancer features suggest mechanistic avenues to pursue in translation.
Macrophage activation induces widespread change in the protein‐coding gene transcriptome,18 but lincRNA modulation during macrophage activation is largely unexplored. During M1 activation, 96 lincRNAs were DE (fold change >2, FDR‐adjusted P<0.01), with 73 up‐ and 23 downregulated, of which 22 were newly annotated (Figure 5A and Table S11). In contrast to M1 activation, only 5 lincRNAs were DE (all upregulated) during M2 activation (Figure 5B; Table S12), consistent with the modest difference in mRNAs between IL‐4–derived M2 HMDM and their macrophage colony‐stimulating factor–differentiated M0 HMDM.18 Relative to all macrophage lincRNAs, those that were DE during macrophage activation were more likely to be syntenic and annotated as expressed in mouse (15% versus 27%; P=0.0013), suggesting that synteny may be a feature of some physiologically relevant macrophage lincRNAs. In contrast, primary sequence conservation with mouse was not associated with DE lincRNAs in macrophage activation (≈25% in both groups). We focused further on lincRNAs that were upregulated in M1 or M2 activation and also macrophage‐enriched, because mRNAs with such features have been shown to contain many important protein‐coding genes with functional roles in macrophage biology.18, 43 Indeed, relative to noninduced lincRNAs, activation‐induced lincRNAs (73 M1‐induced and 5 M2‐induced DE lincRNAs) were more likely to be macrophage‐enriched (55 of 78, P<2.2×10−16) and to overlap PU.1 and C/EBPβ binding (P=3.94×10−3 and P=3.39×10−3, respectively, in TSS and P=9.04×10−5 and P=2.96×10−5, respectively, in gene body). These data highlight a promising subset of human macrophage lincRNAs for follow‐up.
Based on abundance, extent of induction, tissue enrichment, and TF binding, we selected 10 lincRNAs (8 most upregulated in M1 activation and 2 most upregulated in M2 activation) for qRT‐PCR validation and translational exploration. Using a set of independent macrophage samples (n=8 subjects), qRT‐PCR analysis replicated the pattern of activation‐induced lincRNA expression identified at RNA‐seq for all lincRNAs (Figure 5C for M1‐induced, and Figure 5D for M2‐induced). Of these, MIR155HG is nonsyntenic; RP11‐10J5.1, RP11‐701P16.5, CTB‐41I6.2, and RP5‐836N10.1 are syntenic but not annotated as expressed in mouse; and linc‐HEATR6‐2, linc‐SLC39A10‐10, MIR146A, RP4‐794H19.4 and RP11‐184M15.1 are syntenic and annotated in mouse (Figure 5C and and5D).5D). None of these lincRNAs showed significant sequence conservation in mouse. Of these 10 lincRNAs, 6 had enhancer‐like histone signatures (see Table S13 for a summary). As an example, we showed the qRT‐PCR validation of CUFF.15750, one of the most abundantly expressed de novo annotated lincRNAs, which was suppressed during both M1 and M2 activation and has PU.1 and C/EBPβ binding and enhancer‐like features. Public cap analysis of gene expression peak data were consistent with the apparent TSS for CUFF.15750 revealed by our RNA‐seq (Figure S6).
The majority of genetic variants associated with complex diseases are found within noncoding regions of the genome, where the functional consequences of the variation are largely unknown. Consequently, we explored the overlap of macrophage lincRNAs with disease‐associated genetic variations in public data sets. First, to probe broadly whether macrophage lincRNAs may underlie disease associations, we explored genomic loci for the 595 macrophage‐enriched lincRNAs that contained SNP data within the comprehensive GWAS catalog27 with trait associations of P<1×10−5. We identified 66 macrophage‐enriched lincRNAs containing trait‐associated SNPs and highlighted those traits for which macrophages have been implicated, including metabolic (eg, obesity‐related traits, visceral fat, and waist–hip ratio) and immune disorders (eg, Crohn disease, multiple sclerosis, and celiac disease; boldface in Table S14).
Second, because of the central role of macrophage activation in multiple CMDs, we interrogated SNPs within all macrophage‐expressed lincRNAs for their specific association with 13 cardiometabolic traits (Table S2). Of the 2766 macrophage‐expressed lincRNAs, 2340 lincRNAs contained SNPs that were tested in at least 1 of the 13 GWAS data sets. Using our published pipelines,20, 30 lincRNAs containing significant trait‐associated SNPs were filtered stringently to include only those that contained the strongest and independent (r 2<0.3; based on 1000Genomes CEU data44) SNP‐level P value in the region (±500 kb of the lincRNA; see Methods for details). By further filtering for the most prominently expressed lincRNAs (FPKM >0.1, corresponding to top ≈35% expressed macrophage lincRNAs), we identified 3 independent trait‐associated SNPs—for waist‐to‐hip ratio adjusted for body mass index, plasma triglycerides, and plasma low‐density lipoprotein‐cholesterol—that fall within ±1 kb of highly expressed macrophage lincRNAs (Table).
A top trait association is at a lincRNA annotated as RP11‐472N13.3, which we named MacORIS. MacORIS overlaps rs7081678, an SNP associated with central obesity (waist–hip ratio adjusted for body mass index); maps to the chromosome 10p11.22 locus (Figure 6A through through6C);6C); and is a macrophage‐enriched lincRNA that is syntenic but not annotated in mouse. MacORIS is expressed predominantly in M0 HMDM (fractional expression value: 0.44), is barely detectable in human primary adipocytes, and is found at low levels in human adipose cells (Figure 6D) and T cells (Figure 6E). A genome browser view of MacORIS shows abundant PU.1 and C/EBPβ binding (Figure 6F) but no annotation in Genome Reference Consortium Mouse Build 38 (GRC38/mm10) and no expression in published high‐quality RNA‐seq of murine bone marrow–derived macrophages (Figure 6G and Table S15). MacORIS does not contain a conserved open reading frame, and in vitro transcription and translation of MacORIS did not produce any detectable peptides (Figure S7A). The qRT‐PCR of cell fractions revealed that MacORIS is predominantly located in cytoplasm (Figure 6H), suggesting potential posttranscriptional regulatory roles. M1, but not M2, stimulation suppressed MacORIS expression (Figure 6I and and6J).6J). To examine the functional impact of MacORIS on M1 activation, we used GapmeR antisense oligonucleotide to knock down MacORIS in THP‐1Φ and found enhanced expression of IFN‐γ–induced negative regulators SOCS1 and SOCS3 but no effect on lipopolysaccharide‐induced inflammatory genes such as TNF, TNFAIP3, and IL1B (Figure 6K). Cytoplasmic localization suggests that MacORIS modulates cytoplasmic activation, rather than nuclear expression, of IFN‐γ–signaling molecules. IFN‐γ activates IFNGR1 and IFNGR2, via transphosphorylation of JAK1 and JAK2 and with downstream phosphorylation of STAT1 leading to oxidative burst and expression of IFN‐γ–inducible genes, including IL12.45 Knockdown of MacORIS in THP‐1Φ enhanced the phosphorylation of STAT1 without altering total protein levels of IFNGR1, JAK2, or STAT1 (Figure 6L and and6M).6M). Independent validation in M1 THP‐1Φ with knockdown of MacORIS by small interfering RNA showed generally consistent results except that knockdown of MacORIS enhanced JAK2 as well as STAT1 phosphorylation (Figure S7B and S7C); this difference may be attributable to the differential activity and mechanisms of antisense oligonucleotide– versus small interfering RNA–mediated knockdown for nuclear relative to cytoplasmic targets. Overall, these data suggest that cytoplasmic MacORIS serves as a repressor of macrophage IFN‐γ signal transduction by modulating, via as‐yet‐unknown mechanisms, JAK2/STAT1 phosphorylation, thus regulating downstream IFN‐γ–responsive gene expression. Whether the central obesity–associated SNPs at this locus modulate MacORIS expression and function remains to be determined, but macrophage IFN‐γ signaling by MacORIS is a very plausible mechanism for modulation of central obesity and related metabolic disorders at this locus.46
It is important to consider human‐relevant strategies and to develop tools for functional interrogation of human lincRNAs not expressed in mouse. IPSDMs are a renewable source of subject‐specific macrophages and provide a powerful functional genomic tool to address human macrophage biology. We reported previously that IPSDMs had comparable phenotypes, protein‐coding transcriptomes, and functional characteristics as HMDMs and can be used for functional genomic modeling of protein‐coding genes.18 In this article, we extended our IPSDM model for our current lincRNA perspective by examining DE lincRNAs between induced pluripotent stem cells and IPSDMs and comparing resting and activation profiles of lincRNAs in IPSDMs versus HMDMs (Figure 7A).18 A multidimensional scaling plot based on expression of lincRNAs (Figure 7B) revealed that HMDMs and IPSDMs cluster together and are completely distinct from induced pluripotent stem cells; M1 HMDMs and M1 IPSDMs also cluster together and separately from M0 or M2 HMDMs. Differentiation of induced pluripotent stem cells to IPSDMs induced marked lincRNA transcriptome changes with 313 DE lincRNAs. Compared with all other IPSDM lincRNAs, the 153 lincRNAs upregulated during differentiation of induced pluripotent stem cells to IPSDMs (Figure 7C and Table S16) had higher expression (Figure 7D) and had enriched PU.1 and C/EBPβ TF binding (Figure 7E).
The vast majority (>90%) of the M0 HMDM lincRNAs were also present in M0 IPSDMs, and their expression was moderately correlated (r=0.51; Figure 7F). Remarkably, for ≈95% of lincRNAs, there was a similar pattern of activation‐related change in expression in both HMDMs and IPSDMs with strong correlations (eg, r=0.81 between IPSDMs and HMDMs for M1‐activation–induced fold change of lincRNAs; Figure 7G and and7H;7H; Tables S11 and S12). Indeed, only very few lincRNAs were DE between HMDMs and IPSDMs (Figure 7F and Table S17). For the very small number of lincRNAs that were expressed at lower levels in M0 IPSDMs than in M0 HMDMs (eg, linc‐SLC39A10‐10), on activation, their expression in M1 or M2 IPSDMs was comparable to that in M1 or M2 HMDMs (Table S18). As a relevant example, a genome browser view of MacORIS shows consistent expression patterns for HMDMs and IPSDMs at rest and during M1 activation (Figure 6F). Overall, IPSDM lincRNA expression and activation profiles resemble those of HMDMs, supporting the utility of the IPSDM system for functional modeling of lincRNAs in human macrophage genomics.
Macrophages modulate many human pathophysiologies and have emerged as potential therapeutic targets in complex diseases.5 Although a recent microarray‐based study has characterized lncRNAs in M1‐ and M2‐activated HMDMs,14 there is a lack of RNA‐seq–based, unbiased cataloging of the human macrophage lncRNA transcriptome. By exploiting de novo transcriptome reconstruction of deep RNA‐seq data, we provide the most comprehensive inventory and genomic profile, to our knowledge, of polyadenylated human macrophage lincRNAs. We identified 2766 macrophage‐expressed lincRNAs, 861 of which are newly annotated. Most (85%) macrophage lincRNAs are nonsyntenic or are syntenic but not annotated as expressed in mouse. Many lincRNAs are enriched in macrophages, overlap PU.1 and C/EBPβ transcription factor binding sites, and display enhancer‐like chromatin signatures, and multiple macrophage‐enriched lincRNAs were also found to overlap GWAS loci for CMD traits. Macrophage activation, particularly to the M1 phenotype, markedly alters the lincRNA expression profiles, suggesting a role for lincRNAs in macrophage functional activation. MacORIS, a human macrophage‐specific cytoplasmic lincRNA that contains SNPs associated with central obesity, functions as a brake on macrophage IFN‐γ signaling and inflammatory responses. Finally, because many human macrophage lincRNAs are not conserved in mouse, our efficient and scalable human IPSDM system provides a valuable cellular model for functional assessment of lincRNAs in human macrophage biology.
Although reductionist relative to in vivo phenotype complexity, in vitro activation to M1 or M2 macrophage phenotypes has proven useful in defining functional states toward which macrophages can be driven in distinct inflammatory milieu.6, 7 For example, multiple macrophage protein‐coding genes (eg, IL6, TNF, IL1B) of functional importance are markedly induced during M1 activation in vitro and in vivo.18 Consistent with the pattern for mRNAs,7, 8, 9, 18 M1 activation induces profound changes in lincRNA expression with induction of dozens of lincRNAs. Correlation of activation‐dependent change in enhancer lincRNA expression with that of the nearest protein‐coding genes maps to regulation of immune system processes and suggests an integrative regulatory role for some lincRNAs during macrophage activation. Indeed, through our prioritization strategy, we identified that lincRNAs reported previously to modulate myeloid cell functions (eg, linc‐HEATR6‐2, also named lnc‐DC) were recently reported to regulate dendritic cell maturation and function.47 Furthermore, 2 prioritized lincRNAs, MIR155HG and MIR146A, are microRNA host genes for miR‐15548 and miR‐146a,49 2 well‐characterized microRNAs that regulate macrophage inflammatory responses.48, 49 This strategy identified multiple other lincRNAs as promising cis‐regulatory candidates for functional and translational interrogation (eg, AC002480.2 proximal to IL6 and CUFF.135177 proximal to CCL8). A recent in vitro RNA‐seq study of lipopolysaccharide‐stimulated (4 hours) monocytes discovered DE lncRNAs that modulate monocyte response to lipopolysaccharide.39 We identified 49 DE lincRNAs induced in lipopolysaccharide‐treated monocytes, and 14 of these overlapped lincRNAs also upregulated during M1 HMDM activation. Nevertheless, a much larger proportion (59 of 73) of lincRNAs induced during M1 activation are not identified in monocyte activation, suggesting specific macrophage induction and function (Table S19). M1 and M2 activation in vitro, however, provides a relatively narrow window into the diversity of macrophage activation states observed in vivo; future transcriptional profiling of resident macrophages across diverse tissues and settings will provide deeper insight into the in vivo complexity of the human macrophage noncoding transcriptome.
Recent GWASs have revealed novel functional lncRNAs in disease, for example, ANRIL at the 9p21.3 locus for coronary heart disease50 and Lnc13 at 2q12.1 for celiac disease.51 These human genetic studies suggest that lincRNAs may play important modulatory roles in human diseases. Indeed, we identified hundreds of macrophage lincRNAs that reside within intergenic loci previously identified by GWASs for complex traits. We performed a deeper interrogation of lincRNAs in 13 CMD data sets and identified several promising candidates including MacORIS, which we found to act as a repressor of IFN‐γ signaling by regulating phosphorylation of JAK2 and STAT1. Notably, IFN‐γ deficiency protects mice from high fat diet–induced white adipose tissue inflammatory cell accumulation and glucose tolerance.46 Thus, MacORIS modulation of IFN‐γ signaling in macrophages is a plausible mechanism underlying the 10p11.22 locus for central obesity. However, the causal variant at MacORIS and the precise genetic and cellular mechanisms of action of MacORIS require further investigation.
MacORIS is one of many human lincRNAs not present in mice. This lack of conservation combined with historical limitations of human macrophage models presents a specific challenge to functional studies of lincRNAs in human macrophage biology. RNA interference and antisense oligonucleotide–based knockdown approaches in primary monocytes and macrophages are challenging, given low transfection rates and heterogeneity between experiments. THP‐1 monocyte and macrophage lines, although useful, as we demonstrate for MacORIS, are karyotypically abnormal and phenotypically immature,15 thus also not an ideal model for human functional genomics. We developed a high‐fidelity model for human macrophage functional genomics studies.18 Our results in this study reveal comparable lincRNA transcriptome profiles and dynamic regulation during activation in isogenic IPSDMs and HMDMs. Coupled to CRISPR/Cas9 gene editing that precisely introduces targeted mutations and deletions,52 IPSDM provides a powerful tool to decipher the genomic and molecular regulation of human macrophage lincRNAs in human physiology and disease.
Recently, the FANTOM CAT (CAGE‐associated transcriptome)—a human transcriptome meta‐assembly based on cap analysis of gene expression data across 1829 samples from major human primary cell types and tissues as well as transcript models from GENCODE V19, the Cabili set, miTranscriptome, and ENCODE—has defined 27 919 lncRNAs, of which 13 105 were lincRNAs.53 We found 901 of 2766 of our macrophage lincRNAs overlapped FANTOM CAT lincRNAs within ±250 bp of the TSS (Table S20). Because the FANTOM5 CAT included human macrophages from only 3 donors, additional macrophage lincRNAs will be added to such public resources as sample size and sequencing depth increase, as in our study. Nonetheless, the precise 5′‐end transcript mapping in FANTOM5 CAT lincRNAs is complementary to but less comprehensive than our deep RNA‐seq–based human macrophage lincRNA catalog.
In the current work, we focused on lincRNAs for both technical and translational reasons. LncRNAs that either overlap (ie, antisense) or share a TSS interval with protein‐coding genes confound simple interpretation of regulatory features in the region and complicate genetic manipulation in functional studies. There are also analytic challenges in dissecting the contribution of GWAS disease‐associated SNPs residing in lncRNAs that overlap protein‐coding genes. Consequently, lncRNAs excluded from our analysis, including lincRNAs proximal to the coding genes, antisense lncRNAs shown to regulate THP‐1Φ function,54 and single‐exon transcripts are likely to provide additional layers of information about the macrophage noncoding transcriptome.
Our study has many strengths but limitations too. Our lincRNA catalog derived from poly(A) capture RNA‐seq fails to include nonpolyadenylated lncRNAs and short noncoding RNAs. It has been reported that 84.2% and 74.2% of the annotated expressed lncRNAs are poly(A)+ in H9 and Hela cells, respectively; 13.1% and 23.3% are bimorphic, found in both the poly(A)+ and poly(A)− populations, respectively; and 2.7% and 2.5% are poly(A)−, suggesting the majority of the lncRNAs are poly(A)+ or bimorphic, respectively.55 The classification, however, has not been performed in human macrophage. Coding potential was assessed by computational prediction using iSeeRNA33 and Pfam34 with validation by PhyloCSF35 but not with experimental approaches. A fully comprehensive macrophage lncRNA catalog derived from RNA‐seq of ribosomal RNA–depleted samples combined with both bioinformatic and experimental approaches for coding potential assessment will further refine the human macrophage lncRNA catalog for future study. In the meantime, a large number of prioritized lincRNAs in our study remains to be functionally validated to gain deeper mechanistic insights into lincRNA modulation of human macrophage biology and their role in human diseases.
Our work underscores the importance of lincRNA discovery studies, using deep RNA‐seq and de novo assembly, in a species‐ and tissue‐specific manner. It also provides a resource to parse the polyadenylated lincRNA circuitry of macrophage activation and to identify specific lincRNAs for functional studies in macrophage activation and macrophage‐related human diseases, as we have explored for MacORIS. Our IPSDM model provides a unique framework with which to pursue the human macrophage–specific functions of novel lincRNAs in macrophage biology and related diseases and for gene‐editing strategies to advance mechanism‐based clinical and therapeutic translation of human genomic discoveries.
This work was supported by National Institutes of Health (NIH) grants R01‐HL‐113147 (to Reilly and Li), R01‐HL‐132561 (to Reilly), K24‐HL‐107643 (to Reilly), and U01‐HG006398 (to Rader and Morrisey). Reilly is also supported by NIH R01‐HL‐111694. Li is also supported by NIH R01‐GM‐108600. H. Zhang is supported by NIH K99‐HL‐130574. Yang and the University of Pennsylvania Induced Pluripotent Stem Cell Core Facility are supported by the University of Pennsylvania Institute for Regenerative Medicine. Research reported in this publication also used the resources of Columbia University's Cancer Center Flow Core Facility funded in part through NIH Center Grant P30CA013696, and Columbia's CCTI Flow Cytometry Core, supported in part by the Office of the Director, NIH under awards S10RR027050. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.
Data S1. Detailed Methods.
Table S1. Subject Demographics of RNA‐seq (RNA Sequencing) Studies
Table S2. Cardiometabolic Traits Evaluated in Genetics Consortia
Table S3. Primers for qRT‐PCR (Quantitative Reverse Transcriptase–Polymerase Chain Reaction) Validation
Table S4. Coding Potential Assessment of Macrophage lincRNAs (Long Intergenic Noncoding RNAs) With PhyloCSF Using Score Cutoff of 100
Table S5. Annotation, Expression, Synteny, Conservation, Tissue Enrichment, and Chromatin Signature of Macrophage lincRNAs (Long Intergenic Noncoding RNAs)
Table S6. DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs); Monocyte Versus M0 HMDM (Human Peripheral Blood Mononuclear Cell–Derived Macrophages)
Table S7. Top Enriched Gene Ontology Terms for the Nearest Protein Coding Genes to the Upregulated Enhancer‐Associated lincRNAs (Long Intergenic Noncoding RNAs) in M1 Activation
Table S8. Top Canonical Pathways for the Nearest Protein Coding Genes to the Upregulated Enhancer‐Associated lincRNAs (Long Intergenic Noncoding RNAs) in M1 Activation
Table S9. Top Diseases and Biological Functions for the Nearest Coding Genes to the Upregulated Enhancer‐Associated lincRNAs (Long Intergenic Noncoding RNAs) in M1 Activation
Table S10. Fold Change of Enhancer‐Associated lincRNAs (Long Intergenic Noncoding RNAs) and Their Nearest Coding Genes Upregulated in M1 Activation
Table S11. DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs); M0 HMDM (Human Peripheral Blood Mononuclear Cell–Derived Macrophages) Versus M1 HMDM and M0 IPSDM (Induced Pluripotent Stem Cell–Derived Macrophages) Versus M1 IPSDM
Table S12. DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs); M0 HMDM (Human Peripheral Blood Mononuclear Cell–Derived Macrophages) Versus M2 HMDM and M0 IPSDM (Induced Pluripotent Stem Cell–Derived Macrophages) Versus M2 IPSDM
Table S13. Prioritized Differentially Expressed lincRNAs (Long Intergenic Noncoding RNAs) in M1 and M2 Activation
Table S14. Macrophage lincRNAs (Long Intergenic Noncoding RNAs) Harbor Genetic Variants Associated With Traits in GWASs (Genomewide Association Studies)
Table S15. QC Summary of Public Data Sets of RNA‐seq (RNA Sequencing) of Murine Bone Marrow–Derived Macrophages
Table S16. DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs); Induced Pluripotent Stem Cells Versus M0 IPSDM (Induced Pluripotent Stem Cell–Derived Macrophages)
Table S17. DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs); M0 HMDM (Human Peripheral Blood Mononuclear Cell–Derived Macrophages) Versus M0 IPSDM (Induced Pluripotent Stem Cell–Derived Macrophages)
Table S18. Expression of DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs) Between M0 HMDM (Human Peripheral Blood Mononuclear Cell–Derived Macrophages) and M0 IPSDM (Induced Pluripotent Stem Cell–Derived Macrophages) on M1 and M2 Activation
Table S19. Comparison of DE lincRNAs (Differentially Expressed Long Intergenic Noncoding RNAs) in Lipopolysaccharide‐Treated Monocytes and M1 HMDM (Human Peripheral Blood Mononuclear Cell–Derived Macrophages)
Table S20. Macrophage lincRNAs (Long Intergenic Noncoding RNAs) With FANTOM CAT lncRNA (Long Noncoding RNA) Catalog Annotation
Figure S1. Correlation between biological replicates.
Figure S2. The “activation state”‐specific long intergenic noncoding RNAs (lincRNAs) were more likely to be previously unannotated lincRNAs.
Figure S3. Characteristics of long intergenic noncoding RNAs (lincRNAs) vs protein coding genes.
Figure S4. Profound differentiation‐induced long intergenic noncoding RNA (lincRNA) profile change from CD14+ monocytes to macrophages.
Figure S5. Characteristics of enhancer‐associated long intergenic noncoding RNAs (lincRNAs) and gene ontology and Ingenuity Pathway Analysis of the nearest protein coding genes to the upregulated enhancer‐associated lincRNAs during M1 activation.
Figure S6. Quantitative reverse transcriptase–polymerase chain reaction (qRT‐PCR) validation and genome browser view of a newly annotated long intergenic noncoding RNA CUFF.15710.
Figure S7. In vitro transcription and translation of MacORIS (macrophage‐enriched obesity‐associated long intergenic noncoding RNA serving as a repressor of IFN‐γ [interferon γ] signaling) and knock down of MacORIS by small interfering RNA.
(J Am Heart Assoc. 2017;6:e007431 DOI: 10.1161/JAHA.117.007431.)
The abstract of this article was presented at the American Heart Association Scientific Sessions, November 11–15, 2017, in Anaheim, CA.
Hanrui Zhang, Email: ude.aibmuloc.cmuc@8142zh.
Muredach P. Reilly, Email: ude.aibmuloc.cmuc@4412rpm.