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.


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][3][4][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.

Ivyspring International Publisher
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][18][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][24][25][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.
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.

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 proteincoding 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 miRNAgene 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].

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 mockelectroporated 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 peerreviewed 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).  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).
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 miRNAgene 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).
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). 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. 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.
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). 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. 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. 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.
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 proteincoding 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][84][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][95][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 proteincoding 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.