Single-cell profiling reveals GPCR heterogeneity and functional patterning during neuroinflammation

Multiple sclerosis (MS) is a neuroinflammatory disease affecting the CNS. It is characterized by chronic focal demyelination and neurodegeneration due to immune-mediated destruction of myelin sheaths and inflammatory attack toward the neuronal compartment (1, 2). With a prevalence of 2.5 million patients worldwide, MS is the most frequent autoimmune disorder of the CNS. Studies in humans and mice demonstrated a prominent role of CD4 T cells in the pathogenesis of autoimmune demyelinating diseases, in particular of T-bet1/IFNγ–positive Th1 cells and RORγt/IL-17–positive Th17 cells (3–6). Current therapeutic concepts for the treatment of MS include i.v. corticoids during relapse and β IFNs, glatiramer acetate, teriflunomide mitoxantrone, natalizumab, or alemtuzumab during remission (1, 7). Relatively new to the treatment of MS are drugs modulating the function of GPCRs, for example the sphingosine 1 phosphate (S1P) receptor modulator fingolimod (8) or dimethyl fumarate, which activates the hydroxycarboxylic acid receptor HCA₂ through its metabolite monomethyl fumarate (9). Other GPCR-based therapeutic strategies currently in clinical and preclinical investigation include the modulation of chemokine (10) or endocannabinoid receptor signaling (11). However, given the fact that the human genome encodes approximately 800 GPCRs, among them 367 nonolfactory receptors (12), it seems likely that the therapeutic potential of GPCRs in MS is not yet fully exploited (13).

GPCRs are the largest family of transmembrane receptors in eukaryotes; they transduce signals as diverse as neurotransmitters, hormones, local mediators, metabolic signals, olfactory cues, and light (14). Various studies addressed the role of individual GPCR families in immune cells — most prominently among them, the chemokine receptors (15), lysophospholipid receptors (16, 17), or prostanoid receptors (18, 19). A systematic analysis of GPCR expression in murine lymphoid organs revealed that approximately 255 GPCRs are expressed in spleen, thymus, and BM (20). However, since these studies were performed on pooled RNA/cDNA, they only report an average expression for a given population. This is a major limitation, since immune cell populations are not homogenous, but consist of various subpopulations of varying activation states. There is ample evidence that the functional behavior of immune cell subpopulations is influenced by their GPCR repertoire; in circulating naive T cells, for example, egress from lymph nodes is critically regulated by the expression of the S1P₁ receptor (17, 21). Furthermore, T cell activation is characterized by downregulation of CCR7 and upregulation of various inflammatory GPCRs, the identity of which will determine tropism of the effector cells (22, 23). Also, GPCRs on endothelial cells (ECs) importantly contribute to inflammatory processes by regulating endothelial permeability and inflammatory gene expression (24, 25). It seems crucial for the development of new therapeutic strategies to identify GPCRs that are associated with specific functional states within these various cell types. However, systematic analyses of the correlation between functional state and GPCR expression were, so far, not possible, since single-cell–based techniques such as flow cytometry or in situ hybridization only allowed the targeting of a limited number of proteins or genes at a time and were often hampered by the availability of reliable antibodies, probes, or sufficient amount of cells. In order to overcome this limitation, protocols for single-cell expression analysis have been developed; they rely on the reverse transcriptase PCR–based (RT-PCR–based) amplification of mRNAs of cells isolated by laser capture microdissection, fluorescence-activated cell sorting, or microfluidic approaches (26, 27). These approaches allow for the analysis of tens to hundreds of genes in single cells and may be adapted to perform whole transcriptome analyses either by using microarrays or next-generation sequencing (28). Here, we report a microfluidic-based GPCR expression analysis of different cell types implicated in the pathogenesis of human MS and its mouse model experimental autoimmune encephalomyelitis (EAE) and show examples of how this information may be used to modulate cell functions.

Array design and quality control. To reduce the number of GPCRs to be investigated on the single-cell level, we determined GPCR expression in bulk RNA of freshly isolated murine T cells from naive lymph nodes or inflamed spinal cord, as well as from ECs. Of 424 tested GPCRs, 146 were expressed in at least one of these 3 cell types (Supplemental Figure 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.95063DS1). For these 138 receptors, as well as for 6 additional GPCRs, a primer array for single-cell expression analysis was designed (Table 1). In addition, 11 genes identifying individual cell types and 39 function-defining genes — including the 3 reference genes β-actin (Actb), glyceraldehyde 3-phosphate dehydrogenase (Gapdh), and hypoxanthine phosphoribosyltransferase (Hprt) — were included in the array (Table 1). GPCRs that did not allow intron-spanning design had to be excluded from the analysis because they gave positive results in all tested cell types, indicative of contamination with genomic DNA (strikethrough in Table 1). The expression of the remaining 175 genes was systematically analyzed in freshly isolated individual CD4 T cells, macrophages, and ECs from blood or lymph nodes of healthy mice as well as draining lymph nodes at EAE onset or spinal cords at peak of disease. Despite stringent sorting criteria, single-cell analysis of individual cells showed a high degree of contamination with other cell types (Supplemental Figure 2A). Only cells with clear identity and positive expression of reference genes as cDNA quality control were used for further analyses (Supplemental Figure 2B).

Table 1

Overview over genes included in the array

Single-cell GPCR expression in naive CD4 T cells. In a first set of experiments, we studied GPCR expression in individual CD4 T cells from lymph nodes of naive mice (CD4ln) and found it to be highly heterogeneous (Figure 1A and Supplemental Figure 3). On average, individual cells expressed 8.6 ± 0.5 of the 125 analyzed GPCRs, though the individual values varied between 1 and 18 GPCRs per cell (Figure 1B). Though all mice were C57BL/6 and kept under similar pathogen-free conditions, GPCR frequencies varied to some degree between individual mice (Figure 1C). High heterogeneity was also observed in naive murine CD4 T cells from peripheral blood (CD4bl), as well as in human blood CD4 T cells (Figure 1D). To further validate this high degree of GPCR heterogeneity, we compared single-cell expression data with data obtained in GPCR reporter mice or antibody-mediated GPCR staining. Flow cytometric analysis of the percentage of Ccr7- or Ccr6-expressing CD4ln cells resulted in values closely matching the percentages determined by single-cell RT-PCR, and the same was true for the flow cytometric analysis of CD4ln from knockin mice that express EGFP under control of the Gpr183 promoter (Figure 1E).

Figure 1

Single-cell GPCR expression in freshly isolated murine and human CD4 T cells from lymph node (CD4ln) and blood (CD4bl). (A) Heat map of GPCR expression in murine CD4ln. GPCRs are sorted by frequency; horizontal bars on the right side visualize expression frequency (%). GPCRs expressed in less than 5% of cell are not shown (each column one cell; for complete data set, see Supplemental Figure 3). Ptprc (encoding CD45), Cd4, Actb, and Gapdh are shown as quality controls (50 cells from 3 mice). (B) Number of GPCRs expressed in individual CD4ln. (C) Frequency of GPCR expression in CD4ln from different donor mice (16–17 cells per donor). (D) Frequency of GPCR expression in CD4bl from different murine (32 cells from 5 mice) and human donors (140 cells from 4 donors) (ranked by human expression; cutoff 5%). (E) Comparison of GPCR expression frequency as judged by single-cell RT-PCR (sc RT-PCR) or flow cytometric analysis of the percentage of antibody-stained (Ccr7, Ccr6) or GFP-expressing (Gpr183) CD4ln and 2D2 CD4ln (only for GPR183). Expression data are calculated as 2^{(Lod Ct – sample Ct)}; LoD Ct set to 24.

Activation-dependent changes in GPCR heterogeneity in CD4 T cells. We next investigated GPCR heterogeneity in in vivo–activated CD4 T cells. To do so, we used the EAE model and compared CD4ln from naive mice with CD4 T cells isolated from draining lymph nodes (CD4dr) during disease induction (days 8–12) and spinal cord–infiltrating CD4 T cells (CD4sc) harvested at peak of disease (days 14–17) (Figure 2A). The percentage of cells expressing a given GPCR was, in most cases, greatly increased in CD4sc (Figure 2, A and B); this was particularly true for chemokine receptors, but also for receptors of inflammatory mediators such as thrombin receptors F2r and F2rl2 (encoding protease activated receptor subtypes PAR₁ and PAR₃), leukotriene receptor Cysltr1, prostaglandin E2 receptor subtypes Ptger2 and Ptger4 (encoding EP₂, EP₄), or oxysterol receptor Gpr183 (encoding EBI2) (Figure 2A and Supplemental Figure 3). In addition, orphan receptors such as Cd97, Gpr65, Gpr107, Gpr108, or P2ry10 were upregulated in frequency and/or intensity (Figure 2A). Other receptor mRNAs were downregulated — for example, those of chemokine receptor Ccr7, calcitonin receptor-like receptor Calcrl, or vasoactive intestinal polypeptide receptor 1 Vipr1 (Figure 2A). Changes in CD4dr were, in most cases, less pronounced than in CD4sc (Figure 2, A and B, and Supplemental Figure 4A). A significant increase in GPCR expression was also observed after Th1- or Th17-directed in vitro MOG_35-55 stimulation of splenic CD4 T cells from mice carrying the MOG_35-55-specifc T cell receptor 2D2 (Figure 2C), but k-means cluster analysis showed that in vitro– and in vivo–differentiated cells were clearly distinct with respect to their GPCR profile (Figure 2D and Supplemental Figures 4, B–D). We also investigated whether the GPCR expression patterns observed in CD4sc of actively MOG_35-55-immunized mice were conserved in other murine models of MS, such as adoptive transfer EAE (CD4scAT), and found that the majority of GPCRs showed similar regulation in both EAE models (Supplemental Figure 4E).

Figure 2

Single-cell GPCR expression in CD4 T cells after MOG_35-55-dependent activation in vivo or in vitro. (A) Heat map of GPCR expression in CD4ln, CD4dr, and CD4sc (50, 43, 115 cells from 3, 3, 12 mice, respectively); horizontal bars on the right side visualize expression frequency (%). (B) Number of GPCRs expressed in individual CD4ln, CD4dr, CD4sc. (C) Number of GPCRs expressed in individual splenic CD4 cells (CD4spn) from 2D2 TCR transgenic mice in the naive state or after in vitro differentiation toward Th1 or Th17, respectively. (D) T-SNE plot showing the degree of similarity between individual in vitro–differentiated (Th1, 26 cells; Th17, 43 cells) and in vivo–differentiated CD4 T cells. The closer together cells are plotted, the more similar they are; k-means cluster assignment is indicated by color; cell type is indicated by symbols. Expression data are calculated as 2^{(Lod Ct – sample Ct)}; LoD Ct was set to 24. ***P < 0.001 (B, C) by unpaired t test.

Subgroups within CD4sc. To identify GPCRs that are enriched in rare subpopulations, we performed k-means cluster analysis within CD4sc. This analysis identified various subgroups that showed characteristic differences in the expression of function-defining genes and GPCR repertoire (Figure 3A). Cluster 6 cells were characterized by high expression of Foxp3, indicating they were spinal cord–infiltrating Treg (Figure 3B). In cluster 4 cells, expression of Tbx21 (encoding T-bet1) was overrepresented, indicating these cells were Th1 cells (Figure 3C). Cluster 1 cells were characterized by high expression of granulocyte-macrophage CSF (GM-CSF, encoded by Csf2), IL-17 (Il17), tumor necrosis factor α (Tnf), and transcription factor Rorγt (Rorc), indicative of a highly pathogenic Th17 subtype (Th17^high_path, Figure 3D). Cluster 3 cells also showed characteristics of Th17 cells (expression of Il17 and Rorc) but low expression of pathogenicity markers such as Csf2 and Tnf (Figure 3E), suggesting they represent a Th17 population with low pathogenicity (Th17^low_path). A direct comparison of Rorc-positive cells from clusters 1 and 3 showed significantly increased expression of GPCRs Gpr18, Gpr132, S1pr1, P2ry10, Cxcr3, and Cxcr4 in Th17^low_path (Figure 3F), whereas expression of Ccr8, Gpr65, and Cxcr6 was increased in Th17^high_path. Spearman’s rank analysis of the correlation between individual genes confirmed that genes Klf2, P2ry10, S1pr1, Gpr132, Cxcr3, and Cxcr4 on the one hand and genes Il17, Rorc, Csf2, CD25, Tnf, and Ccr8 on the other hand show correlated expression (Figure 3G).

Figure 3

Functional subgroups within spinal cord–infiltrating CD4 T cells (CD4sc). (A) T-SNE plot showing similarities between individual CD4sc and results of k-means cluster analysis (the closer together cells are plotted, the more similar they are; k-means cluster assignment is indicated by color; 136 cells). (B–E) Genes differentially expressed in cluster 6 (B), cluster 4 (C), cluster 1 (D), or cluster 3 (E) compared with all cells. Only genes with fold enrichment > 1.5 or < 0.7 and P < 0.05 are shown. (F) Heat map of CD4sc cells from cluster 1 (top) and cluster 3 (bottom) (only Rorc-positive cells are shown). (G) Representation of Spearman’s rank coefficients for the correlation of expression between individual genes (width of connecting lines indicates strength of correlation). (H and I) Flow cytometric analysis of TNFα and GM-CSF expression in CCR8-negative and CCR8-positive Th17-CD4sc. (H) Gating strategy for cell sort and cytokine expression analysis in CD4⁺IL-17A⁺CCR8^– (left) and CD4⁺IL-17A⁺CCR8⁺ T cells. (I) Statistical evaluation (n = 3). (J) Statistical analysis of the percentage of TNFα/GMCSF–double positive cells in CXCR4-negative and CXCR4-positive Th17-CD4sc (n = 4). Expression data are calculated as 2^{(Lod Ct – sample Ct)}; LoD Ct was set to 24. Function-defining genes are shown in blue.*P < 0.05; **P < 0.01; ***P < 0.001 (F, unpaired t test; I–J, paired t test).

To investigate whether these subgroups of Th17 cells also existed on the protein level, we analyzed expression of TNFα and GM-CSF in CCR8-positive and CCR8-negative Th17 cells isolated from the inflamed spinal cord. To do so, CD4sc were sorted into a CCR8-negative and a CCR8-positive population, and levels of intracellular cytokines were determined in IL-17A–positive cells by flow cytometry (Figure 3H). In line with our transcriptional analysis, we found that CCR8-positive Th17-CD4sc contained a significantly higher percentage of TNFα/GM-CSF–double positive cells than CCR8-negative Th17-CD4sc (Figure 3, H and I), indicating that the subgroups observed on the transcriptional level also exist on the protein level. In line with this notion, we found that CXCR4-positive spinal cord–infiltrating Th17 cells — which, according to transcriptional analysis, show reduced expression of Tnf and Csf2 — indeed showed a reduced percentage of TNFα/GM-CSF–double positive Th17 cells (Figure 3J).

Pharmacological modulation of Th17 pathogenicity. Th17 cells are known to be a highly plastic population (29, 30), and it seems likely that, depending on the local micromilieu, Th17^low_path have the potential to develop into Th17^high_path and vice versa. We therefore investigated whether GPCRs enriched in Th17^low_path have the potential to modulate the transition into Th17^high_path. One of the receptors negatively correlated with expression of Csf2 and Tnf was Cxcr4 (Figure 4A), a chemokine receptor that has been suggested to regulate not only cell migration (31, 32), but also T cell receptor signaling (33). It is unclear whether CXCR4 is able to modulate the functional state of differentiated Th17 cells. We found that in vitro–differentiated Th17 cells that had been cultured for 5 days in the presence of CXCR4 antagonist AMD3100 showed a reduced percentage of IL-17–expressing cells (Figure 4, B and C). The secretion of IL-17, CCL3, CCL4, and TNFα was correspondingly diminished in AMD3100-treated cells (Figure 4, D–G), while other cytokines such as IL-2 or IL-10 were not altered (data not shown). The inhibitory effect of AMD3100 on the murine Th17 population was also observed in human Th17 cells (Figure 4H).

Figure 4

Modulation of Th17 pathogenicity by the CXCR4 antagonist AMD3100 and S1P₁ agonist SEW2871. (A) Negative correlation between expression of Cxr4 and Csf2 or Tnf in CD4sc (r indicates Spearman rank correlation coefficient). (B and C) Flow cytometric analysis of the effect of CXCR4 antagonist AMD3100 on the percentage of IL-17–expressing cells in in vitro–differentiated murine Th17 cells (B, exemplary dot plots [gate: CD4]; C, statistical evaluation) (n = 6). (D–G) Concentration of cytokines as determined by multiplex ELISA in supernatants of AMD3100- or vehicle-treated in vitro–differentiated murine Th17 cells (n = 6). (H) Effect of AMD3100 on the percentage of IL-17–expressing cells in in vitro–differentiated human Th17 cells (gate: CD4; n = 3). (I) Negative correlation between expression of S1pr1 and Rorc or Tnf in CD4sc (r indicates Spearman rank correlation coefficient). (J and K) Flow cytometric analysis of IL-17A and IFNγ expression in in vitro–differentiated murine 2D2 Th17 cells treated for the last 24 hours during induction with vehicle or S1P₁ agonist SEW2871 (gate: CD4 T cells; n = 4). (L–N) Concentration of cytokines as determined by multiplex ELISA in supernatants of SEW2871- or DMSO-treated in vitro–differentiated Th17 cells (n = 4). (O) Motor score determined in RAG1^–/– mice for 30 days after adoptive transfer of Th17 cells that were cultured for 24 hours in the presence of SEW2871 or vehicle DMSO (n = 8–9). Expression data are calculated as 2^{(Lod Ct – sample Ct)}; LoD Ct was set to 24. *P < 0.05,**P < 0.01; ***P < 0.001 (C–H and K–N, paired t test; O, multiple t tests).

Also, the S1P receptor subtype S1pr1 was negatively correlated to expression of pathogenic markers such as Rorc or Tnfa (Figure 4I). In splenic CD4 T cells, S1P₁ was suggested to facilitate Th17 differentiation (34, 35), but whether S1P₁ signaling modulates the fate of already-differentiated Th17 cells is unclear. We tested the effect of S1P₁ agonist SEW2871 on in vitro–differentiated Th17 cells and found that it significantly reduced the proportion of IL-17A^posIFNg^neg cells (Figure 4, J and K). This was accompanied by reduced production of chemokines such as CCL3 and CCL4, whereas IL-4 secretion was, in trend, increased (Figure 4, L–N). In line with a reduced pathogenicity of these cells, we found that adoptive transfer of SEW2871-treated 2D2-trangenic Th17 cells into RAG1^–/–-deficient recipients resulted in delayed development of EAE (Figure 4O). Together, these data suggest that subgroups of Th17 cells are characterized by specific GPCR expression patterns and that differentially expressed GPCRs may be targeted to modulate Th17 pathogenicity.

GPCR heterogeneity in myeloid cells. To investigate activation-dependent changes in GPCR heterogeneity in individual myeloid cells, we compared naive splenic Ly6C^high monocytes with CD11b-positive cells obtained at peak of disease from spinal cords of MOG_35-55-immunized mice (CD11b_sc). Some receptors showed reduced expression frequency in CD11b_sc, such as Emr4, Cd97, Ptgir, Adora2a, or Gpr18, but the majority of receptors showed increased expression (Figure 5, A and B). Also, in vitro differentiation of BM-derived macrophages resulted in a clear increase in GPCR expression (Figure 5C), but k-means cluster analysis showed that both M1- and M2-differentiated macrophages differed significantly from spinal cord–infiltrating CD11b_sc (Figure 5D and Supplemental Figure 5, A–D).

Figure 5

Single-cell GPCR expression in splenic monocytes versus spinal cord–infiltrating CD11b-positive cells (CD11b_sc). (A) Heat map of GPCR expression in CD11b⁺Ly6C^high monocytes and CD11b_sc cells (monocytes: 42 cells from 2 mice; CD11b_sc: 98 cells from 6 mice); horizontal bars on the right side visualize expression frequency (%). (B) Number of GPCRs expressed in individual monocytes and CD11b_sc cells. (C) Number of GPCRs expressed in monocytes and in vitro–differentiated M1- or M2-polarized BM-derived macrophages (BMDM_M1: 31 cells; BMDM_M2: 38 cells). (D) T-SNE map representation of transcriptome similarities of in vitro–differentiated macrophages (BMDM_M1, BMDM_M2) or in vivo–differentiated myeloid cells (Cd11b_sc). K-means cluster assignment is color-coded; cell type is indicated by symbol.

Further cluster analysis showed that the CD11b_sc population consisted of 2 subpopulations (Figure 6, A and B). Cluster 2 cells were positive for macrophage-associated genes such as Ccr2, H2-Ab1, and Il1b and showed an increased expression of GPCRs Gpr132, Fpr2, Cd97, Fpr1, and Cxcr4. Cluster 1 cells, in contrast, were positive for Cx3cr1, P2ry13, P2ry12, Gpr56, and Gpr34 (Figure 6B), indicative of a microglial population (36). Spearman’s rank correlation analysis confirmed the existence of a Cx3cr1-associated gene cluster predominantly expressed in microglia and a Ccr2-associated gene cluster predominantly expressed in macrophages (Figure 6C). However, k-means analyses also indicated that intermediate cell types exist: Ccr2^negCx3cr1^pos microglia contained a subpopulation of approximately 20% (cluster 2 in Figure 6D) that was characterized by an increased expression of GPCRs predominantly found in macrophages, such as Cxcr4, Cd97, Ccr1, or Il1b (Figure 6E).

Figure 6

Subgroups within spinal cord–infiltrating CD11b-positive cells (CD11b_sc). (A) T-SNE plot showing similarities between individual CD11b_sc and results of k-means cluster analysis (k-means cluster assignment indicated by color). (B) Heat map showing expression of selected genes in CD11b_sc after grouping according to clusters identified in A. (C) Representation of Spearman coefficients for the correlation of gene expression (width of connecting lines indicates strength of correlation). (D) T-SNE plot showing similarities between individual Ccr2^negCx3cr1^pos putative microglia and results of k-means cluster analysis (k-means cluster assignment indicated by color). (E) Expression of selected genes in Ccr2^negCx3cr1^pos after grouping according to clusters identified in D. Expression data are calculated as 2^{(Lod Ct – sample Ct)}; LoD Ct = 4. Function-defining genes are shown in blue. *P < 0.05; **P < 0.01; ***P < 0.001 (unpaired t test).

GPCR expression in individual resting and activated spinal cord ECs. Finally, we investigated activation-dependent changes in GPCR heterogeneity in individual spinal cord ECs from healthy mice (ECscN) and mice at peak of EAE (ECscEAE). Neuroinflammation resulted in an increased expression of mRNAs encoding adhesion molecules such as E-Selectin (Sele), Icam1, and vascular adhesion molecule 1 (Vcam1), as well as transcription factors Hif1a, Fos, and Egr1 (Figure 7A). In addition, several GPCRs showed upregulation, for example Lphn2, F2r, and Darc, while other GPCRs were downregulated in frequency and intensity, such as Gabbr1, Cd97, Gpr124, Cxcr7, Gpr107, and Lphn1 (Figure 7A). The total number of GPCRs expressed per individual cell did not differ between naive and inflamed ECs (ECscN: 14.2 ± 0.5; ECscEAE: 14.1 ± 0.5). K-means cluster analysis within the group of ECscEAE identified a subgroup of approximately 25% of cells that were characterized by high expression of Sele, Icam1, Vcam1, Fos, Egr1, Darc, and F2r, indicative of a highly activated EC population (Figure 7, B and C). These cells showed reduced expression of Gpr107, Cxcr7, Gabbr1, Gpr124, and others (Figure 7, B and C). Spearman’s rank correlation analysis confirmed the association between Vcam1, Icam1, Sele, Fos, Egr1, and Darc on the one hand and Klf2, Gpr56, Cd97, Gabbr1, and Lphn2 on the other hand (Figure 7D).

Figure 7

Single-cell GPCR expression in spinal cord endothelial cells from naive mice (ECscN) and mice at peak EAE disease (ECscEAE). (A) Heat map of GPCR expression in ECscN (33 from 5 mice) and ECscEAE (125 cells from 15 mice); horizontal bars on the right side visualize expression frequency (%). (B) T-SNE plot showing similarities between individual CD4scEAE and results of k-means cluster analysis (k-means cluster assignment is indicated by color). (C) Genes differentially expressed in cluster 1 cells. (D) Representation of Spearman coefficients for the correlation of gene expression (width of connecting lines indicates strength of correlation). (E and F) Transcriptional changes observed after knockdown of LPHN2 with 2 independent siRNAs (n = 4 experiments). (G and H) Adhesion of calcein-labeled THP1 monocytes to resting (G) or TNFα-activated (H) HUVEC. The proportion of calcein-positive adherent cells was determined after washing and trypsinization by flow cytometry and is expressed as percent of all trypsinized cells. Expression data are calculated as 2^{(Lod Ct – sample Ct)}; LoD Ct = 24.*P < 0.05; **P < 0.01; ***P < 0.001 (E and F, one sample t test; G and H, paired t test).

We finally investigated whether data from single-cell expression analysis can be used to develop new strategies to modulate EC activation. Lphn2, an orphan receptor with unknown function in EC, was the only GPCR upregulated in ECscEAE that was not associated with the inflamed gene cluster (Figure 7, B–D), suggesting it might play a role in maintaining/reestablishing quiescence in EC. We therefore investigated the consequences of LPHN2 knockdown in human umbilical vein endothelial cells (HUVEC). mRNA sequencing revealed that knockdown of LPHN2 resulted in an inflammatory activation of ECs that included upregulation of mRNAs of proinflammatory molecules such as IL-8, various chemokines and chemokine-binding proteins such as syndecans (encoded by SDC1 and -2, and CD44 (37) (Figure 7E). Also, production of proadhesive matrix protein fibronectin was increased, as well as expression of enzymes involved in the synthesis of proinflammatory mediators such as prostanoid-synthesizing cyclooxygenases 1 and 2 (encoded by PTGS1 and -2) or S1P-synthesizing sphingosine kinase (SPHK1) (Figure 7E). S1P degrading enzyme S1P phosphatase 1 (SGPP1), in contrast, was reduced, as was expression of endothelial NO synthase (encoded by NOS3) (Figure 7E). However, there was also transcriptional evidence for reduced adhesiveness — for example, a reduction of mRNAs for E- and P-selectin (SELE, SELP), as well as increased expression of endothelial inhibitors of leukocyte adhesion such as GDF15 (38) and KLF4 (39) (Figure 7F). To investigate the net effect of these changes on leukocyte-endothelial interaction, we studied adhesion of human THP1 monocytes to HUVEC after knockdown of LPHN2. We found that both under basal conditions and in TNFα-stimulated HUVEC, THP1 adhesion was increased (Figure 7, G and H), suggesting that orphan GPCR LPHN2 mediates antiinflammatory and antiadhesive effects in EC.

We report here the first single-cell GPCR expression analysis, to our knowledge, in primary immune cells and spinal cord ECs during neuroinflammation. Central goals of this study were to determine GPCR expression in individual cells of a given type under resting conditions and after activation, and to determine whether specific GPCR repertoires can be assigned to functional subgroups. We found that in spinal cord–infiltrating CD4 T cells, a subpopulation of quiescent Th17 cells exists that is characterized by an increased expression of receptors such as Gpr18, Gpr132, S1pr1, P2ry10, Cxcr3, Cxcr4, and that antagonism at CXCR4 or agonism at S1P₁ reduces pathogenicity of Th17 cells in vitro and in vivo. We furthermore found that orphan GPCR Lphn2 was in inflamed spinal cord ECs associated with a less activated phenotype, and we show that knockdown of this receptor enhances endothelial adhesiveness.

All cell types analyzed in this study showed a surprisingly high heterogeneity in GPCR expression, and a comparable degree of heterogeneity was found in a recent study employing mRNA sequencing in murine Th17 cells (40). Studies in other fields, mainly from developmental biology, suggest that this is not a specialty of GPCRs, but is also observed in other gene families (41, 42). However, due to the limited amount of RNA obtained from an individual cell, it cannot be excluded that — due to the technical limitations of single-cell expression analysis — very low abundance transcripts remain undetected, in particular if the RNA quality is compromised. Since GPCRs are, in many cases, expressed at very low levels (43), it is in this context important to consider the sensitivity of the detection system used for single-cell expression analysis. A recent study directly compared the 2 major readout systems for single-cell expression analysis, RT-PCR and RNA sequencing (RNA-Seq), with respect to their sensitivity for GPCR detection in murine aortic smooth muscle cells. Single-cell RT-PCR detected most GPCRs more efficiently than RNA-Seq, whereas mRNAs encoding cytoskeletal proteins were detected equally well by both methods (44). A possible explanation for this difference in sensitivity lies in the fact that single-cell RT-PCR uses target-specific preamplification of mRNAs, while preamplification for RNA-Seq is unbiased. We therefore conclude that, for projects addressing transcriptional differences in individual cells in a more general manner, RNA-Seq is superior because of the breadth of data obtained. Projects specifically addressing genes that are known to be expressed at low levels, in particular when dealing with transcriptionally rather quiescent cell types, benefit from the higher sensitivity of single-cell RT-PCR. In line with the notion that single-cell RT-PCR faithfully reports expression of low-abundance transcripts, we found that the expression frequencies observed on the transcriptional levels largely matched the frequencies determined in flow cytometric reporter analyses.

A comparison of murine and human peripheral blood CD4 T cells showed that the high degree of GPCR heterogeneity was present in both species, though murine CD4bl expressed less inflammatory chemokine receptors and showed higher expression of receptors associated with the naive T cell state (Ccr7, Vipr1, S1p1), most likely due to the specific pathogen–free husbandry conditions of experimental mice. However, even under specific pathogen–free conditions, certain fluctuations in the health status are to be expected, and this most likely underlies the observed variation in the GPCR repertoire of naive CD4 T cells from individual mice.

Inflammatory activation resulted in all cell types except ECs in an increased number of GPCRs per cell. This was not only true for chemokine receptors or receptors of inflammatory mediators such as thrombin, but also for receptors that are less well studied in T cells, such as Gpr183 or Gabbr1. The oxysterol receptor EBI2 (encoded by Gpr183) was previously shown to guide B cells during the germinal center reaction and to promote follicular T helper cell differentiation (45, 46). With respect to Gabbr1, various studies suggested immunomodulatory functions of GABA, but the inhibitory effect of GABA on T cells was, so far, mainly attributed to the GABA_A receptor subtype, a ligand-gated ion channel (47–49). The upregulation of Gabbr1 in CD4sc indicates that GABA_B may also contribute to these effects, suggesting that GABA_B agonist baclofen might have beneficial effects in MS treatment beyond management of spasticity. Despite the general upregulation of GPCR expression, 63 of 125 tested GPCRs remained undetectable in CD4sc, raising the question whether they are truly not expressed or just below the detection threshold in individual cells. Sixty of these GPCRs were also negative in NanoString bulk RNA analysis, indicating that these receptors are truly not expressed in CD4sc. The remaining 3 GPCRs were detected in bulk RNA analysis of CD4sc, but according to our single-cell analysis, they are only present in spinal cord myeloid cells: P2ry12, P2ry13, and C3ar1. Since our data showed that more than 25% of FACS-isolated spinal cord CD4 T cells are positive for myeloid marker Itgam, indicating they are actually doublets, we concluded that the detection of these GPCRs in bulk RNA analysis is an artifact due to myeloid contamination.

An important aim of this study was to investigate whether functionally distinct subgroups can be assigned to specific GPCR patterns, thereby allowing a more selective pharmacological modulation of these cells. In line with a previous single-cell expression study (40), we found that the population of spinal cord–infiltrating CD4 T cells consists not only of Treg, Th1, and Th17 cells, but also functionally diverse subgroups of Th17 cells. Our data show that Th17 cells with high expression of pathogenic marker genes such as TNF or GM-CSF express a GPCR repertoire that significantly differs from that of Th17 cells with low pathogenicity, and we provide evidence that pharmacological targeting of these GPCRs is able to modulate the pathogenicity of these cells.

Both our transcriptional and flow cytometric data show, for example, that expression of CXCR4, the receptor for chemokine SDF-1α/CXCL12, is associated with reduced pathogenicity in spinal cord–infiltrating CD4 cells. Previous studies suggested that CXCR4 regulates gene expression and migration in a number of cell types, including leukocytes, hematopoietic progenitor cells, and tumor cells (31, 32). In addition, CXCR4 impacts on thymic β-selection in T cells (50). Our data show that blockade of CXCR4 by AMD3100 reduces pathogenicity of Th17 cells, resulting in reduced production of proinflammatory mediators. As to how CXCR4 regulates Th17 activation, CXCR4 was shown to physically interact with the T cell receptor resulting in ERK phosphorylation, costimulation, and enhanced cytokine expression (33).

Another receptor associated with the more quiescent Th17 population was S1P₁, a well-known modulator of T cell trafficking (21) and efficient Th cell production (51). Activation of S1P₁ was shown to facilitate Th17 development during differentiation (34, 35), but our data show that, in mature Th17 cells, S1P₁ negatively regulates production of IL-17, CCL3/MIP-1a, and CCL4/MIP-1b, resulting in delayed onset of adoptive transfer EAE.

Although both are of myeloid origin, macrophages and microglia show characteristic differences in their GPCR repertoire. GPCRs that were predominantly or exclusively expressed in individual microglia are Cx3cr1, P2ry12, P2ry13, Gpr56, or Gpr34; these data are in line with a recent expression analysis in pooled cDNA from microglia (36). Within the activated microglial population, we identified a subpopulation of cells that were characterized by expression of markers normally associated with macrophage fate — for example Il1b, Ccr1, Cd97, or Cxcr4. These cells continued to express the GPCR repertoire usually observed in microglia, suggesting they were not simply freshly entered patrolling monocytes (CCR2⁻GR1⁻CX3CR1^hiLY6C⁻) (52). Instead, these cells might be either BM-derived monocytes/macrophages that undergo microglial transdifferentiation (53) or microglia that assumed a macrophage-like phenotype.

Inflammatory activation of ECs plays a crucial role in leukocyte transmigration across the blood-brain barrier (54). Spinal cord ECs harvested at peak of disease contained a subpopulation of approximately 30% that showed clear signs of inflammatory activation, and these cells were characterized by increased expression of GPCRs F2r and Darc. Chemokine GPCR DARC has previously been shown to be upregulated in brain ECs during neuroinflammation and to shuttle inflammatory chemokines across the blood-brain barrier (55). Activation of thrombin receptor PAR₁, encoded by F2r, has been show to increase endothelial permeability and to enhance surface expression of adhesion molecules, thereby promoting platelet and leukocyte margination (25). Furthermore, PAR₁ antagonism was shown to ameliorate EAE (56). Also, orphan receptor latrophilin-2 (encoded by Lphn2) was upregulated in ECsc derived from EAE mice, but in contrast to F2r and Darc, it was rather associated with a less-activated endothelial phenotype. The latrophilin family of adhesion receptors has been implicated in regulation of cell-cell and cell-matrix interactions, but data regarding specific cellular functions of LPHN2 are, except for a putative role in EMT, lacking (57, 58). Our knockdown experiments in human brain EC showed that loss of LPHN2 results in increased expression of mediators implicated in leukocyte adhesion and transmigration through the endothelial barrier, such as IL-8 (59), CXCL11 (60), CXCL12 (61), or chemokine binding proteins syndecan and CD44 (37). Also mRNA levels of synthesizing enzymes for inflammatory lipid mediators such as prostaglandins or S1P were increased, whereas endothelial NO synthase expression was reduced. In line with these transcriptional changes, endothelial adhesiveness for monocytes was increased in LPHN2-knockdown HUVEC. These findings indicate that inflammation-induced upregulation of LPHN2 is a counter-regulatory mechanism helping to restrict inflammatory activation.

Taken together, we show here that GPCR expression is surprisingly heterogeneous in the different cell populations implicated in neuroinflammation and that, within each cell type, functionally relevant subpopulations with specific GPCR patterns can be identified. We show in 3 examples how the knowledge about GPCR expression in functionally relevant subpopulations can be exploited to modulate the functional state of these cells; our findings, therefore, provide an ideal basis for the development of GPCR-based therapeutic strategies that selectively target pathogenic cell subpopulations.

Mouse lines. Female C57BL/6J and 2D2 (62) mice were purchased from The Jackson Laboratory and bred at the Max Planck Institute for Heart and Lung Research. RAG1^–/–, EBI2-FL-EGFP (63), and EBI2-FL-EGFP 2D2 mice were bred at the Institute for Molecular Medicine at the University Medical Center of the Johannes Gutenberg-University Mainz. All mouse lines were on a pure C57BL/6 background.

Induction of active EAE. For active EAE induction, mice were immunized s.c. with 250 μg MOG_35-55 myelin oligodendrocyte glycoprotein peptide (GenScript) emulsified in Complete Freund’s adjuvant (Becton Dickinson) containing 4 mg/ml of heat-inactivated Mycobacterium tuberculosis (H37Ra, Becton Dickinson). On days 0 and 2, mice received 500 ng pertussis toxin (Sigma-Aldrich) in PBS i.p. Clinical scoring of EAE was conducted as follows: 0, no clinical disease; 1, limp tail; 2, impaired righting reflex and gait; 3, partial hind limb paralysis; 4, tetraparesis; 5, dead. Draining lymph nodes were analyzed after mice developed first clinical symptoms (day 8–12 after immunization, clinical sore ≤ 1); spinal cord infiltrate was examined at peak of disease (day 14–17 after immunization, clinical score 2–3.5).

Induction of adoptive-transfer EAE. For adoptive-transfer EAE, mice were sacrificed on day 10 after active immunization without pertussis toxin, and cells from spleen and lymph nodes were isolated and cultured in cell culture medium containing 50 μg/ml MOG_35-55, 10 μg/ml anti-IFNγ antibodies (BioXCell) and 10 ng/ml IL-23 (Miltenyi) for 4 days at 37°C with 5% CO₂. Afterward, 2 × 10⁵ blasting Th17 cells were injected i.v. into RAG1^–/– mice, with injection of 200 ng pertussis toxin along with transfer and on day 2 after transfer. For adoptive transfer of in vitro–differentiated Th17 cells, differentiation was performed as described above, and 7 × 10⁶ cells were injected i.v. into RAG1^–/– mice; 200 ng pertussis toxin were injected i.p. on days 0 and 2.

Isolation of spinal cord leukocytes and ECs. Spinal cords were homogenized with a glass tissue homogenizer in PBS containing 1% glucose and 0.1% BSA. After centrifugation, spinal cords were resuspended in 6 ml of 30% Percoll (Sigma-Aldrich) and layered on a gradient consisting of 4 ml 45% Percoll and 2 ml 70% Percoll. Gradients were spun for 20 minutes (970 g, room temperature, without break), and interphases between the layers were harvested. After washing, cells were stained with fluorochrome-labeled antibodies, and desired leukocyte populations were purified by flow cytometry.

For the isolation of ECs, spinal cord homogenates were digested in the presence of 3.8 mg/ml collagenase II (Worthington), 0.6 U/ml dispase (Sigma-Aldrich), 1 U/ml DNase I (New England Biolabs) for 1 hour at 37°C under shaking. After centrifugation, digested spinal cords were resuspended in 25% BSA and centrifuged (2000 g, 20 minutes, 4°C), and the myelin layer on top of the suspension was sucked off. After washing once in PBS containing 0.5% BSA and 2 mM EDTA, ECs were stained with antibodies directed against CD45 (clone 30-F11, catalog number 12-0451-81, Thermo Fisher Scientific) and CD31 (clone 390, catalog number 25-0311-81, Thermo Fisher Scientific); CD45^-CD31⁺ ECs were purified by flow cytometry.

Isolation of human PBMC. For peripheral blood mononuclear cell (PBMC), isolation heparinized whole blood was layered over Ficoll (GE Healthcare) and centrifuged at 400 g for 30 minutes at room temperature without break. The PBMC containing interphase was washed 3 times with PBS and prepared for cell sorting.

Cell sorting and flow cytometric analysis. The following fluorochrome-labeled antibodies were used for cell sorting and flow cytometric analysis: rat anti–mouse CD4 (clone GK1.5, catalogs 25-0041-81 and 11-0041-82), mouse anti–human CD4 (clone OKT4, catalog 12-0048-41), hamster anti–mouse CD11c (clone N418, catalog 12-0114-81), rat anti–mouse CD197 (clone 4B12, catalog 48-1971-80), rat anti–mouse CD31 (clone 390, catalog 25-0311-81), rat anti–mouse CD45 (clone 30-F11, catalog 12-0451-81), mouse anti–human CD45RA (clone HI100, catalog 11-0458-41), mouse anti–human CD45RO (clone UCHL1, catalog 15-0457-41), rat anti–mouse CXCR4 (clone 2B11, catalog 12-9991-81), rat anti–mouse IgG2b κ (clone eB149/10H5, catalog 12-4031-81), mouse anti–human IFNγ (clone 4S.B3, catalog 48-7319-41), rat anti–mouse anti–IL-17A (clone eBio17B7, catalog 48-7177-80), mouse anti–human IL-17A (clone eBio64DEC17, catalog 17-7179-41), rat anti–mouse RORγt (clone B2D, catalog 17-6981-80), and rat anti–mouse TNFα (clone MP6-XT22, catalogs 25-7321-80 and 12-7321-81) (all eBioscience); mouse anti–human CD4 (clone OKT4, catalog 317413), hamster anti–mouse CD196 (clone 29-2L17, catalog 129815), rat anti–mouse CD198 (clone SA214G2, catalog 150313), rat anti–mouse F4/80 (clone BM8, catalog 123111), rat anti–mouse GM-CSF (clone MP1-22E9, catalog 505409), rat anti–mouse IFNγ (clone XMG1.2, catalog 505807), and rat anti–mouse Ly6C (clone HK1.4, catalog 128017) (all BioLegend); and rat anti–mouse CD11b (clone M1/70, catalog 561688 Becton Dickinson).

Single-cell suspensions were incubated with rat anti–mouse CD16/32 (clone 93, catalog 101301, BioLegend) for 10 minutes at room temperature to block Fc receptors, followed by staining with fluorochrome-labeled antibodies for 20 minutes at 4°C in the dark. After washing with staining buffer (PBS, 0.5% BSA, 2 mM), cells were either purified by cell sorting using a JSAN cell sorter (Bay Bioscience Corporation) or analyzed using a FACS Canto II (Becton Dickinson). Intracellular staining was performed as previously described (64). Flow cytometric data were analyzed by FlowJo software (Tree Star Inc.).

Cell culture. Murine BM-derived macrophages (BMDM) were generated by culturing BM-derived cells for 7 days in DMEM-F12 medium supplemented with 10% FCS and 100 U/ml penicillin/streptomycin, 2 mM glutamine (Gibco, Thermo Fisher Scientific) supplemented with murine recombinant 100 U/ml M-CSF (Peprotech). On day 7, M1-type BMDM were differentiated by adding 150 U/ml murine recombinant IFNγ (Peprotech) for 18 hours, together with 10 ng/ml LPS (InvivoGen) during the last 4 hours. For generation of M2 BMDM, 10 U/ml murine recombinant IL-4 (Peprotech) was added for 18 hours.

Murine primary T cells were cultured in RPMI1640 (Gibco, Thermo Fisher Scientific) containing 10% FCS, 100 U/ml penicillin/streptomycin, 2 mM glutamine, 1 mM pyruvate, nonessential amino acids, and 50 μM β-mercaptoethanol. To generate in vitro–differentiated Th17 cells for single-cell expression analysis, 2D2 splenocytes were stimulated for 3 days with 20 μg/ml MOG_35-55 (GenScript) in the presence of 10 ng/ml recombinant murine IL-12 (R&D Systems) for Th1 or with 2 ng/ml human recombinant TGF-β1, 5 ng/ml recombinant IL-6 (both Peprotech), and 20 ng/ml murine recombinanti–IL-23 (R&D Systems) for Th17 differentiation.

For Th17 in vitro differentiation in the presence of S1P₁ agonist SEW2871 (Sigma-Aldrich) or the CXCR4 antagonist AMD3100 (Sigma-Aldrich), naive CD4 T cells from 2D2 or C57BL6/J mice were isolated from spleen by MACS purification (Miltenyi Biotec) according to the manufacturer protocol. Cells were cultured at 1 × 10⁵ cells/well in 200 μl T cell medium (TCM) (RPMI with 10% FCS, 2 mM L-Gln, 100 U/ml penicillin/streptomycin, 1 mM sodium pyruvate, 50 μM 2-mercaptoethanol, 10 mM HEPES, and 1% nonessential amino acids [MEM]) in 96-well plates. For Th17 differentiation, 1 μg/ml anti-CD3 (clone 2C11, catalog 100313, BioLegend), 6 ng/ml anti-CD28 hamster anti–mouse CD28 (clone 37.51, catalog 102111, BioLegend) and 10 μg/ml rat anti–mouse IFNγ (clone XMG1.2, catalog 505811, BioLegend) with 2 ng/ml recombinant human TGF-β1 (R&D Systems), 5 ng/ml recombinant murine IL-6 (Promocell), and 20 ng/ml recombinant murine IL-23 (R&D Systems) were used. Cells were cultured at 37°C and 5% CO₂ for 5 days for Th17 differentiation. Th17 cell culture was supplemented with 5 μg/ml AMD3100 on day 0 or 10 μM SEW2871 (Sigma-Aldrich) on day 4.

For human Th17 cell induction, naive CD4⁺CD45RA⁺CD45RO^– were isolated from human PBMCs by flow cytometry. Cells (5 × 10⁵/ml) were cultured in X-VIVO 15 medium (Lonza) in the presence of 10 μg/ml plate-bound mouse anti–human CD3 (clone OKT3, catalog 317303, BioLegend), 1 μg/ml soluble mouse anti–human CD28 (clone CD28.2, catalog 302913, BioLegend), 5 ng/ml recombinant human TGF-β1, 10 ng/ml recombinant human IL-1β, 10 ng/ml recombinant human IL-6, and 10 ng/ml recombinant human IL-23 (all Peprotech) in the absence or presence of 5 μg/ml AMD3100 for 7 days.

For intracellular cytokine staining, cells were activated for 4 hours in TCM or X-VIVO 15 medium with 50 ng/ml PMA (Sigma-Aldrich), 500 ng/ml Ionomycin (Invitrogen), and 1× Monensin (eBioscience) at 37°C and 5% CO₂.

HUVEC were obtained from Lonza and cultured on dishes coated with 30 μg/ml Collagen (Corning Inc.) in EBM-2 medium supplemented with EGM-2 SingleQuot Kit Supplemental & Growth Factors (both Lonza). For endothelial activation, HUVEC were stimulated with 10 ng/ml human recombinant TNFα (Peprotech) for 10 hours.

THP1 cells were obtained from Sigma-Aldrich and cultured in RPMI1640 containing 10% FCS, 100 U/ml penicillin/streptomycin, and nonessential amino acids.

Single-cell quantitative PCR. For single-cell expression analysis, lymph nodes of healthy mice or mice at disease onset of MOG_35-55-induced active EAE were dissected, minced, filtered, and FAC-sorted for the corresponding lineage marker. Cells from spinal cords and human PBMCs were isolated as described above. After sorting, cells were loaded onto the microfluidic-C1 Single-Cell Auto Prep System (5–10 μm mRNA Arrays, Fluidigm) followed by RNA isolation and cDNA synthesis according to the manufacturer protocol. cDNA derived from visually empty chambers or chambers containing more than one cell was excluded from further analysis; also, cDNAs that showed poor expression of reference gene Hprt in conventional RT-PCR were excluded. High-throughput quantitative PCR (qPCR) of harvested cDNA from single cells was performed on 96.96 Dynamic Array IFC with a BioMark system (Fluidigm) using SsoFast EvaGreen Supermix low ROX (Bio-Rad) and Delta Gene primer assays (Fluidigm) as listed in Supplemental Tables 1 and 2 (murine and human, respectively). Only single-cell cDNAs negative for lineage markers Cd8 (CD8 T cells), Cd19 (B cells), Ly6g (neutrophils), Cdh1 (epithelial cells), and Myh11 (smooth muscle cells) but positive for reference gene Hprt and markers Cd4 (for CD4 T cells), Itgam (for CD11b-positive myeloid cells), or Cdh5 (ECs) were included in statistical analyses. Because C1-generated cDNA samples are potentially contaminated with genomic DNA, intron-spanning primer design was used for all genes. The limit of detection for the BioMark HD System has been estimated to be at a Ct value of 24 cycles (limit of detection [LoD] Ct); all sample Ct values were therefore subtracted from the LoD Ct using the formula: gene expression = 2^{(LoD Ct – sample Ct)} (65).

Bioinformatic analyses of single-cell RT-PCR. Heat maps of single-cell RT-PCR expression data (calculated using equation 2^{(Lod Ct – sample Ct)}) were generated in Perseus software (66). Further bioinformatic analyses were performed using customized functions from the RaceID R package (42). For k-means clustering, we utilized the clustexp and clustheatmap functions; cluster-specific genes were exported by a customized clustdiffgenes function. t-SNE plots were generated by the plottsne function of RaceID R package (42).

Graphical representation of bioinformatic data. In T-SNE plots, t-distributed stochastic neighbor embedding (t-SNE) is a nonlinear dimensionality reduction algorithm that maps multidimensional data (i.e., all expression data available for a given cell) to 2 (or more) dimensions suitable for human observation. In these plots, the similarity of individual cells (each dot represents 1 cell) is displayed as proximity — the closer together 2 cells are, the more similar they are. Clusters identified by k-means clustering are indicated by color; some graphs, furthermore, indicate different cell types by symbols.

In bar graphs showing fold enrichment/reduction in specific clusters, cluster specific genes (P < 0.05 according to expected transcript count probability from binomial testing; ref. 42) were exported by a customized clustdiffgenes function. Graphs show only fold changes (>1.5 or <0.7).

In Spearman’s rank correlation analysis, the correlation of expression between individual genes was determined by calculating the Spearman’s rank correlation coefficient in Perseus software (66). To visualize the strength of correlation between individual genes, the ExpressionCorrelation plugin in Cytoscape (67) was used. This plugin computes a similarity network for the tested genes in which the width of the line connecting 2 genes indicates the strength of the expression correlation (the thicker the connector, the higher the correlation). Only correlation coefficients > 0.3 are displayed.

siRNA transfection. siRNA for human LPHN2 (Hs_LPHN2_2, target sequence: AAGCTGGTTGACACTAATAAA; Hs_LPHN2_4, target sequence: TCGCCAAACAACAACATATAA) and AllStars negative control siRNA were obtained from Quiagen. HUVECs (8 × 10⁴/ml) were cultured in 24- or 96-well plates and transfected with 28.5 nM of the corresponding siRNA using RNAiMax lipofectamine (Invitrogen). Gene expression analysis and adhesion assays were conducted 82 hours after transfection.

Adhesion assays. To determine THP1 adhesion, THP1 cells were stained with 1 μM Calcein-AM (AAT-Bioquest) for 30 minutes at 37°C in the dark. Calcein-labelled THP1 cells (2 × 10⁵/ml) were added to siRNA-transfected HUVEC monolayers at 37°C for 15 minutes. In some cases, HUVEC were pretreated with 10 ng/ml TNFα for 10 hours. Subsequently, not-attached THP1 cells were removed by pipetting and washed 3 times with PBS, followed by trypsinization. The amount of Calcein-labeled THP1 cells was assessed by flow cytometry and normalized to a defined number of 10.1 μm AccuCount fluorescent particles (ACFP-100-3, Spherotech).

Multiplex bead array assay. Cytokines and chemokines in supernatant of in vitro–differentiated Th17 cells were determined according to the manufacturer protocol using a mouse cytokine/chemokine multiplex bead array (Millipore) and a Luminex MAGPIX analyzer.

RNA isolation for cell pool analysis. RNA for cell pool Nanostring analysis and RNA-Seq was isolated with the RNeasy Micro Kit (Qiagen) according to the corresponding manufacturer`s instructions. For RNA-Seq, DNA was digested on-column with DNaseI digestion (DNase-Free DNase Set, Qiagen) to avoid contamination by genomic DNA.

mRNA sequencing and expression analysis. RNA and library preparation integrity were verified with a BioAnalyzer 2100 (Agilent) or LabChip Gx Touch 24 (Perkin Elmer). Total RNA (4 μg) was used as input for Truseq Stranded mRNA Library preparation following the low sample protocol (Illumina). Sequencing was performed on the NextSeq500 instrument (Illumina) using v2 chemistry, resulting in minimum of 32 M reads per library with 2 × 75 bp paired-end setup. The resulting raw reads were assessed for quality, adapter content, and duplication rates with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Trimmomatic version 0.33 was employed to trim reads after a quality drop below a mean of Q18 in a window of 5 nucleotides (68). Only reads above 30 nucleotides were cleared for further analyses. Trimmed and filtered reads were aligned versus the Ensembl human genome version hg19 (GRCh37.p5) using STAR 2.4.0a with the parameter “--outFilterMismatchNoverLmax 0.1” to increase the maximum ratio of mismatches to mapped length to 10% (69). The number of reads aligning to genes was counted with featureCounts 1.4.5-p1 tool from the Subread package (70). Only reads mapping at least partially inside exons were admitted and aggregated per gene. Reads overlapping multiple genes or aligning to multiple regions were excluded. Differentially expressed genes were identified using DESeq2 version 1.62 (71). Only genes with a minimum fold change of ± 2 (log₂ ±1), a maximum Benjamini-Hochberg corrected P value of 0.05, and a minimum combined mean of 5 reads were deemed to be significantly differentially expressed. The Ensemble annotation was enriched with UniProt data (release 06.06.2014) based on Ensembl gene identifiers (Activities at the Universal Protein Resource [UniProt]).

NanoString analysis of GPCR expression in bulk RNA. NanoString analyses were performed as described previously (72, 73). In brief, 250–500 ng RNA from sorted cells was applied in a total volume of 30 μl in the assay. Barcodes were counted for about 1,150 fields of view per sample. Counts were first normalized to the geometric mean of the positive control spike counts; then, a background correction was done by subtracting the mean plus 2 SDs of the 8 negative control counts for each lane. Data were not normalized to reference genes because none of the reference genes showed sufficiently stable expression in all cell types according to the geNorm algorithm. Values that were < 20 were fixed to background level.

Statistics. Data are presented as means ± SEM if not otherwise indicated. Comparisons between 2 groups were performed using 2-tailed 2-sample t test or paired t test as indicated; normalized data (control group set to 1) were analyzed by one-sample t test. P values are indicated as follows: *P < 0.05, **P < 0.01, ***P < 0.001. Single-cell RNA-Seq data have been deposited in GEO database (https://www.ncbi.nlm.nih.gov/geo/), under the accession code GSE99778.

Study approval. All animal experiments were conducted in accordance with the corresponding institutional guidelines and permission of the states Hesse and Rhineland-Palatinate. Analyses in blood samples of healthy human donors were performed after obtaining written consent and approved by the local ethics committee (University of Frankfurt).