A group of 15 children with ALL (8 B-ALL and 7 T-ALL) aged between 0 and 18 years were considered for the study of differentially expressed microRNA (DE-miRNA) in bone marrow cells according to the marrow stroma sHLA-G levels. For real-time quantitative PCR validation experiments, we compared the levels of DE-miRNA in another group of ALL patients (23 B-ALL and 11 T-ALL). To demonstrate that the effect observed was related to the lymphoid cell type, we also evaluated samples from 31 children with acute myeloid leukemia (AML). We also included a control group with 14 samples from children whose myelogram confirmed the absence of leukemia. The expressions of the HLA-G5 and RREB1 target genes were evaluated in the bone marrow cells of 19 children with B-ALL, 8 with T-ALL, and 28 with AML. All patients were referred, diagnosed, and treated at the IMIP Hospital, Recife, Brazil. Bone marrow aspirates were obtained from each patient at admission and submitted for the isolation of mononuclear cell fractioning for leukemia diagnosis confirmation, performed as previously described (Marques et al., 2011). The samples were stored under –80 °C conditions provided by a laboratory deep freezer, which was protected against power outage by an uninterruptible power supply (UPS) and emergency line. All the patients with leukemia presented at least 70% of blasts in the bone marrow. The samples were obtained after the children’s legal guardians provided informed consent, approving their participation in the study. The study protocol was previously approved by the local ethics committee (CAAE: #13296913.3.0000.5190 and #0073.0.095.000-10). The patients’ (age and sex) and blast (immunophenotype and genetic alterations) features are shown in Table 1.

A sandwich ELISA assay was used to measure the soluble HLA-G (shredded HLA-G1 and HLA-G5 isoforms) levels, following the manufacturer’s instructions (BioVendor Laboratory Medicine, Inc., Czech Republic), with the limit of detection of 0.6 Units/mL. Our previous study detected an average of 200 U/mL ± 25 SD (standard deviation) of sHLA-G levels in the bone marrow stroma of healthy children (Almeida et al., 2018). Based on this, patients presenting between 150 and 250 U/mL of sHLA-G levels in the bone marrow stroma, that is, 200 U/mL plus two standard deviations above or below, were defined as the intermediate producers, those presenting with more than 250 U/mL of sHLA-G were high producers, while those who produced less than 150 U/mL were low producers.

We used the miRNA sequencing database to evaluate the miRNA expression related to the sHLA-G levels in the marrow stroma. Total RNA extraction, quality assessment, library construction, and miRNA sequencing were performed as described previously (Almeida et al., 2019). miRNA sequencing data have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under the accession number E-MTAB-11621. The sequencing analysis included read quality control and contamination assessment using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Cutadapt (Martin, 2011) programs considering a Q-score ≥30 and reads with a length ≥ 17 nucleotides. We used Bowtie (http://bowtie-bio.sourceforge.net/index.shtml) for indexing of human reference genome hg38 version, which is deposited in the UCSC Genome Browser (https://genome.ucsc.edu/). The miRDeep2 2.0.0.8 software (Friedländer et al., 2008) was applied for sequence alignment and miRNA identification, considering miRBase release 21 (http://www.mirbase.org/) (Griffiths-Jones et al., 2006; Kozomara and Griffiths-Jones, 2010). Differentially expressed (DE) miRNA profiles were obtained using the edgeR package (Robinson et al., 2010) and the standard analysis and quantile normalization parameters in the R software (https://cran.r-project.org/), considering at least 20 reads in a minimum of 1 sample, a false discovery rate (FDR) ≤0.05, and a log fold change (logFC) cutoff point of 1 or −1. A comparison of the bone marrow miRNA levels between lower versus higher sHLA-G producers was performed. Target prediction of DE-miRNAs was performed using the miRWalk 2.0 (Dweep et al., 2015), and functional annotation was determined by DAVID tools v.6.7 (Huang et al., 2009; Huang et al., 2009), considering the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms: biological process and FAT level, both with Benjamini–Hochberg (BH)–corrected p-values ≤ 0.05. The DE-miRNA alignment with the HLA-G gene (NG_029039.1) and mRNA sequences (NM_002127.5) was performed using the RNAhybrid v.2.2 tool (Kruger and Rehmsmeier, 2006), considering the essential features for the interaction of the two molecules, that is, Watson and Crick base pairing, few gaps in the interaction, especially on the seed sequence, seed (2–8 miRNA nucleotide), low free energy (≤−20 Kcal), and interaction with target 3′UTR, coding sequence, and promoter region (Castelli et al., 2014).

A search for genes encoding proteins related to HLA-G transcription’s positive and negative regulation was performed, considering previous studies that describe or review the action of such molecules (Moreau et al., 1999; Gobin et al., 2002; Flajollet et al., 2009; Castelli et al., 2014; Yaghi et al., 2016). The positive regulators that were considered were CREB1, CREBBP, JUN, ATF2, IRF1, HIF1A, and IL10. The negative regulators that were sought were RREB1, HDAC1, CTBP1, and CTBP2. We also considered REST, EHMT1, ZEB1, ZEB2, ZNF217, and LSD1 genes since the proteins are members of the CTBP core complex (Shi et al., 2003; Shi et al., 2004) and may exert an indirect influence on HLA-G expression. All the positive and negative regulators of HLA-G those were considered were analyzed for their ability to interact with the differentially expressed miRNAs in this study, according to the miRTarBase v. 8.0, a database of experimentally validated interactions (Chou et al., 2018).

For validation experiments, we selected the four miRNAs most likely to target the HLA-G gene (NG_029039.1) and messenger RNA (NM_002127.5) based on the sequence alignment analysis (RNAhybrid v.2.2) (Krûger, Rehmsmeier, 2006). The representative scheme showing the interaction site between HLA-G and the four miRNAs selected for validation was constructed using the ApE v2.0.61 software (https://jorgensen.biology.utah.edu/wayned/ape/). The TaqMan Advanced miRNA cDNA Synthesis Kit (Life Technologies, Foster City, California, USA), TaqMan Advanced miRNA Assay (reference: miR-191-5p; targets: miR-5096, miR-4516, miR-4488, miR-486-5p; Life Technologies), and TaqMan Fast Advanced Master Mix (Life Technologies) were used according to the manufacturer’s instructions to evaluate the miRNA expression. Reverse transcription PCR (RT-PCR) assays were performed in a SimpliAmp Thermal Cycler (Applied Biosystems, Foster City, California, USA) and quantitative PCR (q-PCR) in a QuantStudio 5 Real-Time System (Applied Biosystems) and 7500 Real-Time System (Applied Biosystems) according to the manufacturer’s instructions.

To study the relative expression of HLA-G5 and RREB1, cDNA synthesis was performed from total RNA using the enzyme M-MLV-RT 200 U/µL (Invitrogen, Carlsbad, California, USA) and SimpliAmp Thermal Cycler equipment (Applied Biosystems). HLA-G5 primers have been described in Gomes et al. (2018), and they were designed to target all RREB1 isoforms (RREB-1F: 5′-AAA​GAT​GGT​AGA​AGA​CGG​G-3′ and RREB-1R: 5′-GTG​GGT​TAT​CTG​AAT​GGG​TC-3′). Expression was performed using the SYBR Green DNA intercalator (Applied Biosystems).

The normality distribution of the miRNA–mRNA expression was determined using the Shapiro–Wilk and Kolmogorov–Smirnov tests. For comparison between two or three groups, the Mann–Whitney U and Kruskal–Wallis tests were used, respectively. A Spearman’s correlation coefficient analysis was performed between miRNA and mRNA expressions. In different experiments, the number of samples may have differed due to the shortage of clinical samples that did not allow all analyses to be performed. The GraphPad Prism V.5.01 (GraphPad Software, Inc.) was used to perform the analyses, considering a significant p-value ≤ 0.05. For miRNA relative expression analysis, ∆Ct (cycle threshold) values were determined based on the following equation: ∆Ct = Ct (target miRNA) − Ct (reference miRNA). The Ct values were the average duplicates with a standard deviation (SD) ≤ 0.5. The same equation and parameters were used to calculate the mRNA expression, considering HLA-G5 or RREB1 as the target gene and GAPDH as the reference gene.

Bone marrow sHLA-G levels in childhood AML, T-ALL, and B-ALL showed no statistical differences (p = 0.3483). There were low, intermediate, and high sHLA-G producers in each leukemia subtype (Figure 1).

The analysis of differentially expressed miRNA profiles in the bone marrow cells of non-treated children with ALL revealed 10 miRNAs upregulated in high sHLA-G producers (logFC >2.0) when compared with low sHLA-G producers (Table 2).

The analysis of target gene prediction with the 10 miRNAs showed 14,518 potential gene targets, of which only the hsa-miR-5096 was predicted as a putative regulator of HLA-G mRNA by three different algorithms and, by less number, the hsa-miR-4516, hsa-miR-4488, and hsa-miR-486-5p miRNAs (Figure 2). All four miRNAs presented several anchor sites at the promoter and coding regions of the HLA-G gene. Some of the binding sites of miRNAs are in transcription factor zones. The cAMP-responsive element (CRE) is a predicted site for hsa-miR-5096 binding; the heat shock element (HSE) for hsa-miR-486-5p; the hypoxia-responsive element (HRE) for hsa-miR-4516 and hsa-miR-4488; the Kappa B1, Kappa B2 (NF-κB responsive element), interferon-stimulated response element (ISRE) module, and the SXY module for hsa-miR-4488, hsa-miR-4516, and hsa-miR-486-5p; and the Ras-responsive element (RRE) and progesterone-responsive element (PRE) for hsa-miRNA-4516 and hsa-miR-4488. The hsa-miR-5096 did not bind to any of these transcription-binding sites. Only the hsa-miR-4516 targets the HLA-G 3′ untranslated region at the position covering the +3035 C/T polymorphic site.

The analysis of the miRNA gene targets for functional annotation revealed several biological pathways involving genes already described in the literature that may be positive or negative regulators for the HLA-G gene. Table 3 shows the genes involved in the induction or repression of HLA-G transcription and their putative miRNA regulators identified in this study.

Considering the 10 most significant KEGG pathways related to all upregulated miRNAs in high sHLA-G producers, 6 included the genes encoding known positive (CREB1, CREBBP, JUN, and IL10) and negative (CTBP1/2 and HDAC1) regulators of HLA-G expression (Figure 3), as well in other pathways associated with leukemogenesis, namely, hsa04310:Wnt, hsa04660:T-cell receptor, hsa04062:chemokine, hsa04662:B-cell receptor, hsa04330:Notch, and hsa04350:TGF-beta signaling pathways. Most of the statistically significant GO biological processes involved in regulating transcription and cell signaling cascade include inducers (CREB1, CREBBP, ATF2, JUN, and IL10) and repressors (RREB1, CTBP1/2, and HDAC1) of the HLA-G expression (Figure 3).

The comparison of the miRNA levels in the bone marrow showed that miR-486-5p, miR-4488, and miR-5096 levels were significantly higher in controls than in ALL, particularly B-ALL, and only miR-486-5p was higher in controls than in T-ALL (p < 0.05). No significant differences in miRNA levels were observed between controls and AML bone marrow samples (p > 0.05). In addition, the AML samples showed higher miR-4488 and miR-5096 levels than did B-ALL and higher miR-486-5p levels than did T-ALL (p < 0.05) (Figure 4).

In T-ALL, the hsa-miR-5096 levels correlated positively with the HLA-G5 mRNA expression (rho = 1; p = 0.0167) (Figure 5). In myeloid leukemia, the hsa-miR-4516 (rho = 0.4638; p = 0.0258) and hsa-miR-4488 (rho = 0.6509, p = 0.0008) levels were also positively correlated with the HLA-G5 mRNA levels. However, the increase in HLA-G5 mRNA expression was translated into a significant decrease in sHLA-G only in myeloid leukemia, with moderate and significant Spearman’s coefficient (rho = 0.475; p = 0.0397), but neither in B-ALL nor T-ALL (Figure 6). This was assumed considering that the delta Ct values are inversely proportional to the mRNA levels.

In addition, increased hsa-miR-5096 (rho = 0.72; p = 0.0144) and hsa-miR-4516 (rho = 0.67; p = 0.0277) levels (low ΔCt) correlated with decreased sHLA-G protein levels in T-ALL, but only hsa-miR-5096 correlated also with the HLA-G mRNA (Figure 7), but only hsa-miR-5096 correlated also with the HLA-G mRNA.

For a detailed analysis, the samples were categorized according to the bone marrow miRNA levels in the low or high miRNA level group, and sHLA-G levels in both groups were compared. In T-ALL, patients with high levels of hsa-miR-5096 and miR-4516 had a median sHLA-G value of 46 U/mL, while patients with low levels of miRNA had a median sHLA-G value of 200 U/mL (p = 0.0519). Overall, high miRNA expressions were associated with homogenous low sHLA-G levels, while low miRNA levels were associated with largely variable sHLA-G levels, which contributed to the borderline significance of the differences. In B-ALL, the groups of low and high miRNA levels were not capable of segregating samples with different sHLA-G levels. In AML, the difference between the median value of sHLA-G between the low- and high-miRNA-level groups was not significant (Figure 8).

Considering that the RREB1 gene is a target for the four studied miRNAs and that the RREB-1 protein has three potential binding sites in the HLA-G gene promoter, the relationship between the RREB1 mRNA levels and each miRNA and HLA-G5 mRNA levels were evaluated. The results revealed that only in B-ALL, the RREB1 and HLA-G5 mRNA expressions were positively correlated (rho = 0.5632, p = 0.0018). In addition, only the hsa-miR-4488 correlated positively with RREB1 mRNA expression (rho = 0.4368, p = 0.0615), but it did not reach significance.