Rheumatoid arthritis reprograms circadian output pathways

Study design

We recruited adults (age 16–80 years old) with a clinical diagnosis of RA and age- and gender-matched controls from The Kellgren Centre for Rheumatology, Manchester University NHS Foundation Trust. All patients were seropositive for rheumatoid factor (RF) and/or anti-citrillunated peptide antibodies (ACPA), fulfilled 2010 ACR/EULAR Criteria for RA [12] and had active disease with a DAS28 > 3.2. Controls who did not suffer from any chronic inflammatory disease were recruited. Exclusion criteria in both groups were as follows: regular or recent (1 month) medication affecting sleep, glucocorticoid use, sleep disorder, BMI > 30, shift working, and recent (< 2 weeks) trans-time zone travel. The medication for each patient is shown in Additional file 1: Supplemental file S1. Two patients were on no disease modifying or biologic therapy at the time of study. One was a recently diagnosed early RA, and the other had stopped methotrexate and hydroxychloroquine and was being worked up to start a new therapy. See Additional file 2: Supplementary methods file for power calculations.

Subjects recorded a sleep diary and wore an actigraph, (PAL Technologies Ltd., UK) [13] for 1 week at home, before being admitted to the NIHR/Wellcome Trust Manchester Clinical Research Facility for a 24-h stay. Subjects were not subjected to circadian conditions, but rather were in a room with access to natural light, and no television from midnight until 7 am. Prescribed inactivity was between 10 pm and 7 am. Meals were provided as per standard MCRF protocol, at 7 pm, 8 am, and 1 pm. All subjects had free access to water. The study protocol was approved by local ethics committee (IRAS 102779, REC ref. 12/NW/0815).

Statistical analysis

Univariate analysis was carried out using a normal data general linear model, repeated measures ANOVA or a Wilcoxon rank sum test. Normality was determined using a Q-Qplot (using qqnorm) where appropriate. Perseus version 1.6 was used to determine mean lipid levels (with a permutation-based false discovery rate). Rhythmic lipids were determined using a Gaussian process regression with a periodic kernel [14], implemented in python (using Enthought solutions). The data are normalised to have zero mean and unit standard deviation per participant and per lipid. Then, the hierarchical statistical model (based on [15]) accounts for structure in the data in the following way. The concentration of each lipid is modelled using a hidden Gaussian process, which is common to all participants in the study. This Gaussian process is given a rhythmic structure using the covariance function described by [14]. Each participant’s lipid level is modelled as a further deviation from this Gaussian process by another Gaussian process, which is lipid- and participant-specific, using a standard Matern 3/2 covariance function. Further variation is modelled as Gaussian noise. A small number of parameters controlling the relative variance of the Gaussian processes in the hierarchy are estimated by maximum likelihood.

The data likelihood ratios were generated for healthy and RA patients; 5% level (chi-squared distribution) was considered significant. Missing values (due to one broken vial) were imputed with the MICE package in R [16]. There were missing values for healthy subject 8 (due to inadequate morning blood collection) and saliva from RA subject 3. RA10 left the study early and was replaced by RA11. Gene ontology analysis was carried out using PANTHER v11 [17] .

RNAseq

RNA integrity was determined using a Tape Station (Thermo); no loss of RNA integrity was observed during storage. Four samples were run per lane on a HiSeq (Illumina). Sequencing data were processed with Bowtie-TopHat (with Hg38). Differential expression was determined using EdgeR/LIMMA in the R environment. EdgeR implements generalised linear models. Using this approach, we were able to perform within- and between-group differences with a strict control over the error rate.

Proteomics and phosphoproteomics

See Additional file 2: Supplementary methods file.

Mediator lipidomics

Mediator lipidomics was performed using targeted quantitative LC-MS/MS assays [18–20]. See Additional file 2: Supplementary methods file.

CIA mouse model

All animal procedures were performed in accordance with the United Kingdom Animals (Scientific Procedure) Act 1986. Animals were housed in 12:12 lighting conditions (LD; 12-h light to 12-h dark) with ad libitum access to normal rodent chow. All times are Zeitgeber (ZT), where ZT0 is lights on.

Male DBA/1 mice (Envigo) age 8–12 weeks were immunised with 1 mg/ml bovine type II collagen (MD Bioscience, Zurich, Switzerland) emulsified with an equal volume of complete Freund’s adjuvant (CFA; MD). The emulsion (50 μL) was injected intradermally, under anaesthesia in two sites just above the base of the tail. As a control, animals were injected with an emulsion of saline and CFA in equal volumes. Twenty-one days later, a booster (200 μL 1 mg/ml collagen in saline) was administered intraperitoneally to CIA animals. Immunisation and boost were routinely administered at zeitgeber time (ZT) 6. Disease incidence and severity and hind paw thickness were monitored from 18 days post immunisation (DPI). Disease severity was scored in each limb on a 4-point scale: (1) one inflamed digit, (2) two or more inflamed digits, (3) swelling of the foot pad and minor ankylosis and (4) severe swelling of the foot pad and joint and severe ankylosis, and the sum was calculated. RNA was extracted using phenol-chloroform (limbs) or a SV RNA extraction kit (liver) (Promega). First strand synthesis was carried out using a GoScript cDNA conversion kit; qPCR was carried out using GoTaq qPCR master mix (see Additional file 2: Supplemental methods for primer sequences).

Assessment of parameters of diurnal activity in RA patients

We recruited patients with seropositive RA and age-/gender-matched healthy controls and examined their sleep profile for 1 week, before inpatient measurements of saliva, serum, and peripheral blood mononuclear leukocytes (PBML) (Fig. 1a, b). Our RA (n = 10) cohort was selected on the basis of active disease, DAS28 score > 3.2, but without recent steroid use (Additional file 1: Supplemental file S1 and Additional file 2: Supplementary methods for power analysis and study design). We did not detect any group differences in sleep metrics, or hypothalamic-pituitary-adrenal axis function or rhythm (Fig. 1c, d, Additional file 3: Figure S1A-E).

Identification of diurnal gene expression changes specific to RA patients

To obtain sufficient cells for analysis, as we were limited in sampling frequency, we selected dawn and dusk (06:00 and 18:00); these time points had previously been found to show the greatest difference in multiple animal models of circadian control of inflammation [21]. Additionally, the peak of clinical RA disease activity is at 06:00 [22] at the end of the rest phase.

Comparing gene expression between RA patients and healthy controls, we found 1547 genes (Fig. 2a) to be differentially expressed at 06:00, but in marked contrast, only 287 genes (Fig. 2b) at 18:00, emphasising the magnitude of the circadian effect. The overlap between these two datasets is shown in Additional file 4: Figure S2E. There was no effect of age, BMI or the salivary concentration of cortisone on differential gene expression at these time points (using linear regression of differentially expressed genes against age, BMI or 06:00/18:00 salivary cortisone concentrations). Genes identified as significant in these comparisons showed no correlation with gene expression level. Within the RA group, we identified 104 genes differing between AM and PM, but only 25 in healthy controls, suggesting a gain in rhythmic gene with disease (Fig. 2c, d), which was not due to differences in circulating immune cell repertoire (Additional file 4: Figure S2A-D).

PER3 is a robust circadian gene in humans [23], and this showed similar time-of-day changes in both RA and control groups, peaking in the morning (Fig. 2e), indicating that the immune cell core clock operation was likely unaffected by disease and that the study conditions did not significantly perturb underlying circadian rhythmicity. Some genes gained a time-of-day variation in expression with RA, such as OCRL (Fig. 2f), and also included IL6ST, SOCS3, TLR2 and HCAR3 (shown in Additional file 4: Figure S2F-I), again supporting gain of rhythmic function in disease. We also found a number of genes whose expression was determined by disease status, but not responsive to time of day, LTBR4 and TNFRSF13B (Additional file 4: Figure S2J, K).

Analysis of the differentially expressed genes at 06:00 revealed enrichment of terms implicating phagocytosis and antibody-mediated complement activation and B cell activation pathways (Fig. 2g) and immune cell activation, particularly TLR pathways and innate immune responses (Fig. 2h and example of TL4 cascade genes shown in Additional file 4: Figure S2L and M). Figure 2i shows the expanded innate immune network from Fig. 2h, which includes MAPK14 (P38a). Likely, transcription factors involved were SP1 and STAT3 (Fig. 2j). Although fewer genes were differentially expressed at 18:00, similar Reactome pathways were identified.

Integrated analysis of the HPA axis and PBML gene expression

The connections between the core cellular circadian pacemaker, the HPA axis (CORT) and the HPA axis biomarker (GILZ) were investigated further using non-parametric analytical approaches (Additional file 5: Figure S5B, C). In all subjects, robust correlations between the core clock genes were seen, as a positive association between GILZ and salivary cortisone, as expected. The relationships between core clock gene expression were all stronger in the RA group compared to the controls (depicted by diameter of the circles in the plot), again strongly supporting the emergence of a more robust, more tightly cross-coupled circadian oscillator in the presence of chronic, active inflammation (Additional file 5: Figure S5B and C).

Identification of diurnal proteome changes in RA

Because we observed such striking differences in gene expression at 06:00, we measured the phosphoproteome in pursuit of upstream causative pathways. The total immune cell proteome was analysed at 06:00 and 18:00, and we found changes between healthy and RA patients with the greatest difference again being at 06:00 (Fig. 3a). Of the 2710 proteins detected, we found only 27 that were increased in RA. At 18:00, only 3 proteins differed in expression between RA and controls.

Protein determination thus identified relatively few changes dependent on disease state. In marked contrast, phosphoproteomic profiling revealed extensive differences, with the striking observation showing clear separation between dawn and dusk specific to RA patients (Fig. 3b). This was a similar gain of temporal effect as was seen with the earlier transcriptome analysis (Fig. 2c, d).

A total of 147 phosphosites varied between dawn and dusk in the RA group (Fig. 3c) from which we identified a consensus arginine and proline directed phosphoserine site as strongly enriched within the temporally controlled RA phosphoproteome (Fig. 3d). Interestingly, some sites such as phospho-T198 of cAMP-dependent protein kinase alpha (PRKACA, PKAa) were higher in the morning, while a closely related site phospho-T196 was higher in the evening. Phospho-sites that were unique to RA included HSP27-S82, MAPK14(p38)-Y182 and cAMP-dependent PKA-T198 (Fig. 3e). The time-of-day-regulated kinases formed a functional network (Fig. 3f), with phosphosites found to differ by time of day and to have known regulatory function marked on the figure, with inferred connecting nodes as white. This analysis revealed multiple members of the MAP kinase family, and PKA, as being regulated in RA cells by time of day.

Immune cells responses to activation by time of day in RA

As we observed major changes in basal gene expression and phospho-proteome in immune cells from RA patients at dawn, we sought responses to cell activation ex vivo. Circulating immune cells were activated (with LPS and anti-CD3/28; from here LPS denotes the mixture of stimuli). This revealed more than 8000 transcripts regulated by LPS (Fig. 4a). We further analysed differentially regulated genes using Edge set enrichment analysis (using healthy AM as the control character and genes activated in RA at both times), showing several well-characterised disease targets as being implicated in the differential pattern of gene expression (Fig. 4b). This, again, included members of the MAP kinase family (MAPK1, MAPK3, MAP3K7). In the PM samples, STAT3 activity was inferred (Fig. 4b: see Additional file 6: Figure S3A and B for full networks). Pathway analysis of genes regulated in only in RA at 06:00 and 18:00 implicated AKT/PTEN (06:00) and RNA polymerase components (18:00) (Additional file 6: Figure S3C,D).

More genes responded to activation in RA (Fig. 4c, d). Analysis of RA-specific responder genes at 06:00 (Fig. 4e) or 18:00 (Fig. 4f) highlighted metabolic pathways and translational control respectively. Analysis of likely transcriptional regulators of 06:00 and 18:00 LPS-regulated genes in RA identified a number of candidate transcription factors, including RXRa (06:00) and CREB1 (18:00) (Additional file 7: Figure S4A-C), both targets for PKA, which we find is time of day regulated in RA (Fig. 3f).

Phosphoproteomic analysis identified 1494 phosphorylated motifs in total, of which 664 were regulated by cellular activation (Fig. 4g, Additional file 8: Supplemental data file). Phosphorylation sites on HSP27 (S15, 65 and 82) and STAT3 (S727 and Y705) were highly responsive to LPS, but in these ex vivo studies, we did not detect any time-of-day difference in response. Motif analysis highlighted a role for MAPK (PxSPx) and AKT/PKA (RxxS), a similar pattern to that previously found (in Fig. 3d); however, the PKAa T198 site that we previously found to be regulated by time of day (Fig. 4g) was unresponsive to LPS.

To quantify more precisely changes in phosphorylation of these and other members of the MAP kinase-signalling network, we moved to a quantitative antibody-based bead array method (Luminex) (Fig. 4h). Phospho-STAT3(S727) was increased by LPS and was found higher in RA (p = 0.005). There was also a disease and time interaction (p = 0.03), with STAT3 phosphorylation increased at 18:00 in the healthy controls. ATF2, a transcription factor and target for ERK1/2 (and other MAPKs), was phosphorylated in response to LPS (T71) (p < 0.001) and showed a time-of-day increase in the RA group (disease state × time, p = 0.03), further adding to the previous finding that MAPK signalling is altered in the RA group by time of day.

Serum lipidomic analysis in RA patients reveals newly rhythmic ceramide species

In additional to regulating mRNA expression, circadian rhythms also regulate the concentration of other biomolecules, including circulating serum components. As metabolite circadian rhythms can act on clock gene circuits to affect their function, and in turn clock output pathways can regulate the metabolome, we also analysed the serum metabolome in the same two groups of subjects [24–26]. In previous studies in healthy human volunteers, the major circulating metabolome components lying under circadian control were lipids, many of which have well-defined signalling functions in inflammation and immunity [27]. Therefore, we focussed our attention on serum lipids capable of transmitting signals to target cells.

We used a targeted LC-MS/MS approach and measured 116 ceramides, eicosanoids and endocannabinoids (Fig. 5a). The ceramides were further classified according to structure, comprising either a sphingosine (S) or dihydrosphingosine (DS) sphingoid base and either a non-hydroxy (N) or alpha-hydroxy (A) fatty acid. Overall, we could not distinguish study groups by PCA analysis of serum lipids (not shown), and there was no effect of gender, age, or BMI (using regression analysis, p > 0.05).

A Gaussian process model was used to identify rhythmic lipids [14, 15]. We found more rhythmic lipids in the RA group, particularly ceramides (Fig. 5b). The peak times for each rhythmic species were also estimated (Fig. 5c, d), with the acrophase of the newly rhythmic ceramides occurring at 23:00. The majority of rhythmic lipids in both healthy and RA subjects peaked at 06:00 or 18:00 (Fig. 5c, d). This suggested that the newly rhythmic ceramides were products of a newly rhythmic enzymatic pathway.

A small number of serum lipids were found to be rhythmic in both RA patients and controls, including 9HOTrE, 12(13)EpOME, 13HOTrE and the ceramides CER[N(18)DS(24)], CER[N(18)DS(26)], CER[N(22)DS(18)], CER[N(24)DS(18)] and CER[N(29)S(18)]. The newly rhythmic ceramide species in RA were particularly of the CER[NDS] class (examples shown in Additional file 5: Figure S5A), with fewer rhythmic alpha-hydroxy ceramides (CER[AS(18)DS(18)] shown as an example in Fig. 5e). Some lipids showed a stronger and significant oscillation only in healthy subjects (9-HODE shown as an example in Fig. 5f), interesting as both 9- and 13-HODE species had previously been identified as rhythmic in a healthy human circadian study [27, 28], usefully providing an external quality control for our clinical protocol.

Serum ceramides similarly acquire a new circadian rhythmicity in mouse experimental arthritis (CIA)

To identify if a similar acquisition of ceramide rhythmicity was a consistent feature of chronic inflammatory arthritis, we assessed a CIA mouse model (Additional file 9: Figure S6). This model results in chronic, destructive arthritis (Additional file 9: Figure S6A,B,C), and the disease expression, measured by paw swelling or cytokine expression, shows a strong diurnal variation in severity, as in human RA [21, 29]. Mice were kept in 12:12 L:D conditions, and times are expressed at Zeitgeber time (ZT) where ZT0 is lights on or 06:00. We identified the gain of circadian oscillation in similar ceramide species in the CIA mouse serum as in RA (Fig. 6a, b). Examples of newly rhythmic ceramide species in the CIA mice are shown in Fig. 6 c.

Ceramide synthesis is complex [30] (Fig. 6d, de novo pathway shown). A number of ceramide synthase genes were found to be expressed in the inflamed limb (Additional file 10: Figure S7A-E), but while some (CERS2) were suppressed by inflammation, and others (CERS4, 5, 6) were stimulated, there was no variation in expression over time. Gene expression of other enzymes regulating ceramide biosynthesis did not vary by time in CIA or control mouse limbs (Additional file 10: Figure S7F).

In contrast, in the liver, we identified cycling ceramide synthase gene expression only in CIA mice. CERS2 and CERS4 reached their nadir at ZT6 (Fig. 6e), while CERS6 did not oscillate but was significantly higher in the arthritic mouse liver tissue (Fig. 6f). Importantly, there was no impact of arthritis on liver core clock gene oscillation (Fig. 6g).

In order to explore mechanisms by which CIA might drive cycling hepatic ceramide synthesis, we investigated whether inflammatory mediators from the inflamed joints were acting on the liver to drive rhythmic gene expression changes over time. In the CIA model, serum IL-1β and IL-6 show circadian variation [21] and these were chosen as candidate mediators. IL-6 did not affect mouse liver cell CERS expression (Fig. 6h), but IL-1β potently suppressed CERS4 and stimulated CERS6 (Fig. 6h). Therefore, circadian production of inflammatory cytokines may convey a timing signature to hepatic ceramide synthase expression. The long-chain ceramides we see varying by time of day are principally regulated by CERS2 and CERS4, both of which cycle in the liver, rather than CERS6 which did not significantly vary across time. It is also notable that serum IL-1β peaks at ZT6 [21], the nadir of CERS4 gene expression, and serum ceramide concentrations (Fig. 6c, e).