Graphical abstract [37]graphic file with name ga1.jpg [38]Open in a new tab Regulatory classification of dexamethasone-responsive genes in macrophages by integration of ribo-seq and RNA-seq data Keywords: Glucocorticoid receptor, Macrophages, Inflammation, Ribosome profiling, RNA-seq, Translational regulation Abstract Glucocorticoids such as dexamethasone (Dex) are widely used to treat both acute and chronic inflammatory conditions. They regulate immune responses by dampening cell-mediated immunity in a glucocorticoid receptor (GR)-dependent manner, by suppressing the expression of pro-inflammatory cytokines and chemokines and by stimulating the expression of anti-inflammatory mediators. Despite its evident clinical benefit, the mechanistic underpinnings of the gene regulatory networks transcriptionally controlled by GR in a context-specific manner remain mysterious. Next generation sequencing methods such mRNA sequencing (RNA-seq) and Ribosome profiling (ribo-seq) provide tools to investigate the transcriptional and post-transcriptional mechanisms that govern gene expression. Here, we integrate matched RNA-seq data with ribo-seq data from human acute monocytic leukemia (THP-1) cells treated with the TLR4 ligand lipopolysaccharide (LPS) and with Dex, to investigate the global transcriptional and translational regulation (translational efficiency, ΔTE) of Dex-responsive genes. We find that the expression of most of the Dex-responsive genes are regulated at both the transcriptional and the post-transcriptional level, with the transcriptional changes intensified on the translational level. Overrepresentation pathway analysis combined with STRING protein network analysis and manual functional exploration, identified these genes to encode immune effectors and immunomodulators that contribute to macrophage-mediated immunity and to the maintenance of macrophage-mediated immune homeostasis. Further research into the translational regulatory network underlying the GR anti-inflammatory response could pave the way for the development of novel immunomodulatory therapeutic regimens with fewer undesirable side effects. 1. Introduction Gene regulatory networks (GRNs) are continually adjusting mRNA and protein levels according to cellular needs, affecting all steps of the expression process, from transcription to mRNA maturation, transport, translation, and protein degradation [39][1]. These adjustments are further customized in the presence of external stimuli that transform the transcriptional and translational landscapes within the cells. Over the last two decades, assessing gene transcription levels using next generation sequencing methods such as RNA-seq, have become central to the understanding of the molecular mechanisms underpinning the biological processes and phenotypes in various biological contexts. RNA-seq analysis can, however, not always capture the full picture of gene expression, as the final fate of mRNA molecules depend on mRNA stability, mRNA degradation, and translation efficiency (TE). ribo-seq offers a quantitative approach to study translational regulation of global transcriptomes by quantifying the capture of ribosome-protected RNA fragments (RPFs) [40][2]. Normalizing ribosome footprint density by mRNA abundance determines a gene's translation efficiency. In ribo-seq, the number of RPFs between different conditions for a given gene is used as a proxy for a change in the translation of the encoded protein. A major caveat with ribo-seq is, however, the mRNA abundance of the transcript, which directly affects the probability of ribosome occupancy, and thus may not necessarily reflect the transcriptional efficiency of each gene. Synthetic glucocorticoids (GCs) such as dexamethasone (Dex), are commonly used in the treatment of a variety of immune and inflammatory diseases due to their potent anti-inflammatory effect. Their immunomodulatory properties stem from a complex mechanism of action that dampens cell-mediated immunity by inhibiting the expression of pro-inflammatory cytokines and chemokines, while activating the expression of anti-inflammatory immune mediators [41][3], [42][4], [43][5], [44][6]. The effects of GCs are mediated through the ubiquitously expressed intracellular, ligand-dependent glucocorticoid receptor (GR). GR classically exerts its transcriptional regulation, following GC binding, by interacting with evolutionarily conserved 15 bp palindromic consensus DNA sequences (AGAACANNNTGTTCT) designated as glucocorticoid response elements (GREs). GREs are primarily present in the enhancer and promoter regions of GR target genes and drive, predominantly through protein–protein interactions with co-activators and histone acetyl transferases (HATs), the transcription of anti-inflammatory factors [45][7], [46][8], [47][9], [48][10]. Transcriptional inhibition by GR is, on the other hand, exerted via recruitment of co-repressors such as Glucocorticoid receptor interacting protein 1 (GRIP1) and histone deacetylases (HDACs), disrupting pro-inflammatory NF-κB/interferon response factor3 (IRF3) complexes, and interaction with classical palindromes or with cryptic GREs within NF-κB and AP-1 motifs [49][11], [50][12], [51][13], [52][14] ([53]Fig. S1). Though regulation by GCs/GR has been extensively studied on the transcriptional level, the impact of GCs on post-transcriptional gene regulatory mechanisms is largely unknown, with the exception of a small subset of pro-inflammatory genes. These include effector molecules such as Tumor necrosis factor alpha (TNFα), interleukin 6 (IL6), and Chemokine (C—C motif) ligands (CCL)2 and CCL7, suggesting a role for GCs/GR in post-transcriptional gene expression modulation [54][15], [55][16], [56][17], [57][18]. In the present study, we aimed at systematically investigating the global impact of Dex treatment on the transcriptional and post-transcriptional regulation of gene expression in an inflammatory setting. For this, a well-established in vitro model of inflammation using macrophage-like cells derived from phorbol 12-myristate 13-acetate (PMA)-differentiated human THP-1 cells was utilized [58][19]. THP-1 is a human monocytic cell line derived from a patient with acute monocytic leukemia, which has been widely used to investigate monocyte/macrophage activities [59][20]. Differentiated THP-1 cells are extremely sensitive to LPS and respond by releasing a plethora of inflammatory cytokines [60][21], [61][22], [62][23]. To investigate the effect of Dex on inflammation, Dex treatment was performed concurrently with LPS stimulation, with the latter inducing a pro-inflammatory M1 macrophage phenotype. Further, using previously published methodology and integrating matched mRNA-seq data with ribo-seq data, the change in translation efficiency (ΔTE) was calculated [63][24]. This approach allows the subcategorization of genes into buffered, forwarded, intensified, or exclusively (exclusive) translationally regulated genes. We show that the expression of the majority of the Dex-responsive (activated and suppressed) and LPS-stimulated genes can be classified as “intensified”. These genes undergo statistically significant transcriptional and translational changes and display a change in translational efficiency. Overrepresentation analysis (Reactome pathways, KEGG pathways, and gene ontology (GO)), combined with STRING protein network analysis and gene-by-gene manual functional exploration, identified these genes to exhibit immune effector and immunomodulatory functions that contribute to maintaining macrophage-mediated immune homeostasis. Overall, our findings show that, in addition to transcriptional control, most Dex-responsive genes are subject to translational regulation in LPS-stimulated macrophages, and that the protein products of these genes are involved in macrophage-mediated immunity. 2. Materials and methods 2.1. Cell propagation, culture, and treatment THP-1 cells were obtained from ATCC and cultured in RPMI 1640 medium (Thermo Fisher Scientific; A1049101) supplemented with 10 % heat inactivated FBS (Sigma Aldrich; F9665), 1 % Pen/Strep (Sigma Aldrich; P4333). The cells were grown at 37 °C and 5 % CO[2] and routinely tested negative for mycoplasma contamination. For THP-1 differentiation, cells were grown in culture media containing the dialyzed FBS (Sigma Aldrich; F0392) and incubated with 10 ng/ml (final concentration) phorbol 12-myristate 13-acetate (PMA, Sigma Aldrich; P1585) for 24 h, followed by 24 h incubation with fresh growth media (without PMA). The differentiated THP-1 cells were treated for 3 h with vehicle (0.1 % EtOH), 100 ng/ml lipopolysaccharide (LPS, Sigma Aldrich; L6529), or with a combination of 100 ng/ml LPS and 1 μM Dex (Sigma Aldrich; D4902). All treatments were performed on five independent biological replicates (n = 5). Following the 3 h incubation, the cells were harvested for ribosome profiling and mRNA sequencing. 2.2. Ribosome profiling Ribosome profiling of THP-1 derived macrophage-like cells was performed according to previously published work by our group and others [64][25], [65][26], [66][27]. In short, 10 million cells were lysed for 10 min on ice in 1 mL lysis buffer consisting of 20 mM Tris-Cl (pH 7.4), 150 mM NaCl, 5 mM MgCl[2], 1 % Triton X-100, 0.1 % NP-40, 1 mM dithiothreitol, 10U/ml DNase I, cycloheximide (0.1 mg/ml) and nuclease-free H[2]O. The lysate was homogenized by immediate repeated pipetting and multiple passes through a syringe fitted with a 21G needle, allowing for dissociation of cell clumps and facilitating quick and equal lysis of the cells. Samples were next centrifuged at 20,000g for 10 min at 4 °C to pellet cell debris. Per sample, 200 or 400 μL of lysate were digested with RNase 1 (N6904K; Biozym), purified with Microspin Sephacryl S-400 HR columns (Sigma-Aldrich; GE27-5140–01) and 1μg footprints were used for removal of the ribosomal RNA according to the ribo-Zero Gold rRNA Removal Kit (h/m/r) (Illumina; MRZG12324). The footprints were purified by excision from a Novex 15 % TBE Urea PAGE gel (Fisher Scientific; EC68852BOX) followed by a treatment of the 3́ends with T4 PNK (Biozym; P0503K) to allow ligation to a pre-adenylated linker. The RNA was then reverse transcribed with Reverse Transcriptase (Biozym; ERT12925K) and the cDNA purified via a Novex 10 % TBE Urea PAGE gel (Fisher Scientific; EC68752BOX). The fragments were circularized using Circligase I (Biozym; CL4115K) followed by a PCR amplification and size selection using a Novex 8 % TBE PAGE gel (Fisher Scientific; EC62152BOX). For all samples, ribosome profiling library size distributions were checked on the Bioanalyzer 2100 using a high sensitivity DNA assay (Agilent; 5067–4626), multiplexed and sequenced on a NovaSeq 6000 Illumina producing single end 1x51nt reads. The samples were processed in one batch to avoid a sample processing bias. THP-1 macrophage ribo-seq libraries were sequenced to an average depth of 80 M (min. 64 M, max. 91 M) raw reads. 2.3. Stranded mRNA sequencing Total RNA was isolated using TRIzol Reagent (Invitrogen; 15596018) using 3 million of the exact same human THP-1 derived macrophage-like cells processed for ribosome profiling. RIN scores were measured on a BioAnalyzer 2100 using the RNA 6000 Nano assay (Agilent; 5067–1511). Poly(A)-purified mRNA-seq libraries were generated from high quality RNA (average RNA Integrity Number (RIN)) of 7.7. RNA-seq library preparation was performed according to the TruSeq Stranded mRNA Reference Guide, using 500 ng of total RNA as input. Libraries were multiplexed and sequenced on an Illumina NovaSeq 6000 producing paired 2x101nt reads. THP-1 macrophage mRNA-seq libraries were sequenced to an average depth of 160 M (min. 81 M, max. 207 M) raw reads. 2.4. RNA-seq and ribo-seq data analysis Ribo-seq and RNA-seq were analyzed as previously described [67][28]. Before genome mapping, Ribo-seq reads were clipped to remove residual adapters using the FASTX toolkit ([68]https://hannonlab.cshl.edu/fastx_toolkit/). Reads mapping to the ribosomal RNA and tRNA sequences were removed with Bowtie2 v2.4.5 [69][29]. The paired mRNA-seq reads (2x101nt) were trimmed to 29-mers (the average length of Ribo-seq reads) to reduce any read length or filtering bias towards mapping or quantification. This was done in order to establish the comparability of the data obtained from Ribo-seq and RNA-seq. Ribo-seq and trimmed RNA-seq reads were then mapped to the human reference genome (GRCh38, Ensembl v87) using STAR v2.7.10a [70][30]. During the alignment step, the following parameters were used: “--limitIObufferSize300000000-- limitOutSJcollapsed10000000--outSAMattributesAll--outFilterMultimapNmax 20--outFilterMismatchNmax 2 --alignSJoverhangMin 1000--twopassMode Basic”. 2.5. Prediction and differential translation of actively translated ORFs Actively translated ORFs were detected with ORFquant 1.00 [71][31] using a pooled set including all mapped ribo-seq samples. We only considered ORFs with a minimum number of 15P-sites in: i) cells treated with lipopolysaccharide (LPS), and ii) cells treated with a combination of lipopolysaccharide (LPS) and Dex. Predicted ORFs were divided into six considered sORF biotypes: CDS (annotated coding sequences, partially or totally in-frame), lncRNA ORFs (lncRNA-ORFs, encoded by presumed long non-coding RNAs), ncRNA ORFs (ncRNA-ORF, encoded by processed transcripts of protein-coding genes), upstream ORFs (uORFs, encoded by 5′ UTR sequences), internal ORFs (intORFs, fully overlapping an annotated CDS in an alternative frame), and downstream ORFs (dORFs, encoded by 3′ UTR sequences)). Next, P-site counts were extracted with RiboseQC ([72]https://doi.org/10.1101/601468). In-frame P-site counts were quantified for each called ORF and used as input to identify differentially translated ORFs with DESeq2 v1.26.0 [73][32]. ORFs with an absolute fold change equal or higher than 1.5 and an adjusted p-value lower than 0.01 were defined as differentially translated. 2.6. Quality control (QC) analysis of ribo-seq data An algorithm called ribo-TIS Hunter (ribo-TISH) was used to perform quality control analysis on the mapped RPFs [74][33]. RPFs distribution around annotated coding genes and grouping by read length were checked using the “quality” function of this toolkit. RPF 5′ end distribution near annotated gene start codons is assessed for each read group before estimating P-site offsets. The RPF count between the 15 bp upstream of the first base of the start codon and the 12 bp upstream of the first base of the stop codon was taken into account when calculating the RPF count distribution in three reading frames as well as the CDS metagene profile. The ratio between the highest RPF count across all three reading frames and the total of all RPF counts was used to calculate the percentage of RPF counts in the dominant frame. Furthermore, RPF counts between −40 and +20 bp of the start/stop codon's first base were summed across all annotated genes to generate an RPF count profile near the start/stop codon. Further QC analysis was performed using the RiboWaltz R package to identify P-sites and calculate the distance of the P-site from the start and stop codons of the coding sequence. Furthermore, the tool returns the data structure used to calculate the distribution of P-site across transcript regions (5′ UTR, CDS, and 3′ UTR) and trinucleotide (3-nt) periodicity covered by the P-site [75][34]. 2.7. Identification of translationally regulated genes The quantification of gene expression was accomplished by using featureCounts v2.0.1 to count the number of reads that mapped to a known genomic and transcriptomic reference [76][35]. DESeq2 v1.26.0 statistical model (design=∼Condition + SeqType + Condition:SeqType) was used to normalize gene counts with respect to size factor or other normalization parameters for both ribo-seq and RNA-seq data [77][32]. The translation efficiency of genes was calculated by combining RPF and mRNA counts, and then the deltaTE (ΔTE) method in R was used to identify genes that were differentially and translationally regulated, as described by Chothani et. al. [78][20]. A false discovery rate (FDR) threshold value of 0.05 was utilized as a cutoff to quantify change in translation efficiency for each gene and categorize genes into various regulatory classes utilizing the ΔTE approach. 2.8. Pathway enrichment and regulatory network analysis The gene list from each regulatory class was analysed using the Enrichr R package to test the probability that annotated gene sets were statistically enriched for Gene Ontology (GO) terms (version 2015) based on their participation in biological processes [79][36]. Associations were considered statistically significant if their adjusted p-values for multiple testing were less than 0.05. We further performed integrative pathway enrichment analysis on RNA-seq, ribo-seq, and ΔTE datasets using ActivePathways R package with the default parameters [80][37]. We provided the p-values of differential changes in mRNA, ribosome footprints (RPF), and translation efficiency (TE) of genes listed in each regulatory class as the first input in this three-step method. As a second input, we used gene sets corresponding to the Reactome Databases Molecular Pathways [81][38] and Gene Ontology Biological Process [82][39] downloaded from the g:Profiler web server [83][40]. Specifically, in the first step of the ActivePathways method, adjusted p-values from different datasets for each gene were merged using a statistical data fusion approach, resulting in an integrated gene list [84][41]. The integrated gene list was then ranked in decreasing order of significance and filtered using the default threshold (unadjusted p < 0.1). Step two involved statistical enrichments of pathways using the ranked hypergeometric test, which determines the pathways significantly enriched in the integrated gene list. Step three generated separate gene lists from individual input datasets and performed similar pathway enrichments (as in step 2) with the ranked hypergeometric test, to find supporting evidence for each pathway from the integrative analysis. To identify key RNA-binding proteins (RBPs) that were differentially expressed in the most prominent “intensified” regulatory class in our data, we used experimental evidence-based databases of RBP-RNA bindings derived from ENCODE [85][41] and POSTAR3 [86][42]. To build the RBP - targets regulatory network, the R packages GENIE3 [87][43], igraph [88][44], and RCy3 [89][45] were used. As input data files, log transformed adjusted p-values for genes that were differentially transcribed (DTG) and translationally efficient (DTEG) from a specific biological process (derived from Enrichr analysis) were used. Cytoscape was employed to visualize the RBP-target network [90][45]. 2.9. Expression analysis Differential expression analysis of genes within each regulatory class was performed using Volcano plots in RStudio (RStudio Team (2022)). RStudio: Integrated Development Environment for R. RStudio, PBC, Boston, MA URL [91]https://www.rstudio.com/.) using the ggplot2 package (Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, [92]https://ggplot2.tidyverse.org). Lists of significant and non-significant DEGs are compiled in [93]Supplementary file 1. Comparative gene analysis using Venn diagrams of significantly up (log2FC ≥ 0.58, p adj ≤ 0.05) and downregulated (log2FC ≤ -0.58, p adj ≤ 0.05) differentially expressed genes (DEGs) was done using the online software tool VENN DIAGRAMS generated by the Van de Peer Lab ([94]https://bioinformatics.psb.ugent.be/webtools/Venn/). Comparisons were done between genes categorized as “intensified” within each platform (RNA-seq and ribo-seq). Output control parameters were set to Non-symmetrical Venn Diagram Shape and Colored Venn Diagram Fill. 2.10. Enrichment and STRING network analysis Enrichment analysis was done on significantly up (log2FC ≥ 0.5849625, p adj ≤ 0.05) and downregulated (log2FC ≤ -0.5849625, p adj ≤ 0.05) differentially expressed genes (DEGs) (LPS + Dex vs LPS) identified in the RNA-seq dataset. Over-Representation Analysis (ORA) of Reactome Pathways was done using the WEB-based Gene SeT AnaLysis Toolkit (WebGestalt, [95]https://www.webgestalt.org) [96][46]. Analysis was done using ensemble gene IDs as input data, “genome” as Selected Reference Set, and the default advanced parameters provided by the website (statistical significance (FDR ≤ 0.05) was calculated using Benjamini-Hochberg Procedure). Plots of the Top 10 ORA Reactome Pathways were generated in RStudio (RStudio Team (2022). RStudio: Integrated Development Environment for R. RStudio, PBC, Boston, MA URL [97]https://www.rstudio.com/.) using the ggplot2 package (Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978–3-319–24277-4, [98]https://ggplot2.tidyverse.org.). STRING network analysis was performed in Cytoscape 3.9.1 [99][47] using the Cytoscape stringAPP [100][48] (medium confidence cutoff, 0.400). To group the proteins in the network based on their interactions from STRING, clusterMaker2 [101][49] was used to run Markov clustering (MCL) [102][50]. The MCL granularity parameter (infiltration value) was set to 4 and array source was set to use the STRING confidence score attribute (stringdb::score) as weights. All other settings were kept at default. Functional enrichment analysis was performed on the top four subnetworks (clusters 1–4) using the stringAPP with an FDR threshold of 5 % (FDR ≤ 0.05), generating lists of terms spanning 14 categories: GO Biological Process, GO Molecular Function, GO Cellular Component, DISEASE, KEGG Pathways, Reactome Pathways, and WikiPathways. Functional enrichment analysis of the networks from the upregulated “intensified” DEGs resulted in 399 statistically significant terms, and analysis of the networks from the downregulated “intensified” DEGs resulted in 1167 statistically significant terms. Each list of terms was filtered to eliminate redundant terms (using the default redundancy cutoff of 0.5), resulting in reduced lists of 108 enriched terms for the networks from the upregulated “intensified” DEGs and 369 enriched terms for the networks from the downregulated “intensified” DEGs. Enrichment analysis of STRING network clusters of “buffered” DEGs identified 44 enriched terms before filtering and 9 enriched terms after filtering, enrichment analysis of “exclusive” DEGs identified 163 terms before filtering and 24 after filtering. Of these, the four most significant terms (FDR ≤ 0.05) from each analysis were chosen. 2.11. Statistical analysis For differential expression analysis with DESeq, the Wald test was used when comparing the two treatment conditions. In addition, a deltaTE equivalent interaction coefficient that is equal to translation change and independent of transcriptional change was introduced into the DESeq statistical model [103][24]. In ActivePathways, Brown's extension of the Fisher's combined probability test was utilized to generate merged p-values from multiple datasets. The regulatory link generated by GENIE3 was ranked by a random forest algorithm [104][43]. 3. Results 3.1. Quality control analysis of ribo-seq data The human monocyte cell line THP-1 was utilized to examine the overall effects of the synthetic glucocorticoid Dexamethasone (Dex) on the transcriptional and post-transcriptional regulation of gene expression in the context of macrophage innate immunity. Phorbol 12-myristate 13-acetate (PMA) was used to differentiate THP-1 cells into macrophages (henceforth referred to as macrophages), and then LPS stimulation was used to establish a pro-inflammatory M1 macrophage phenotype. The modulatory effect of Dex on LPS induced macrophages was analyzed using an integrated approach combining the mRNA sequencing and Ribosome profiling ([105]Fig. 1A). In general, ribo-seq data provides a global snapshot of all the translating ribosomes and predicts the protein abundance in a population of cells at a given time. Fig. 1. [106]Fig. 1 [107]Open in a new tab Experiment summary, computational analysis, and subsequent quality check of ribo-seq data (A) Overview of the ribosome profiling experimental workflow and data analysis. The processed data include five biological replicates for each treatment condition for ribo-seq and four biological replicates for each treatment condition for RNA-seq. (B-D) ribo-TISH analysis for quality control of ribo-seq data of PMA-differentiated THP-1 derived macrophage-like cells treated with control vehicle, LPS, or a combination of LPS and dexamethasone (LPS + Dex). Upper panel: RPF length distribution mapped to annotated protein-coding regions. Lower panel: Different quality profiles/metrics for annotated RPFs. The data for first, second, and third reading frame is represented by the colors red, green, and blue, respectively. Each row displays RPFs with their corresponding lengths. Column 1: Distribution of RPF 5′ ends across all codons and three reading frames, indicating the proportion of RPF counts from the dominant reading frame. Column 2: RPF 5′ end count distribution near annotated TISs. P-site offset and the ratio between RPF counts at annotated TISs and sum of RPF counts near the annotated translation initiation sites (TISs) after correction for P-site offset are shown. Column 3: RPF 5′ end count distribution in close proximity to annotated stop codon. Column 4: RPF count distribution across three reading frames in protein-coding regions (CDS). (For interpretation of the references to