Network- and systems-based re-engineering of dendritic cells with non-coding RNAs for cancer immunotherapy.

Dendritic cells (DCs) are professional antigen-presenting cells that induce and regulate adaptive immunity by presenting antigens to T cells. Due to their coordinative role in adaptive immune responses, DCs have been used as cell-based therapeutic vaccination against cancer. The capacity of DCs to induce a therapeutic immune response can be enhanced by re-wiring of cellular signalling pathways with microRNAs (miRNAs).
Methods: Since the activation and maturation of DCs is controlled by an interconnected signalling network, we deploy an approach that combines RNA sequencing data and systems biology methods to delineate miRNA-based strategies that enhance DC-elicited immune responses.
Results: Through RNA sequencing of IKKβ-matured DCs that are currently being tested in a clinical trial on therapeutic anti-cancer vaccination, we identified 44 differentially expressed miRNAs. According to a network analysis, most of these miRNAs regulate targets that are linked to immune pathways, such as cytokine and interleukin signalling. We employed a network topology-oriented scoring model to rank the miRNAs, analysed their impact on immunogenic potency of DCs, and identified dozens of promising miRNA candidates, with miR-15a and miR-16 as the top ones. The results of our analysis are presented in a database that constitutes a tool to identify DC-relevant miRNA-gene interactions with therapeutic potential (https://www.synmirapy.net/dc-optimization). Conclusions: Our approach enables the systematic analysis and identification of functional miRNA-gene interactions that can be experimentally tested for improving DC immunogenic potency. © The author(s).

Introduction

Dendritic cells (DCs) play an important role in regulating adaptive immunity by presenting antigens to T cells 1. Due to the unique function of DCs in the coordination of adaptive immune responses, they have been tested as cell-based vaccination against tumours 2-5. To obtain immunogenic potency ex vivo, monocyte-derived DCs need to go through a complex maturation process, in which DCs are exposed to a monocyte-conditioned medium or a cocktail of cytokines 6,7. These treatments result in various phenotypic changes in DCs, such as upregulation of co-stimulatory surface markers (e.g., CD80 and CD40) and secretion of pro-inflammatory cytokines (e.g., IL-12 and TNFα). The matured DCs loaded with cancer antigens are infused into patients and trigger a selective immune response by migrating into the peripheral lymphatic tissue, where they encounter and activate tumour-specific T cells 8.

The capacity of DCs to induce an immune response can be improved by molecular engineering. Pfeiffer and co-workers enhanced DCs through electroporation with mRNA encoding a constitutively active variant of IKKβ (caIKK), a kinase upstream of NF-κB that is a key regulator of the immune response 9,10. Specifically, the kinase phosphorylates IκB resulting in desequestration of the transcription factor (TF) NF-κB and its translocation into the nucleus, where it regulates expression of immune-related genes such as cytokines 11,12. The engineered caIKK promotes constant activation of NF-κB signalling, and the cells expressing it (hereafter labelled caIKK-DCs) displayed increased expression of co-stimulatory molecules and pro-inflammatory cytokines 10,13. When those cells were in addition loaded with the Melan-A melanoma antigen, they were able to induce repeated expansion of Melan-A-specific cytotoxic T cells with a memory phenotype 10. Besides, caIKK-DCs displayed an increased ability to activate autologous NK cells 14. They are currently under evaluation as vaccine in a phase-I clinical trial for the treatment of uveal melanoma patients (NCT04335890).

Our hypothesis is that DCs can be further improved using non-coding RNAs, in particular microRNAs (miRNAs), interacting with key regulators of DC activation and maturation. miRNAs are a class of small endogenous non-coding RNAs with a length of 22-25 nucleotides. Through the inhibition and modulation of the transcription and translation of specific protein-coding genes, miRNAs can alter the basal state of cells and the outcome of stimulatory events 15,16. Increasing evidence shows that miRNAs play a crucial role in the development and function of DCs 17-19. They serve as important regulators of complex networks by targeting key signalling genes to regulate proliferation and cell death as well as homeostasis 20. It has also been found that miRNAs are pivotal in both adaptive and innate immunity, e.g., by controlling the differentiation of immune cell subsets and their immunological functions 21. In particular, miRNAs can modulate the immune response by inducing apoptosis, affecting homeostasis, and changing the cytokine profile of DCs 22. Further, one can use miRNAs, alone or in combination, in therapeutic setups to inhibit expression of selected genes in cancer and other targeted cells 23-26.

To facilitate the re-wiring of DCs, it is crucial to understand the intracellular regulatory processes involved in DC maturation and activation. However, the regulatory networks eliciting the activation and maturation of DCs involve multiple interconnected signalling and transcriptional circuits, and their understanding and proper manipulation requires the combined use of high-throughput data and systems biology methods 27,28. We here present a systems biology approach to understanding the role that miRNAs play in regulating the function of DCs in immunotherapy (Figure 1), and exploit this knowledge to enhance their potential to stimulate an immune response using miRNAs.

Figure 1

A systems biology approach to study miRNA regulation in DCs. The study starts with preparing of donor DCs followed by cocktail maturation and subsequent electroporation of mRNA encoding Melan-A (control DCs) or a constitutively active variant of IKKβ (caIKK-DCs). The obtained RNA-seq data are processed and analysed for annotating and quantifying protein-coding genes and miRNAs. The identified differentially expressed genes in DCs are used for pathway enrichment analyses and reconstruction of gene regulatory networks. A network topology-oriented scoring model is employed to prioritize miRNAs in different pathway categories of DCs. Finally, a literature review of the top ranking miRNAs in immune signalling pathways elaborates their potential function for improving the immunogenic potency of caIKK-DCs.

In this study, we chose the caIKK-DCs as an ideal model system to identify miRNAs that are involved in DC activation via NF-κB signalling and can boost pro-inflammatory signals. We think the identified miRNAs can enable the DCs to repetitively stimulate T cell expansion. To this end, we performed RNA sequencing (RNA-seq) to obtain the transcriptomic profile (i.e., protein-coding genes and miRNAs) of caIKK-DCs in relation to standard DCs. Next, we analysed miRNA-gene interactions at the pathway level and reconstructed regulatory networks underlying immunological functions of DCs. We then performed network-based prioritization of miRNAs by integrating their expression profiles and their strength of association with other protein-coding genes.

Our analysis identified dozens of miRNA candidates in the regulation of caIKK-DCs, with miR-15a-5p and miR-16-5p as prominent examples. We showed that both miRNAs may exert a strong regulatory effect on genes involved in NF-κB signalling and also target chemokines and cytokines regulating T-cell responses. Moreover, we delineated molecular mechanisms through which the miRNAs alter the immunogenic potency of caIKK-DCs. The results of our analysis are available in a web database that facilitates their exploration and visualization (https://www.synmirapy.net/dc-optimization), thereby providing researchers with a tool to select functional miRNA-gene interactions with therapeutic potential in DCs for experimental investigation.

Materials and methods

Generation of monocyte-derived dendritic cells from blood samples

Monocyte-derived DCs were generated as previously described 29. In brief, blood samples from seven healthy donors were collected after approval was granted by the responsible institutional review board (Ethikkommission der Friedrich-Alexander-Universität Erlangen-Nürnberg, Ref. no. 4158) and written informed consent was obtained. Peripheral blood mononuclear cells were isolated from whole blood using density gradient centrifugation. Monocytes were extracted from the non-adherent fraction by plastic adherence, were cultured in DC medium (RPMI; Lonza, Verviers, Belgium) containing 1% non-autologous human plasma (Sigma-Aldrich, St Louis, USA), 2 mM L‑glutamine (Lonza), and 20 mg/L gentamicin (Lonza), and were differentiated into DCs by application of 800 U/mL GM‑CSF (Miltenyi Biotec, Bergisch Gladbach, Germany) and 250 IU/mL IL‑4 (Miltenyi Biotec) on days 1, 3, and 5. After six days in culture, DCs were matured for 24 hours using a cytokine cocktail consisting of 200 IU/mL IL‑1β (CellGenix, Freiburg, Germany), 1 000 U/mL IL‑6 (CellGenix), 10 ng/mL TNF (Beromun, Boehringer Ingelheim Pharma, Germany), and 1 μg/mL PGE2 (Pfizer, Zurich, Switzerland).

RNA in vitro transcription and DC electroporation

Generation of in vitro transcribed RNA with mMESSAGE mMACHINE™ T7 ULTRA transcription-kits (Thermofisher scientific, Waltham, USA) and the electroporation of cocktail-matured DCs (cmDCs) was carried out as previously described 13. For transcriptome analyses, cmDCs were electroporated using 30 µg RNA encoding constitutively active IKKβ (caIKK) since its introduction into mature DCs has been shown to improve activation of T cells 10 and natural killer cells 14. DCs electroporated with RNA encoding melanoma antigen recognized by T cells 1 (Melan-A) were used as control, and such DCs were shown to have no influence on the DCs' transcriptome profile 30. After electroporation, DCs were cultured in DC medium containing GM‑CSF and IL‑4 at the concentrations indicated above.

RNA sequencing processing and differential gene expression analyses

Total RNA including small RNAs was extracted four hours after electroporation using the RNeasy Plus Mini Kit (QIAGEN GmbH, Hilden, Germany) and the generated samples were sequenced using an Illumina HiSeq-2500. Demultiplexed reads were filtered for ribosomal RNAs, transfer RNAs, mitochondrial rRNAs, and mitochondrial tRNAs. The reads were aligned to the human reference genome (hg19) using STAR (v2.5.2b) and assigned to genes using Subread (v1.5.2). Only uniquely mapping reads that could unambiguously be assigned to a single gene were considered for analysis (Supplementary Table S1).

For miRNA expression quantification, we performed a quality check of the RNA-seq reads with FastQC 31, mapped the short sequences to the human reference genome (hg19) using BWA 32, and calculated raw read counts of mature miRNAs that are known and annotated in miRBase v21 33 (Supplementary Table S2; See Supplementary Materials for details).

Before differential expression analysis, we aggregated read counts of Ensembl identifiers that represent the same gene and discarded genes with less than 5 read counts in any sample to increase power for detecting differentially expressed genes 34,35. Next, we used DESeq2 36 in R version 3.6.3 37 to assess differential expression for protein-coding genes and miRNAs. Then, we performed independent filtering on the results to remove genes that have no or little chance of showing significant evidence (Supplementary Table S3 and S4). Specifically, the independent filtering uses the mean of normalized counts as a threshold to optimize the number of adjusted p-values ≤ 0.05 36. If the normalized expression of a gene was lower than the threshold, it was discarded. The Benjamini-Hochberg method was then used on the set of remaining genes to correct for multiple comparisons 38. Genes with adjusted p-values ≤ 0.05 were regarded as significantly differentially expressed.

Gene set enrichment analysis

We extracted all curated pathways from the Reactome pathway knowledgebase (release 68) 39 together with their hierarchical and biological classification according to the database developers. We retraced Reactome's pathway hierarchy by assigning every pathway from Homo sapiens to its corresponding root categories, such as signalling transduction and immune system (see Supplementary Materials for details). As a result, we obtained a table of Reactome pathways matched to the 26 root categories (Supplementary Figure S3 and Supplementary Table S5).

We applied a competitive gene set test to perform gene set enrichment analyses for Reactome pathways. The algorithm CAMERA 40 tests whether the genes in the set are lowly or highly ranked in terms of differential expression relative to genes not in the set, with a positive gene set score indicating a shared tendency for upregulation of the corresponding genes, and vice versa (see Supplementary Materials for details). All genes identified as differentially expressed from our RNA-seq data were used as the background gene list for the enrichment analysis. All obtained p-values were corrected using the Benjamini-Hochberg method. Pathways with false discovery rate (FDR) ≤ 0.05 were regarded as significantly up- (positive score) or down-regulated (negative score) in our comparison of caIKK-DCs with controls (Supplementary Table S5). The gene set enrichment analysis was performed using the CAMERA implementation in the package limma 41 in R.

Regulatory network reconstruction

We downloaded functional interactions from the Reactome database (release 68). The collection includes protein-protein interactions, transcriptional regulation, gene co-expression, protein domain interaction, gene ontology annotations and text-mined protein interactions, which cover almost half of the human proteome 42. There are different types of directional molecular interactions including: activation, inhibition or repression, and co-expression or complex formation. The biochemical reactions covered are phosphorylation and ubiquitination. We processed the list to transform bidirectional interactions into their two unidirectional constituents. The result list contained 435,043 unidirectional interactions among 13,852 protein-coding genes.

To derive miRNA-gene interactions, we first obtained conserved and non-conserved miRNA binding sites as predicted by Targetscan version 7.2 43. Then, we filtered the predicted interactions with experimental evidence from miRTarBase version 8.0 44 and starBase version 2.0 45. By doing so, we obtained a list of miRNA-gene interactions (Supplementary Table S6) that contain putative miRNA binding sites with experimental support, such as high-throughput experiments (e.g., RNA-seq, microarray, and AGO CLIP-sequencing) and/or low-throughput experiments (e.g., q-PCR, reporter assay, and Western blot).

The obtained lists of protein-coding genes' functional interactions and miRNA-gene interactions were used to reconstruct gene regulatory networks from Reactome pathways. Specifically, for a pathway of interest, we built a network around its participating genes and calculated the pairwise Pearson correlation coefficients between interaction partners from their normalized count values in caIKK-DCs. The normalized counts were obtained using the regularized logarithm method 36. In the reconstructed networks, we used the Pearson correlation coefficients to filter out interactions that disagree with their regulation type. We assumed that positive interactions (i.e., activation) require positive Pearson correlations and negative interactions (i.e., inhibition) negative Pearson correlations between interacting molecules. Interactions annotated as gene co-expression or formations of protein complexes were kept, assuming that the involved genes can affect each other's expression or activity in both directions.

Furthermore, we added annotation to the reconstructed networks' components in the form of differential expression profiles (i.e., fold-change and FDR), types of genes (e.g., protein-coding gene or miRNA), gene interaction types (e.g., functional interaction or post-transcriptional regulation), gene interaction strengths as denoted by the Pearson correlation coefficients introduced above, and immune categories of genes. Immune categories of protein-coding genes were annotated using curated data from the Immport database 46. Data from the TcoF database were used to identify TFs in our networks 47.

Gene prioritization in regulatory networks

We prioritized genes in a network using SANTA in R 48. The algorithm determines a score of relative importance for each node in a network through a clustering model that accounts for network topology (distances between nodes) and node weights (in our case, a measure of differential expression called perturbation). Briefly, a gene is assigned a high score when itself and its close neighbours in the network have a higher-than-average node weight. The closeness, or distance, between genes is calculated by finding the shortest path through the network.

The node weight is given by the gene's perturbation (i.e., -log10(adj-p) • |log2(fold-change)|). Both adjusted p-value and log2 fold-change of the gene were taken from the differential gene expression analysis. The distances between neighbouring nodes were calculated as 1 - |p|, where p represents the Pearson correlation coefficient between the two interacting genes. Higher correlation coefficients (i.e., higher interaction strengths) correspond to lower edge length and thus shorter distance between the nodes. The calculated score was used to prioritize genes (see Supplementary Materials for details). As our networks contain both miRNAs and protein-coding genes that have different types of interactions, miRNAs and protein-coding genes were ranked separately.

Mapping of microRNA-gene interactions into the curated DC network

To generate the curated DC network, we made use of a previously published network of macrophage pathways 49, since macrophages and DCs are generally considered to be quite similar. We manually added pathways for antigen processing and presentation that were not present in the macrophage map through a comprehensive database search in Reactome. The enriched DC network is a restructuring (see 49 for details on the algorithm) of the curated version that also incorporates miRNA-gene interactions identified in this study. The reconstructed network is accessible at https://vcells.net/dendritic-cell.

MicroRNA cooperativity analysis

We used the TriplexRNA database 50 to identify miRNAs that can cooperate with significantly differentially expressed miRNAs in caIKK-DCs to repress protein-coding genes of interest. The obtained RNA triplexes were further filtered using pre-computed equilibrium concentrations and minimum free energies. We kept the RNA triplexes with equilibrium concentrations ≥ 50 nM and minimum free energies ≤ -25 kcal/mol and regarded the participating miRNAs as cooperative partners to repress protein-coding genes.

Data visualization

The gene regulatory networks for significantly enriched pathways were drawn using ggraph 51 and igraph 52 in R. Heat maps were plotted using Complexheatmap 53 in R. Scatter and bar plots were drawn using ggplot2 54 in R. Sankey diagram was drawn using networkD3 55 in R. The clustered Reactome pathways were visualized using Cytoscape version 3.72 56.

Results

Transcriptome analysis reveals miRNA expression changes in caIKK-DCs

To characterize the gene expression profiles induced by a constitutively active IKKβ in DCs, we collected blood samples from seven healthy donors and generated monocyte-derived DCs (see Materials and Methods). The DCs were matured with a cytokine cocktail of IL-1β, IL-6, TNF, and PGE2. The cmDCs were split for the subsequent experiments. One half was electroporated with mRNA encoding constitutively active IKKβ (caIKK, encoded by an engineered IKBKB). The other half was electroporated with mRNA encoding the melanoma antigen Melan-A (encoded by MLANA). The Melan-A protein was used because it is non-functional in DCs as it is naturally involved in melanocyte-specific pathways 57,58. Additionally, we have previously shown that maturated DCs transfected with Melan-A mRNA did not show any significant changes in their transcriptomic profiles compared to mock-electroporated DCs 30. To confirm that caIKK electroporation of cmDCs induced activation in DCs, while electroporation of cmDCs with the control RNA did not, we analysed their phenotypes and cytokine secretion profiles. In caIKK-DCs, we verified the secretion of pro-inflammatory cytokines such as IL-6, IL-8, IL-12, and TNFα (Supplementary Figure S1A), upregulation of co-stimulatory surface markers such as CD25, CD40, CD70, CD80, CD86, and OX40L (Supplementary Figure S1B), and showed that the DC preparations used for RNA sequencing were successfully transfected as indicated by the expression of CD25 and CD70 in the seven donors (Supplementary Figure S1C).

Four hours after electroporation, RNA was isolated and assessed via bulk RNA sequencing (RNA-seq). We chose this early time point because we were interested in mRNA levels, which are expected to quickly respond to the activation of NF-κB as a result of continuous IKKβ activation. From the RNA-seq data, we identified 63 protein-coding genes and 44 miRNAs that were significantly differentially expressed (DE) between caIKK-DCs and controls (Supplementary Table S3 and S4; see Materials and Methods). Among the protein-coding genes, MLANA (encoding Melan-A) and IKBKB (encoding IKKβ) were the most down- and upregulated in caIKK-DCs, respectively (Supplementary Figure S2A). This is in line with the fact that the mRNA content of the two genes was artificially altered in the respective populations and can be considered a quality control for the experimental results. For the miRNAs, miR-146a/b and miR-155 were upregulated in caIKK-DCs, in consistence with them being transcriptional targets of NF-κB 59 and being upregulated in mature DCs 60. By performing principal component analysis, we assessed the clustering tendency in the RNA-seq data. Controls and caIKK-DCs showed better separation when restricting the input to the measured miRNAs rather than the whole transcriptome (Supplementary Figure S2B). In addition, the DE miRNAs unequivocally separated the caIKK-DCs from the controls in hierarchical clustering (Supplementary Figure S2C). These results suggested that caIKK-DCs harbour a distinct miRNA expression profile.

The gene signature induced in caIKK-DCs is associated with NF-κB activation

To understand the molecular function of the identified DE genes in the caIKK-DCs, we performed gene set enrichment analysis using the Reactome pathway database. The database contains more than 2,000 cellular pathways curated from 30,721 peer-reviewed publications and classified into 26 root categories 39, thereby enabling a systematic and comprehensive analysis of DE genes. The 26 categories consist of a set of pathways that are annotated to be hierarchically and functionally linked. We calculated enrichment scores for each pathway which reflect the degree to which its corresponding gene set tends to be up- or downregulated in caIKK-DCs (Figure 2A; see Materials and Methods).

Figure 2

Gene set enrichment analysis on the differential expression analysis between caIKK-enhanced and control DCs. (A) Enrichment plots of two exemplary gene sets. The genes' weighted log2 fold-change (i.e., log2 fold-change divided by the standard error of the log2 fold-change) obtained from the differentially gene expression analysis are sorted from the smallest to the largest (the barcode plot). Genes from TNF signalling (red bars) accumulate at the upregulated end while genes from ligand-receptor interaction (blue bars) in Hedgehog signalling accumulate at the downregulated end. The curves show the local enrichment score of the vertical bars in the barcode plot. For the red curve, parts above the dashed line signify enrichment while parts below the line signify depletion. For the blue one, parts below the dashed line signify enrichment while parts above the line signify depletion. (B) Network visualization of the gene set enrichment results for Reactome's 26 root categories (nodes annotated with texts). Categories are connected when they have shared genes (edges). The size of a node denotes the number of pathways that belongs to the respective category. The node colour represents the number of significantly enriched pathways in a category. A grey node means that no significant pathways were identified in the category. The width and colour of an edge represent the number of shared genes between the two connected categories. A detailed map containing all pathways and their corresponding enrichment scores can be found in Supplementary Figure S3.

We found that the DE genes in caIKK-DCs are significantly enriched in 195 Reactome pathways, most of which belonged to the Reactome categories signal transduction and immune system (Figure 2B; Supplementary Figure S3 and Table S5). In the category of immune system, 65 out of 182 pathways were identified as significantly enriched, including cytokine signalling pathways and pathways associated with innate and adaptive immune response. This suggested that the continuous activation of IKKβ in DCs has a generalized effect on DC-mediated immune responses. In addition, we identified 12 enriched pathways that are directly associated with NF-κB activation and signalling (see Supplementary Table S5), such as NF-κB activation by IκB kinase complexes. This was consistent with the current model of the canonical NF-κB activation pathway, in which IKKs phosphorylate IκB resulting in desequestration of NF-κB and its translocation into the nucleus, where it regulates expression of immune-related cytokine genes and others 11,12. All these NF-κB pathways had positive enrichment scores, indicating that the genes involved are more likely to be up-regulated in the caIKK-DCs. The results were in line with our expectation that the caIKK-DCs can trigger a stronger immune response as a result of NF-κB desequestration by constitutive activation of IKKβ.

Significantly differentially expressed miRNAs in caIKK-DCs regulate an abundance of enriched immune pathways by targeting hundreds of their protein-coding genes

To identify potential miRNA-gene interactions regulating the immunogenic potency of DCs, we first obtained putative miRNA-gene interactions for the significantly DE miRNAs (see Materials and Methods). We kept the putative interactions that are validated by experiments. For each identified miRNA-gene interaction, we then computed the Pearson correlation coefficient between miRNA and target gene expression. The interactions with negative correlation were regarded as reliable and functional, as miRNAs canonically repress translation initiation or stimulate mRNA degradation 61 and miRNA-mediated gene activation usually results from indirect regulation mechanisms 62. The data showed that 36 out of the 44 miRNAs are involved in the regulation of protein-coding genes belonging to 195 enriched pathways of the 26 Reactome root categories (Figure 3A).

Figure 3

miRNA targeting profiles in Reactome pathway categories. (A) Overview of the 36 DE miRNAs (in columns) targeting the 26 Reactome pathway categories (in rows). On the heat map grid, the number of protein-coding genes targeted by the respective miRNA is given, and the colour represents the category's number of significantly enriched pathways that are regulated by the miRNA. For example, if an entry shows n with a colour corresponding to m on the figure legend, it means a miRNA regulates n targets in m significantly enriched pathways of a category. The bar plot on the left indicates the per-category percentage of the protein-coding genes (black bars) that were found in our RNA-seq data, with the total number of molecules in the category given next to it. The bar plot on the right indicates the percentage of enriched pathways (black bars) per category, with the total number of pathways in the category given next to it. The figures at the bottom tabulate how many categories a miRNA regulates. The top annotation shows statistics from the differential expression analysis for the miRNAs (i.e., fold-change in log2 scale and FDR). (B) The network shows miRNA-gene interactions in the category signal transduction. The four miRNAs (miR-16, miR-15a, miR-20a, and miR-424) that have the largest number of targets were selected. The node size is proportional to the node degree. The node colour represents a gene's fold-change in log2 scale. The node shape denotes the type of a gene, including protein-coding (square) and miRNA (circle), with their names shown in blue and black labels, respectively. TFs are drawn as diamonds with their names shown in purple font. The colour of node borders represents different categories of annotated immune genes, with the gene names given as labels. The edge colour shows Pearson correlation coefficients between the expression of miRNAs and that of their targets.

In individual pathway categories, the number of molecules (i.e., protein-coding genes, DNA/RNA, drugs, and chemical compounds) ranged from 27 to 2727, and genes identified by our RNA-seq data covered between 51% and 95% of the molecules in the respective category. For all categories except digestion and absorption, the DE miRNAs were identified to regulate 1 to 103 protein-coding genes (denoted by the numbers on the heat map grid cells in Figure 3A). Some miRNAs were found to regulate the expression of dozens of protein-coding genes in more than ten categories, suggesting that they act as regulatory hubs in caIKK-DCs. For example, miR-15a-5p, miR-16-5p, miR-20a-5p, and miR-424-5p can potentially regulate more than 80 protein-coding genes in the category signal transduction (Figure 3B), and they also have more than 60 targets in immune system. In contrast, miR-15a-3p and miR-9-3p target only BCL2 and MAPK1 in immune system, suggesting a specific role for them in regulating cytokine signalling in caIKK-DCs (Supplementary Figure S4). Furthermore, we found that some miRNAs target a high fraction of enriched pathways belonging to specific Reactome root categories (denoted by the colour of heat map grid cells in Figure 3A). For instance, miR-16-5p regulates 61 out of 64 enriched immune system pathways and 34 out of 36 enriched pathways in signal transduction. This suggested that it plays a vital role in regulating the immunogenic potency of caIKK-DCs. On the other hand, some categories included abundant enriched pathways that are regulated by multiple DE miRNAs. Interesting cases were pathways associated with protein metabolism, RNA metabolism, programmed cell death, and cell cycle. This result suggests that the DE miRNAs in caIKK-DCs are involved in regulating synthesis, processing, and modification of mRNAs and proteins and can also participate in other biological processes, such as cell cycle and cell apoptosis. Taken together, the DE miRNAs in caIKK-DCs target and potentially coordinate the activity of immune-relevant pathways in a pleiotropic fashion.

Network-based prioritization of miRNAs in caIKK-DCs

The ubiquitous, pleiotropic, and concerted gene regulation by miRNAs makes it challenging to quantify the relative impact of each individual miRNA. To prioritize the DE miRNAs according to their potential to act synergistically with NF-κB in DC activation, we applied a network-based method that integrates their expression and interaction profiles.

First, we reconstructed one gene regulatory network for each of the 26 pathway categories. The reconstructed networks were composed of miRNA-gene interactions and functional interactions among protein-coding genes. Interactions were discarded when the sign of their Pearson correlation coefficient of expression disagreed with their regulation type, such as inhibition or activation (see Materials and Methods). Depending on the category, the size of the corresponding networks varied from 1,915 genes and 57,520 interactions (for signal transduction) to 30 genes and 153 interactions (for mitophagy). To prioritize the network components involved in regulating the immunogenic potency of DCs, we used a clustering model 48 to calculate a node score (Figure 4A; see Materials and Methods).

Figure 4

Ranking of miRNA relevance based on expression and regulatory network neighbourhood. (A) Network nodes were scored with an algorithm that uses the guilt-by-association principle to rank genes. In other words, a gene inside of or close to a cluster of important genes is potentially more important than a gene that is further away. In our case, the importance of a gene for a phenotype is quantified by their perturbation in expression (denoted by node colours). In a network, the distance of the gene in question (red or blue border) to other genes is calculated via the weighted shortest-path method. The range of observed distances is then subdivided into discrete bins (denoted by circles in the figure) and an estimate of neighbourhood importance calculated for each bin, i.e., for all nodes up to the respective distance. The area under the curve (AUC) in the plot of bins (d) vs neighbourhood importance is used as the gene's score. In the example, gene k is much closer to genes with high node weights (i.e., their perturbation in caIKK-DCs is large) than gene j. As a result, the blue area is bigger than the red area, and thus gene k ranks higher than j. (B) Heat map of miRNA ranking in pathway categories. The columns of the matrix indicate the 36 DE miRNAs sorted by perturbation (i.e., node weights), and the rows of the matrix indicate 25 categories of Reactome pathways. The category digestion and absorption is not shown, as we did not identify functional and reliable interactions among its genes. On the heat map grid, the rank of a miRNA in a category is given as a number, and the colour represents the number of protein-coding genes targeted by it. For example, if an entry shows k with a colour corresponding to n on the figure legend, it means that the miRNA ranks kth (1st is the highest ranking) and regulates n targets in the category. A white grid cell means that the miRNA has no targets in the category and thus no ranking. The top annotation shows node weights of the miRNAs. The numbers in parentheses on the left side list how many genes and edges the reconstructed regulatory network of the category possessed. The box plots on the left show the distribution of edge weights (denoted by Pearson correlation coefficients between genes) in the networks. The bar plots on the right show the Pearson correlation coefficients between a miRNA's perturbation and its score. The numbers at the bottom show the number of times miRNAs ranked 1st in the pathway categories.

The score ranked IKBKB, whose expression was greatly increased by mRNA electroporation, as the top protein-coding gene in 51 out of 59 networks in which it is involved (Supplementary Table S7; Supplementary Figure S5). This result is consistent with our expectation that the intentionally modulated gene in experiments is prioritized, and thus demonstrating the ability of the model to identify crucial regulatory genes in the experiments. In the two prominent categories signal transduction and immune system, the NF-κB family and genes related with immune signalling or antigen processing and presentation tended to rank higher than other genes (Supplementary Figure S6). This result again justified the ability of the model to prioritize important genes in networks, as members of the NF-κB family are downstream targets of IKBKB while signalling and antigen presentation genes are supposed to be crucial regulating immune function of DCs.

Furthermore, we analysed the data to identify crucial miRNAs for each Reactome root category. As shown in Figure 4B, miRNAs with higher node weights (i.e., stronger perturbation) generally ranked higher in a category, as miRNA scores and node weights showed a positive correlation, ranging from 0.19 to 0.98. Specifically, miR-503-5p, miR-503-3p, and miR-146-5p had the highest perturbation in the DE miRNAs, and they ranked top in 22 out of 26 categories. However, the interaction profile also plays a role, as for example in signal transduction, the three top-ranking miRNAs miR-101-3p, miR-16-5p, and miR-15a-5p had lower perturbation than miR-146-5p but interacted with more protein-coding genes. In addition, the three above miRNAs and miR-144-3p ranked top in immune system, most probably due to the reason that they regulate a large number of protein-coding genes associated with immune signalling pathways. To facilitate the visualization of our results, we integrated the data and the identified miRNA-gene interactions into a comprehensive, manually curated regulatory network including key pathways in DC priming and activation according to the literature (see Materials and Methods; https://vcells.net/dendritic-cell).

Taken together, the reconstructed regulatory networks underlying different cell functions allowed us to identify important miRNA regulators based on their expression and interaction profiles. The miRNAs with the highest scores possibly exert regulatory functions, and manipulation of their expression levels may enhance the immunogenic potency of DCs.

Potential miRNA-gene interactions to improve caIKK-DCs

To characterize the functional role that miRNAs play in caIKK-DCs, we delineated landscapes of miRNA-gene interactions in the significantly enriched pathways that were found in corresponding categories (see Figure 5 and https://www.synmirapy.net/DC-optimization). The interaction landscapes are a way of systematically mapping relevant gene interactions, and in our case they served as a tool for identifying functional miRNA-gene interactions in DC priming and activation. As we were particularly interested in identifying miRNAs that can enhance the caIKK-DCs' immunogenic potency, we focused on analysing miRNA-mediated gene regulation in the category immune system. In this category, we identified hundreds of miRNA-gene interactions in significantly enriched pathways, including toll-like receptor, cytokine, and interleukin signalling as well as MHC processing and presentation. All of these pathways had positive enrichment scores, indicating that the involved genes tended to be upregulated in caIKK-DCs according to our analysis. Most protein-coding genes used as indicators of DC activation and maturation 5,60,63,64 were found to be upregulated in the enriched pathways (Supplementary Table S8). The activation of NF-κB signalling led to upregulation of surface proteins that can prime T cells (e.g., CD40, CD70, CD80, and CD86), chemokines (e.g., CCL3 and CXCL10) that are necessary for T-cell migration, TNF superfamily members that can induce crosstalk between T cells and DCs (e.g., TNF, TNFRSF4, and TNFSF9), and cytokines that are responsible for stimulating proliferation and activation of T cells (e.g., CXCL8, IL6, IL12A, and IL12B).

Figure 5

Landscape of miRNA-mediated DC gene regulation in immune signalling pathways. The heat map has two components that share a set of columns corresponding to 98 DE protein-coding genes that are targeted by the DE miRNAs. In the upper component, rows represent pathways from the category immune system, and grid cell colours indicate whether a protein-coding gene is involved in the pathway (grey), involved in the pathway and an immune gene (grey grid cells with red borders), or not involved in the pathway (white). The figures next to a pathway name indicate how many DE immune genes (left) and how many protein-coding genes found in our RNA-seq data (right) belong to it. The top annotation highlights genes with different characteristic immune function using a colour code. The annotation on the right side shows the statistics of the gene set enrichment analysis including the enrichment score and the FDR. The bar plot between the heat map components shows the log2 fold-change of the genes in caIKK-DCs (blue: downregulated; red: upregulated). In the lower component, the rows represent the ranking miRNAs in immune system (from high to low) and the grid cells show the regulative influence of a protein-coding gene by a miRNA, which is estimated by the Pearson correlation coefficients between their expression profiles. If a gene is a known immune gene, the corresponding grid cell has a red border. The numbers in the parentheses next to the miRNA names show the number of DE immune genes and the number of DE protein-coding genes that are regulated by a miRNA. The right annotation shows the results of the differential expression analysis including the log2 fold-change of miRNA expressions and their FDRs. For lack of space, we show only enriched pathways with more than 30 protein-coding genes picked up in the RNA-seq data, and in each pathway, only a subset of protein-coding genes that are estimated to be strongly influenced by the miRNAs (Pearson correlation ≤ -0.3) are shown. The complete landscape of miRNA-gene interactions in immune system is shown in Supplementary Figure S7.

Furthermore, our data showed that the identified DE miRNAs have a regulative influence (represented by Pearson correlation ≤ -0.3) on protein-coding genes associated with NF-κB activation, cytokines, chemokines, and TFs that are associated with immunophenotypes of DCs (Figure 6). Some of the DE miRNAs were found to cooperate with other miRNAs to regulate the expression of a protein-coding gene (see Materials and Methods). This mechanism, known as miRNA cooperativity, is characterized by more efficient inhibitory effects on the target's expression compared to the regulation by individual miRNAs 23,25,26. Moreover, for most of the identified miRNAs, our analysis proposed specific modulation of their expression levels to improve immunogenic potency of DCs. However, in some cases, such as miR- 34a-5p and miR-20a-5p, up- or downregulating their expression levels may result in contradictory effects in DC-mediated immune response, thereby requiring further analysis and experimental investigation. In the following paragraphs, we illustrate and discuss specific functions of the miRNAs in caIKK-DCs.

Figure 6

Potential miRNAs for improving immunogenic potency of caIKK-DCs. The Sankey diagram contains three columns made up of nodes representing the DE miRNAs from our RNA-seq data and their cooperating miRNAs, protein-coding genes targeted by the miRNAs, and DC phenotypes associated with the protein-coding genes, respectively. miRNA pairs that were identified to cooperatively repress a protein-coding gene are connected by a brace. Colours of miRNAs and protein-coding genes indicate whether or not they were significantly DE in caIKK-DCs. Arrows in miRNA nodes indicate how the expression of miRNAs should be manipulated to obtain DCs with higher immunogenic potency: upregulation (↑), downregulation (↓), and no suggestion due to potentially conflicting effects of the miRNA (-). Connections between miRNAs and protein-coding genes show regulative influence of protein-coding genes by miRNAs (strong: Pearson correlation ≤ -0.5; weak: -0.3 ≤ Pearson correlation < -0.5). The connections between protein-coding genes and phenotypes denote how a gene regulates a phenotype according to literature. For instance, miR-424-3p and miR-224-5p target IRF4 that is known to positively regulate differentiation of DCs. The two miRNAs cooperatively repress the protein-coding gene, but the observed downregulation of miR-424-3p results in a decreased inhibitory effect on the expression of IRF4. A detailed discussion of the results can be found in the main text. The corresponding miRNA-gene interactions in immune system as well as annotated gene-phenotype associations can be found in Supplementary Table S9 and Table S10, respectively.

Chemokines direct cell migration via induction of chemotaxis. For example, CCL5 and CXCL10 improve CD8+ T-cell infiltration 65 and CCL20 plays a role in recruiting regulatory T cells and T helper (Th) 17 cells 66,67. CXCL10 is repressed by miR-16-5p and CCL20 is targeted by miR-21a-5p with its cooperating miR-25-3p. This suggested that the miRNAs exert an inhibitory function in the recruitment of T cells. miR-21a-5p also targets IL12A, a subunit of the inflammatory cytokine IL-12 that is necessary for CD8+ T-cell clonal expansion, function and memory 63,68.

Control of DC survival is necessary for maintaining their homeostasis 69,70. We showed that miR-15a-5p with its cooperating miRNAs (i.e., miR-156-5p and miR-876-5p) and miR-20a-5p target BCL2, and miR-101-3p targets MCL1. The repression of the BCL2 family of anti-apoptosis genes by these miRNAs suggested their ability to undermine the survival mechanism of DCs.

miR-424-5p and miR-224-5p can co-repress IRF4, which is a member of the interferon-regulatory family and can regulate differentiation of specific DCs that can induce Th 2 cell responses 71. miR-20a-5p and miR-144-3p regulate the MAPK signalling pathway by targeting MAPK1 and MAP3K8 respectively, and these MAP kinases have been found to activate the IKK complex that triggers NF-κB activation 72,73 and also to regulate release of TNFα by DCs 74. miR-34a-5p has a strong regulative influence on CD44 whose presence is important for the immune synapse between DCs and T cells that subsequently regulates T-cell activation 75 and apoptosis 76. miR-34a-5p also targets TNFAIP8 whose knockdown in DCs has been found to promote DC maturation and activation followed by increased proliferation and differentiation of T cells 77. miR-9-5p can cooperate with miR-139-5p to repress CXCR4 that is required for DC migration into the skin's draining lymph nodes 78. miR-142-3p with its cooperating miR-429 and miR-142-5p target the small GTP-binding protein RAC1 that controls the formation of dendrites in mature DCs and their migration toward T cells 79.

Some identified DE miRNAs target protein-coding genes involved in regulating the DC-mediated secretion of cytokines that are important for the T-cell response. The repression of NFATC3 by miR-424-5p and its cooperating miR-370-3p suggested a regulating influence on the production of IL-2 that is involved in T-cell priming 80,81. STX3 has been shown to play a role in trafficking of IL-6 or MIP-1α in DCs and thus regulating their secretion 82 and is targeted by let-7e-5p, miR-146a-5p, and miR-146b-5p with its cooperative partner miR-519d-3p. The deficiency of IRAK1 in plasmacytoid DCs abrogates IFNα production, leading to a remodulation of T cell function 83-85, and IRAK1 is a target of miR-142-3p.

Finally, miR-16-5p and miR-15a-5p can cooperate with miR-203a-3p to repress IL-15, an interleukin which can induce T-cell proliferation, enhance cytolytic effector cells including natural killer and cytotoxic T cells, and reinforce B-cell stimulation 86. A recent in vivo study has shown that an IL15-enhanced DC vaccine is a potent delayer of tumour growth, improves mouse survival, and induces a stronger Th1-skewed T-cell response 87. The two miRNAs also target the receptor TNFSF9 (also known as CD137), whose stimulation in DCs by its ligand CD137L can lead to secretion of IL-6 and IL-12 and induce T-cell proliferation 63,88,89. In addition, miR-16-5p and miR-15a-5p were identified to strongly repress IKBKB itself. Since both miRNAs were found to be downregulated in caIKK-DCs, this implied a positive feedback loop in NF-κB signalling as the miR-15/16 cluster is a transcriptional target of NFKB1 90. The results suggested both miRNAs as promising candidates for improving the immunogenic potency of caIKK-DCs, as they not only have the ability to strengthen NF-κB activation but also to improve DC-induced immune responses through regulating cytokines and chemokines.

Discussion and Conclusions

We applied a systems biology approach to investigate the regulatory functions of miRNAs in caIKK-DCs. Due to the promiscuous binding of miRNAs, it is challenging to identify relevant miRNA-gene interactions for experimental validation and cell re-engineering 91,92. Our approach, which integrates transcriptomic profiling, networks of curated signalling pathways, and a prioritisation score, allows the systematic identification of condition-specific miRNA-gene interactions.

Through RNA sequencing of monocyte-derived DCs matured with a cytokine cocktail and electroporated with caIKK mRNA, we identified DE protein-coding genes and miRNAs in the caIKK-DCs. The identified DE miRNAs correctly separated the caIKK-DCs from the control, suggesting a well-defined transcriptional response to caIKK that is consistent with our understanding that miRNAs act as post-transcriptional regulators of expression in DC differentiation and function 93. Among the identified miRNAs there were several, such as miR-15a-5p, miR-16-5p, miR-20a-5p, and miR-424-5p, which target a considerable number of genes. Such hubs have been shown to be important regulators, as they represent sites of signalling convergence in gene regulatory networks and coordinate cell development and function 94-96. In contrast, other DE miRNAs, such as miR-15a-3p and miR-9-3p, exert a narrow function by regulating the expression of specific protein-coding genes in the caIKK-DCs.

Integration of the transcriptomic response into the curated pathways from Reactome provides an understanding of the functional changes at the pathway level. The gene set enrichment analysis highlighted cytokine, interleukin, and toll-like receptor signalling pathways that are involved in regulating various aspects of innate and adaptive immune responses 97. Such results may be compromised when other pathway databases such as KEGG 98 and WikiPathways 99 are employed, as the relevant pathways and molecular interactions in the pathways are different from Reactome 100,101. To circumvent this issue, one possibility could be to extract the overlapped networks between the different databases of pathways; however, this is not always possible due to the differences in annotation of genes and interactions. An alternative option is to integrate the data and the detected miRNA-gene interactions into comprehensive, manually curated regulatory networks based on the current literature on DC regulation. This way, one can put the newly discovered relevant interactions into the context of the existing knowledge and facilitate the mining and interpretation of the omics data 102,103. However, when used inappropriately, knowledge-based networks mainly rediscover existing knowledge but may overlook insights gained from the evaluation of an all-encompassing network.

We reconstructed regulatory networks from Reactome pathways and used them to rank genes and miRNAs according to their predicted impact on DC function. Systematic computation of such a ranking supports and facilitates experimental efforts, allowing them to focus on the most promising candidates. Gene prioritization algorithms have been widely used in recent times to rank genes in networks 104,105. For instance, the PageRank algorithm designed to analyse the relative importance of websites was adapted to identify crucial genes in biological networks 106, and diffusion-based methods were used on dense networks to prioritize genes 107. We used a gene prioritization algorithm that utilizes the guilt-by-association principle to rank genes based on their own perturbation, i.e., differential expression profile, and their weighted distances to other perturbed genes in a network. The algorithm prioritized dozens of miRNAs, of which miR-16-5p and miR-15a-5p are the top candidates to regulate the immunogenic potency of DCs.

Finally, an in-depth analysis of the identified miRNA-gene interactions in immune signalling pathways showed diverse roles of the DE miRNAs in regulating DC-mediated immune response. For instance, miR-16-5p and miR-15a-5p may have strong regulatory influence on IKBKB that activates NF-κB and on TNFSF9 that controls cytokine secretion of DCs; both miRNAs could have weak inhibitory effects on BCL2 that maintains DC homeostasis and on cytokines (such as CXCL10 and IL15) that regulates T-cell response but may cooperate with other miRNAs to more efficiently repress the protein-coding genes. While these predictions were made based on validated miRNA binding sites and negative correlation between the expression levels of miRNAs and their targets, they cannot quantify the strength of individual repression effects 108. The results suggested both miRNAs as potential candidates for improving immunogenic potency of caIKK-DCs through strengthening NF-κB stimulation and also synergistically regulating other genes related with immunogenic potency.

For most identified crucial miRNAs, our analysis suggested up- or downregulation of their expression levels to improve immunogenic potency of DCs, but in some cases the pleiotropic nature of miRNAs in regulating gene expression makes it difficult to decide how to experimentally modulate their expression. In addition, it is worth noting that the results reflect the early transcriptional response that may differ from that in the long-term. From an experimental perspective, the next step would be to analyse the kinetics of the expressions of miRNA and mRNA after the activation of NF-κB. Further, it remains to be tested how co-electroporating the selected miRNAs, or artificial antagonists thereof, with caIKK or introducing them into the cells after a delay will alter the DCs' phenotype and immunogenic potency.

Taken together, our approach enables the systematic analysis and identification of functional miRNA-gene interactions that can be experimentally tested for improving DC immunogenic potency. Since the results were produced in a computationally reproducible manner and were stored in a public database, experimental tests of the predictions can be performed in the future not only by our group but also by other researchers working on DC-based cancer immunotherapy. Additionally, since the approach is not specific for DCs, it can be adapted to study miRNAs in other immune cells and relevant immunotherapies.

Supplementary figures.

Click here for additional data file.

Supplementary tables.

Click here for additional data file.