Data analysis

2.2.11.1 Quantification of transcript levels from RT-qPCR data

Raw RT-qPCR data was attained using the CFX Manager software. Quantification of transcript levels was done using the 2^-DCq method. Cycle quantification values (Cq) of target genes were normalized to housekeeping control (GAPDH) and resulting DCq values used to calculate relative expression ratios 2^-DCq. Fold change expression ratios were determined by dividing the relative expression ratios of samples by that of their respective control.

2.2.11.2 Quality control, differential gene expression, and GO analysis of RNA sequencing data

Raw sequencing reads of Illumnia RNA sequencing were provided as FASTQ files by the BIH Genomics Core Facility. Quality control of raw reads, as well as read mapping and counting was performed in Python according to [345]–[347]. Ensemble Genome Reference Consortium human build 38 (GRCh38.all.fa) was used as reference genome with its respective annotation library (Homo_sampiens_GRCh38.93.gft). The genome index was generated, read files merged and raw reads aligned using STAR algorithm (as described in [345], [348]). STAR alignments were stored in SAM/BAM format and results were visualized in Rstudio (as described in [349]). Alignment success was controlled using the SAMtools (samtools view and samtools flagstat scripts), and RseQC (bam_stat script) Python based software packages according to the developer’s instructions and [345]. Biases in read distributions, in silico calculated RIN (RNA Integrity Number), and similarity between replicate samples were assessed using the anaconda3 (read_distribution and tin scripts), as well as RseQC (geneBody_coverage script) Python based software packages according to the developer’s instructions and [345]. Following successful quality assessment, read counting per gene was performed using the featureCounts script of the Python subread package according to [345]. Hits were called if overlaps (1 bp or more) were found between reads and a single genomic feature (no multiple overlaps allowed). Read normalization, correlation and differential gene expression analysis, as well as visualization of results was performed in Rstudio using the DESeq2 and associated

than 10 read counts across samples were excluded from the differential gene expression analysis. Gene ontology analysis of differentially expressed genes was conducted with the web-based application GOrilla. Significantly regulated genes were provided as target and all expressed genes in U-2 OS cells as background list to determine enrichment of GO terms (biological function).

2.2.11.3 Analysis of raw bioluminescence counts

Raw bioluminescence signals were analyzed using the in-house developed software ChronoStar. Data was trend-eliminated by dividing raw counts by their 24 hour running average. Subsequently, a sinusoidal function (Equation 3) was fitted to the detrended data and circadian parameters were extracted.

𝑥(𝑡) = 𝑒^!"#∗ 𝐴 ∗ cos (𝜔𝑡 ∗ 24 − 𝜔 ∗ 𝜙) (3)

A= amplitude d = damping constant

w = 2p/period [hs]

f = phase [hs]

t = time [hs]

2.2.11.4 Determination of Dphase and Dperiod

Stimulation of U-2 OS Bmal1:Luc or Per2:Luc reporter cells or PER2::LUC murine tissue explants were performed as described. Raw time series were trimmed to start at the timepoint of stimulation and circadian parameters were extracted using ChronoStar. Phases were normalized to periods. Absolute phase differences between samples and their respective controls were determined and transformed into “circular Dphase” values, i.e. to range between ±12 hours.

Determination of Dperiod was done by analyzing raw time series in ChronoStar, extracting period values and calculating the absolute difference between periods of samples and their respective controls.

Phase-pulling experiments

Peak phases of the U-2 OS Per2:Luc reporter cells were determined for the first, second, and third circadian cycle of bioluminescence oscillations post-synchronization.

To do so raw time series were trimmed to start 12 hours (first peak), 36 hours (second peak), or 60 hours (third peak) post-synchronization. Circadian parameters were determined in ChronoStar and phases were normalized to periods.

2.2.11.5 Determination of AUC and fold change AUC

Stimulation of U-2 OS 7xCRE:Luc, 7xmutCRE:Luc, or 7xSRE:Luc was performed as described. Raw time series were trimmed to start at the timepoint of stimulation and area under the curve was determined in GraphPad PRISM (baseline: 0 cps, peak threshold: < 10% of the distance form minimum to maximum cps). Fold AUC was determined by dividing AUC values of samples by their respective controls.

2.2.11.6 Dose-responses to pharmacological TGF-b receptor inhibitor AUC of 7xCRE:Luc reporter cells following stimulation with conditioned or control medium containing TGF-b receptor inhibitor, as well as circadian parameters of Per2:Luc report cells following the addition of TGF-b receptor inhibitor were determined as described. Relative fold change AUC was determined as follows: relative AUCs were calculated by dividing sample AUCs by their respective solvent controls, i.e. 0.0 µM LY2109761 in conditioned or control medium, relative AUCs of conditioned medium were normalized to respective relative AUCs of control medium.

Relative amplitude and damping parameters were determined by dividing amplitude and damping parameters of by their respective solvent controls, i.e. 0.0 µM LY2109761. Dperiod values were determined as described above and relative to solvent control, i.e. 0.0 µM LY2109761. EC50 values were determined in GraphPad PRISM by non-linear regression fitting of an asymmetric sigmoidal curve to relative fold change AUC, relative amplitude/damping, or Dperiod values.

2.2.11.7 Calculation of % recovery and selection of active chromatography fractions

Percent recovery of conditioned medium activity following stimulation of U-2 OS 7xCRE:Luc reporter cells was determined as follows:

(i) “assay AUC” was determined in GraphPad PRSIM as described above (for 1:5 dilution of input or fractions in reporter medium)

(ii) “5-fold AUC” was calculated by extrapolating the assay AUC (based on a previously determined standard curve of a conditioned medium dilution series) (iii) Total activity of input and chromatography fractions was calculated by

multiplying the 5-fold AUC with the absolute input or fraction volume

(iv) % recovery was calculated by dividing total activity of the fractions by total activity of the input (multiplied by 100)

Active fractions were defined as those fractions with % recovery > mean ± SD % recovery of all fractions, as well as with protein content < mean protein content of all fractions (excluding negative absorbance values at 280 nm).

2.2.11.8 Identification of secreted protein hits

Identification of protein hits from peptide sequences was done by our collaboration partners at the Protein Purification and Analysis Unit of the Max Planck Institute for Infectious Biology (Berlin, GER). In brief, observed mass spectrometry spectra were compared to a contaminant, as well as the SwissProt (release 2018_11, taxonomy:

homo sapiens) primary sequence databases. Following parameters were used during the Mascot search:

Table 2-6: Mascot search parameters

Type of search MS/MS Ion search

Enzyme Trypsin/P

Variable modification Acetyl (Protein N-term),

Carbamidomethyl (C), Gln à pyro-Glu (N-term Q), Oxidation (M)

Mass values Monoisotopic

Protein mass Unrestricted

Peptide mass tolerance ± 5 ppm

Fragment mass tolerance ± 0.03 D

Max missed cleavage 2

Instrument type ESI-FTICR

FDR 1%

Resulting Mascot protein hits of active and inactive fractions were filtered for human secreted proteins predicted by MDSEC using Rstudio (human protein atlas:

https://www.proteinatlas.org/search/protein_class:Secreted+proteins+predicted+by+

MDSEC, accessed Feb 2019). Subsequently, secreted proteins hits identified in inactive fractions were removed from hits of the active fractions.

2.2.11.9 Statistical analysis

All statistical analyses were carried out in GraphPad PRISM.