Mapping Rora expression in resting and activated CD4+ T cells

The transcription factor Rora has been shown to be important for the development of ILC2 and the regulation of ILC3, macrophages and Treg cells. Here we investigate the role of Rora across CD4+ T cells in general, but with an emphasis on Th2 cells, both in vitro as well as in the context of several in vivo type 2 infection models. We dissect the function of Rora using overexpression and a CD4-conditional Rora-knockout mouse, as well as a RORA-reporter mouse. We establish the importance of Rora in CD4+ T cells for controlling lung inflammation induced by Nippostrongylus brasiliensis infection, and have measured the effect on downstream genes using RNA-seq. Using a systematic stimulation screen of CD4+ T cells, coupled with RNA-seq, we identify upstream regulators of Rora, most importantly IL-33 and CCL7. Our data suggest that Rora is a negative regulator of the immune system, possibly through several downstream pathways, and is under control of the local microenvironment.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 IL-33 and CCL7. Our data suggest that Rora is a negative regulator of the immune system, possibly through several downstream pathways, and is under control of the local microenvironment.
EMBL Interdisciplinary Postdoctoral fellowship, supported by H2020 Marie Skłodowska Curie Actions and awarded to LHV, the Swedish Research Council in the form of a grant awarded to JH (#2016-06598), the European Research Council in the form of a grant awarded to SAT (ThDEFINE), FUV -Fondazione Umberto Veronesi in the form of funds awarded VP, UK Medical Research Council in the form of a grant awarded to JW and ANJM (U105178805), Wellcome Trust in the form of a grant awarded to JW and ANJM (100963/Z/13/Z) and a Single Cell Gene Expression Atlas grant awarded to ZM (108437/Z/15/Z), and the Open Targets Grant awarded to ZM (OTAR2067). Wellcome Sanger Institute core facilities are supported by a grant (WT206194). Microsoft Research Cambridge, a subsidiary of Microsoft, provided support in the form of a salary for JF. The specific roles of these authors are articulated in the 'author contributions' section. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests:
The authors have read the journal's policy and have the following competing interests: JF was an employee of Microsoft Research Cambridge at the time the study was conducted, but is no longer an employee of the company. EN is an employee of Aleph Labs, but was not employed by the company at the time the study was conducted. GK is an employee of AstraZeneca, but was not employed by the company at the time the study was conducted. This does not alter our adherence to PLOS ONE policies on sharing data and materials. There are no patents, products in development or marketed products associated with this research to declare. health were killed by a Schedule 1 method. Humane endpoints were determined utilising a tailored score sheet that incorporated observations of an animal's extremities, abdomen, coat condition, behavior and weight. No adverse events were experienced.
Cell sorting. Cells were sorted with a BD Influx Cell Sorter or SY3200 Synergy cell sorter (iCyt) and analysed on Fortessa or LSRII BD cell analyzers using FlowJo. Surface and intracellular stainings were carried out according to the eBioscience protocols. A full list of the antibodies used is provided in S3 File.
Histology. Lung lobes were fixed for 48 hours in 10% neutral-buffered formalin (Sigma), washed in PBS and transported to AML Laboratories (Florida) for embedding, sectioning and staining with haematoxylin and eosin.
Microscopy and image analysis. Images were captured using a VHX keyence microscope, using objective #2D1510064 at magnification "200" (2.1495um/px). The 2D stitching feature was used to capture the entirety of each lung section. Two sections were imaged from each mouse, with focus automatically calibrated for each section. All images were acquired in one seating, with the same settings for brightness.
Image analysis was done in R using the EBimage package [20]. In brief, pixels corresponding to background, cytoplasm and nuclei were manually annotated in one image. A classifier was set up to separate pixels into the 3 classes based on manual inspection of the RGB space.
To estimate emphysema, the distribution of inner holes was calculated by first detecting non-background area. A binary dilation operation was applied to filter out 1-2 pixel holes caused by improper tissue detection. Segmentation was then performed to find disjunct regions. The area of each region was calculated.
From the histogram H of hole sizes, we estimate the unevenness (indicating emphysema, compensating for overall compression of slides due to fixation) as S = H 99% /H 80% . Other ranges were tested and gave a similar outcome. Veins, surrounding space, and other artifacts are filtered out above the 1% top percentile, while the lower percentile reflects tissue shrinkage. Alternatively, the score can be interpreted as a measure of hole size unevenness/variance. To test for emphysema in the Rora KO we fit a linear model S~replicate + genotype, with genotype = 1 for Cre and +/Flox, 2 for Cre and Flox/Flox, and 0 otherwise (p = 0.018). A simpler ttest of +/+ vs Flox/Flox, not using any of the +/Flox samples, and ignoring replicate batch effects, gives p = 0.17. This is not significant but in the direction of the more inclusive linear model.
In vivo bulk RNA-seq. RNA was extracted with SPRI beads, followed by Smart-seq2. Read counts were estimated with Kallisto [21]. DE genes were estimated by DESeq2 [16] after counts had been rounded to the nearest integer count. Simple control vs treatment linear models were used. Counts and condition matrices are provided in S2 File.

Quantitative single-cell gene expression analysis
Single cell RT-qPCR measurement and analysis. Single-cell gene expression analysis was performed using BioMark 96.96 Dynamic Array platform (Fluidigm, San Francisco, CA) and TaqMan Gene Expression Assays (Applied Biosystems, Carlsbad, CA). Single cells were sorted into 5μl of CellsDirect reaction mix and immediately stored in -80C. Control wells containing no cells were included. On thawing, a mix containing 2.5μL gene specific 0.2x TaqMan gene expression assays (Applied Biosystems), 1.2 μL CellsDirect RT/Taq mix, and 0.3 μL TE buffer were added to each well. RT-PCR pre-amplification cycling conditions were: 50˚C, 15min; 95˚C, 2min; 22x(95˚C, 15s; 60˚C, 4min). Samples were diluted 1:5 in TE buffer and 6% were mixed with TaqMan Universal PCR Master Mix (Applied Biosystems). The sample mix and TaqMan assays were loaded separately into the wells of 96.96 Gene expression Dynamic Arrays (Fluidigm) in presence of appropriate loading reagents. The arrays were read in a Biomark analysis system (Fluidigm). ΔCt values were calculated in reference to the average of Atp5a1, Hprt1 and Ubc.
The expression of each gene was fit by a bivariate normal distribution (R mixtools package [22]) and a cut-off set at the average position between the lower and upper gaussian midpoints. Genes above this were considered expressed if above this cut-off. qPCR levels and metadata are provided in S1 File.
S. mansoni C1 single cell RNA-seq and analysis. For S. mansoni, three small (5-10 μm) C1 Single-Cell Auto Prep IFC chips (Fluidigm) were primed and 5000 cells were sorted directly into the chip. To allow estimation of technical variability, 1 μl of a 1:4000 dilution of ERCC (External RNA Controls Consortium) spike-in mix (Ambion, Life Technologies) was added to the lysis reagent. Cell capture sites were visually inspected one by one using a microscope. The capture sites that did not contain single cells were noted and were removed from downstream analysis. Reverse transcription and cDNA preamplification were performed using the SMARTer Ultra Low RNA kit (Clontech) and the Advantage 2 PCR kit according to the manufacturer's instructions on the C1 device. cDNA was harvested and diluted to 0.1-0.3 ng/μl and libraries were prepared in 96-well plates using a Nextera XT DNA Sample Preparation kit (Illumina) according to the protocol supplied by Fluidigm. Libraries were pooled and sequenced on an Illumina HiSeq2500 using paired-end 75-bp reads.
Salmon was used to estimate gene expression counts [23]. Poor quality libraries were eliminated using Scater [24] based on exonic and mitochondrial read counts. For all queries, a gene was considered expressed if the Log10 (1+normalized count) was above 0.5. Gene overlap was tested using Fisher's method. Counts and condition matrices are provided in S2 File.
N. brasiliensis C1 single cell RNA-seq. T cells from N. brasiliensis infected mice were captured with C1 analogously to the S. mansoni mice. The cells were sorted for CD3+CD4+SELL-, from mesLN, medLN, lung and spleens. This data is available but was not used due to low quality.
Single cell transcriptomes were generated by using the Smartseq-2 protocol [25]. Single cells were sorted into 96 well plate that contained 5μl of Triton-X lysis buffer, 1 μl of 10 μM olido-dT30-VN, 1 μl dNTP mix (25 mM each) and ERCC controls at a final dilution of 1:64m and immediately stored in -80˚C. cDNA was submitted for the amplification for 25 cycles on an Alpha Cycler 4 thermal cycler. Amplified cDNA went through two subsequent rounds of cleaning by using Agencourt AMPure XP beads (Beckman Coulter UK Ltd, High Wycombe, UK) at a 1.0x ratio on a liquid handler Zephyr G3 NGS Workstation (PerkinElmer) and eluted in 20ul of RNAse free water. Purified cDNA was subjected to quality control using 1 μl of cDNA on an Agilent 2100 BioAnalyser (Agilent Technologies, Santa Clara, CA, USA) using the Agilent High Sensitivity DNA kit. Samples were normalised to a concentration of 0.3 ng/ μl. Nextera libraries were prepared using Nextera XT DNA Sample Preparation kit (Illumina) according to the protocol supplied by Fluidigm. The library preparation has been done in combination with benchtop liquid handlers: Zephyr G3 NGS Workstation (PerkinElmer) and Mantis (Formulatrix). Purified pool of Nextera libraries was subjected to quality control using 1 μl of cDNA on an Agilent 2100 BioAnalyser using the Agilent High Sensitivity DNA kit.
Salmon v0.8.2 [23] was used to map the reads to ensembl cDNA reference GRCm38 to quantify the gene expression. Counts were used for all the analyses. Low quality cells with >15% mitochondrial genes or <1500 total genes were filtered out. Seurat [26] was used for single-cell analysis: Highly variable genes were calculated (x.low.cutoff = 0, x.high.cutoff = 5, y. cutoff = 0.1), PCA performed on these (pcs.compute = 50), and t-SNE on the first 25 PCA components. Finally clustering was done (reduction.type = "pca", dims.use = 1:25, resolution = 0.5). Common marker genes were plotted on top of the clusters and the cluster identities were assigned manually (Gata3, Tbx21, Foxp3, Sell, and Rorc). The TCR sequences for each single T cell were assembled using TraCeR [27] which allowed the reconstruction of the TCRs from scRNA-seq data and their expression abundance, in addition to identification of the size, diversity and lineage relation of clonal subpopulations. Cells for which more than two alpha or beta chains were identified were excluded from further analysis. iNKT cells were detected by their characteristic TCRA gene segments (TRAV11-TRAJ18). Only CD4+ T cells were kept in the final plot; other cells remain annotated in the ArrayExpress submission. Counts and condition matrices are provided in S2 File.
For the repeated screen, culturing was done in the same manner. The cells were then stained with PI and anti-CD4 eFluor 660 (eBioscience 50-0041-82, clone GK1.5) for 30 minutes in IMDM media. The cells were then FACS sorted to retain live CD4+ T cells, on average 5000 cells per well.
During Th0 conditions, it is possible for the cells to become somewhat Th1 polarized; we have looked at markers in the data but can neither confirm nor disprove that this is the case for these experiments.
RNA-sequencing. For the first screen, RNA was extracted using Agencourt AMPure XP beads (Catalog #A63881) at 1.5:1 ratio into 10ul elution volume. RNA-seq was done following Smart-seq2 [25] but with the following changes: The input was increased to 7ul and Mix 1 consisted of 2ul 10uM dNTP and 2.5ul 10uM oligo-dT. Mix 2 consisted of 5ul SmartScribe buffer, 5ul 5M Betaine, 1.25ul SmartScribe, 1.25ul (100uM) DTT, 0.15ul 1M MgCl2, 0.38ul 100uM TSO, and 0.47ul SUPERase inhibitor. SPRI cleanup was done at ratio 0.8:1 following firststrand synthesis, eluted into 12ul NFW. Mix 3 consisted of 12.5ul KAPA HiFi 2x master mix and 0.5ul 10uM ISPCR. PCR was done with 10 cycles. In the second screen, RNA-seq was done similarly with the following modification: No SPRI-cleanup was done after first-strand synthesis. Instead a larger amount of mix 3 was added directly.
On average 0.2ng cDNA was then used for Nextera XT library preparation, done at 25% of volume recommended by manufacturers protocol. Sequencing was done with 2 lanes of Illumina 50bp PE HiSeq 2500.
Upstream analysis. Salmon was used to estimate gene expression counts. Poor quality libraries were eliminated using Scater based on exonic and mitochondrial read counts. Genes with average normalized counts below 0.3 were removed. The increase/decrease in Rora expression level was computed by DEseq2 [16] using a linear model: normalized countst reatment + mouseReplicate, where treatment is a factorial variable. The correlation heatmap between treatments was calculated as cor(cor(FC)). The first correlation is used as a dimensionality reduction while the second correlation is used to remove the technical bias toward positive correlation. Clustering was done using R hclust, standard parameters.
The CellTrace FACS data was analyzed with FACSanadu [28] to obtain cells in each division stage. The cell counts were analyzed across replicates and compared to controls using R.

Metabolite LC-MS
Metabolites were extracted using a modification of a previous method [29]. Briefly, roughly 200k cells were resuspended in chloroform/methanol (2:1). After mixing samples, water was added (300 μL) and centrifuged at 13,000 x g for 20 min. The organic layer (lower fraction) was collected and dried under a stream of nitrogen. The organic fraction was reconstituted in chloroform/methanol (1:1; 100 μL) and 10 μL were mixed in 90 μL IPA/acetonitrile/water (2:1:1). The samples were analysed by LC-MS using an LTQ Orbitrap Elite Mass Spectrometer (Thermo Scientific). Sample (5 μL) was injected onto a C18 packed-tip column (75μm x 100mm, Thermo Scientific) at 55˚C. The mobile phase A was acetonitrile/water 60:40, 10 mM ammonium formate whilst B was LC-MS-grade acetonitrile/isopropanol 10:90, 10 mmol/L ammonium formate. In negative ion mode, ammonium acetate was used to aid ionisation. In both positive and negative ion mode, the gradient used (S5 File) was run at a flow rate 0.5 mL/ min. Data were acquired with a mass range of 100-2000 m/z. Chromatograms were converted to mzML format and annotated using an in-house R script. The annotated metabolites were size factor normalized to enable comparison. PCA was performed on the in vitro and in vivo samples separately (S1 Fig in S1 File). Lungs and spleens from in vivo samples appear to cluster together, suggesting no need to correct for tissue differences. In vitro samples were compared with Limma [30], with all naive cells (negatively selected as CD8-CD11b-CD11c-CD19-CD24-CD25-CD44-CD45R-CD49b-TCRγ/δ-TER119-) compared to all day 6 samples Th0/1/2. For the in vivo samples, we used the simplest model of~1+isko. Bayesian error correction was performed.

Rora is expressed in activated CD4+ T helper cells in vivo
We wanted to get an overview of Rora expression in CD4+ T helper cells during worm infection. As most of the mouse models used in this study are type 2 immunity, we focus more specifically in Th2 cells. Mice were infected with Nippostrongylus brasiliensis (N. brasiliensis), a tissue migrating parasitic roundworm of rodents. N. brasiliensis initiates infection at the skin, the site of invasion, and then migrates via the circulatory system to the lungs. Through exploitation of host lung clearance, worms are transported to the intestine, from where they are expelled. In total, this infection is resolved within 7-10 days (Fig 1).
To investigate Rora activity in activated Th cells, we examined Rora expression using a Rora +/teal reporter mouse, in which the gene encoding Teal fluorescent protein has been inserted at the 3' end of the Rora coding sequence and expressed from a self-cleaving T2A site [14]. Flow cytometric analysis of reporter gene expression indicated that Rora was expressed in a proportion of CD4+ T cells isolated from medLN, spleen, lung and siLP (Fig 2C). We noted that all Rora expressing cells possessed an activated phenotype, as evidenced by their expression of CD44, ICOS and CD38 and down-regulation of SELL ( Fig 2D). Accordingly, Rora +/teal T cells were therefore significantly enriched within the activated T cell subset, comprising 40-60% of this population, but absent from the pool of naive cells (CD4+-CD25-CD44-SELL+) (Fig 2E). This correlation of Rora expression with T cell activation is in contrast with RNA-seq analysis made during the first 72h of in vitro culture of Th0 and Th2 cells, as previously reported (Fig 2F) [15]. There, in both mouse and human, Rora is expressed in naive cells, with a sudden increase around 2-4 hours of activation, before a gradual decrease in expression (overall slightly higher expression in Th2 over Th0). This agrees with another in vitro activation dataset where the Rora expression level is unchanged and low, in both Th17 and Th1 cells, following activation [31]. To compare with the in vivo situation, we generated RNA-seq data from CD4+ T cells, isolated from the spleens of N. brasiliensis infected wildtype mice 1, 3 and 10 days after infection, and overlapped with the in vitro data (Fig 2F). Here Rora is seen to increase, reaching a high and fairly sustained level around day 3. This is in agreement with a previous single-cell analysis of Th1/Tfh cells during malaria infection where Rora expression is seen to increase over time in Th1 [32] (S2A, S2B Fig in S1 File). It is also consistent with a study on atopic dermatitis [10], a study on Tregs [33] (S2C Fig in S1 File), as well as on T cells in a melanoma model [34] (S2D Fig in S1 File). Taken together, as Rora is consistently present in vivo but not in vitro, these 6 datasets suggest that Rora is induced not primarily by activation, but possibly by an external cue from the micro-environment.
Because the particular (micro-)environments might influence Rora expression, we examined our scRT-qPCR data across different organs (Fig 2B). Overall, the expression of Rora was relatively constant in CD3+CD4+ cells 3-7 days after infection (Fig 2F). However, in terms of tissue, it was expressed more frequently in the siLP and lungs, where 37% of cells expressed Rora, compared to 20% of the cells in lymph nodes. This pattern was confirmed using flow cytometry, where 10.9±1.9% of CD4+ cells in the lung were Rora +/teal , 10.1±1.1% from the spleen, 27.9±5.8% from the small intestine lamina propria and only 4.0±0.6% from medLN ( Fig 2D). This tissue distribution likely reflects the relative frequency of activated Th cells in each location.
Rora has so far been seen in activated cells, partly overlapping with primary markers of Th2 and Treg (Gata3 and Foxp3). To fully confirm the presence of Rora among different T helper cell types, as defined by the full transcriptional program, we used scRNA-seq on lymphoid and non-lymphoid tissue, isolating CD3+CD4+SELL-cells from the lungs 30 days after infection (N. brasiliensis infected wild-type mice 1, 3 and 10 days after infection, and overlapped with the in vitro data ( Fig 2G). As reference points we included CD3+CD4+Rora teal + and CD3 +CD4+Rora teal -cells from N. brasiliensis infected Rora +/teal reporter mice. Lungs, spleens, medLN and mesLN were also taken from infected and uninfected Rora +/teal mice, 7 days after infection. The cells were clustered and each cluster annotated using common markers (Sell, Ifng, Gata3, Foxp3) as Naive, Th1, Th2 and Treg cells (Fig 2H, all cells are shown in S2E Fig in  S1 File). Cells from the lymph nodes and spleen are mainly Naive or Treg, otherwise the cell origin and T cell fate appear fairly uncorrelated. The distribution of Rora expression level across the different CD4+ T cell types was investigated (Fig 2H and 2I). Again, we note that Rora is expressed among all activated T cells (however, less induced in Th1).
To confirm that this distribution of Rora expression is not specific to the response to N. brasiliensis infection we compared it to another worm model, Schistosoma mansoni (S. mansoni). CD3+CD4+SELL-cells were collected 6 weeks after infection from spleens and mesLN (S2F Fig in S1 File). Single cell RNA-seq was performed and analysis shows that Rora is expressed in 59% of the activated cells (expressed defined as greater than zero counts). Approximately all (98%) cells that expressed Rora also expressed any of the key regulators Gata3, Foxp3 or Tbx21.
To conclude, we have performed FACS, scRT-qPCR and scRNA-seq, and confirmed that Rora is generally expressed in a proportion of activated T cells during N. brasiliensis infection. Rora induction is likely to depend on a local environmental cue not present in the usual CD3/ CD28 in vitro activation system. Thus an in vivo system is required to study the physiological role of Rora.

Rora is expressed in activated T helper cells during immune responses
To further corroborate these findings, we also investigated a non-worm immunity model, Ragweed pollen (RWP), a common allergen which causes lung inflammation [35]. The Rora +/teal reporter mice were exposed to intranasal RWP or PBS (Phosphate Buffered Saline, a negative control administration), administered every 2-3 days for 2 or 6 weeks ( Fig 3A). As anticipated, RWP exposure resulted in increased CD4+ T cell counts in the lung and mediastinal lymph node (medLN) after 2 weeks as compared to PBS controls, and this was maintained after 6 weeks of continuous RWP administration (Fig 3A). In lung, but not medLN, this increase in CD4+ T cell frequency was accompanied by an elevated proportion of T cells expressing Rora, increasing from 30% to 40% between 2 and 6 weeks of RWP administration (Fig 3B and S3 Fig  in S1 File). Consistent with our findings in untreated mice, Rora +/teal T cells elicited by RWP stimulation had an activated phenotype, characterised by a CD44+SELL-surface profile (S3 Fig in S1 File). Notably, the Rora +/teal T cell subset was enriched for cells that had acquired effector function, identified by their expression of the cytokines IFNg, IL-13 and IL-17A, or the Treg-associated transcription factor, FoxP3 (Fig 3C). Here again we note that Rora expressing cells include diverse activated cell phenotypes, that is, Th1, Th2 and Th17.
To understand the involvement of Rora in effector cells, we examined its expression in cytokine-expressing cells. From the N. brasiliensis model single-cell RT-qPCR, 83% of the cells that expressed either IL-4, IL-13 or IL-10, also expressed Rora (Fig 2B). The corresponding number in the S. mansoni-infected mice is 71% (S2F Fig in S1 File, scRNA-seq).
To investigate whether Rora was expressed by Th cells responding to secondary immune challenge, we challenged Rora +/teal and control mice with 2W1S peptide intranasally, in conjunction with the protease allergen papain, which promotes a type-2 immune response ( Fig  3D). After secondary immunisation with 2W1S/papain 3 weeks after the initial stimulation, we were able to identify 2W1S-responsive CD4+ T cells by tetramer staining [36] and FACS analysis. Tetramer+ CD4+ T cells in both the lung and medLN were largely CD44+ and Rora teal + (Fig 3E), indicating that recently activated Th cells express Rora. As with the N. brasiliensis infection, in the S. mansoni model we found significant overlap of Foxp3 and Rora in activated T cells (S2F Fig in S1 File, Fisher test, p = 2 � 10 −4 ). A higher fraction of FOXP3+ cells was also found in sorted Rora +/teal relative to RORA-cells, in the RWP model (Fig 3C). The gap is increased when the exposure to the pollen was longer. The same could also be seen after mice were exposed to 2W1S+papain. The Rora +/teal population contained 20% of the FOXP3+ cells, while the RORA-, only 4% of the cells. Overall our data show that Rora is expressed in Th cells in worm immune models as well as in allergen type 2 immunity. We find Rora transcripts present in activated Th cells, not only in Tregs, but also in Th1 and Th2, as well as activated CD4+ T cells responding to a secondary immune challenge. Rora is also more correlated with cytokine-releasing cells, which suggests a function in effector cells.

Rora KO affects inflammation severity
To isolate the in vivo physiological function of Rora in Th cells, we created a conditional knockout mouse line, where exon 4 of the Rora gene was deleted only in CD4+ cells (CD4 cre Rora fl/fl , verified by RNA-seq, S3A Fig in S1 File), resulting in a dysfunctional RORA protein.
We noted no change in Th cell number when comparing CD4 cre Rora fl/fl to control (S3B Fig in  S1 File). We then infected CD4 cre Rora fl/fl mice with N. brasiliensis and examined lung inflammation after 30 days. An unbiased image analysis algorithm was developed to detect emphysema and score stained lung sections (Fig 4A). Generally, uninfected samples score the lowest for emphysema, as expected (Fig 4B and 4C). During infection, single-allele Rora KO (CD4 cre Rora +/fl ) had a higher score of emphysema, and double-allele KO (CD4 cre Rora fl/fl ) even higher score (linear model fit, p = 0.018). We detected no differences in our negative controls, no-Cre KO mice (Rora fl/fl ), and WT with Cre KO mice (CD4 cre ). As emphysema can be caused by eosinophil infiltration, in analogy with a previous model [13] we analyzed eosinophil infiltration rate and noted an increase but not statistically significant difference (S3C Fig in S1  File).
We conclude that Rora expression in CD4+ T cells modestly promotes lung inflammation associated with N. brasiliensis infection. This may reflect a role for Rora in regulating the activity of T helper cells, and may also result from the modulation of Treg activity, as has previously been reported for Rora+ Tregs in skin [13].

Rora overexpression affects several key immune regulatory genes
To probe the role of Rora in T cells, we next proceeded to identify the molecular targets of Rora in Th2 cells by means of retroviral overexpression. Th2 cells were induced in vitro, infected the day after induction, and on day 5, gene expression was compared by RNA-seq (empty viral vector vs overexpression vector, Fig 5A and 5B). The in vitro differentially expressed (DE) genes were compared to the genes that are differentially expressed between skin Tregs from wt vs CD4 cre Rora fl/fl mice (Fig 5C). Consistent with a previous report [13], we see a strong effect on tumor necrosis factor receptor superfamily member 25 (Tnfrsf25), which has been shown to be required for Th2 effector function, e.g., allergic lung inflammation [37]. The receptors C-C chemokine receptor type 2 (Ccr2), Ccr5, and Tnfrsf23 are in agreement, although with lower fold changes. We also find Ninjurin 2 (Ninj2) which has recently been implicated in endothelial inflammation [38].
We also looked at a previous dataset of Rora small interfering RNA (siRNA)-treated human Th17 cells [39] and find that Tnfrsf25 is DE (p = 2.2 � 10 −4 , DESeq2). In another study of Tregs during steady-state in mouse [40], we find that Tnfrsf25 is higher in colon Tregs than in lymphoid tissue Tregs (p = 2 � 10 −2 , DESeq2). Thus four datasets, in different T helper cell types, confirm Tnfrsf25 as a downstream gene.
Overall the two in vitro systems strengthen the claim of Rora as regulator of Tnfrsf25. However the expression of Tnfrsf25 across Th2, Treg and Naive, but not Th1 (Fig 2I), shows that the function of Rora is not limited to Treg cells.

Rora affects activation in vivo
To see if Rora has additional functions in vivo, we generated bulk RNA-seq from activated (SELL-/CD44+) CD4+ T cells from control and CD4 cre Rora fl/fl mice from lung, spleen, siLP and colon (30 days after N. brasiliensis infection), and compared the DE genes between the different tissues, the overexpression DE genes, and the previous Treg Rora KO data. Volcano plots are shown in Fig 5E (full list of fold changes in S2 File). From manual comparisons of these conditions, we have tried to find DE genes in common between the tissues, and if these in turn can explain the function of Rora. We generally do not see any DE cytokines of high statistical significance, but highlight some other types of genes here. We find the transcription factor Aryl hydrocarbon receptor nuclear translocator-like (Arntl) to be differentially expressed in activated CD4+ T cells from siLP as well as in the overexpression analysis and in Tregs in skin (Fig 5B and 5C). Arntl is involved in the circadian responsiveness, and is expressed in CD4+ T cells [41]. The Arntl promoter contains ROR elements (RORE), and was found to be regulated by Rora in mice [42]. RORA promotes Arntl expression and thereby maintains the circadian rhythm. We find it to be expressed in all cell types (Fig 2I). The S100 calcium-binding protein A4 (S100a4) is another gene which is differentially expressed in activated T cells from lung, colon, and to lesser extent in Tregs in skin. In the single cell dataset, it is expressed mostly in Tregs but generally in Th1/Th2 (Fig 2H). It exists in both intra-and extracellular forms and can activate NF-κB [43]. Using an antibody, S100a4 has been shown to mediate T cell accumulation, and Th1/Th2 polarization, at tumour sites [44]. Thus S100a4 might act downstream of Rora to control T cell activation.
C-X-C chemokine receptor type 6 (Cxcr6) is differentially expressed in siLP, spleen and in Tregs in skin, and is consistently regulated in other tissues. Cxcr6 is expressed in all Th cells but mostly in Th2 cells (Fig 2I). It might be involved in Th cell activation [45,46].
The sialyltransferase St6galnac3 is downregulated in colon and overexpression (Fig 5B and  5F). This enzyme transfers sialic acids to glycolipids and glycoproteins. Along with the previously DE gene Alox8 (Fig 2I), we wondered if Rora may have any function on lipid metabolism. RORA has already been found to bind cholesterol [47] and is evolutionarily related to the nuclear receptor and lipid metabolism gene Peroxisome proliferator-activated receptor gamma (Pparg) [48]. We performed lipid LC-MS on in vivo CD4 cre Rora fl/fl and WT cells from N. brasiliensis infected mice, as well as on in vitro WT naive and mature Th0/Th1/Th2 cells for reference (S4 Fig in S1 File, cells from negative magnetic bead selection as described in methods). The CD4 cre Rora fl/fl presented small differences, and little overlap with the 85 lipids changed in vitro the first 6 days (cutoff p = 0.05). However, all the top DE lipids are phosphocholines or phospoethanolamines, and some phosphocholines can induce a type 2 immune response [49]. All the differentially expressed lipids are listed in S5 File.
In the spleen, Sell is downregulated in the KO (however it is still expressed as other cells are included as well), and opposite in the previous overexpression. This fits the idea of Rora driving activation and exit from the spleen, but why is it not DE in the other organs? We speculate that once a T cell has left the spleen or a lymph node, Rora might have already filled its purpose, and thus the effect on Sell will be smaller in any peripheral tissue.
We further noted the gene Emb (Embigin) in which is significantly DE in siLP, and in the same direction in overexpression, lung and spleen. It was also strongly DE in the previous skin Treg dataset. It has been shown to regulate cell motility in pancreatic cancer, and be controlled through the TGF-β pathway [50]. It is possible that it could have a similar role in T cells.
Overall we see that Rora affects several genes in vivo, that are differentially expressed between different T helper cell types. Several of the downstreams genes may explain the effect of Rora, including for example Arntl, Cxcr6, St6galnac3, Emb and Sell. However, Rora abbrogration does not abolish T cells. It is possible that Rora acts as a helper factor, or tissue residence checkpoint, for T cell activation and migration.

Rora expression can be extrinsically regulated
A remaining question is which genes regulate Rora, and whether regulation is solely due to the T cell itself or from an external signal. More insights into when Rora is induced would provide understanding of its role. As previously seen, Rora is generally not expressed during in vitro activation, but is upregulated during in vivo activation. This discrepancy may be explained by an unknown signal from the environment that is present in vivo but not in vitro. To find candidate upregulators we screened for in vitro upregulated genes during 39 treatments, including cytokines and chemokines, with 15 unstimulated controls, in 3 biological replicates. In short, naive cells were taken from wild type mouse spleens, plated on CD3/CD28 coated plates, with addition of test cytokines/chemokines 24 hours later. 5 days after plating the cells were measured by bulk RNA-seq (Fig 6A).
The stimulants and their effects are shown in Fig 6B (DE genes in S2 File). To validate our perturbation assay we looked at known marker genes. Of the genes having the strongest effect we find the expected genes: IL4 for Gata3 (Th2), Transforming growth factor beta 1 (TGFb1) for Foxp3 (Treg), IL6 and partially TGFb1 and IL21 for Rorc (Th17) [51] (S5a- S5c Fig in S1 File). We then looked for the effect on Rora expression. The most upregulating stimulants are CCL5, IL15 and IL33.
CCL5 is chemotactic for T cells, eosinophils, and basophils. CCL5 signaling through CCR3 has been reported to regulate Th2 (IL4+CD4+) cellular responses to promote metastasis of luminal breast cancers [52]. Interestingly also CCR5, another CCL5 receptor, was upregulated during overexpression, implicating a potential autocrine role for this factor. In vitro, CCR3 and CCR5 are expressed primarily in Th1, Th2 and Th17 [31].
IL15 is produced by non-lymphoid cells and important for the survival of several lymphoid subsets [53]. IL15 KO mice are markedly lymphopenic due to decreased proliferation and decreasing homing to peripheral lymph nodes [54]. Increase in IL15 thus fits our observation of Rora expression primarily in activated Th cells outside the lymph nodes. Il15ra is expressed in all murine CD4+ T cell subsets [31].
IL6 resulted in upregulation of Rora (statistical test significant if both batches of IL6 are pooled), consistent with previous IL6 in vitro culture data [31]. IL6 has several roles, beyond generating Th17 cells in vitro, both pro-and anti-inflammatory [55]. Il6ra is expressed in all CD4+ T cell subsets [31], though in our scRNA-seq data less in Th1 (Fig 2I).
IL33 also induces Rora, and is mainly associated with Th2 or ILC2 cells [56]. This agrees with our in vivo scRNA-seq where Rora is slightly elevated in Th2 cells compared to other Th subtypes (Fig 2I). Our scRNA-seq shows that the corresponding receptor Il1rl1 is expressed in all CD4+ T cell subsets but primarily in Th2 and Treg (Fig 2I).
Ccl7 upregulates Rora and can do so through several receptors; these include for example CCR1, CCR2, CCR3, CCR5, and CCR10 [57]. All of these receptors are present in CD4+ T cells to some extent, but Ccr5 is somewhat higher expressed than the rest [31].
The most downregulating stimulants are C-C motif chemokine 22 (CCL22) and stromal cell-derived factor 1 (SDF1a). CCL22 is secreted by dendritic cells and macrophages, and is suggested to induce Treg cell infiltration into the pleural space in patients with malignant effusion [58,59]. In addition, its receptor CCR4 is required for CD4+ T cell migration to the skin [60]. In bulk in vitro data it is present in all CD4+ T cell subsets, especially after activation [31].
SDF1a (also known as CXCL12), a ligand for the chemokine receptor CXCR4, promotes bone marrow homing for T cells. It has been shown that Naive T cells downregulate CXCR4 to avoid this [61]. A further look at the Human Protein Atlas shows the tissues with expression of SDF1a to be particularly high in spleen, cervix/endometrium, and adipose tissue. In our study we compared with lung and gut, which have less Cxcl12 expressed. This is consistent with our screening data, suggesting SDF1a may be a negative regulator of Rora. Cxcr4 is expressed in all CD4+ T cell subsets at roughly equal level [31] according to bulk and our scRNA-seq data ( Fig  2I).
To further validate if a treatment induces Rora, we also checked if a treatment induces the same downstream genes as in the Rora overexpression experiment. To do this we took genes from the Rora overexpression with a Log2 fold change > 2, and correlated the fold change with the corresponding fold change in a treatment (Fig 6C). Many of the most Rora influencing treatments also score similarly by this score. In particular, IL33 and CCL7 have a high level of agreement, followed by CCL22 (replicate a). Other treatments appear to have a less direct effect on Rora, for example IL6.
Based on our and previous data we conclude that Rora expression is mainly associated with an inflammatory state. Several cytokines may influence the expression of Rora. The effect on other genes can be studied using a similar screening approach, based on our supplementary data.

Conclusions
We have here investigated the impact of Rora on different subsets of CD4+ T cells, with an emphasis on Th2 cells, in several in vivo immune models. Rora is already known to be important for Th17 identity [12], and more recently it was discovered to regulate the function of Tregs in the skin [13]. Here we show that Rora is actually expressed in all activated T helper cells during different types of inflammation, including worm infections and allergy. We expand on the previous literature and suggest that Rora plays a wider role in the regulation of the immune response, in many T helper cell types-Th1, Th2 and secondary activated cells as well. Its function is likely to also involve further key downstream genes than just the previously shown Tnfrsf25. Our model of the Rora regulatory network is depicted in Fig 7. Rora deletion did not alter the overall CD4+ T cell numbers, and therefore does not seem to be involved in T cell proliferation. This was also shown before in a study where Th2 cells were examined in RORα-deficient Staggerer mice [62]. In this study we found that Th2 cells from the RORα-deficient mice were still able to differentiate and proliferate. However, we noticed higher inflammation in the lungs of CD4 cre Rora fl/fl mice, after N. brasiliensis infection, indicating a functional role in CD4+ T cells. A similar effect was shown in skin, when Rora was deleted in Tregs only [13]. The effect was demonstrated to be through Treg Tnfrsf25 signaling. Tnfrsf25 has been shown to be required also for Th2 cell function in lung inflammation [37]. Our analysis of the DE genes of activated Th cells, in the lung and colon of worm infected mice, as well as overexpression, and reanalysis of previously published data, confirms the importance of Tnfrsf25 in Tregs, but that its origin may include also other activated T helper cell subsets.
A higher fraction of activated T helper cells (both Th and Tregs) contained Rora, when located in non-lymphoid tissues (lung and siLP) compared to lymph node tissues. The same pattern was previously shown for skin-located versus LN-located Tregs [13]. The enrichment of Rora +/teal cells in peripheral tissues may reflect its correlation with T cell activation and effector function, or alternatively may relate to a role in influencing T cell migration. Our data is also unable to conclusively show if Rora is expressed before or after cytokines are expressed. Further research is required to see how Rora is involved in these regulatory programs.
One theory is that Rora has a complex interaction with circadian rhythm genes. T cells migrate in and out of the lymph node throughout the day [63], and Rora is already known to be involved in circadian rhythm. Our data strengthens this claim as the circadian gene Arntl is DE in both Tregs and other activated CD4+ T cells, in agreement with previous studies [42,64]. It is possible that Rora is regulated through external signalling. We investigated this with a small screen and found, for example, SDF1a as a negative regulator. This could be interpreted in the context of leukocyte migration during circadian rhythms, as blocking SDF1a disables the circadian migration of leukocytes [63]. A caveat with the circadian rhythm influence is that, if the Rora RNA cycles over the day, so will the Rora protein-but with a delay induced by translation speed. This significantly complicates the interpretation of our data, and possibly helps explain any discrepancies in our RNA-seq vs FACS analysis. Further detailed work is required to account for any circadian effects.
Other theories can be brought forward based on DE genes downstream of Rora. If Sell is a direct downstream gene, it would explain most of the activation. However, other downstream genes such as S100a4, Emb and Cxcr6 have also been linked to activation and cell migration. As Rora is a transcription factor which may have thousands of direct targets, the effect may be through more than one downstream gene, making it hard to pinpoint a single molecular mechanism.
The previous Rora KO sg mice display a range of phenotypes, e.g., reduced body weight gain and neuronal disorders. Assuming the gene regulatory network is preserved in other cell types, our work provides hypothetical mechanisms to explore further in future experiments. Since Rora controls PPARa [65], it may be a possible intermediate, with our data suggesting effects on phosphocholines. In this work we did not investigate the effect of potential Rora ligands but these could further contribute to a feedback loop. As for the neuronal disorders, knowing that the SDF1a-CXCR4 interaction is also important for neuronal guidance [66], SDF1a-mediated regulation of Rora may also explain neuronal impacts.
In conclusion, Rora has been found to be important for several immune cell types, in different contexts. It is involved in multiple pathways and appears to link cell migration with cell cycle, but also affects the outcome of inflammation through, but not limited to, Tnfrsf25. In its natural context, Rora is mainly associated with non-lymphoid tissue and activated T cells. In vitro culture suggests that this is due to environmental stimuli, possibly SDF1a. By generating a large systematic transcriptomic characterization of cytokine effects on T helper cells, we have found several cytokines that may explain be responsible for Rora induction in other tissues. With the large number of conditions tested in this study we are one step closer to understanding the multiple functions of this gene.