miR-15a and miR-15b modulate natural killer and CD8+T-cell activation and anti-tumor immune response by targeting PD-L1 in neuroblastoma

Neuroblastoma (NB) is an enigmatic and deadliest pediatric cancer to treat. The major obstacles to the effective immunotherapy treatments in NB are defective immune cells and the immune evasion tactics deployed by the tumor cells and the stromal microenvironment. Nervous system development during embryonic and pediatric stages is critically mediated by non-coding RNAs such as micro RNAs (miR). Hence, we explored the role of miRs in anti-tumor immune response via a range of data-driven workflows and in vitro & in vivo experiments. Using the TARGET, NB patient dataset (n=249), we applied the robust bioinformatic workflows incorporating differential expression, co-expression, survival, heatmaps, and box plots. We initially demonstrated the role of miR-15a-5p (miR-15a) and miR-15b-5p (miR-15b) as tumor suppressors, followed by their negative association with stromal cell percentages and a statistically significant negative regulation of T and natural killer (NK) cell signature genes, especially CD274 (PD-L1) in stromal-low patient subsets. The NB phase-specific expression of the miR-15a/miR-15b-PD-L1 axis was further corroborated using the PDX (n=24) dataset. We demonstrated miR-15a/miR-15b mediated degradation of PD-L1 mRNA through its interaction with the 3'-untranslated region and the RNA-induced silencing complex using sequence-specific luciferase activity and Ago2 RNA immunoprecipitation assays. In addition, we established miR-15a/miR-15b induced CD8+T and NK cell activation and cytotoxicity against NB in vitro. Moreover, injection of murine cells expressing miR-15a reduced tumor size, tumor vasculature and enhanced the activation and infiltration of CD8+T and NK cells into the tumors in vivo. We further established that blocking the surface PD-L1 using an anti-PD-L1 antibody rescued miR-15a/miR-15b induced CD8+T and NK cell-mediated anti-tumor responses. These findings demonstrate that miR-15a and miR-15b induce an anti-tumor immune response by targeting PD-L1 in NB.

INTRODUCTION
Neuroblastoma (NB) is the most common pediatric cancer affecting children younger than 5 years. 1-3 NB develops from immature nerve cells, most commonly in adrenal glands situated above the kidney. NB accounts for 15% of childhood cancer-related mortality, and approximately 50% of children treated for high-risk NB have more aggressive tumor relapse with less than 20% 5-year overall survival. [4][5][6][7] During embryogenesis, multipotent neural crest cells differentiate into multiple cell types, including sympathetic neurons, adrenal medullary cells, and modified postganglionic sympathetic neurons. 8,9 The deregulation of signaling pathways involved in the differentiation of neural crest cells leads to NB development in the sympathetic nervous system or medullar region of the adrenal glands. NBs are highly heterogeneous tumors consisting of various sub-cell types, and they display various genomic alterations, including MYCN amplification, 10 different DNA ploidy patterns, 11 deletion in the short arm of chromosome 1, 12 a gain of chromosome 17q, 13 and chromosome 11q deletion. 14 These genomic alterations, along with International NB Staging System, stage histopathology and age are used to classify NB patients into risk groups. 15,16 The high-risk NB patients are challenging to treat and require high doses of chemotherapy and radiotherapy in the clinic. [17][18][19] The chances of tumor regression (>50%) and death rate (41.7% compared with 0% in low-risk and 5.2% in intermediate-risk patients) in high-risk patients are higher than other groups. 19,20 Moreover, the use of intense treatment regimens has significant adverse effects on patients' quality of life.
Several investigators have tried different therapeutic approaches like high-dose chemotherapy, surgery, radiation therapy, stem cell transplantation, retinoid therapy, and immunotherapy to improve clinical outcomes in aggressive high-risk NB patients. 17,19,21,22 Immunotherapy has recently gained traction in high-risk NB treatment and improved survival rates in patients. 22,23 Immunotherapy treatment harnesses the CD8 + T and natural killer (NK) cells to recognize, target, and eradicate malignant cells. [24][25][26] The immunotherapy approach of using chimeric monoclonal anti-disialoganglioside (GD 2 ) antibody dinutuximab has shown some success in NB clinical trials but is associated with severe side effects including neuropathic pain, infection, infusion-related reactions, capillary leak syndrome, decreased sensation, and paresthesia. 23,27 In NB, the levels of infiltrating CD8 + T and NK cells correlate with therapy response. 28,29 Elevated levels of programmed death-ligand 1 (PD-L1 or CD274) have been found in NB tumors. 28,30,31 Two recent reports suggested that higher PD-L1 expression positively correlates with worse outcomes in NB patients. 30,31 Patients with high PD-L1 expression had an increased risk of NB relapse and less overall survival compared with patients with low PD-L1 expression. [30][31][32][33] Therefore, therapies targeting PD-L1 could strengthen anti-tumor immunity and offer a promising strategy to treat highrisk NB.
Checkpoint inhibitor drugs that inhibit the PD-L1/PD-1 pathway have shown promising results in boosting immune responses and improving clinical outcomes in cancer patients. Food and Drug Administration (FDA) approval of PD-L1/PD-1 targeting drugs to treat cancer in recent years has attracted considerable interest in the discovery and development of improved PD-L1/PD-1 signaling pathway inhibitors. However, many patients do not respond well to these therapies and develop resistance and relapse. 30,31,34 Therefore, there is a need for novel approaches that target PD-L1 upstream signaling pathways to overcome these challenges.
Micro RNAs (miRs) are key upstream regulators in NB. [35][36][37] Several miRs have been reported to have roles in NB progression recently. [38][39][40] For example, suppression of MYCN targeting tumorsuppressor miRNA let-7 by LIN28B upregulates MYCN and promotes proliferation in non-malignant neuroblasts. 41 Amplification and overexpression of LIN28B or chromosomal loss of let-7 drives NB development and growth. 42 Moreover, downregulation of miR-34a, 43 miR-34b, 44 miR-124b, 45 miR-145, 46 miR-542, 47 miR-204, 48 and miR-27b, 49 is strongly associated with NB progression and therapy resistance. We recently showed that exosomal miR-155 induces chemotherapy resistance by targeting TERF1, an inhibitor of telomerase complex activation within the tumor microenvironment. 50 Several investigators have reported that a distinct set of miRs activate/inhibit the immune system at multiple levels, including regulation of the checkpoint molecule PD-L1. [50][51][52] For instance, miR-424(322) inhibits PD-L1 and CD80 expression, which promotes CD8 + T cells activation and reverses chemoresistance. The suppression of miR-21 in tumor-associated macrophages (TAMs) reprograms TAMs to proinflammatory anti-tumor type, improves CD8 + T cell tumor infiltration and inhibits angiogenesis. Our previous study shows that NB cells transfer miR-21 in exosomes to monocytes, which upregulates oncogenic miR-155 expression in monocytes. Monocytes transfer this miR-155 to NB cells via exosomes and induce chemotherapy resistance. 50,53,54 Several miR-based therapies, for example, tumor-suppressor miRs, miR-34 and miR-122, are in clinical development to treat cancer or other diseases. 55 Nervous system development during embryonic and pediatric stages requires precise epigenetic regulation of proliferation and differentiation pathways that are critically mediated by non-coding RNAs (ncRNAs). These ncRNAs, primarily miRs, also contribute to several tumor suppressive and oncogenic pathways within the tumor microenvironment. However, how miRs influence T and NK cell function in NB remains poorly understood. In this Figure 1. High-risk NB patients and PDX tumors show higher tumor-infiltrating lymphocyte exhaustion signatures together with reduced miR-15a and miR-15b (A) The Volcano plot depicting top-ranked up-(red) and down-regulated (blue) miRs (miR genes) in deceased versus living cohorts. The log-transformed TARGET dataset (n=249) was analyzed to prioritize miRs with statistically significant fold changes. The vertical red lines depict the fold change threshold of 1.77, corresponding to the log 2-fold change of 0.25. The horizontal red line indicates the P-value threshold of 0.03. The study prioritized miR-15a and miR-15b as tumor suppressor miRs for further characterization and validation. (B) The heatmap depicting the differential expression of miR-15a, miR-15b, and the T-and NK-cell signature genes across the TARGET dataset patient cohorts. The 249 TARGET dataset samples were classified into six categories with different aggressiveness based on a multi-attribute criterion depicted in the top annotation. The mean miR-15a and miR-15b log-transformed & scaled expression values for stromal-high and stromal-less tumor categories were -0.28 & -0.669. Two of the six categories were significantly enriched with stromal-high tumors. The values in stromal-high tumors were significantly lower than the corresponding mean values across the four other categories. This indicates potential antitumor immune response functions of miR-15a, and miR-15b. (C) Heat map of Pearson correlation Z-values between the mean miR-15a & miR-15b expression values and the eight genes of T-and NK-cell signatures. CD274 showed a strong signal with significant negative correlation profiles in stromal-low tumors versus the stromal-high tumors. This is consistent with the previous observations that miR-15a, miR-15b are tumor suppressors whose expression is higher in stromal-low tumors leading to better mR inhibitory functions captured as the significant negative correlations with target mRNA in these cohorts. (D) Heatmap depicting the differential expression of miR-15a, miR-15b, and the T-and NK-cell signature genes across the PDX dataset patient cohorts. The 24 PDX dataset samples were classified into six categories with different aggressiveness based on the Phase and Age depicted in the heatmap's top annotation. The mean miR-15a scaled expression values for the two less aggressive categories were high (0.5679 & 0.463), while its target mRNA CD274 demonstrated suppressed expression (-0.399 & -0.3206), corroborating our hypotheses that miR-15a/miR-15b-CD274 axis might contribute to immune evasion and tumor aggressiveness. (E) Kaplan-Meier curves showing eventfree and overall survival probability rates with different levels of PD-L1 in NB patients. (F) Immunofluorescence images of PD-L1 on NB patient tumor and normal tissue microarrays photographed at 5X. (G) Western blotting analysis of PD-L1 in whole lysates of PDX tumors of NB patients collected at the diagnosis and progression stages. Images were captured at low and higher exposures (exp). (H) IHC staining of PD-L1 in PDX tumors of NB patients at the diagnosis and progression stages. study, we found that infiltration of fewer T and NK cells, diminished miR-15a-5p, miR-15b-5p (hereafter miR-15a, miR-15b) levels, and enhanced PD-L1 expression were positively associated with poor survival in NB patients and patient-derived xenograft (PDX) tumors. The miR-15a and miR-15b activate CD8 + T and NK cell mediated anti-tumor immune response against NB by PD-L1 mRNA degradation through in vitro and in vivo approaches.

RESULTS
Higher PD-L1 and reduced miR-15a and miR-15b associate with higher tumor-infiltrating lymphocyte exhaustion signatures in NB patients and PDX tumors To prioritize tumor suppressor like miRs, we identified statistically significant downregulated miR (precursor miRs) genes using the survival status in the TARGET (Therapeutically Applicable Research to Generate Effective Treatments), NB patient (n=249) dataset as depicted in the Volcano plot Figure 1A. MiR-15a and miR-15b were selected to investigate further their tumor suppression mechanisms, especially their role in immune evasion. Immune-(NK cells, T lymphocytes), and stromal-[extracellular matrix (ECM), cancer-associated fibroblasts (CAFs), and mesenchymal stromal cells (MSCs)] cells are an area of current interest in the context of immune evasion mediated by tumor microenvironment in NB. The stromal cells comprising ECM, CAFs, and MSCs suppress effector immune cell activation and infiltration at the tumor site leading to immune evasion. Towards this, the 249 TARGET patient samples were classified into six categories using a multi-attribute criterion based on (1) overall survival (OS)-(2) event-free survival (EFS)-status (deceased/living), (3) INSS stage, (4) histology of tumor samples (favorable/unfavorable), (5) risk group (high/low), and (6) stroma cell percentages (low/high) as depicted in the heatmap top annotation ( Figure 1B). Furthermore, to investigate the regulatory accepts of miR-15a and miR-15b against immune effectors like T-& NK-cells, we performed differential-and co-expression analyses against the Tumor Inflammation Signatures (TIS) like T-cell exhaustion and T-& NK-cell abundance, which was presented as heatmaps ( Figure 1B-D). The TIS signature derived using 9000 tumor samples in TCGA comprised an 18-gene signature associated with four immune response processes 108,33 . TIS serves as an immune-phenotyping tool and can be a proxy for critical immune mediators. T cell/ NK cell abundance (NKG7, CD8A & HLA-E) and T cell exhaustion (CD274, CD276, TIGIT, PDCD1LG2 & LAG3), which constitute two of the four immune processes, were deemed particularly relevant for investigation of the mechanistic aspects of miR-15a and miR-15b mediated antitumor immune responses. Our preliminary observations from the TARGET dataset (n=249) differential-and co-expression analyses demonstrated significant downregulation of tumor suppressor miR-15a and miR-15b (mean scaled expression = -0.28 & -0.669) among the two stromal-high cohorts. The immune evasion functions were further delineated by assessing the correlations with the eight TIS T-cell and NK-cell function signature genes across the six tumor categories based their aggressiveness ( Figure 1B). Among the six clinical attributes used for patient stratification, the stromal cell percentage measure is the most relevant and conceptually related to immune evasion. It has demonstrated that tumors with lower stromal cell percentages (primary tumors) have potent antitumor immune responses. Higher stromal cell percentages (aggressive tumors) were associated with pro-tumor immune evasion phenotypes. Further, we determined the correlation between the mean expression values of miR-15a, miR-15b, and T-and NK-cell gene signatures in TARGET, NB patient cohorts. The expression of miR-15a and miR-15b is lower in stromal-high & aggressive tumors, as shown in Figure 1C. The coexpression analyses further corroborated these findings with stromallow & primary tumors demonstrating high miR-15a and miR-15b expression and high mRNA inhibitory activity observed as a significant negative correlation (mean correlation z-value: -1.943) with CD274, the TIS signature genes with the most significant differential co-expression and a clear signal. Stromal-high & aggressive tumors with little or no miR-15a, and miR-15b expression/activity demonstrated relatively insignificant negative correlations with CD274 (mean correlation z-value: -1.461) ( Figure 1C). We further performed a similar parallel analysis using PDX tumor data available from the dataset of genomic profiling of childhood tumor PDX models 109 . The 24 sample PDX dataset classified into six categories with varying aggressiveness using clinical attributes Phase and Age ( Figure 1D). The miR-15a demonstrated high expression for the two least aggressive categories (mean: 0.5679 & 0.463), while its most significant target, CD274, showed suppressed expressions (mean: -0.399 & -0.3206) in the same two least aggressive categories. These observations from the TARGET dataset further corroborated with the PDX dataset reinforce our hypotheses that miR-15a/miR-15b-CD274 axis might contribute to immune evasion and tumor aggressiveness.
From the TIS signature genes, we observed a remarkable, highly statistically significant correlation between increased PD-L1 expression and poor survival in the TARGET, NB patients (Figures 1E and S1A). We further screened the relationship between miR-15a and . SK-N-BE(2) cells were transfected with miR-15a, miR-15b, or ctrl miRs for 48 h followed by immunoprecipitation (IP) with a-Ago2 antibody. Ago2 bound RNA complexes were eluted, purified, and quantified for Ago2 bound mRNA and miRs by RT-qPCR using TaqMan assays. GAPDH mRNA and U6, a non-coding small nuclear RNA, were used as non-specific controls. The data were compared with the control IgG-bound mRNA or miR and set to one for normalization. Data represent the mean ± standard error of three to four independent biological experiments. Statistical analyses were performed using a two-sided unpaired t-test. *p < 0.001. survival status from other NB patient datasets. We observed that the expression of DLEU2, a host gene of miR-15a (GSE16476, n=88), is associated with relapse-free and overall survival in NB patients (GSE16476, n = 88), ( Figure S1B). Similarly, high miR-15b expression is significantly associated with better progression-free survival (p=0.048) and overall survival (p=0.056) in NB patients (Tumor NB ALT-WestermannÀ144-tpm-gencode19 in R2 database) ( Figure S1C). Further, as given in Figure 1D-G, higher levels of miR-15a found in low-risk vs. high-risk. Later, we performed a tissue microarray (TMA) analysis to examine the PD-L1 protein expression level in NB and normal tissues. As shown in Figures 1F and S1J, PD-L1 protein expressed much higher in NB tumor tissues than in nonneoplastic normal tissues. Moreover, we checked PD-L1 expression in NB PDX models (derived from patient tumor cells collected at the diagnosis and relapse stages) that we established by expansion in nude mice, as illustrated in Figure S1K. In western blotting analysis, we saw a trend of PD-L1 upregulation in relapse-specific PDX tumors versus those taken at diagnosis ( Figure 1G, lanes 3, 6, 7 versus 1, 2, 4, 5). Next, immunohistochemical (IHC) staining of PD-L1 showed similar results in PDX tumors ( Figures 1H and S1L). Furthermore, high levels of PD-L1 expression were detected on the surface of GD2+ tumor cells isolated from PDX tumor tissues from NB patients ( Figures S1M and S1N). These observations demonstrate that higher PD-L1 and diminished miR-15a and miR-15b are associated with adverse clinical features and poor survival in NB patients.

MiR-15a and miR-15b target PD-L1 in NB cells
To evaluate the role of miR-15a and miR-15b in PD-L1 regulation if any, we first transfected SK-N-AS, NB-19, and SK-N-BE(2) cell lines with precursors of miR-15a, miR-15b, and non-targeting control miRs for 48 h. The quantification of miRs is shown in Figure S2A. Cells showed a significant reduction of PD-L1 mRNA upon receiving miR mimics (Figure 2A). A similar decrease was also seen in PD-L1 total protein levels in miRs transfected cells ( Figures 2B and S2B), and the downregulation was more pronounced in SK-N-AS and NB-19 cells than SK-N-BE(2). Hence, we chose these two cell lines for further studies. To explore the mechanistic aspects of PD-L1 regulation by miRs, we checked if any complementary sequence exists between PD-L1 mRNA and miRs using in silico analysis tools including TargetScan. Figure S2C shows the presence of complementary binding sites between miR-15a, miR-15b, and 3 0 UTR of PD-L1 mRNA. To investigate whether these miRs regulate PD-L1 through 3 0 UTR, we cloned PD-L1-3 0 UTR (NM_014,143.4) wild-type region and seven base mutations at nucleotides 567-573 in the miR binding sites of the 3 0 UTR region referred to as 3 0 UTR mutant into a commercially available luciferase reporter vector pEZX-MT06 (Cat. #HmiT117860-MT06, Cat. #CS-HmiT117860-MT06-01; GeneCopoeia, Rockville, MD). SK-N-AS cells were transfected with precursors of miR-15a, miR-15b, or control miR oligonucleotides together with PD-L1-3 0 UTR wild-type or mutatant constructs for 48h. MiR-15a and miR-15b showed a reduction in luciferase activity with PD-L1-3 0 UTR wild type but not with PD-L1 mutant ( Figure 2C). These observations indicate that miR-15a and miR-15b target PD-L1 through 3 0 UTR (Figure 2A-C). Further, we evaluated whether surface PD-L1 expression is also regulated by these miRs. We analysed surface PD-L1 expression in miR-15a and miR-15b transfected SK-N-AS, NB-19 cells, and an SK-N-B(E)2 cell line stably expressing miR-15a by flow cytometry. As shown in Figures 2D-2F, reduced surface PD-L1 expression was found, suggesting that miR-15a and miR-15b target PD-L1 in NB. The representative flow cytometric plots showing fluorescence minus one (FMO) control to omit background signals are shown in Figure S2D. The quantification of miR-15a in stable miR-15a expressing SK-N-B(E)2 cell line is shown in Figure S2E. To further confirm that the effect is not off-target, we silenced endogenous miR-15a and miR-15b in NB cells by transfecting antisense oligonucleotides (miR inhibitors). Conversely, with the introduction of miR inhibitors, cells exhibited increased PD-L1 expression ( Figures 2G and S2F). The quantification of miR-15a and miR-15b in NB cells after transfection with miR inhibitors is shown in Figure S2G. Next, we performed Ago2 immunoprecipitation to see whether these miRs bind to PD-L1 mRNA in NB. A schematic presentation of Ago2-RNA-IP is depicted in Figure S2H. As shown in Figure 2H, anti-Ago2 antibodies successfully immunoprecipitated PD-L1 mRNA but not GAPDH mRNA. Further, Ago2 also immunoprecipitated miR-15a and miR-15b but not U6, a small nucleolar RNA, suggesting that miR-15a and miR-15b regulate PD-L1 via direct interaction with PD-L1 mRNA in NB cells.

MiR-15a and miR-15b promote CD8 + T cell proliferation and induce CD8 + T cell-mediated NB cytotoxicity
We tested the effect of miR-15a and miR-15b on CD8 + T cell-mediated cytotoxicity against NB cells. NB cells were treated with miR-15a and miR-15b mimics for 24 h and cocultured with activated CD8 + T cells for an additional 24 h. As shown in Figures 3A and 3B, a higher Untouched CD8 + T cells were isolated from PBMCs of healthy human blood donors by negative selection using the MojoSort human CD8 + T cell Isolation Kit (BioLegend). For 7 days, CD8 + T cells were activated and expanded using human T-activator CD3/CD28 Dynabeads and used for coculture experiments with NB cells. NB cells were transfected with miRs for 24 h and used for coculture experiments with activated CD8 + T cells for an additional 24 h. Non-targeting control (ctrl) miR cells served as a control group. Cells were permeabilized (intracellular), fixed, and stained for respective antibodies described in the materials and methods section and analyzed by flow cytometry. The percentage of Granzyme B + (A), Perforin + (B), CD3/CD8 + (C), ki67 + (D), and cleaved caspase-3 + cells (E) are shown in each of their respective plots. Data represent the mean ± standard error of three to five independent biological experiments. Statistical analyses were performed using a two-sided unpaired t-test.
percentage of CD8 + T cells showed intracellular granzyme B and perforin expression after coculture (E:T=1:1) with miR-15a and miR-15b expressing SK-N-AS ( Figures 3A, 3B top panels, and S3G), and NB-19 ( Figures 3A, 3B, bottom panels, and S3H) cells. Next, we checked proliferation capacity and Ki-67 (a nuclear protein associated with cell proliferation) staining in cytotoxic CD8 + T cells after coculture (E:T=1:1) with miR-15a and miR-15b expressing NB cells. Interestingly, CD3 + CD8 + T cells showed more proliferation (61% in miR-15a and miR-15b versus 33% control miRs, Figures 3C, S3I1, and S3I2) and higher Ki-67 expression ( Figures 3D and S3J) in the presence of miR expressing SK-N-AS ( Figures 3C and 3D, top panels), and NB-19 ( Figures 3C and 3D, bottom panels) cells than control group cells. These results showed that NB cells expressing miR-15a and miR-15b induce proliferation in CD8 + T cells. In the same experiments, we further checked the functional effect of activated CD8 + T cells on NB cytotoxicity by measuring cleaved caspase-3, an apoptosis marker in target NB cells. As shown in Figures 3E and S3K, caspase-3 cleavage increased from 6% to 16.2% and 22.8% (control miR versus miR-15a and miR-15b) in SK-N-AS cells ( Figure 3E, top panel), and 19% to 24.7% and 25.5% (control versus miR-15a and miR-15b) in NB-19 cells ( Figure 3E, bottom panel). These results suggest that miRs expressing NB cells were more prone to CD8 + T cellmediated apoptosis. Further, to confirm these findings, we cocultured activated CD8 + T cells (E:T=5:1) with stable luciferase gene expressing miR-15a and miR15b treated SK-N-BE(2) cells for 24 h and performed the luciferase reporter assay as described in the materials and methods section. We found significantly reduced luciferase activity in miR expressing cells than in their control group, suggesting that miR-15a and miR-15b induce CD8 + T-cell-mediated NB cytotoxicity ( Figure 3F).

MiR-15a and miR-15b promote NK cell activation and induce NK cell-mediated NB cytotoxicity
We tested the effect of miR-15a and miR-15b on NK cell-mediated anti-tumor immune response against NB cells. Representative flow cytometric pseudo color plots showing human NK cells purity before and after expansion ex vivo using irradiated K562-mbIL21 feeder cells are depicted in Figures S4A and S4B. The FMOs related to the experiments are shown in Figures S4C-S4E. First, NB cells were treated with miR-15a and miR-15b mimics for 48 h and cocultured with activated NK cells for an additional 5 h. NK cells displayed remarkably higher perforin expression upon coculture with miR expressing NB cells than control-miRs ( Figures 4A, 4B, S4F, and S4G). Further, fluores-cence intensity analysis showed increased accumulation of CD107a molecules on the NK cells placed on NB cells expressing miRs (17% to27.2% and 42.5%, SK-N-AS; 21.8% to31.3% and 37.7%, NB-19) (Figures 4C and S4H). These results show that miR-15a and miR-15b in NB cells might enhance the activation of NK cells.
To understand the functional consequences of NK cell activation on tumor cell cytotoxicity, we assayed for cleaved caspase-3, a marker for cell death, and the number of dead cells in miR-15a and miR-15b expressing SK-N-AS and NB-19 cells after coculture (E:T=1:1 for SK-N-AS and E:T=0.25:1 for NB-19) with or without NK cells for 5 h. NK cells significantly induced caspase-3 cleavage from 27.1% to 46.9% and 44.6% in SK-N-AS and 49% to 62% and 61% in NB-19 cells in the presence of miRs compared with control-miRs ( Figures 4D and  4E). Quantification of cleaved caspase-3 is shown in Figures S4I  and S4J. Next, we detected dead NB cells using Zombie Aqua, an amine-reactive, a water-soluble fluorescent dye that stains dead cells due to their compromised membranes. As shown in Figure 4F, NK cells coculture with SK-N-AS cells induced a dose-dependent increase of dead NB cells from 3.2% to 25.4% and 47.1% at an E&T ratio of 0.5: 1 and 1:1. The overexpression of miR-15a and miR-15b further sensitized SK-N-AS cells to NK cell-mediated cytotoxicity. A similar trend of NK cell-mediated cytotoxicity was further noticed in miR expressing NB-19 (E:T=0.25:1 and 0.5:1) and SK-N-B(E)2 (E:T=1:1) cells ( Figures 4G and 4H). Quantified NK cell-mediated NB cytotoxicity levels are given in Figures S4K-S4M. These results were further confirmed in phase-contrast images of miR-15a stably expressing SK-N-B(E)2 cells showing higher frequency of dead cells upon coculture with activated NK cells ( Figures 4I, S4N, and S4O). These results collectively show that miR-15a and miR-15b promote NK cell-mediated release of cytotoxic granule perforin and cleaved caspase-3 activation, enhancing the efficacy of NK cell cytotoxicity against NB cells.

PD-L1 blockade rescues miR-15a and miR-15b mediated immune activation and cytokine secretion
We determined to test the role of PD-L1 in miR-15a and miR-15b mediated regulation of CD8 + T and NK cell functions in NB. First, SK-N-AS cells were blocked with anti-PD-L1 or IgG control antibody for 24 h. followed by transfection with miR-15a and miR-15b mimics for an additional 24 h. As shown in Figure 5A, miR-15a and miR-15b effectively silenced (47% to 26% and 28.4%) surface PD-L1 in cells  treated with control antibody but not PD-L1 antibody. These results confirmed an effective blockage of surface PD-L1 in SK-N-AS cells. These cells were cocultured (E:T=1:1) with CD8 + T cells for 24 h, NK cells for 5 h and analyzed for activation markers granzyme B on CD8 + T cells and CD107a on NK cells. As shown in Figures 5B, 5C, S5A, and S5B, blocking cell surface PD-L1 significantly enhanced granzyme B secretion from 40% to 54% in CD8 + T cells and degranulation marker CD107a from 34.2% to 63.7% in NK cells. Consistent with previous results, SK-N-AS cells treated with miR-15a and miR-15b significantly increased granzyme B secretion (40% to46.7% and 49.3%, p=0.0034 and p=0.0008) in cocultured CD8 + T cells. Further, we noticed a similar trend with CD107a (34.2% to 49.7% and 56%, p=0.0025 and p=0.0006) in cocultured NK cells in the control IgG antibody treatment group, but no significant difference in PD-L1 antibody treated group.
Next, to confirm the role of PD-L1 in miR-15a and miR-15b mediated anti-tumor response, we followed another approach of silencing PD-L1 using small interfering RNAs (siRNA). As shown in Figure 5D, PD-L1 siRNAs effectively silenced (47.3% to31.5%, p<0.05) surface PD-L1 in NB-19 cells compared with control siRNAs. These cells were used for coculture experiments with NK cells. NK cells showed a significant increase of CD107a expression upon coculture with miR-15a and miR-15b expressing NB-19 cells that are treated with control siRNAs but no significant difference in cells treated with PD-L1 siRNA ( Figures 5E and S5C).These results further confirmed that PD-L1 suppression plays an important role in miR-15a and miR-15b-mediated immune cell activation. In addition, we analyzed the secretion of soluble factor interferon (IFN)-g by ELISA in cell culture supernatants of effector NK cells coincubated (E:T=1:1) with miR-15a and miR-15b expressing target SK-N-AS cells and treated with control IgG or anti-PD-L1 antibody. Interestingly, we did not see a significant difference in INF-g secretion in control miR versus miR-15a, miR-15b treated SK-N-AS cells and cocultured with NK cells. However, anti-PD-L1 antibody treatment increased INF-g concentration in the media upon coculture with NK cells, which was further decreasing after miR-15a and miR-15b overexpression, implying that these miRs suppressed INF-g secretion from NK cells ( Figure 5F). Altogether, these results indicate that suppression of PD-L1 by miR-15a and miR-15b mediate immune cell activation and cytokine secretion response in NB.

MiR-15a regulates CD8 + T and NK cell-mediated anti-tumor immune response in NB in vivo
Based on these in vitro findings, we next sought to investigate the role of miR-15a and miR-15b in modulating immune response using an im-mune-competent C57BL/6 mice model. We chose miR-15a for our in vivo experiments as both miR-15a and miR-15b share a common complementary sequence with PD-L1 mRNA ( Figure S2C). We subcutaneously injected C57/BL6 mice with murine stable NB975 cells expressing either GFP-control miR or GFP-miR-15a for 30 days and evaluated tumor growth. A significant (p<0.05) repression of tumor growth, weight, and volume were observed in the mice that received miR-15a cells compared with control miR cells (Figures 6A, 6B, and S6A). Tumor tissues were harvested, freshly prepared single-cell suspensions, and analyzed PD-L1 on GFP+ve tumor cells by flow cytometry. As shown in Figures 6C and 6D, tumor cells of the mice that received miR-15a cells showed a significant (p<0.05) reduction of surface PD-L1 (64.9%-33.4%) compared to control miR cells, suggesting that miR-15a targets PD-L1 in vivo. Next, CD8 + T and NK cells infiltration was evaluated using a flow cytometer. As depicted in Figures 6E,  6F, S6B, and S6C, miR-15a tumors showed a significant increase in CD8 + T cell infiltration (30.8% versus 45.8%, p=0.002), and Ki-67, a CD8 + T cell proliferation marker from 6.73% to 9.48% (p=0.018) when compared with the control miR group. Further, the infiltrated CD8 + T cells in miR-15a-specific tumors produced higher levels of INF-g (0.72% versus 1.72%, p=0.038) (Figures 6G and S6D) than the control miR group. These results suggest that miR-15a promoted T cell priming and infiltration into NB tumors.
Tumor-infiltrating CD25 + CD4 + regulatory T cells (Tregs) with Foxp3 expression play a key role in suppressing anti-tumor activity and promoting tumor growth. 56 Thus, we examined the effect of miR-15a on CD25 + CD4 + Foxp3 + Tregs in NB tumors. As shown in Figures S6E and S6F, tumors of mice that received miR-15a cells had a significantly lower (p=0.0016) Treg population than the control group, indicating that miR-15a repressed Tregs in the tumor microenvironment to prevent immune escape. We next evaluated the activation status of tumor-reactive NK cells in tumors by analyzing perforin on CD3 À (to exclude T cells), NK1.1 + (a marker for mouse NK cells) gated population from the single-cell suspension of tumor tissues. We found higher perforin expression (16.5% versus 27.1%, p=0.013) in the miR-15 group tumors than the control miR group tumors ( Figures 6H and S6G). Last, we examined microvessel infiltration by IHC staining using an anti-mouse CD34 antibody, which stains endothelial cells. As shown in Figures 6I and S6H, tumors of mice that received control miR cells showed intense neovascularization as indicated by high CD34 staining, which was significantly decreased in tumors of mice that received miR-15a cells . The mean microvessel count/field in control miR group tumors was 233±44 (SD) and 101±52 (SD) in miR-15a group tumors ( Figure S6I). These observations indicate that miR-15a decreases tumor microvessel density in vivo. Overall, our data strongly support the role of miR-15a and miR-15b in anti-tumor immune response through PD-L1 degradation.

DISCUSSION
Tumor-immune interactions are intricate, and efforts to capture this complexity via a single marker such as PD-L1 expression as measured by IHC, or tumor mutation load as a surrogate of potential tumor antigenicity, yields limited and incomplete information. This is a significant concern for the analyses of pediatric tumors' immunogenicity since their tumor mutation burden is negligible. The complex and dynamic nature of the tumor-immune microenvironment also involves stromal cells consisting of extracellular matrix, cancer-associated fibroblasts, and mesenchymal stromal cells. These cells prevent the infiltration of immune mediators while enabling drug resistance and metastasis via mesenchymal cancer cell permeation. This data-driven study represents one of the earliest attempts to use stromal cell percentages to interrogate non-coding RNA molecules with tumor antigenic function. Using differential-and co-expression studies on TARGET (n=249) and PDX (n=24) datasets, the stromal-cell specific expression and co-expression of miR-15a, miR-15b, and its target gene PD-L1 were established. These findings were experimentally validated in great detail using in vitro and in vivo approaches. Various cell surface and other immunophenotyping markers were quantified, reinforcing the crucial role of the miR-15a/miR-15b & PD-L1 axis in regulating various immune mediators involved in NB tumor immunogenicity.
PD-L1 surface expression on cancer cells is primarily involved in immune cell dysfunctions, and many studies have reported its prognostic significance in several cancer types. 57-60 PD-L1 expression in tumor cells attenuate T cell activation and functions. 60,61 Currently, the FDA has approved the use of three PD-L1 inhibitors: (1) atezolizumab for non-small cell lung cancer (NSCLC), 62 (2) durvalumab for NSCLC and urothelial carcinoma (UC), 63 and (3) avelumab for UC and renal cell carcinoma. 64,65 Moreover, several other inhibitors are in different phases of clinical trials for various cancer types. 66 However, one of the challenges limiting the effectiveness of anti-PD-L1 therapy is that checkpoint molecules acquire several mutations during tumor progression, 67 suggesting the importance of targeting upstream regulators of checkpoint molecules. 68,69 MiRs have been shown to regulate various oncogenes in NB. 70,71 Hence miR-mediated targeting of PD-L1 could have great translational potential to stage 4 NB patients. Hence, we have developed workflows for elucidating regulatory networks of immune evasion to prioritize novel miRs upstream of PD-L1 regulation. The present study explored a miR-based approach to target immune checkpoint molecule PD-L1, which sensitized NB cells to CD8 + T and NK cellmediated cytotoxicity. Multiple studies have demonstrated the association of poor survival with fewer tumor-infiltrating lymphocytes, including T cells and NK cells in NB. 28,33,[72][73][74][75] Moreover, NB tumors employ several immune evasion strategies, including the expression of immune checkpoint molecules. [76][77][78][79][80] Our bioinformatic co-expression analysis from the TARGET data highlighted the presence of exhaustive immune signatures, including PD-L1, in NB patients and PDX tumors. We further found that lower expression of PD-L1 is associated with better progression-free and overall survival in NB patients. Moreover, PD-L1 mRNA levels are higher in patients with tumor relapse. In contrast, high expression levels of miR-15a were positively correlated with good progression-free and overall survival of NB patients. Patients who had a high risk or died from the disease expressed lower miR-15a levels. Besides, miR-15a expressed differentially across different NB stages, with the lowest expression found in most aggressive and therapy-resistant stage 4 patients. Further, we found an inverse correlation between PD-L1 and miR-15a expression in NB patients. These analyses implied that higher PD-L1 and lower miR-15a expressions are associated with poor clinical features in NB.
Our experimental results showed that miR-15 family in particular miR-15a and miR-15b degraded PD-L1 mRNA and protein in several NB cell lines, including SK-N-AS, NB-19, and SK-N-BE(2). Further, enrichment of miR-15a, miR-15b, and PD-L1 mRNA with Ago2, a core component of RNA-induced silencing complex (RISC), revealed a specific interaction between PD-L1, miR-15a, miR-15b within the RISC complex. RISC plays a major role in regulating gene expression through Argonaute RISC Catalytic Component2 (Ago2)-mediated interaction with miRs and the complementary sequence of the target mRNA, thus serving as an indicator of target gene regulation. 81 PD-L1, a membrane-anchored protein, binds to its receptor protein PD-1 on the surface of immune cells. Thus, surface expression of PD-L1 on cancer cells is a critical component of determining the anti-tumor immune response. CD8 + T cells are the main components of an adaptive immune response against tumor cells. 82,83 NK cells are the main effector immune cells that play essential roles in tumor surveillance and demonstrate anti-tumor cell cytotoxicity without the requirement of the major histocompatibility complex. 84 Both CD8 + T and NK cells express PD-1 receptors, making them more susceptible to PD-L1-mediated dys-function by tumor cells. 85,86 MiR-15a treatment significantly reduced cell surface PD-L1 in NB cell lines suggesting that miR-15a could promote CD8 + T and NK cell-mediated anti-tumor immune responses in NB cells. 52 NB cells upon receiving miR-15a and miR-15b induced CD8 + T and NK cell activation markers such as granzyme B, perforin, and CD107A, as well as CD8 + T cell clonal expansion. CD107A, also known as lysosomal-associated membrane protein-1 (LAMP-1), is a marker for NK cell activation and protects NK cells from degranulation-associated damage. 87,88 Cytotoxicity experimental results revealed that NK cells were more efficient in lysing NB cells upon miRs treatment. Enhanced granzyme B and perforin secretions and activation of CD8 + T and NK cells by miR-15a and miR-15b are associated with a more prominent cleavage of caspase-3 in target NB cells. Caspase-3 is one of the endpoint executioner caspases, which cleaves various target proteins, including PARP-1, leading to an apoptotic breakdown of the cell, 89 suggesting that miR-15a and miR-15b play an essential role in executing immune cell-mediated NB cell apoptosis.
Further, mechanistic experiments interfering PD-L1 by siRNAs or blocking the surface PD-L1 using anti-PD-L1 antibody revealed efficient CD8 + T and NK cell-mediated anti-tumor immune responses against NB cells. Studies have also shown that tumor necrosis factor (TNF)-a is critical for NK cells in tumor response as it promotes cell recruitment at the tumor site. 90 TNF-a binds to TNF receptor-2 (TNFR2) on the surface of NK cells and promotes INF-g production in NK cells. NK cells control tumor growth by releasing TNF-a and IFN-g. 90-92 IFN-g and TNF-a synergistically enhance NK cell-induced tumor cytotoxicity through the NF-kB-dependent pathway. 91 However, we did not see a significant difference in the induction of IFN-g in NK cells upon treatment with miRs. Interestingly, NB cells treatment with miR-15a, miR-15b in the presence of anti-PD-L1 antibody significantly decreased IFN-g secretion from cocultured NK cells. MiR-15a, miR-15b have complementary binding sequences to 3 0 UTR of IFN-g mRNA, and studies have shown that miR-15a, miR-15b inhibit IFN-g production in NK cells. 93 Since miR-15a, miR-15b are not good for NB cells, we believe NB cells can secrete out these miRs in the medium. The uptake of miR-15a, miR-15b by NK cells could be responsible for decreased IFN-g production in NK cells after coculture with miRs overexpressing NB cells. IFN-g is known to induce surface PD-L1 expression on cancer cells, protecting cancer cells from anti-tumor immunity and playing a role in chemoresistance. [94][95][96][97] MiR-15 mediated inhibition of IFN-g secretion might act as a double-edged sword that prevents IFN-g-induced PD-L1 upregulation and promotes CD8 + T and NK cell functions by directly targeting PD-L1. In vivo experiments provided evidence that miR-15a promotes the activation and infiltration of CD8 + T and NK cells and reduces tumor vasculature, tumor size, and tumor-promoting CD4 + CD25 + FoxP3 + Tregs in tumors. These results further suggest that miR-15a mediates immune regulation and, therefore, should be considered in treating patients together with approved drugs .
In summary, we have provided evidence that diminished miR-15a, miR-15b, and higher PD-L1 expression are associated with immune cell exhaustion signatures and tumor relapse in NB patients and PDX tumors. PD-L1 mediates immune evasion, whereas targeting PD-L1 using miR-15a and miR-15b oligonucleotides leads to immune cell activation, cytokine secretion, and immune cell-mediated NB cytotoxicity and anti-tumor immune response in vitro. Injection of miR-15a expressing cells in immune-competent C57/BL6 mouse reduced tumor growth and tumor vessel formation while increasing the activation and infiltration of immune cells, suggesting that degradation of PD-L1 by miR-15a and miR-15b activate anti-tumor immune responses in high-risk, low immunogenic NB patients.

Patient stratification using clinical attributes and tumor-stromal cell percentages
A systematic meta-analysis of tumor stroma quantified in 16 solid cancer types from 2732 patients revealed that the stroma fraction was highly variable within and across the tumor types. Kidney cancer showed the lowest while pancreato-biliary type periampullary cancer www.moleculartherapy.org demonstrated the highest stroma proportion (median 19% and 73% respectively) 110 . Despite the wide variation, no significant association between stromal abundance and prognosis was observed. Given this context, we have chosen a stromal cell percent cut-off that yielded the approximately 50-50 stratification of the 249 patients in concert with survival status.

Missing value imputation using multivariate imputation by chained equations
Multivariate imputation by chained equations (MICE) is an imputation method based on Fully Conditional Specification, where different models impute incomplete attributes. Hence, MICE can impute missing values in datasets with continuous, binary, and categorical attributes by using either logistic or linear regression for each attribute. While expression/count data mostly had continuous data, the tumor sample attribute file had continuous and categorical variables. The integrity of the variable types was verified before imputation was attempted. The number of multiple imputations (m) and the number of iterations (maxit) varied from their default value of 5. The best parametric settings were determined after preliminary explorations and assessments of the imputation quality. Both TARGET and PDX datasets had multiple missing values in the expression and tumor sample attribute files.

Representative probe prioritization
Both TARGET and PDX datasets had multiple probes per gene. Furthermore, both these datasets were derived from different array platforms. Any meaningful meta-analyses and corroboration assessments require a standardized approach for prioritizing the best representative probe of the gene. Prioritization of the best representative probe prior to the application of statistical tests can, in some instances, be detrimental to the consistency between the obtained results and the proposed hypotheses. However, such prioritization represents a best practices protocol with an unbiased, objective, and data-driven approach with robust statistical significance and predictive utility.

Human tissue samples and peripheral blood mononuclear cells
Deidentified formalin-fixed and embedded paraffin blocks of NB patient tissue samples were obtained as a tissue microarray (#NB642a) from US Biomax Inc (26 patient samples, in duplicates). All procedures with NB patient tissue microarray followed our institutional ethical standards and were approved by the human ethics committee. These blocks were used for the Immunofluorescence (IF) analysis. Freshly obtained human peripheral blood mononuclear cells (PBMCs) from the Elutriation Core Facility of our institute were used for the isolation of CD8 + T and NK cells.

PDX tumors
We have established NB-specific PDX tumors in collaboration with C. Patrick Reynolds, MD, PhD, from Texas Tech University Health Sciences Center through the Alex's Lemonade Stand Foundation (ALSF)/Children's Oncology Group (COG) Childhood Cancer Repository (https://www.cccells.org) and the Childhood Solid Tumor Network at St. Jude Children's Research Hospital. 70 In brief, patient tumor samples were surgically removed, cut into multiple pieces in sterile Hanks' balanced salt solution supplemented with antibiotics. The deidentified patient tumor cells were subcutaneously injected into 4-to 6-week-old nude mice. Mice were observed for palpable tumors, and the tumor tissues were harvested once the tumor reached end volume, then a single-cell suspension was prepared for re-implantation in the next set of mice, with multiple passages for expansion. These tumors were used for the western blotting analysis, P0 = Passage 0, P1 = Passage 1 ( Figure S1K). Our institutional ethics committee approved the study methodologies. The PDX tumor tissues COG-N-564x, COG-N-453x, COG-N-496x, COG-N-603x, COG-N-470x, COG-N-424x, and SJNBL013763_x1 were used for the Western blotting and IHC analysis.  (2), NB-19, K562, and NB975 cells were cultured in RPMI and SK-N-AS cells in DMEM media, each containing 10% heat-inactivated fetal bovine serum (FBS) supplemented with penicillin/streptomycin (50 U/mL), L-glutamine, sodium pyruvate, and non-essential amino acids, as described earlier. Cell lines were cultured in a 5% CO 2 humidified atmosphere at 37 C. Cell lines were authenticated by short tandem repeats analysis and were frequently tested to ensure they were free from mycoplasma.

CD8 + T cell isolation, activation, and culture
PBMCs from healthy human donors were processed using Histopaque-1077 reagent and CD8 + T cells were isolated using MojoSort Human CD8 + T cell isolation kit (Biolegend, #480012) according to the manufacturer's instructions. The CD8 + T cells were stimulated by Dynabeads Human T-Activator CD3/CD28 (Thermo Fisher Scientific, #11131D) at a 1:1 ratio of T cells to beads in the growth medium of 10% FBS in RPMI-1640 for 48 h. The activated CD8 + T cells were harvested, and the purity was analyzed via flow cytometry using anti-CD3, anti-CD8, and anti-CD4 antibodies. The CD8 + T cells used in the study were CD3 + CD8 + CD4 À . The gating strategy is given in Figure S3A

NK cell isolation, expansion, and culture
Untouched NK cells were isolated from PBMCs using the MojoSort Human NK Cell isolation kit (BioLegend, #480054) per the manufacturer's instructions. In brief, PBMCs (1Â10 7 /100 mL MojoSort buffer) were incubated with Biotin-Antibody Cocktail (10 mL) against lineage-specific targets for 15 min on ice. The cells were washed with MojoSort buffer, followed by the addition of streptavidin magnetic nanobeads (10 mL) for 15 min on ice, and resuspended in MojoSort buffer (2.5 mL). The cells were placed in the magnet for 5 min at room temperature (RT), and the unbound NK cells were separated by pouring the liquid into a fresh tube. The cells were cultured in RPMI-1640 supplemented with 10% FBS and 1% penicillin/streptomycin together with recombinant (r) human (h) IL-2 (20 ng/mL), and rhIL-15 (50 ng/mL) for 4 days and were used for expansion. The obtained NK cells were cultured with irradiated K562 CSTX002 mbIL21.41bbL cells at a ratio of 1:1 in the presence of rhIL-2 (20 ng/mL) and expanded for 2 weeks. NK cells were replenished with fresh K562 cells once a week and rhIL-2 twice a week. The purity of NK cells was determined by staining with anti-CD56-APC/ Fire750, and CD16-PE antibodies, whereas contaminating T cells were labeled by anti-CD3-AF700 antibodies, using flow cytometry. The NK cells used in the study were CD56 + CD16 + CD3 À (Figures S4A and  S4B) and were used for coculture experiments with NB cells or frozen in liquid nitrogen in a freezing medium containing 95% FBS for future use. These activated NK cells were used in coculture experiments with SK-N-AS cells (E:T=1:1), NB-19 cells (E:T=0.25:1), and SK-N-B(E)2 cells (E:T=1:1) with stable miR-15a expression for 5 h.

PD-L1 cell surface blocking
Surface PD-L1 on NB cells was blocked by treatment with control immunoglobulin (Ig) G or anti-PD-L1 antibody (15 mg/mL, BioXCell, Lebanon, NH) for 24 h followed by transfection with miR-15a and miR-15b oligonucleotides for an additional 24 h. After wash out the unbound antibody, these cells were used for coculture experiments with CD8 + T and NK cells.

Western blotting and antibodies
Cells or tumor tissue samples were lysed in radioimmunoprecipitation assay buffer containing freshly added EDTA-free protease and phosphatase inhibitor tablets for 30 min on ice. The protein samples were resolved on 10% gradient SDS-PAGE and transferred to a nitrocellulose membrane, followed by blocking with 5% BSA for 1 h at RT. The membranes were incubated with primary antibodies, including monoclonal mouse anti-human PD-L1/CD274 (Proteintech, #66248-1-Ig, 1:8,000, overnight at 4 C), horseradish peroxidase (HRP)-conjugated anti-b-actin (Proteintech, #HRP60008, 1:12000, for 2 h at RT). After secondary antibody incubation at RT for 1 h, membranes were developed using an enhanced chemiluminescent reagent. The images were captured on X-ray film and quantified using ImageJ analysis software. 104

Flow cytometry analysis
The assessment of PD-L1, CD8 + T, and NK cell surface markers, cleaved caspase-3, perforin, granzyme B, Ki-67, Foxp3, INF-g and cell survival analysis was performed by immunofluorescence staining followed by flow cytometric analysis. In brief, cells were collected at 350 Â g for 5 min, washed with, and suspended in PBS (1.0 Â 10 5 / 100 mL), and all steps were performed in the dark unless stated. The dead cells were excluded by incubation with Zombie Aqua fixable viability dye (BioLegend, # 423102) at RT for 30 min. Before immunostaining, cells were incubated with Fc receptor blocking solution (TruStain FcX for human, #422302 or TruStain FcX PLUS for mouse #156604, BioLegend) diluted in ice-cold cell staining buffer (BD Biosciences) on ice for 15 min to reduce non-specific antibody binding. For the analysis of surface GD 2 , PD-L1, CD3, CD4, CD8, CD16, CD25, CD56, CD107a, and NK1.1 proteins, the cells were incubated with fluorochrome-conjugated anti-GD 2 , anti-PD-L1, anti-CD3, anti-CD4, anti-CD8, anti-CD16, anti-CD25, anti-CD56, anti-CD107a, and anti-NK1.1 antibodies on ice for 20 min. For measurement of intracellular cleaved caspase-3, perforin, granzyme B, Ki-67, Foxp3, and INF-g, cells were washed with staining buffer (2x), fixed and permeabilized with Cyto-Fast Fix/Perm Buffer (BioLegend) at RT for 20 min followed by staining with the fluorochrome-conjugated anti-perforin, anti-granzyme B, anti-Ki67, anti-Foxp3, anti-INF-g, or unconjugated anti-cleaved caspase-3 antibodies at RT for 25 min. The cells were further incubated with the fluorochrome-conjugated secondary antibody at RT for another 40 min if stained with an unconjugated primary antibody. After incubation, cells were washed with Cyto-Fast Perm buffer (2x) and cell staining buffer (1x), then suspended in cell staining buffer. The specificity of the antibodies was verified by staining with respective isotype control antibodies and FMOs. The samples were analyzed freshly or within 24 h using the Attune Acoustic Focusing Flow Cytometer (Thermo Fischer Scientific) or LSRII (Becton Dickinson). The data evaluation was performed using the FlowJo v10.8.1 program.

Enzyme-linked immunosorbent assay
The concentration of INF-g in the culture supernatants was assayed by Enzyme-Linked ImmunoSorbent Assay (ELISA) using LEGEND MAX Human IFN-g ELISA (BioLegend, #430107) kit per the manufacturer's instructions. In brief, the cell culture supernatants were collected and human INF-g standards at 500, 250, 125, 62.5, 31.3 and 15.6 pg/mL concentrations were freshly prepared. The ELISA plates were washed (4x) with PBS containing 0.05% Tween 20, and blocked (200 mL/well) with PBS containing 10% FBS at RT on a shaker for 1 h. The cell culture supernatant samples (3x diluted), or human IFN-g standards (100 mL/well) were added at RT for 2 h on a shaker. The plates were washed (4x), then the human INF-g antibody was added (100 mL/well) and incubated for 1 h at RT with shaking. The plates were washed (5x), and the avidin-HRP solution was added (100 mL/well) and incubated at RT for 30 min with shaking. The plates were washed (5x) before adding tetramethylbenzidine substrate solution (100 mL/well) at RT for 5 min in the dark or until the development of blue color. The reaction was stopped by a stop solution (100 mL/well), and the optical density was measured at 450 nm and 570 nm on a plate reader. The final readings were obtained by subtracting 570 nm absorbance from 450 nm and compared with IFN-g standard curve. Each sample was assayed in triplicate and repeated more than three times.

Direct tumor cell killing quantification by luminescence assay
Luciferase gene expressing stable SK-N-B(E)2 cells were used in CD8 + T cell-mediated direct tumor cell killing assay. The cells were seeded (0.9 Â 10 5 /well) into a 12-well plate for 48 h, followed by transfection with miR-15a, miR-15b, or non-targeting control mimics for 9 h. Next, activated CD8 + T cells (effector [E]) were cocultured with SK-N-B(E)2 (target [T]) cells at a 5:1 E and T ratio for 48 h. The supernatants containing the effector CD8 + T cells were removed, and the remaining dead detached target cancer cells were carefully washed away by rinsing twice with PBS (1 mL). The target cancer cells were collected, washed with PBS, and lysed (100 mL) in reporter lysis buffer (Promega Corporation, Madison, WI; #E3971) with a freeze-thaw cycle to achieve complete cell lysis. The cleared cell lysates (80 mL) were transferred to a Corning/Costar white, flat-bottom 96-well plates, followed by the addition of substrate (80 mL/well) D-luciferin (Promega Corporation, #E151A). The luminescence was measured using a luminescence plate reader (Infinite 200 PRO, Tecan). The luminescence units were normalized to control (tumor target cells without any effector CD8 + T cells) set to 100. The results were reported as the percentage of relative luminescence units.

PD-L1 3 0 UTR constructs and luciferase reporter assays
Publicly available online bioinformatics databases such as TargetScan (www.targetscan.org), miRNA Target Prediction Database (mirdb. org), and microRNA target prediction (www.microrna.orgg) were used to predict the potential miR-15a, miR-15b binding sites in the 3 0 UTR of PD-L1 mRNA. The predicted miR binding sites in the 3 0 UTR region of PD-L1 mRNA referred to as 3 0 UTR wild type (WT) and mutations in the miRNA binding sites (seven bases) of the 3 0 UTR region referred to as 3 0 UTR mutant were cloned into a luciferase vector pEZX-MT06 (Cat. #HmiT117860-MT06, Cat. #CS-HmiT117860-MT06-01; GeneCopoeia, Rockville, MD). NB cells were transfected with these luciferase plasmid constructs (1 mg) with miR-15a, miR-15b, or control miR oligonucleotides (100 nM), using lipofectamine reagent. After 48 h, the cells were lysed and the luciferase activity was measured by luciferase assay reagent (Promega) using the luminometer following manufacturer's instructions. The luciferase activity was expressed relative to control miR, which was set to 100.
Ago2 protein-associated RNA immunoprecipitation Immunoprecipitation of Ago2 protein-bound RNA-protein complexes was performed as described previously. 70 NB cells were transfected with miR-15a, miR-15b, and control miRs for 48 h. After transfection, cells were collected at 300Âg, washed with PBS, and lysed in polysome lysis buffer (5 mM MgCl2, 100 mM KCl, 10 mM HEPES [pH 7.0], 0.5% Nonidet P-40, 1 mM DTT, 100 U of RNase inhibitor/mL) supplemented with 20 mM EDTA and protease inhibitors on ice for 20 min. The lysates were centrifuged to remove cell debris, and clear lysates were collected. The lysates were diluted to 1:10 in freshly made NT2 buffer (50 mM Tris [pH 7.4], 150 mM NaCl, 1 mM MgCl2, 0.05% Nonidet P-40, 1 mM DTT, 100 U of RNase inhibitor/mL) supplemented with 20 mM EDTA and protease inhibitors and were incubated with a-Ago2 primary antibody at 4 C for 6 h, followed by incubation with protein A beads for an additional 2 h. The beads were then washed with NT2 buffer five times and treated with DNase I and proteinase K. The beads were eluted with NT2 buffer containing 0.1% SDS, and RNAs were extracted using phenol-chloroform and ethanol precipitation. The Ago2 bound immunoprecipitated RNA was used for the analysis of PD-L1 mRNA and miR-15a, miR-15b by RT-qPCR assay and normalized to IgG-immunoprecipitations as described. 50 Immunohistochemistry and immunofluorescence IHC and IF staining for PD-L1 and CD34 proteins was performed through our institutional core facility on a Discovery Ultra staining system (Roche Tissue Diagnostics, Ventana Medical Systems, Inc.). In brief, tissue samples were processed by fixing (70% ethanol), embedded in paraffin, and cut into 5-mm sections. The sections were deparaffinized in EZ prep, a mild detergent solution at 69 C for 24 min (Roche Diagnostics, #950-102), and antigen retrieval was achieved by incubation with cell-conditioning solution (Roche Diagnostics, #950-124), a Tris EDTA buffer (pH 8.2), for 32 min at 95 C. To prevent unspecific signal, we blocked endogenous peroxidase activity by treating with Inhibitor CM, available within the DISCOVERY ChromoMap DAB Kit (Roche Diagnostics, #760-159) for 8 min at RT. The sections were incubated with either rabbit monoclonal PD-L1 antibodies (clone 28-8, Abcam, #ab205921) or anti-mouse CD34 for endothelial cells (clone MEC14.7, BioLegend, #119302) at 37 C for 30 min using the Discovery Benchmark ULTRA advanced staining system. The HRP-conjugated secondary antibody (Roche Diagnostics, #760-4815) incubation was performed at 37 C for 16 min followed by diaminobenzidine (DAB) signal detection using DISCOVERY ChromoMap DAB Kit. The sections were counterstained with hematoxylin II for microscopic examination. The images were acquired using a fluorescence microscope (Leica Microsystems Inc, USA). Tumor microvessel density was assessed according to the criteria described by Weidner and colleagues. 106,107 First, the entire tumor section was scanned at low magnification to find areas having high vessel density. Then, images of four such fields were digitally acquired at high magnification (20Â). Microvessel density was quantified by manual counting using Photoshop (version 22.0.1), and an average of four fields from each tumor was determined. Four mice from each group were included in the analysis.

NB mouse xenograft model
All animal studies were conducted in accordance with and approved by the UNMC's Institutional Animal Care and Use Committee. In vivo immune response was evaluated using NB cell line-derived xenografts. We purchased 5-to 6-week-old C57/BL6 mice from The Jackson Laboratory (Bar Harbor, ME). We used only female mice for in vivo experiments to avoid biological variables. NB975, a murine NB cell line that stably expresses murine GFP-miR-15a, or non-targeting GFP-control miR, suspended in Matrigel (Corning Corporation, #CB-40230C), were subcutaneously inoculated into the right flank (3Â10 6 /100 mL) of C57/BL6 mice within the animal care facility at UNMC. Mice were monitored daily to ensure that the injection sites were healthy. On day 30, the mice were sacrificed, and the tumors were excised, photographed, and the final tumor volumes were measured using the equation: V=LÂW2/2 (L: tumor length; W: tumor width). Tumor tissues were lysed using standard protocols, prepared single-cell suspensions and then processed for immunophenotyping by flow cytometry using the methods described above.

Statistical analysis
All the statistical analyses were performed using GraphPad Prism software. The data were presented as the mean ± standard error or individual values derived from independent experiments and/or independent biological replicates. The Kruskal-Wallis test was used for multiple comparisons between groups. The survival analysis was calculated by the Kaplan-Meier method and the log-rank test. Statistical significance was calculated using a two-sided unpaired t-test. P values smaller than 0.05 were considered statistically significant.

DATA AND CODE AVAILABILITY
Data were generated by the authors and are available upon request from the corresponding author.


Inducible Resistance to Amikacin in Mycobacterium abscessus Isolated in Beijing, China

Abstract We aimed to determine the prevalence of amikacin (AMK) resistance of clinical Mycobacterium abscessus (MAB) isolates and to investigate if AMK resistance was induced by AMK exposure. A total of 75 MAB isolates underwent susceptibility testing for AMK after 3 and 14 days of incubation, respectively. The partial fragment of the rrs gene conferring AMK resistance was sequenced. The MIC values for AMK ranged from 0.5 to 128 μg/mL, with MIC50 and MIC90 values of 2 and 32 μg/mL, respectively. In addition, 9.3% of isolates (7/75) were resistant to AMK, all of which harbored a mutation within the rrs locus, including six with A1408G mutation and one with a C1409T mutation. Of note, the MICs of three isolates were significantly increased from 2 μg/mL to 64 μg/mL (one isolate) and 2 μg/mL to 128 μg/mL (two isolates), suggesting that three of the MAB isolates had inducible resistance to AMK. In conclusion, our data demonstrate that approximately one-tenth of clinical MAB isolates in Beijing harbored AMK resistance due to the acquisition of rrs mutations. Additionally, we firstly identified that intrinsic AMK resistance is inducible in MAB isolates, highlighting the urgent need to establish a proper method for the in vitro detection of AMK susceptibility in MAB.

Introduction
Nontuberculous mycobacteria (NTM) are increasingly recognized as important pathogens causing infections in both humans and animals. 1 Generally, NTM are classified as either slow or rapid growers based on their growth rate. 2 Of the rapidly growing mycobacteria (RGM) species, Mycobacterium abscessus (MAB) has been shown to be the most frequently encountered causative agent in human infections, accounting for approximately 65%-80% of RGM lung diseases. 1,3 The leading threat posed by this species is due to its resistance to many antibiotics, thereby resulting in unsatisfactory clinical outcomes. 4 Amikacin (AMK), a member of the aminoglycoside family of antibiotics, exhibits antimicrobial effects by targeting 16S rRNA and inhibiting protein synthesis. 5 It is a cornerstone of antimicrobial chemotherapy treatment for MAB infections in view of its promising efficacy in vitro and in vivo. 6 Multiple molecular mechanisms contribute to AMK resistance in MAB isolates, including target mutations and increased efflux transporters. 4 The major mechanism conferring AMK resistance is spontaneous single mutations within the rrs gene encoding 16S rRNA, yielding a high level of AMK resistance in MAB clinical isolates. 7 Most acquired amikacin-resistant MAB isolates had A1408G or C1409T mutations and other mutations in rrs. [8][9][10] Previous experimental studies have demonstrated that the aminoglycoside-modifying enzyme gene eis2 and the multidrug efflux transporter gene (tap) are also involved in AMK resistance in mycobacteria. 11,12 Another important intrinsic mechanism for macrolides in MAB has been described as inducible resistance conferred by a inducible erythromycin ribosomal methylase erm(41) gene. 13 Considering that the 16S and 23S rRNAs are the major RNA components in the ribosomal subunits, we hypothesized whether this mechanism was involved in AMK resistance. Thus, the primary objectives of this study were to determine the prevalence of AMK resistance of clinical MAB isolates and to investigate if AMK resistance was induced by AMK exposure.

Materials and Methods
For this study, a set of clinical MAB isolates and the reference MAB strain ATCC 19977 were included in our analysis; the AMK susceptibilities of these organisms were tested by using the broth microdilution assay following the guidelines from the Clinical and Laboratory Standards Institute (CLSI). 14 Briefly, freshly cultured colonies were harvested from the surface of Löwenstein-Jensen (L-J) medium and then a bacteria suspension equal to the turbidity of a 0.5 McFarland standard was prepared in normal saline. This suspension was diluted 1:20 in cation-adjusted Mueller-Hinton broth (CAMHB) and incubated in a range of AMK concentrations (0.13 to 128 μg/mL). The inoculated plates were incubated at 30°C for 72 hours. To assess phenotypic induction, the plates were submitted to an extended incubation with reading after 14 days of incubation, respectively. The MIC was defined as the lowest concentration of antibiotics that inhibited visible growth of mycobacteria.
Mycobacterial DNA was extracted using the MasterPure TM Complete DNA and RNA Purification Kit (Lucigen corporation, Wisconsin, USA) and was used as the template for PCR amplification of the rrs gene conferring acquired resistance to AMK in mycobacteria. The 351 bp fragment of the rrs gene was amplified using the primers rrs-F (5′-ATGACGTCAAGTCATCATGCC-3′) and rrs-R (5′-AGGTGATCCAGCCGCACCTTC-3′) as previously reported. 10 The PCR products were sent to the Tsingke Company (Beijing, China) and sequenced with the rrs-F primer by Sanger sequencing method. DNA sequences were aligned with the homologous sequences of the reference MAB strain ATCC 19977 using multiple sequence alignments (http://www.ncbi.nlm.nih.gov/BLAST). The mutations were described corresponding to the positions found in the 16S RNA of Escherichia coli.

Results
All 75 MAB isolates underwent susceptibility testing for AMK, the results of which are summarized in Table 1. The MIC values for AMK ranged from 0.5 to 128 μg/mL, with MIC50 and MIC90 values of 2 and 32 μg/mL, respectively. Using the CLSI-endorsed 64 μg/mL as the breakpoint, 9.3% (7/75) were resistant to AMK, comprising of one with an MIC value of 64 μg/mL and six with MIC values of 128 μg/mL. Out of 68 AMK-susceptible MAB isolates on day 3, eight showed at least an 8-fold increase in MIC of AMK on day 14. Of note, the MICs of three isolates were significantly increased from 2 μg/mL to 64 μg/mL (one isolate) and 2 μg/mL to 128 μg/mL (two isolates), suggesting that three of the MAB isolates had inducible resistance to AMK (Table 2).
We further analyzed the rrs polymorphism conferring AMK resistance. All seven isolates with acquired resistance to AMK harbored a mutation within the rrs locus. Of these isolates, six had a A1408G mutation and one had a C1409T mutation.

Discussion
Due to the intrinsic resistance to a diverse array of antimicrobial agents, clinical management of individuals infected with MAB is challenging since it is difficult to formulate effective therapy regimens; 4 however, amikacin holds great promise for the treatment of MAB infections. 6 In this study, our data demonstrate that approximately one-tenth of clinical MAB isolates harbored AMK resistance due to the acquisition of rrs mutations. This rate was significantly higher than those from the USA (4.8%), 15 South Korea (3.4%), 16 France (4.8%), 17 and Shanghai (3.9%). 18 Several potential explanations could be responsible for the high AMK resistance rate observed in our study. On the one hand, there was a great range in the literature regarding the drug susceptibilities of MAB and M. massiliense. A recent study by Guo and colleagues demonstrated that the rate of AMK resistance for MAB (4.1%) was dramatically higher than that for M. massiliense in Shanghai (0.0%). 18 In contrast, a systematic analysis of a global genome dataset found that a slightly higher proportion of M. massiliense genomes contained the canonical rrs mutation conferring high-level AMK resistance than MAB. 19 Hence, we speculate that the diverse results across studies may be due to the differences in subspecies compositions. Secondly, the majority of microbiology laboratories do not have the capability to differentiate Mycobacterium tuberculosis from NTM. Therefore, patients infected with NTM could be misdiagnosed as drug-resistant TB, resulting in exposure to inappropriate medications. As the former cornerstone drug for drug-resistant TB, AMK is widely used for the treatment of multidrug-resistant and extensively drug-resistant TB.
Although it provides promising efficacy against MAB, the ineffective background regimens would facilitate the emergence of AMK-resistant mycobacteria, which could partly explain the high ratio of AMK resistance in MAB observed in the present study. In view of this potential threat, species identification of mycobacteria is essential to initiate treatment promptly for patients with symptoms suggestive of tuberculosis in this country with a high burden of TB. We also firstly identified that intrinsic AMK resistance is inducible in MAB isolates. Indeed, there is increasing evidence that the methylation of 16S rRNA is an important mechanism of resistance against aminoglycosides among gram-negative bacteria and actinomycetes. 20,21 In the latter case, multiple actinomycetes are intrinsically resistant to aminoglycosides via the inactivation of aminoglycosides by persistent methylation of 16S rRNA. 20 In addition to mutations in the rrs allele, a series of genes, including eis (Rv2416c), tap (Rv1258c), and whiB7 (Rv3197), are reported to confer low-level resistance to AMK in mycobacteria. However, our primary data demonstrate that a proportion of isolates had MICs of 128 μg/mL. Another potential mechanism found in clinical settings is acetylation by the aminoglycoside 6′-N-acetyltransferase type Ib, which prevents the interaction between amikacin and the ribosome; this has been noted in several gram-positive and gram-negative bacteria. 22,23 However, the expression of this gene was not inducible under AMK exposure, which may be not directly associated with inducible AMK resistance as observed in MAB isolates. As a consequence, we hypothesize that an inducible 16S ribosome

2289
methylase may serve as the underlying mechanism associated with inducible AMK resistance, presenting an efficient means to avoid inhibition of its own protein synthesis without additional fitness cost that makes it untenable in vivo. Further experimental research is required to identify the core 16S rRNA methylase involved in this important phenomenon. Although the exact mechanism remains unclear, our findings provide an important hint for the clinical management of MAB patients. Inducible AMK resistance may manifest in diverse clinical outcomes in MAB cases after treatment with AMKcontaining regimens. An extended incubation period is required to accurately identify AMK resistance of MAB isolates. We also acknowledged several limitations to our study. First, although we observed an increase in MICs after AMK exposure, the mechanism for induced AMK resistance was not clearly illustrated in MAB isolates. Further investigations will be conducted to elucidate its underlying mechanisms. Second, given the fact that all the NTM patients enrolled in our study were outpatients, the lack of treatment outcomes from the outpatients inhibits us from verifying the relationship between induced AMK resistance and treatment outcome. Third, the high proportion of AMK resistance was attributed to previous irregular adherence to treatment as well as antibiotic abuse, but the previous exposure history to aminoglycosides was not clear. Finally, recent studies have declared a significant difference in induced CLA resistance between MAB and Mycobacterium massiliense. 24 In the present study, we only included MATB isolates. Thus, it is difficult to clarify whether induced AMK resistance differs significantly among subspecies.
In conclusion, our data demonstrate that approximately one-tenth of clinical MAB isolates in Beijing harbored AMK resistance due to the acquisition of rrs mutations. Additionally, we firstly identified that intrinsic AMK resistance is inducible in MAB isolates, highlighting the urgent need to establish a proper method for the in vitro detection of AMK susceptibility in MAB. Further experimental research is required to identify the core 16S rRNA methylase involved in AMK-inducible resistance in MAB.

Ethics Approval and Consent to Participate
The study was conducted according to the principles found in the Declaration of Helsinki. This study was approved by the Ethics Committee of the Beijing Chest Hospital, affiliated with Capital Medical University. All study participants signed an informed consent that agreed to the anonymous use of clinical samples.


s2fieldsofstudy:
Biology
Environmental Science
Medicine

text:
Species-Specific Patterns of Gut Metabolic Modules in Dutch Individuals with Different Dietary Habits

ABSTRACT Diet is an important determinant of the human gut microbiome. Here, we analyzed fecal metagenomes of Dutch adults following omnivorous, pescatarian, vegan, and vegetarian diets. We compared the taxonomic composition of individuals from our study with publicly available gut metagenomes from westernized and non-westernized societies. We observed that, despite long-term transition to diets rich in plant fibers (vegan or vegetarian), the microbiomes of these were typical of westernized populations, and similar in composition to omnivores. Although there were no major differences in metabolic modules, we identified differences in the species that contributed to particular functions, such as carbohydrate degradation and short-chain fatty acid metabolism. Overall, this study shows functional redundancy of the microbiomes among westernized populations, which is independent of long-term individual dietary habits. IMPORTANCE Diet is an important modulator of the human gut microbiome, which is susceptible to increased consumption of plant fibers in vegan or vegetarian lifestyles. To investigate this, we compared the gut microbiome of Dutch adults following omnivorous, pescatarian, vegan and vegetarian diets. We did not observe major differences in the gut microbiome composition and function between individuals with different dietary habits. However, we observed differences in the species that contribute to the core functions of the gut microbiome. Our study thus emphasizes the need to better understand the species-specific functional changes associated with dietary habits in the human gut microbiome.

I n westernized societies traditionally known to have an omnivorous diet, some individuals have transitioned toward vegetarian and vegan diets. These diets are enriched in microbiota-accessible carbohydrates (MACs) that are not digested and absorbed by the host but support a complex microbiome with supposedly beneficial effect on host health (1). MACs like resistant starch, inulin, xylan, pectin and similar complex substrates favor the growth of specialized taxa like Ruminococcus bromii, Roseburia intestinalis, and Faecalibacterium prausnitzii (2). High-fiber diets rich in MACs are reported to have beneficial effects on host physiology via enhanced butyrate production and promoting higher bacterial diversity (3). Short-term (4-day) dietary interventions have reported reversible changes in the microbiome, indicating rapid adaptation of the microbiota to the changing diet (3). Individuals consuming high-fiber diets for longer periods are shown to have a Prevotella enriched microbiome, which is stable at least up to 10 days, indicating that long-term dietary habits can select for specific taxa that form a stable community (4). However, in our recent study comparing gut microbiome of individuals with different long-term dietary habits ($6 months), we did not observe major differences in microbiome diversity and resistome (5). As high functional redundancy is a hallmark of the gut microbiome, where multiple species can carry out similar functional roles (6, 7), linking abundances of specific bacterial taxa with their contribution to functional potential may provide a better understanding of species-specific associations with dietary habits.
Here, we investigated whether different long-term dietary habits ($6 months) result in a gut microbiome composition that deviates from that observed in westernized populations and identify differences in species-specific contributions to key metabolic pathways.
The participants in our study (NLD-VEGA) were classified into vegetarian, vegan or pescatarian if they followed these diets for at least 6 months, and omnivores if they eat meat at least three time per week (5,8). Details about the NLD-VEGA study population are described in Table S1, and metagenomic data analysis is provided in the supplementary information (Text S1). Information on public data used in this study are provided in Table S2. Taxonomic composition and functional annotations were done using the biobakery tools (9). The fecal microbiomes of participants in the NLD-VEGA study were similar to westernized populations (Wilcox test, P , 2.2e 2 16) and showed high similarity with individuals from another independent Dutch study, Netherlands LifeLines cohort (NLD-LifeLines; n = 1039) (Fig. 1a, b and c) (10). The NLD-VEGA and NLD-LifeLines cohorts shared 39 (61%) of the core microbes (min. relative abundance 0.0001% and 75% prevalence) (Fig. S1). This common core (NLD-VEGA 1 NLD-LifeLines) consists of species involved in polysaccharide degradation and short-chain fatty acid production (e.g., Ruminococcus bromii, Eubacterium rectale, Faecalibacterium prausnitzii, Eubacterium hallii, Coprococcus catus, Anaerostipes hadrus) (Fig. 1d) (11)(12)(13).
We investigated the effect of long-term dietary habits (i.e., omnivore, vegetarian, vegan and pescatarian) on so-called VANISH (volatile and/or associated negatively with industrialized societies of humans) taxa in our cohort (14). VANISH taxa, i.e., Succinovibrionaceae, Paraprevotellaceae, Prevotellaceae, and Spirochaetaceae, are considered to be lost with industrialization, and therefore more prevalent in non-westernized microbiomes. As previously observed in westernized diets (14,15), the NLD-VEGA population showed loss of certain VANISH taxa and overall low abundance of Prevotellaceae (Fig. S2a), with no significant differences between the diet groups. The BloSSUM taxa (bloom or selected in societies of urbanization/modernization) were observed in all diet groups in the NLD-VEGA (Fig. S2b). Furthermore, western-style diet is often associated with a higher Bacteroidetes, while traditional high-fiber diet rich in vegetables and fish is associated with a high Firmicutes abundance (1,16). Here, we observed a significantly higher Firmicutes/Bacteroidetes ratio in the pescatarian and vegan group compared to the omnivore and vegetarian group (Fig. 1e). While diets rich in plant-based fibers were also reported to be associated with a higher Prevotella/Bacteroides (P/B) ratio (17), we did not find these associations in our data set (Fig. 1f). These results highlight the potential non-universality of previously reported associations of higher P/B ratio with a plant-based diet (18). In addition to diet, possibly other factors (e.g., environmental exposures) play a significant role in determining the P/B ratio in non-westernized individuals. Therefore, increasing dietary fibers or becoming vegetarian/ vegan alone might not be enough to increase the P/B ratio in westernized individuals.
Metabolic potential of the gut microbiome is key for maintenance of host health. Here, we investigated the functional repertoire in gut metagenomes and species-specific contributions to key metabolic processes, and its relation to the different long-term dietary habits. We used a curated version of the previously described Gut Metabolic Modules (GMMs) to group KEGG orthologues (KOs) into metabolic modules (see supplementary data) (18). To these GMMs, we added modules related to amino acid metabolism and refined GMMs based on a synthetic microbiome used in a previous study (19). Except for lipid degradation, abundances of amino acid, carbohydrate and glycoprotein degradation modules were not significantly different between diet groups (Fig. S3). However, species-specific differences in abundances of the bacteria involved in the different processes were observed (Fig. 2). We observed nine species that differed in their contribution to amino acid degradation module (Fig. 1a, Kruskal-Wallis test, adj. P , 0.05, followed by Dunn's Test adj. P , 0.05). For carbohydrate degradation module, eight species differed in their contribution between groups Long-Term Dietary Habits and Species-Specific Differences mSphere (Fig. 1b). No species were observed to differ in relative contribution to amino acid and carbohydrate degradation modules between vegans and vegetarians. For glycoprotein and lipid degradation, Bacteroides ovatus was observed to be significantly different between omnivore and vegetarian (Fig. 1c). Notably, no statistically significant species-specific differences in contribution to lipid metabolism was observed in our study.

Long-Term Dietary Habits and Species-Specific Differences mSphere
Except for butyrate metabolism, the GMMs for acetate, propionate and lactate metabolisms did not differ between diet groups (Fig. S3). Similar with higher level GMMs, speciesspecific differences in contribution to these modules were observed. Ruminococcus torques had a significantly higher contribution in omnivores compared to vegan group. Three species showed differences in contribution to lactate metabolism, Streptococcus thermophilus and S. thermophilus CAG236 in omnivores compared to vegan and vegetarian groups. Roseburia hominis was higher contributor in pescatarians compared to omnivores. Ruminococcus torques was higher contributor to propionate metabolism module in omnivores. Although R. torques is not yet reported to produce propionate, it does have  degradation (a, b, and c). Species significantly different between diet groups and their relative contribution to gut metabolic modules for short-chain fatty acid metabolism (d, e, and f). Significance in differences were calculated using Kruskal-Wallis test followed by Dunn's test for pairwise comparisons and adjusted for multiple comparisons with the Benjamini & Hochberg method. The differences were considered significant when adjusted P-value was less than , 0.05.
Long-Term Dietary Habits and Species-Specific Differences mSphere key genes for producing propionate via the propanediol pathway (20). R. torques is known to produce lactate, an important precursor for propionate production via the acrylate pathway. R. torques is reported to be associated with proinflammation (21). Flavonifractor plautii contributed higher to propionate metabolism module in omnivores compared to pescatarians. This bacterium is known to utilize simple sugars and amino acids to produce butyrate and was recently shown to be associated with young-onset colorectal cancer (22,23). No statistically significant species-specific differences in contribution to butyrate metabolism between the four diet group was observed in our study.
Overall, we see that despite showing modest differences in abundances of key GMMs, species-specific differences in contribution to key metabolic pathways were observed between the diet groups. This further highlights the widely accepted functional redundancy of the gut microbiome. Future research investigating species specific activity using a combination of multi-omics and in vitro/in vivo investigations can help to improve our understanding of diet mediated effects on human health.
Data availability. The raw sequencing data are available at EMBL-ENA under accession number PRJEB45944. The codes to reproduce the analysis are available from GitHub (https://github.com/RIVM-IIV-Microbiome/VEGA-2021).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, DOCX file, 0.03 MB.

ACKNOWLEDGMENTS
This work was supported by the Netherlands Centre of One Health (NCOH metagenome call), the Dutch Ministry of Health, Welfare and Sport (grant number V/050305/AA and the Strategic Program RIVM [SPR]), and Utrecht Exposome Hub of Utrecht Life Sciences (www .uu.nl/exposome), funded by the Executive Board of Utrecht University. We thank Jeroen Frank and Jeroen Laros for support with bioinformatic resources.
The authors have no conflicts of interest to declare.