MiR-150-5p Overexpression in Triple-Negative Breast Cancer Contributes to the In Vitro Aggressiveness of This Breast Cancer Subtype

Simple Summary Triple-negative breast cancer (TNBC) is a clinically aggressive type of breast cancer. MicroRNAs (miRNAs) are small molecules that regulate the expression of genes involved in tumor cell signaling. The miR-150-5p is frequently deregulated in cancer, with expression and mode of action varying according to the cancer type. In this study, we investigated the expression levels of miR-150-5p in TNBC, its association with clinical and pathological features of patients, and its role in modulating TNBC cell proliferation, migration, and drug resistance. Our results suggest that miR-150-5p is highly expressed in TNBC and that miR-150-5p expression levels are associated with tumor grade, patient survival, and ethnicity. Our findings also indicate that miR-150-5p contributes to the aggressive phenotypes of TNBC cells in vitro. Abstract MiR-150-5p is frequently deregulated in cancer, with expression and mode of action varying according to the tumor type. Here, we investigated the expression levels and role of miR-150-5p in the aggressive breast cancer subtype triple-negative breast cancer (TNBC). MiR-150-5p expression levels were analyzed in tissue samples from 113 patients with invasive breast cancer (56 TNBC and 57 non-TNBC) and 41 adjacent non-tumor tissues (ANT). Overexpression of miR-150-5p was observed in tumor tissues compared with ANT tissues and in TNBC compared with non-TNBC tissues. MiR-150-5p expression levels were significantly associated with high tumor grades and the Caucasian ethnicity. Interestingly, high miR-150-5p levels were associated with prolonged overall survival. Manipulation of miR-150-5p expression in TNBC cells modulated cell proliferation, clonogenicity, migration, and drug resistance. Manipulation of miR-150-5p expression also resulted in altered expression of its mRNA targets, including epithelial-to-mesenchymal transition markers, MYB, and members of the SRC pathway. These findings suggest that miR-150-5p is overexpressed in TNBC and contributes to the aggressiveness of TNBC cells in vitro.

Introduction
Triple-negative breast cancer (TNBC) is a clinically aggressive breast cancer subtype that contributes to the patient's unfavorable prognosis, resistance to therapy, and a high risk to develop metastatic lesions [1]. TNBCs exhibit the classical "hallmarks" of invasive tumors, including high levels of genomic instability, impaired function of BRCA1 and estrogen signaling, and high frequency of mitotic and angiogenic events [2][3][4].
MicroRNAs (miRNAs) are a class of endogenous non-coding RNA molecules that regulate key biological processes, including cell development and malignant transformation. In cancer, miRNAs can act as oncogenes or tumor suppressor genes by regulating the expression of genes that control cell proliferation, differentiation, apoptosis, invasion, and metastasis [5][6][7]. Therefore, miRNAs have emerged as novel and promising therapeutic targets in cancer [8,9].
In vitro studies on TNBC cell lines have demonstrated that miR-150-5p overexpression is present in cells with more aggressive characteristics, such as the ones with high proliferation rates, the ability to form colonies, migrate, invade, and develop distant metastasis [28][29][30]. However, other studies have shown that miR-150-5p overexpression inhibits breast tumorigenesis [29,31]. Distinct mechanisms have been proposed to underlie the varied roles of miR-150-5p in modulating these phenotypes. These mechanisms include the direct or indirect regulation of its mRNA targets, such as MYB, HGMA2, TP53, SRC1N1, and ZEB1 [25,26,28], as well as its interaction with long non-coding RNAs [25,26,31,32].
In this study, we investigated the expression levels of miR-150-5p in clinically wellannotated TNBC specimens and its role in conferring aggressive tumor phenotypes in TNBC cells.

Clinical Samples
Non-identified breast cancer patients' samples of formalin-fixed paraffin-embedded (FFPE) material was obtained from 113 invasive breast carcinoma tissues and 41 adjacent non-tumor (ANT) tissues. The samples were retrospectively collected from the Histopathology Shared Resources (HSTR) at the Lombardi Comprehensive Cancer Center, Georgetown University, under IRB approval (IRB#1992-048). The tumors were classified as TNBC (n = 56) or non-TNBC (n = 57) subtype by the expression of ER, PR, and HER2 protein markers according to international guidelines [33,34]. Paired ANT tissues were obtained from TNBC (n = 21) and non-TNBC (n = 20) FFPE sections after pathological evaluation.
The patients' clinical data and tumors' histopathological parameters were retrieved by the HSTR personnel. This information included age at diagnosis, ethnicity (self-reported), size (cm), histology, and tumor grade, and lymph node and distant metastasis involvement.
The mean age at diagnosis was 51.7 ± 11.64 (range, 30-70 years) in the TNBC group and 53.3 ± 12.24 (range, 32-88 years) in the non-TNBC group. Tumors sized ≤2 cm were present in 58.7% of patients in the TNBC group and 64.7% in the non-TNBC group; tumors sized between 2 cm and 5 cm were present in 39.1% and 31.4% of patients, respectively. Tumors >5 cm were observed in only one patient in the TNBC group and two patients in the non-TNBC group. All patients in the TNBC group had invasive ductal carcinomas and tumor grade 2/3. In contrast, patients in the non-TNBC group presented both ductal (96%) and lobular invasive carcinomas (12.3%), and only 39.6% had grade 2/3 tumors. Lymph node metastasis was observed in 40.5% of patients in the TNBC group and 50% of those in the non-TNBC group; distant metastases were observed in 6.9% and 18.2% of patients, respectively. Most patients in both groups were African Americans. Clinical outcomes (alive/deceased) were known only for 24 patients with TNBC. The median follow-up time was 128.4 ± 57.98 months (range, 23.9-252.8 months); 20 patients were alive, and four were deceased. Tumor histology and grade differed significantly between patients with TNBC and those with non-TNBC (Table 1).

Breast Cell Lines
One ER+ breast cancer cell line (MCF-7), six TNBC cell lines (BT549, HCC1806, MDA-MB-157, MDA-MB-231, MDA-MB-453, and MDA-MB-468), and one non-tumorigenic human breast epithelial cell line (MCF-10A) were acquired from the Tissue Culture Shared Resource (TCSR) at the Lombardi Comprehensive Cancer Center, Georgetown. The cell lines were authenticated by STR allelotype analysis, following the International Cell Line Authentication Committee (ICLAC) guidelines [35] (Title: Guide to Human Cell Line Authentication. Available online: https://iclac.org/resources/human-cell-line-authentication/, accessed on 1 June 2018). Cells were cultivated in Roswell Park Memorial Institute Medium (RPMI) 1640 or Dulbecco's Modified Eagle's Medium (DMEM) media supplemented with 10% Fetal Bovine Serum (FBS), 5% Penicillin/Streptomycin (P/S), and 1% Fungizone in T75 culture flasks at 37 • C and in 5% CO 2 conditions. Cells were harvested for the functional assays in a confluence of 70% to 80%. Considering the baseline levels of the HCC1806 and MDA-MB-231 cells, their transfection efficiency using our reverse transfection protocol for both miR-150-5p mimic and inhibitor, and their high metastatic potential, these cells were selected for the functional assays performed.

RNA Isolation
Areas of tumor and adjacent normal tissue (ANT) were microdissected from 10-µm tissue sections of FFPE tumor blocks, after pathology inspection, according to Torresan et al. [36]. The RecoverAll Total Nucleic Acid Isolation kit (Thermo Fisher Scientific, Waltham, MA, USA), was used to isolate RNA and the 2100 Bioanalyzer (Agilent Tech- nologies Inc., Santa Clara, CA, USA), and a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) were used to assess RNA quality and quantity, respectively. RNA was isolated from breast cell lines using the same procedure.

Real-Time Quantitative PCR (RT-qPCR)
The expression levels of miR-150-5p were assessed in clinical samples and in the cell lines by RT-qPCR using a Taqman MicroRNA Assay for hsa-miR-150-5p (Assay#000473, Life Technologies, Carlsbad, CA, USA), following the manufacturer's protocol. RNU48 was used as endogenous control. MiRNA expression levels were calculated using the 2 −∆∆Ct method [37].

Reverse Transfection Assays
The HCC1806 and MDA-MB-231 cells were subjected to the manipulation of miR-150-5p expression levels using mirVanaTM miRNA mimic (miR-150-5p-mimic) and inhibitor (miR-150-5p-inh) assays (Life Technologies, Carlsbad, CA, USA). In all the transfection assays a negative control (NC) was included. The viability of the cells to different concentrations of lipofectamine was assessed prior to the transfection assays. Reverse transfection was performed in triplicate, according to our previous protocol [38]; inhibition and overexpression of miR-150-5p were achieved using a final concentration of 50 nM after 48 h of transfection. The efficiency of transfection was determined by the expression levels of miR-150-5p by RT-qPCR in the TNBC transfected cells in relation to the NC.

Cell Proliferation Assay
Cell proliferation was assessed in the TNBC transfected cells using the Cell Titer Aqueous Solution Cell Proliferation Assay (MTS; Promega, Madison, WI, USA). HCC1806 and MDA-MB-231 transfected cells (4 × 10 3 ) were plated in 96-well culture plates (20 mL/well) and exposed to the Cell Titer Solution. Proliferation rate was assessed at 24, 48, 72, and 96 h using an ELISA reader by measuring optical absorbance at 490 nm. The period of 48 h after transfection was considered Day 0. Experiments were performed in triplicate.

Clonogenic Assay
A clonogenic assay was performed to assess colony formation in the TNBC transfected cells according to Franken et al. (2006) [39] with modifications. HCC1806 cells (1 × 10 3 cells) were trypsinized 48 h after transfection, seeded in 6-well plates and placed at 37 • C for colony formation. The medium was replaced at seven days and after 14 days the cells were fixed with methanol and acetic acid (3:1 ratio) and stained with 0.5% crystal violet (Sigma Aldrich, St. Louis, MO, USA). The colonies (considered of clusters with >50 cells) of triplicate experiments were counted under a light microscope.

Wound Healing Assay
HCC1806 and MDA-MB-231 cells (3 × 10 4 cells/well) were grown in medium with low concentration of fetal bovine serum and in the presence of mitomycin C (MMC) a mitotic inhibitor, as recommended to suppress proliferation in the wound healing assays [40][41][42][43]. The cells were then seeded at 3 × 10 4 cells/well and reverse transfected with miR-150-5p mimic, miR-150-5p inhibitor, and NC directly into the 2-well culture-inserts in µ-Dish 35 mm (Ibidi GmbH, Gräfelfing, Germany), according to our previous protocol [38]. After 48 h of transfection, the inserts were removed, and the cells were washed with Phosphate Buffered Saline (PBS). Migration activity was evaluated by the gap distance between the inserts, measured in the images captured at 0 h and after 24 h and analyzed using ImageJ software [44].

Cytotoxic Assays
HCC1806 transfected cells were plated (3 × 10 4 cells/well) and exposed to 500 nM of doxorubicin. The MTT assay was used per standard protocols. Cell viability was assessed at

The Cancer Genome Atlas Data Processing and Analysis
MiR-150-5p expression was assessed in The Cancer Genome Atlas (TCGA)-Breast Carcinoma (BRCA) dataset (2016_01_28 data version) obtained using the GDAC FireBrowse TCGA portal [45] ("Broad GDAC Firehouse. Available online: http://firebrowse.org/ accessed on 17 February 2022"). Log2-normalized mature miR-150-5p counts were obtained for 755 primary tumors (PT) and 87 normal tissues (NT). Two hundred thirty-eight PT samples had information on ER, PR, and HER2 expression and were classified as TNBC (n = 46) or non-TNBC (n = 192). The Student's t-test was used to compare miR-150-5p expression levels, using p < 0.05 as a threshold of significance. Association with survival was obtained using the aggregated breast cancer clinical studies extracted from TCGA and Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) databases (selected for the TNBC subtype) available at the KM Plotter [46] ("Kaplan Meier Plotter. Available online: http://kmplot.com/analysis/ accessed on 1 February 2022"). The data was displayed considering hazard ratios (HRs), confidence intervals, and log-rank p values.

Statistical Analysis
MiRNA and mRNA expression data were obtained using SDS 2.4 software (Applied Biosystems, Waltham, MA, USA), and the relative expression analysis was performed by the 2 −∆∆Ct method [37]. Data were presented as the mean and standard error of the mean (SEM) from three independent experiments. A Fisher's exact test was used to compare age at diagnosis (≤53, >53), tumor size (≤2 cm, >2 and <5 cm, >5 cm), tumor histology (ductal, lobular), tumor grade (1/2, 3/4), ethnicity (Caucasian, African American), and Cancers 2022, 14, 2156 6 of 16 lymph node and distant metastasis status (positive, negative) between groups. MiR-150-5p expression levels between TNBC, non-TNBC, and ANT groups were compared using the Wilcoxon signed-rank test. The relationship between miR-150-5p levels and clinical and histological parameters was assessed using the Mann-Whitney U test. Survival analysis was performed using Kaplan-Meier and log-rank tests. Statistical analyses were carried out using GraphPad Prism version 8 (San Diego, CA, USA), using p values < 0.05 as statistically significant.

Statistical Analysis
MiRNA and mRNA expression data were obtained using SDS 2.4 software (Applied Biosystems, Waltham, MA, USA), and the relative expression analysis was performed by the 2 −∆∆Ct method [37]. Data were presented as the mean and standard error of the mean (SEM) from three independent experiments. A Fisher's exact test was used to compare age at diagnosis (≤53, >53), tumor size (≤2 cm, >2 and <5 cm, >5 cm), tumor histology (ductal, lobular), tumor grade (1/2, 3/4), ethnicity (Caucasian, African American), and lymph node and distant metastasis status (positive, negative) between groups. MiR-150-5p expression levels between TNBC, non-TNBC, and ANT groups were compared using the Wilcoxon signed-rank test. The relationship between miR-150-5p levels and clinical and histological parameters was assessed using the Mann-Whitney U test. Survival analysis was performed using Kaplan-Meier and log-rank tests. Statistical analyses were carried out using GraphPad Prism version 8 (San Diego, CA, USA), using p values < 0.05 as statistically significant.

MiR-150-5p Expression Levels Are Associated with Tumor Grade, Ethnicity, and Overall Survival in TNBC
The levels of miR-150-5p were significantly associated with tumor grade in the entire cohort (p = 0.0039) and in the TNBC subgroup (p = 0.0453), with higher expression levels in grade 3/4 tumors (Figure 2A). A significant association was also observed with the patients' ethnicity in the entire cohort (p < 0.0001), as well as in TNBC (p = 0.0107) and non-TNBC subgroups (p < 0.0001), with Caucasians presenting higher miR-150-5p expression levels than African Americans ( Figure 2B). No significant associations were observed between miR-150-5p levels and the median age at diagnosis (>53, ≤53), tumor size, tumor histology (ductal or lobular), lymph node metastasis, and distant metastasis.

MiR-150-5p Expression Levels Are Associated with Tumor Grade, Ethnicity, and Overall Survival in TNBC
The levels of miR-150-5p were significantly associated with tumor grade in the entire cohort (p = 0.0039) and in the TNBC subgroup (p = 0.0453), with higher expression levels in grade 3/4 tumors (Figure 2A). A significant association was also observed with the patients' ethnicity in the entire cohort (p < 0.0001), as well as in TNBC (p = 0.0107) and non-TNBC subgroups (p < 0.0001), with Caucasians presenting higher miR-150-5p expression levels than African Americans ( Figure 2B). No significant associations were observed between miR-150-5p levels and the median age at diagnosis (>53, ≤53), tumor size, tumor histology (ductal or lobular), lymph node metastasis, and distant metastasis.

Inhibition of MiR-150-5p Expression Reduces Cell Proliferation and Clonogenic Growth in TNBC Cells
A significant 34.6% decrease in cell proliferation was observed in HCC1806 cells after 72 h and 96 h of transfection with the miR-150-5p inhibitor ( Figure 3A). Consistently, a decrease by 63.75% in the number of colonies was observed with the inhibition of miR-150-5p expression ( Figure 3C). Cell proliferation and clonogenicity of the HCC1806 cells were, however, not affected with the ectopic expression of miR-150-5p ( Figure 3A,C). In the MDA-MB-231 cells, manipulation of miR-150-5p expression did not affect cell proliferation after 72 h or 96 h ( Figure 3B).

Inhibition of MiR-150-5p Expression Reduces Cell Migration and Resistance to Doxorubicin in TNBC Cells
Wound healing assays were conducted to determine the effects of manipulation of miR-150-5p expression on cell migration in HCC1806 and MDA-MB-231 cell lines. MiR-150-5p expression inhibition significantly decreased the percentage of wound closure by

Inhibition of MiR-150-5p Expression Reduces Cell Migration and Resistance to Doxorubicin in TNBC Cells
Wound healing assays were conducted to determine the effects of manipulation of miR-150-5p expression on cell migration in HCC1806 and MDA-MB-231 cell lines. MiR-150-5p expression inhibition significantly decreased the percentage of wound closure by 30.68% and 26.12% in HCC1806 and MDA-MB-231 cells, respectively ( Figure 4A-D). It is interesting to note that in the MDA-MB-231 cells, the changes in this phenotype were not associated with correspondent changes in their cell proliferation rate. Cell migration was not affected by the ectopic expression of miR-150-5p.

Inhibition of MiR-150-5p Expression Reduces the Expression of MYB and Members of the ERK1/2 Pathway
To determine the effects of miR-150-5p manipulation on the expression of direct mRNA targets (Figures 4F and 5A), we investigated the expression levels of MYB, TP53, ZEB1, and SRC proteins in cells transfected with miR-150-5p inhibitor or mimic. Inhibition of miR-150-5p expression led to increased levels of MYB in the MDA-MB-231 cells ( Figure  5B), but not in the HCC1806 cells (Supplementary Figure S2). No alterations were observed in the expression levels of TP53 and ZEB1 in MDA-MB-231 cells ( Figure 5B). The cytotoxic effects of doxorubicin were evaluated in the miR-150-5p transfected TNBC cells (HCC1806). Cell viability was measured 24 h, 48 h, and 72 h after treatment with 50 nM doxorubicin. At all treatment time points, cell viability was significantly lower in cells transfected with miR-150-5p inhibitor than in NC cells ( Figure 4E).

Inhibition of MiR-150-5p Expression Reduces the Expression of MYB and Members of the ERK1/2 Pathway
To determine the effects of miR-150-5p manipulation on the expression of direct mRNA targets ( Figures 4F and 5A) Next, we investigated the expression levels of members of the SRC pathway. In the MDA-MB-231 cells, the manipulation of miR-150-5p expression did not affect the levels of active SRC (p-Src Y416) or inactive SRC (p-Scr Y527). However, the expression levels of p-ERK1/2, which is downstream of the SRC pathway, were reduced upon miR-150-5p inhibition. In addition, increased p-ERK1/2 protein levels were observed with miR-150-5p ectopic expression, suggesting that miR-150-5p plays a role in the activation of the ERK1/2 pathway ( Figure 5B). The same changes in expression levels of inactive (p-Scr Y527) and pERK1/2 were also observed for HCC1806 after manipulation of miR-150-5p expression (Supplementary Figure S2).

Discussion
TNBC is characterized by aggressive phenotypes, poor outcomes, and limited treatment options. MiRNAs act as post-transcriptional gene expression regulators, suppressing or activating the expression of genes implicated in TNBC development and progression [50][51][52][53]. Mounting evidence suggests that miRNAs can be used to differentiate tumor subtypes [54][55], including breast tumors [56][57][58]. We have previously shown that TNBC and non-TNBC tumors from different patient populations exhibit distinct patterns of miRNA expression [59,60]. Notably, miR-150-5p was one of the miRNAs that were highly expressed in TNBC compared with non-TNBC tumors [59]. In this study, we analyzed breast cancer samples from an independent cohort of 113 patients and confirmed that miR-150-5p expression levels were significantly higher in TNBC tissues than in non-TNBC tissues. ROC curve analysis further supported the ability of miR-150-5p levels to discriminate between TNBC and non-TNBC tumors. In addition, we showed that miR-150-5p expression presented a high tumor specificity, as miR-150-5p expression levels were significantly higher in tumor tissues than in ANT tissues. Previous studies yielded conflicting data regarding the expression of miR-150-5p in TNBC. Our results are in agreement with Lu et al. [28], but contradict reports of lower expression of miR-150-5p in TNBC tissues Next, we investigated the expression levels of members of the SRC pathway. In the MDA-MB-231 cells, the manipulation of miR-150-5p expression did not affect the levels of active SRC (p-Src Y416) or inactive SRC (p-Scr Y527). However, the expression levels of p-ERK1/2, which is downstream of the SRC pathway, were reduced upon miR-150-5p inhibition. In addition, increased p-ERK1/2 protein levels were observed with miR-150-5p ectopic expression, suggesting that miR-150-5p plays a role in the activation of the ERK1/2 pathway ( Figure 5B). The same changes in expression levels of inactive (p-Scr Y527) and pERK1/2 were also observed for HCC1806 after manipulation of miR-150-5p expression (Supplementary Figure S2).

Discussion
TNBC is characterized by aggressive phenotypes, poor outcomes, and limited treatment options. MiRNAs act as post-transcriptional gene expression regulators, suppressing or activating the expression of genes implicated in TNBC development and progression [50][51][52][53]. Mounting evidence suggests that miRNAs can be used to differentiate tumor subtypes [54,55], including breast tumors [56][57][58]. We have previously shown that TNBC and non-TNBC tumors from different patient populations exhibit distinct patterns of miRNA expression [59,60]. Notably, miR-150-5p was one of the miRNAs that were highly expressed in TNBC compared with non-TNBC tumors [59]. In this study, we analyzed breast cancer samples from an independent cohort of 113 patients and confirmed that miR-150-5p expression levels were significantly higher in TNBC tissues than in non-TNBC tissues. ROC curve analysis further supported the ability of miR-150-5p levels to discriminate between TNBC and non-TNBC tumors. In addition, we showed that miR-150-5p expression presented a high tumor specificity, as miR-150-5p expression levels were significantly higher in tumor tissues than in ANT tissues. Previous studies yielded conflicting data regarding the expression of miR-150-5p in TNBC. Our results are in agreement with Lu et al. [28], but contradict reports of lower expression of miR-150-5p in TNBC tissues compared with non-tumor tissues [27,29]. These opposite expression patterns may reflect the ability of miR-150-5p to regulate multiple gene targets, some of which have opposing oncogenic or tumor-suppressive functions [61]. The combined suppression or activation of other miR-NAs regulating the same targets or functional pathways can result in a more pronounced effect of a single miRNA with either an oncogenic or tumor-suppressive function in a given type of tumor [62].
Analysis of the relationship between miR-150-5p expression levels and the clinical and histopathological characteristics of patients revealed that miR-150-5p levels were significantly associated with grade 3/4 tumors (TNBC and non-TNBC combined and TNBC subtype only). In contrast to the findings of other studies, no significant association was observed between miR-150-5p expression and other poor prognostic factors, such as lymph node metastasis and distant metastasis [29,63]. Interestingly, survival analysis in the TNBC subgroup showed that miR-150-5p upregulation was associated with prolonged overall survival. These findings were confirmed in patients with TNBC from the TCGA and METABRIC databases and highlight the inconsistent implications of miR-150-5p expression in prognosis, as reported by Wang et al. [14].
MiR-150-5p expression levels were also significantly associated with ethnicity in the entire cohort and in the TNBC and non-TNBC subgroups, with Caucasians showing higher miR-150-5p levels than African Americans. In contrast, our previous global miRNA profiling study [59] showed that African American patients with TNBC presented higher expression levels of miR-150-5p than non-Hispanic White patients. It is worth pointing out that in this study, ethnicity was self-reported and not determined by ancestral genomic profiling; therefore, the ethnicity of the patients may have affected miRNA expression in a way that could not be determined [64,65]. These results support the hypothesis that biological differences exist in miRNA expression patterns between populations and that these differences may influence prognosis and treatment outcomes.
The expression levels of miR-150-5p also varied among cell lines. Most TNBC cell lines had higher miR-150-5p levels than non-tumorigenic breast cells (MCF-10A) and non-TNBC cells (MCF-7). MiR-150-5p expression was manipulated in the highly metastatic HCC1806 and MDA-MB-231 cells [66,67], to evaluate its potential role in modulating TNBC aggressive phenotypes. MiR-150-5p expression inhibition led to a decrease of the TNBC cells tumorigenicity, as significant decreases in cell proliferation, capacity of colonies formation, cell migration, and resistance to doxorubicin were observed. However, ectopic expression of miR-150-5p did not affect these phenotypes. These results suggest that additional levels of miR-150-5p in TNBC cell lines are not required to confer further aggressiveness. Furthermore, a given miRNA, such as miR-150-5p, very likely does not entirely control a given tumorigenic phenotype, requiring the cooperation of other miRNAs that act cooperatively in the same biological process and/or signaling pathway they are involved, by regulating common mRNA targets [62]. Several miRNAs from common clusters or families share similar functions with miR-150-5p and regulate the same mRNA targets. For example, miR-769-5p is also mapped at 19q13 and regulates epithelial-mesenchymal transition (EMT), invasion, and metastasis [68][69][70]. In addition, long non-coding RNAs (lncRNAs) can interact with miRNAs, allowing them to exert their full cellular action. Several studies have demonstrated a relationship between lncRNAs and miR-150-5p in breast tumorigenesis; these lncRNAs include FOXD2 AS1 [31], RNA 01121 [25], and MAFG AS1 [26].
It is also of interest to note that the miR-150-5p impact on tumor aggressiveness was distinct in the patients' TNBC tissue samples and in the TNBC cell line models. In the patients' tumor tissues, high miR-150-5p expression levels was associated with higher survival rates, which as mentioned above, was also observed in the TNBC patients from the TCGA and METABRIC databases. At the other hand, in the TNBC cell lines, miR-150-5p behaved most compatible with an oncogenic mode of action, where its inhibition was associated with a decrease in tumor aggressiveness. Several reasons can account for these discrepancies, including the fact that the TNBC clinical samples were derived from patients' primary tumors and the TNBC cell lines from highly metastatic tumors. However, even within metastatic samples and metastatic cell lines, genotypic and phenotypic discrepancies can occur. In a study conducted by Liu et al. (2019) [71], a comparative analysis of publicly available genomic and transcriptomic signatures of clinical metastatic breast tumors and metastatic breast cell lines models demonstrated substantial differences among these intrinsic signatures, which can affect distinct cancer associated genes' networks and also affect the classical tumor phenotypes.
Cell adhesion and adherent junctions were the process mostly observed for miR-150-5p functional enrichment analysis. These are critical processes of EMT, which is essential for cell migration and invasion [72,73]. EMT is regulated by several transcription factors and proteins, including E-cadherin, SLUG, Vimentin, SNAIL, ZEB1, and Twist [74]. Alterations in these proteins directly affect the cells tumorigenic behavior by promoting alterations in the cytoskeleton, adhesion, and cell polarity [75]. In this study we evaluated the role of miR-150-5p in the expression of the EMT markers SNAIL, SLUG, Vimentin, and ZEB1. We observed that the expression levels of the early EMT markers SNAIL and SLUG were reduced with the inhibition of miR-150-5p expression. Additionally, significantly higher levels of the oncoprotein MYB (a direct miR-150-5p target) was observed with the inhibition of miR-150-5p expression. Ectopic expression of miR-150-5p, however, did not change MYB expression levels, which were similar to those observed in the NC cells.
Previous studies on lung, gastric, cervical, and breast cancers have shown that miR-150-5p promotes cell proliferation and migration by targeting the SRC kinase signaling inhibitor 1 (SRCIN1), a negative regulator of SRC and an experimentally validated direct target of miR-150-5p [23,28,76]. SRCIN1 inhibits SRC via phosphorylation (tyrosine 597) and the subsequent activation of CSK. Interestingly, a relationship between pERK and SNAIL has been previously reported in breast cancer, with high levels of pERK and SNAIL being associated with increased migration [77]. Therefore, we hypothesized that inhibition of miR-150-5p, which reduced the expression levels of SNAIL in TNBC cells, would negatively affect the transcriptional activity of SRCIN1 and other proteins of the SRC pathway, such as pERK and pFak. These proteins are regulated by the phosphorylation of SRC and promote cell migration and invasion. Lower expression levels of pERK1/2 were observed with the inhibition of miR-150-5p; conversely, ectopic expression of miR-150-5p led to an increase in pERK1/2 levels. In contrast, the levels of active (p-Src Y416) or inactive (p-Scr Y527) SRC or pFak were not affected with the inhibition and ectopic expression of miR-150-5p. Manipulation of the expression levels of these proteins would be necessary to show the reverse activity of miR-150-5p function in these cells and the reversal of the tumorigenic phenotypes.

Conclusions
Our findings suggest that miR-150-5p is involved in TNBC aggressiveness, presenting differential expression in these tumors compared with non-tumor and non-TNBC tissues samples. Our results also indicate that miR-150-5p regulates cell proliferation, clonogenic growth, drug resistance, and migration as observed in the two TNBC cell lines evaluated, HCC1806 and MDA-MB-231. Although several mechanisms can underlie the tumorigenic functions of miR-150-5p in these cells, changes in the expression of its target genes MYB, ZEB1, and members of the SRC pathway can contribute to its oncogenic role in TNBC.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14092156/s1, Figure S1: Expression levels of miR-150-5p by RT-qPCR in seven breast cancer cell lines when compared to the MCF-10A; Figure S2. Western blotting analysis and quantitative readings of MYB, p-Scr Y527, and pERK1/2 expression levels for HCC1806 after manipulation of miR-150-5p expression. Funding: This research was partially funded by the Georgetown University Center of Excellence in Regulatory Science and Innovation (CERSI; U01FD004319), a collaborative effort between the University and the U.S. Food and Drug Administration to promote regulatory science through innovative research and education. This research does not necessarily reflect the views of the FDA.

Institutional Review Board Statement:
The study was conducted in accordance with the guidelines of the Declaration of Helsinki and was approved by the Institutional Review Board (or Ethics Committee) of Georgetown University Medical Center (IRB#1992-048) to the HSTR.
Informed Consent Statement: Patient consent was waived, considering that no patient identifiers were provided to the authors.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: https://app.box.com/s/h5bc82qa1t25xe3z6cstx876mexy4u2j, accessed on 10 March 2022.


Comprehensive Analysis of N6-Methyladenosine-Related Long Noncoding RNA Prognosis of Acute Myeloid Leukemia and Immune Cell Infiltration

N6-Methyladenosine-related long noncoding RNAs play an essential role in many cancers’ development. However, the relationship between m6A-related lncRNAs and acute myelogenous leukemia (AML) prognosis remains unclear. We systematically analyzed the association of m6A-related lncRNAs with the prognosis and tumor immune microenvironment (TME) features using the therapeutically applicable research to generate effective treatment (TARGET) database. We screened 315 lncRNAs associated with AML prognosis and identified nine key lncRNAs associated with m6A by the LASSO Cox analysis. A model was established based on these nine lncRNAs and the predictive power was explored in The Cancer Genome Atlas (TCGA) database. The areas under the ROC curve of TARGET and TCGA databases for ROC at 1, 3, and 5 years are 0.701, 0.704, and 0.696, and 0.587, 0.639, and 0.685, respectively. The nomogram and decision curve analysis (DCA) showed that the risk score was more accurate than other clinical indicators in evaluating patients’ prognoses. The clusters with a better prognosis enrich the AML pathways and immune-related pathways. We also found a close correlation between prognostic m6A-related lncRNAs and tumor immune cell infiltration. LAG3 expression at the immune checkpoint was lower in the worse prognostic cluster. In conclusion, m6A-related lncRNAs partly affected AML prognosis by remodeling the TME and affecting the anticarcinogenic ability of immune checkpoints, especially LAG3 inhibitors. The prognostic model constructed with nine key m6A-related lncRNAs can provide a method to assess the prognosis of AML patients in both adults and children.

INTRODUCTION
Acute myelogenous leukemia (AML) is a hematopoietic malignancy characterized by numerous cytogenetic and molecular aberrations. It accounts for 15%-20% of childhood leukemia, approximately 33% of adolescents and young adults, and 80% of adult leukemia (Creutzig et al., 2008;Molica et al., 2019). Despite the considerable improvements achieved in intensified treatment strategy, supportive care, and riskadapted patient stratification, the overall survival (OS) does not exceed 70% and relapse rates range between 25% and 35% (Lonetti et al., 2019). AML represents a clinical challenge because of its poor prognosis, highlighting an urgent need to explore novel therapeutic targets and new biomarkers for diagnosis and prognosis (Jiang et al., 2021).
The discovery of epigenetic regulation has opened a new realm of gene regulation in eukaryotes. There is increasing evidence showing that N6-methyladenosine (m6A), the most prevalent RNA modification, plays a critical role in RNA regulation (Deng et al., 2018). The reversibility and dynamics features of RNA m6A modification indicate its unique contribution to tumorigenesis (Dai et al., 2018). Several studies have reported that m6A modification is involved in the occurrence and development of AML (citation). The overexpression of m6A methyltransferase METTL14 and METTL3 in AML cells promoted the self-renewal of leukemic stem cells and the development and maintenance of AML (Barbieri et al., 2017;Vu et al., 2017;Weng et al., 2018). WTAP, an oncogenic protein in AML, functions as a regulatory subunit of the m6A methyltransferase complex (Bansal et al., 2014). The m6A demethylase FTO also plays an oncogenic role in AML (Li et al., 2017).
LncRNAs, a kind of noncoding RNAs with approximately 200 nt-100 kb nucleotides in length, regulate the translation, shearing, and degradation of target mRNAs (Silva et al., 2010). LncRNAs also serve as the preferred biomarkers in the diagnosis, prognosis, and therapeutic approach to various types of diseases (Silva et al., 2010). Previous studies have reported that lncRNA dysregulation is associated with the prognosis of AML patients (Sun et al., 2018;Feng et al., 2020). However, there is still a lack of systematic evaluation of the prognostic value of lncRNAs in AML in a large sample. It is known that m6A modifications play an essential role in the dysregulation of lncRNAs in AML. However, the modifying regulatory role of m6A regulators on lncRNAs and the mechanism of association between m6A-related lncRNAs and AML prognosis are still far from being elucidated, and an indepth understanding of the mechanisms of m6A modifications of lncRNAs in AML progression may help to identify relevant prognostic biomarkers . In this study, we analyzed the m6A-related lncRNA expression of AML patients in the therapeutically applicable research to generate effective treatment (TARGET) database and explored the relationship between m6A-related lncRNAs and AML prognosis. Furthermore, we constructed the prognostic model using the key m6A-related lncRNAs identified in this process and evaluated the predictive power of the model in The Cancer Genome Atlas (TCGA) database.
Progression of AML is highly correlated with the physiological state of the tumor microenvironment (TME) (Roma-Rodrigues et al., 2019). It was reported that the immunosuppressive environment favors the immune escape of AML cells (Nahas et al., 2019). Immune checkpoint proteins are highly relevant to the initiation of immunocyte signaling pathways, which could be manipulated by tumor cells to escape immune response and form a TME that is beneficial to neoplastic development (Topalian et al., 2016;Osipov et al., 2019;Bonavita et al., 2020). Here, we found that m6A-related lncRNAs might affect AML prognosis through the immunocyte signaling pathway. Thus, we further explored the association between m6A-related lncRNAs and immune checkpoint proteins, tried to explain the mechanism of m6A-related lncRNAs on the TME heterogeneity, and searched for biomarkers that could serve as potential immunotherapy targets.

Acquisition of Gene Expression and Clinical Data
The RNA sequencing data and the associated clinical data from 358 AML patient samples were downloaded from the TARGET database (https://ocg.cancer.gov/programs/target), the download date was March 20, 2021. After excluding 63 samples without survival information, 295 AML patient samples data were used for screening m6A-related lncRNAs, analyzing the relationship between AML prognosis and immune cell infiltration, and the construction of a prognostic model.
In addition, we downloaded RNA sequencing data of 151 AML patient samples and the clinical data of 200 AML patient samples from the TCGA (The Cancer Genome Atlas) database (https://www.cancer.gov/about-nci/organization/ccg/research/ structural-genomics/tcga) to verify the predictive effect of the prognostic model. The download date was March 20, 2021.

Annotation and
Ribonucleic acid sequencing data were merged into an RNA matrix file using the programming language Perl (http://www.perl. org/). LncRNAs were identified by recognizing the ensemble IDs of the genes, including lincRNA, antisense, processed transcript, sense intronic, 3 prime overlapping ncRNAs, and sense overlapping. We also extracted the expression matrix of 23 m6A regulators from AML patients from the TARGET database for the following analysis.
The Pearson correlation analysis was implemented to identify m6A-related lncRNAs (with the |Pearson R| > 0.4 and p < 0.001). Taking the overall survival (OS) as the prognostic outcome endpoints, we used the univariate Cox regression analysis to screen m6A-related lncRNAs associated with AML prognosis.

Subgroup Identification Based on Consensus Clustering and GSEA Between Different Clusters
"ConsensusClusterPlus" and "limma" packages (Wilkerson and Hayes, 2010;Ritchie et al., 2015) in R software were used to classify AML patients into subtypes (cluster 1, cluster 2, and so on) based on the expression patterns of lncRNAs screened by the univariate Cox regression analysis. Kaplan-Meier curves and the logrank test were used to compare the OS between the different clusters to determine the differential expression of prognosis-related lncRNAs. Genome Set Enrichment Analysis (GSEA) was applied to the different clusters. Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets and phenotype tags with high and low expression files were loaded into the GSEA (v4.0.3; Broad Cambridge University Institute, MA, https://www.gsea-msigdb.org/gsea/index.jsp) software.

Role of Immune Cell Infiltration and the Tumor Microenvironment
After calculating the stromal score and immune scores by the ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) (version 2.15.3) algorithm (https://sourceforge.net/projects/estimateproject/), we obtained the ESTIMATE score by combining the two scores. Our final index was to calculate the tumor purity based on the ESTIMATE score (Yoshihara et al., 2013). According to the study conducted by Aran et al., in 2015, a tumor purity over the minimum threshold of 60% means that the clustering results are reliable (Aran et al., 2015). The Cell-type Identification by Estimating Relative Subpopulations of RNA Transcripts (CIBERSORT) method (Newman et al., 2015) was applied to assess the proportion of 22 immune cell subtypes in AML patient samples to explore the differences in immune cell subtypes among three clusters. Similarly, we compared the differences between five immune checkpoints of three different subtypes. The expression of immune checkpoints is closely related to immunotherapy. Five immune checkpoints were derived from previous studies, including programmed death 1 (PD-1) (Sharpe and Pauken, 2018) and its ligand 1 (PD-L1) (Daassi et al., 2020), cytotoxic T-lymphocyte antigen 4 (CTLA-4) (Agdashian et al., 2019), mucin domain-containing molecule-3 (TIM-3) (Wolf et al., 2020), and lymphocyte-activation gene 3 (LAG3).

Construction and Validation of Prognostic Models
LASSO Cox regression was used to further identify biomarkers associated with AML prognosis (Friedman et al., 2010). Finally, nine m6A-related AML prognosis-related lncRNAs were obtained and the risk score for the prognosis of each AML patient was calculated based on the expression of the nine lncRNAs. The calculating formula is: where coef i means the coefficients and x i is the FPKM value of each m6A-related lncRNA.
AML patients in the TARGET database were divided into high-risk and low-risk groups based on the median risk score. Kaplan-Meier curves and log-rank tests were used to compare the prognostic difference between the low-risk group and the highrisk group. The time-dependent receiver operating characteristic (ROC) curves and the area under curve (AUC) were measured by the package "survivalROC" in R software, which was used to evaluate the prognostic prediction accuracy of the model and 1/3/ 5-year OS (Heagerty and Zheng, 2005). The 1/3/5-year OS of the model was also calculated in the TCGA database to assess the model's predictive ability in adult AML patients.

Constructing a Predictive Nomogram and Decision Curve
A nomogram is widely used to predict the prognosis of cancer. All independent prognostic factors including gender, age, race, FAB category, WBC at diagnosis (WBC), bone marrow leukemic blast percentage (%) (BM), and peripheral blasts (%) (PB) were used to build a nomogram and calculate a nomogram score to investigate the probability of 1, 3, and 5 years OS (overall survival) of AML. Finally, calibration curves were plotted to estimate the calibration capability of the nomogram and a decision curve analysis was used to evaluate the clinical usefulness of the nomogram.

Statistical Analysis
Most analyses were performed with R software (version 4.0.5, http://www.R-project.org). Unless otherwise noted, p < 0.05 was considered statistically significant.

Characteristics of Participants in This Study
Tables 1, 2 show the characteristics of the participants from the TARGET and TCGA databases, respectively. Figure 1 is the flowchart of this study. From a total of 11,535 annotated lncRNAs and the expression matrix of 23 m6Aassociated genes, we identified 1,442 lncRNAs as significant m6A-associated genes and found 315 lncRNAs significantly associated with AML prognosis (p < 0.05), using the univariate Cox regression analysis (Supplementary Table S1).

Subgroup Identification Based on Consensus Clustering and GSEA Between Different Clusters
We performed consensus clustering on 295 AML samples and divided the samples into three subgroups based on the maximum AUC increment of CDF and the expression pattern of m6Arelated lncRNAs (Figures 2A-C). The Kaplan-Meier method and the log-rank test were used to compare the prognosis of the three clusters and the result showed that there are statistical differences in the OS ( Figure 2D) (p = 0.007). The survival curves of cluster 1 and cluster 2 are almost the same, but there are significant differences between cluster 3 and cluster 1 (p = 0.003) or cluster 2 (p = 0.006). Thus, cluster 1 and cluster 2 were considered better prognostic groups, whereas cluster 3 was identified as the poor prognostic group.
The chi-square test was used between clusters and other clinical parameters, including gender, age, Race, FAB category, WBC at diagnosis (WBC), bone marrow leukemic blast The interquartile spacing of age is 10.39 (3.74, 14.64). The minimum value of age is 0.03 and the maximum value of age is 23.51. 7 variables were taken as covariates, including gender, age at diagnosis in days, race, FAB category, ethnicity, WBC at diagnosis, bone marrow leukemic blast percentage (%), and peripheral blasts (%). "Gender" is taken as a categorical variable; "Race" is taken as a categorical variable; "Age" is taken as a categorical variable; "FAB category" is taken as a categorical variable; "WBC at Diagnosis" is taken as a categorical variable; "Bone marrow leukemic blast percentage (%)" is taken as a categorical variable; and "Peripheral blasts (%)" is taken as a categorical variable. The median age as well as the interquartile spacing is 57.50 (44.75, 67.00). The minimum value of age is 18 and the maximum value of age is 88. Five variables were taken as covariates, including gender, age, race, FAB category, ethnicity. "Gender" was taken as categorical variable; "Race" was taken as categorical variable; "Age" was taken as categorical variable; "FAB category" was taken as categorical variable; and "Ethnicity" was taken as categorical variable.
Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 888173 percentage (%) (BM), peripheral blasts (%) (PB), and the cluster to test for statistically significant differences in the composition ratios between populations. The results were plotted in a heatmap ( Figure 2E). Red represents high lncRNA expression and blue represents low lncRNA expression. The horizontal axis represents 295 AML patients and the vertical axis represents 315 lncRNAs. We found that expressions of lncRNAs were low in all of cluster 1 samples and most of cluster 2, and expressions were high in a small part of cluster 2 and all of cluster 3. Significant differences in the FAB categories of AML were also found between the three clusters. Comparison between three different clustering-related clinical factors, the result of the FAB categories in the installment in significant differences between the three clusters, and other factors did not show statistical differences. The results indicated that high m6A-related lncRNA expression is associated with poor AML prognosis. The GSEA result between clusters 1 and 3 showed that the B-cell receptor signaling pathway, acute myeloid leukemia pathway, and T-cell receptor signaling pathway were enriched in cluster 1 (Figures 3A-C). The GSEA result between clusters 2 and 3 is similar to that between clusters 1 and 3. Because of that, most of those signaling pathways were involved in the regulation of immune checkpoint expression, leukemia-related pathways, and cytokine-related regulation pathways. Therefore, we further analyzed immunity, including immune checkpoint expression, tumor immune cell abundance profiles, and TME scores to investigate the differences between the three clusters ( Figures 3D-F).

Analysis of the Immune Status
The observation of enriched immune-related pathways between three clusters drove us to further explore the relationship between m6A-related lncRNAs and cancer immunity. The ESTIMATE algorithm was used to calculate immune scores and stromal scores. The estimate scores and tumor purity were calculated by immune scores and stromal scores for all AML samples. The results showed that estimate scores and tumor purity were not statistically significant between the three clusters whereas the tumor purity of all three clusters was greater than 60% ( Figures 4G,H). In addition to that, the CIBERSORT algorithm was used to calculate the contents of 22 different immune cell types between three clusters. The contents of B cells naive, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells follicular helper, NK cells resting, NK cells activated, macrophage M0, macrophage M2, dendritic cells activated, and eosinophils were significantly different between the three clusters. Monocytes, mast cells activated, eosinophils, T cells CD4 memory resting, and B cells naive accounted for a large proportion of immune infiltrating cells. B cells naive, T cells CD8, NK cells activated, and macrophages M2 were higher in the population with a better prognosis than in those with a poor prognosis. The levels of T cells CD4 naive, NK cells resting, and eosinophils were lower in those with excellent prognosis than in those with poor prognosis ( Figure 4A). In the comparison between cluster 1 and cluster 3, and cluster 2 and cluster 3, we found statistically significant differences in B cells naive, T cells CD4 naive, T cells CD4 memory resting, and macrophages M2 (Figures 4B,C).
Immunotherapy targeting immune checkpoints holds great promise for the clinical treatment of human cancers (Wan et al., 2021). Therefore, we investigated the expression differences of several genes encoding well-known immune checkpoint proteins between three clusters. The results showed that the expression of the LAG3 gene was statistically different among the three clusters (p < 0.05) ( Figure 4D): it was high in cluster 1 and cluster 2, which had a better prognosis (p < 0.05) ( Figures 4E,F). Lymphocyte-activation gene-3 (LAG3, also named CD223) is a cell surface molecule expressed on activated T cells. It plays a key role in the immune checkpoint, which makes it a potential cancer immunotherapeutic target. The results suggest that patients in clusters 1 and 2 may be more sensitive to the immunotherapy targeting LAG3.

Construction and Validation of the Prognostic Model
The LASSO Cox regression analysis revealed that nine out of the 315 lncRNAs were with closely associated with AML prognosis ( Figures 5A,B), which were LINC00852 (p < 0.001, coefficient = 0.0256), AL157392.3 (p < 0.001, coefficient = 0.0278), AC127459.1 (p < 0.001, coefficient = 0.0946), AC106820.3 (p < 0.001, coefficient = 0.0024), AC092757.2 (p = 0.004, coefficient = 1.6278), AC124248.1 (p < 0.001, coefficient = 0.0054), AC145207.5 (p < 0.001, coefficient = 0.0008), DNAAF4-CCPG1 (p < 0.001, coefficient = 0.0070), and AC023908.3 (p < 0.001, coefficient = 0.1977) with all HR greater than 1, as shown in Table 3. We included the risk score, age, race, gender, WBC at diagnosis, bone marrow leukemic blast percentage (%), peripheral blasts (%), and FAB categories in the univariate and multivariate Cox regression analyses to find correlates that might affect the prognosis of AML in children. The results of univariate and multivariate Cox regression analyses showed that only the risk score was an associated factor affecting the prognosis of AML in children (p < 0.05) ( Figures 5E,F). Patients were divided into high-risk and low-risk groups based on the median risk score. Kaplan-Meier survival curves showed that low-risk patients in the TARGET database had better OS than high-risk patients ( Figure 5C). We analyzed whether clinically relevant factors differed between the high-risk and low-risk groups. We observed significant differences between the high-risk and low-risk groups of patients regarding the FAB categories (p < 0.05) and clusters (p < 0.05) ( Figure 6A). We also compared differences in risk scores between different clinical subgroups. The results showed statistically significant differences in risk scores between clusters 1 and 3, clusters 2 and 3, WBC at diagnosis, bone marrow leukemic blast percentage (%), and FAB categories ( Figures 6B-D) (Supplementary Figures S1A-D). Next, we evaluated the risk score of each AML case in the TARGET database and results showed that AML patients in the lowrisk group have a better survival status and a shorter mortality status than the high-risk group (Figures 7A,C). A time-dependent ROC analysis indicated that the AUC of the risk score predicted OS at 1, 3, and 5 years were 0.701, 0.704, and 0.696, respectively ( Figure 7E).
The prognostic model was then used in the TCGA database to evaluate the predictive power for survival in adult AML patients. Patients were also divided into high-risk and low-risk groups based on the median risk score. The Kaplan-Meier survival curves showed that low-risk patients in the TCGA database had better OS than high-risk patients ( Figure 5D). We also found that AML patients in the low-risk group had better survival status and shorter mortality than the high-risk group (Figures 7B,D). A time-dependent ROC analysis showed that the AUCs of the risk score for predicting OS at 1, 3, and 5 years were 0.587, 0.639, and 0.685, respectively ( Figure 7F). Overall, the aforementioned results suggest that the prognostic model had a good predictive power in predicting the survival time of pediatric patients with AML and adult patients.

Constructing a Predictive Nomogram and Decision Curve
Nomograms and decision curves of 1, 3, and 5 years are shown in Figures 8A-E. Based on the clinical information and risk score of each patient, the probability of survival in 1, 3, and 5 years can be calculated through the nomogram. Risk scores and FAB categories were identified as significant factors for prognosis in patients with AML on the nomogram. The calibration curve is shown in Figure 8B. The 1-, 3-, and 5-year decision curves showed that the nomogram score was the optimal predictor of survival for pediatric AML patients, followed by the risk score.

DISCUSSION
Previous studies have found that m6A-regulated lncRNAs are involved in the biological process of cancer progression. Numerous studies have shown that m6A could modify lncRNAs, which contribute to tumorigenesis in multiple cancer types, including proliferation, invasion, and metastasis (Chen et al., 2020). m6A-related lncRNAs may serve as a more promising prognostic marker in cancer. However, the potential role of m6A-related lncRNAs in the prognosis of AML remains unclear.
In this study, we downloaded the expression data of 295 AML patients from the TARGET database and identified 315 m6Arelated prognostic lncRNAs. Three prognostic sample clusters were identified according to the expression pattern of m6A-related  Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 888173 lncRNAs and cluster 3 was associated with poor prognosis. GSEA results showed that several vital immune-related signaling pathways, such as the B-cell receptor signaling pathway and T-cell receptor signaling pathway, were enriched in better prognostic clusters (clusters 1 and 2). The content of 11 out of the 22 immune cells was significantly different between the better prognostic clusters (cluster 1 or cluster 2) and the poorer prognostic cluster (cluster 3). Patients with good prognosis had lower levels of T cells CD4 naive, NK cells quiescent, and eosinophils than poor prognosis patients. In all populations, monocytes were much more abundant than other cells. The large proportion of acute myelomonocytic leukemia (M4) and acute monocytic leukemia (M5) in the studied population may have contributed to this outcome. In addition, we found that there is not statistical difference between the three clusters for immune scores, estimate scores, stromal scores, and tumor purity. When measuring the effectiveness of predictive checkpoint inhibitor drugs, it is only accurate if the extent of immune cell infiltration into the tumor is clearly quantified. When this aspect of tumor purity is unclear, estimates of the success of immunotherapy may be too high or too low. TCGA initially set a mass threshold that the tumor samples included in the cohort consisted of at least 80% tumor cell nuclei, as determined by a visual analysis. However, this threshold was later lowered to 60%. Low tumor purity can have a devastating effect on the results of the cluster analysis (Aran et al., 2015). Our calculations showed that the tumor purity between all three clusters was greater than 60% while there was no statistical difference, hence, we obtained that the confounding effect between the subtypes was small. Different immune cell types and contents can affect the development and prognosis of AML. Jiang's study predicted the prognosis of pediatric AML patients by immune checkpoint, and the degree of immune cell infiltration calculated in the study was similar to that in our study. (Jiang et al., 2021). Tumor-infiltrating immune cells are highly predictive of tumor progression and patient survival. Therefore, the immune cell infiltration results we obtained may provide some theoretical basis for precise immune cell therapy for patients among different clusters. Interestingly, we found that the population with a better prognosis exhibited relatively high expression of the immune checkpoint protein LAG3 and was more sensitive to immunotherapy. Existing studies on immune-related treatments for AML are also increasing. A study showed that tumorinfiltrating CD8 + cytotoxic T-cells in AML had upregulated inhibitory receptors such as programmed cell death 1 (PD-1), cytotoxic T-lymphocyte antigen 4 (CTLA-4), T-cell immunoglobulin and mucin domain-containing protein 3 (TIM-3), and lymphocyte-activation gene 3 (LAG3) (Ozkazanc et al., 2016). We also found statistically significant differences in the numbers of infiltrated cytotoxic T-cells in the three clusters. There are more infiltrated cytotoxic T-cells in cluster 2 than that in cluster 3. This may explain the higher expression of CTLA4 and LAG3 in cluster 2. It was reported that immune checkpoint inhibitors targeting PD1, PD-L1, and CTLA4 are clinically effective in AML. However, there are not clinical reports on LAG3 inhibitor therapy on AML (Dama et al., 2019;Daver et al., 2019). Our study may shed some new perspectives on the study of LAG3 as a new immune checkpoint inhibitor for AML.
The LASSO Cox regression method confirmed that nine out of 315 m6A-related lncRNAs are valuable in constructing a prognostic model for predicting OS in AML patients. AML patients were divided into the low-risk and high-risk subgroups based on the median risk score, and the high-risk group had poor clinical outcomes. The nine m6A-related lncRNAs are associated with AML prognosis. LINC00852, identified with a role in promoting AML cell proliferation, was highly expressed in the AML high-risk group. Several studies reported that LINC00852 plays a role in the proliferative and aggressive nature of osteosarcoma associated with the receptor complex kinase AXL, a founding member of the TAM receptor complex kinase family and that dysregulation of AXL by chemotherapy may be resistant to Gas6 stimulation-induced acute myeloid leukemia in AML cells . Hong et al. found that AXL-mRNA expression was upregulated in relapsed drug-resistant AML specimens (Hong et al., 2008). In the U937 cell line, chemotherapy drugs (doxorubicin, etoposide, and cisplatin) induced the dose-dependent expression of AXL-mRNA by increasing the methylation level of the AXL promoter CCWGG. Ben-Batalla et al. found that AXL-mRNA was expressed in 57% (64/112) of newly diagnosed middle-risk AML cells with normal genetic karyotype and its expression level was an independent prognostic factor for the overall survival of patients (Ben-Batalla et al., 2013). AXL was closely related to drug resistance and the poor prognosis of AML. Targeting AXL could block the activation of the Gas6/AXL pathway, inhibit the proliferation of AML cells, and overcome the resistance of FLT3-ITD + AML cells to FLT3 inhibitors. In addition, our analysis showed that METTL14, WTAP, YTHDC2, and FMR1 could be positively regulated by LINC00852. Evidence from studies suggested that METTL14 regulates myelopoiesis and leukemogenesis through the SPI1-MET-TL14-MYB/MYC signaling axis, and the myeloid transcriptional regulator SPI1 negatively regulates METTL14 to induce the m6A modification of oncogenes MYB and MYC, thereby inhibiting AML cell differentiation and promoting cell self-renewal (Weng et al., 2018). Therefore, LINC00852 is a promising prognostic biomarker in AML patients. AL157392.3 and AC145207.5 were also found to be highly expressed in our study in the AML high-risk group. AML has recently been shown to be regulated by glycolytic regulators, thereby promoting leukemogenesis (Chen et al., 2014). Glycolysis-related genes and transcription factors can mediate the development of diseases in specific AML subtypes (Federzoni et al., 2014;Robinson et al., 2020). AL157392.3 was found to be significantly correlated with glycolysis in bladder cancer, lowgrade glioma, mesothelioma, pancreatic cancer, and uveal melanoma (Ho et al., 2021). However, the mechanism of AL157392.3 in the glucose metabolic pathway affecting the prognosis of AML patients has not been investigated. Therefore, we suggest that AL157392.3 may influence the development of AML through the glycolytic pathway. Moreover, AC145207.5 was only found to play a role in developing hepatocellular carcinoma, but not in AML (Yu and Zhu, 2021;Zhou et al., 2021). The other six lncRNAs were rarely reported in cancers and further research should be conducted to reveal their potential role in AML development. Interestingly, seven lncRNAs were screened in a study to identify m6A-related lncRNAs for prognosis in adult AML (Li et al., 2021). Since the population we used was pediatric AML patients, this may have led to different results from ours. The subsequent analysis in the TCGA adult dataset showed that the prognostic model we constructed in pediatric AML patients also has predictive power in adult AML patients.
Nine lncRNAs were used to construct the AML prognostic model. The areas under the ROC curves at 1, 3, and 5 years in the TARGET database were 0.701, 0.704, and 0.696, respectively. We also verified the model with the data from the TCGA database. The prognostic model has a strong predictive ability in pediatric AML patients. Further analysis of the TCGA data showed that larger the area under the ROC curve with increasing time, the stronger the predictive power of the prognostic model for adult AML patients. The area under the ROC curve at five years was close to 0.7. Similarly, in a study predicting the prognosis of AML in children, the authors constructed a predictive model from the TARGET database and validated it in the GEO of adult patients (Jiang et al., 2021). In this study, heterogeneity of the lncRNAs involved in this model between adults and children has not been reported. In addition, we incorporated clinically relevant factors into the column plots and decision curves and found that nomogram scores and risk scores had a superior predictive power for the survival of pediatric AML patients than any other clinical factor.
Our study has some limitations. First, AML is a complex disease caused by multiple factors. Because of the excessive absence of some factors, for example, cytogenetic indicators were not evaluated in our study. Second, the description of the potential mechanisms of cancer genes was based on experimental evidence only for LINC00852. The study of the other eight lncRNAs needs to be further confirmed in future by in vitro and in vivo experiments. Then, the processing of blood specimens and the collection of population data will vary among different databases, inevitably introducing errors that may affect the analysis results. Finally, our experimental data were obtained from publicly available databases and pediatric databases lack relevant cohort validation, and our study should next include more prospective pediatric cohort studies.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.


Individual and combined effects of the GSTM1, GSTT1, and GSTP1 polymorphisms on type 2 diabetes mellitus risk: A systematic review and meta-analysis

Backgrounds: Compared with previously published meta-analyses, this is the first study to investigate the combined effects of glutathione-S-transferase polymorphisms (GSTM1, GSTT1 and GSTP1 IIe105Val) and type 2 diabetes mellitus (T2DM) risk; moreover, the credibility of statistically significant associations was assessed; furthermore, many new original studies were published. Objectives: To determine the relationship between GSTM1, GSTT1, and GSTP1 polymorphisms with T2DM risk. Methods: PubMed, Embase, Wanfang, and China National Knowledge Infrastructure Databases were searched. We quantify the relationship using crude odds ratios and their 95% confidence intervals Moreover, the Venice criteria, false-positive report probability (FPRP), and Bayesian false discovery probability (BFDP) were used to validate the significance of the results. Results: Overall, significantly increased T2DM risk was found between individual and combined effects of GSTM1, GSTT1, and GSTP1 polymorphisms on T2DM risk, but, combined effects of the GSTT1 and GSTP1 polymorphisms was not statistically significant. GSTT1 gene polymorphism significantly increases the risk of T2DM complications, while GSTM1 and GSTP1 polymorphisms had no statistical significance. The GSTM1 null genotype was linked to a particularly increased risk of T2DM in Caucasians; the GSTT1 null genotype was connected to a significantly higher risk of T2DM in Asians and Indians; and the GSTP1 IIe105Val polymorphism was related to a substantially increased T2DM risk in Indians. Moreover, the GSTM1 and GSTT1 double null genotype was associated with substantially increased T2DM risk in Caucasians and Indians; the combined effects of GSTM1 and GSTP1 polymorphisms was associated with higher T2DM risk in Caucasians. However, all significant results were false when the Venice criteria, FPRP, and BFDP test were used (any FPRP >0.2 and BFDP value >0.8). Conclusion: The current analysis strongly suggests that the individual and combined effects of GSTM1, GSTT1 and GSTP1 polymorphisms might not be connected with elevated T2DM risk.

Results: Overall, significantly increased T2DM risk was found between individual and combined effects of GSTM1, GSTT1, and GSTP1 polymorphisms on T2DM risk, but, combined effects of the GSTT1 and GSTP1 polymorphisms was not statistically significant. GSTT1 gene polymorphism significantly increases the risk of T2DM complications, while GSTM1 and GSTP1 polymorphisms had no statistical significance. The GSTM1 null genotype was linked to a particularly increased risk of T2DM in Caucasians; the GSTT1 null genotype was connected to a significantly higher risk of T2DM in Asians and Indians; and the GSTP1 IIe105Val polymorphism was related to a substantially increased T2DM risk in Indians. Moreover, the GSTM1 and GSTT1 double null genotype was associated with substantially increased T2DM risk in Caucasians and Indians; the combined effects of GSTM1 and GSTP1 polymorphisms was associated with higher T2DM risk in Caucasians. However, all significant results were false when the Venice criteria, FPRP, and BFDP test were used (any FPRP >0.2 and BFDP value >0.8).

Introduction
Insulin resistance, excessive or improper glucagon secretion, and insufficient insulin secretion are symptoms of type 2 diabetes mellitus (T2DM), a chronic endocrine metabolic disorder. According to the newly released International Diabetes Federation Diabetes Atlas (10th edition), there are an estimated 537 million people living with diabetes in the world now, with that number expected to rise to 784 million by 2,045 based on current trends. In recent years, the incidence of diabetic complications, such as diabetic nephropathy and diabetic cardiovascular disease, has gradually increased. With the increase of people's lifestyle such as poor diet, sedentary and inactive, the prevalence of T2DM will continue to rise, which will undoubtedly bring a huge burden on human health and medical expenditure (https://diabetesatlas.org/en/.). T2DM pathogenesis is complicated and diverse, involving both environmental and genetic risk factors. Currently, antioxidant and detoxification gene polymorphisms are believed to play a significant role in the risk of T2DM and related complications (Wu et al., 2014;Pahwa et al., 2017;Dendup et al., 2018).
Oxidative stress is the result of imbalance of reactive oxygen species generation and clearance, which is closely related to insulin resistance and functional impairment of islet β-cells. GST is a superfamily protease encoded by multiple genes with antioxidant and detoxification functions (Tangvarasittichai, 2015;Vats et al., 2015). It catalyzes the formation of thioether bonds between glutathione and electrophiles, thereby catalyzing the reduction of hydrogen peroxide, reducing fatty acids, phospholipids and DNA bases to corresponding alcohols, and scavenging lipid free radicals and hydrogen peroxide, thus reducing body damage (Hayes et al., 2005;Tabatabaei-Malazy et al., 2017). GSTM1 and GSTT1 gene locus has been mapped on chromosome 1p13.3 and 22q11.2, respectively, resulting in the reduced activity of a functional gene product (Hayes and Strange, 2000;Gonul et al., 2012). GSTP1 gene polymorphism is a single nucleotide polymorphism, a codon 105 A-G mutation at exon 5 in GSTP1 polymorphism leads to change in isoleucine (IIe) to valine (Val), is linked to a change in decreased enzymatic activity (Zimniak et al., 1994;Mastana et al., 2013). Thus, these three gene mutations may increase the risk of T2DM based on biological effects.
Independent and combined effects of GSTM1, GSTT1 and GSTP1 IIe105Val polymorphisms with the risk of T2DM have been described in a number of original investigations, and the association of GSTs gene polymorphisms with T2DM complications has also been reported. Of note, these results were contradictory. Five published meta-analyses (Tang et al., 2013;Yi et al., 2013;Zhang et al., 2013;Saadat, 2017;Nath et al., 2019) reported the GSTM1, GSTT1 and GSTP1 IIe105Val polymorphisms with T2DM risk. Three published metaanalyses have investigated the association between GSTs gene polymorphisms and T2DM complications (Orlewski and Orlewska, 2015;Sun et al., 2015;Nath et al., 2019). Yet, their studies had certain limitations. This could be due to the limited sample size or the lack of quality assessments of the literature. More importantly, past meta-analyses did not assess the reliability of positive results. As a result, we conducted a systematic review of the literature and a latest meta-analysis to assess whether GSTM1, GSTT1, and GSTP1 gene contribute to the risk of T2DM and related complications. We are also investigated the effect of gene-gene interactions on T2DM risk.

Search strategy
The present study followed the Preferred Reporting Items for Systematic Review and Meta-Analysis (PRISMA) group's guidelines. The PubMed, Embase, Wan-fang, and China National Knowledge Infrastructure (CNKI) databases were used to conduct the literature search (deadline to 16 May 2022). The following was utilized as a search strategy (polymorphism OR variant OR variation OR mutation OR SNP OR genome-wide association study OR genetic association study OR genotype OR allele) and (diabetes OR mellitus OR diabetes mellitus) AND (glutathione S-transferase M1 OR GSTM1 OR glutathione S-transferase T1 OR GSTT1 OR glutathione S-transferase P1 OR GSTP1). References of the retrieved articles and review articles on this topic were also carefully examined in additional suitable studies. If necessary, the respective authors were contacted by e-mail. Moreover, there were no restrictions or limitations on language.

Inclusion and exclusion criteria
The inclusion criteria were: 1) case-control or cohort design; 2) studies must examine the correlation between the individual and combined effects of GST polymorphisms with T2DM or T2DM complication risk; 3) studies with odds ratio (OR) and their 95% confidence interval (CI), or enough data to calculate these numbers. Exclusion criteria were: 1) case reports, reviews, commentaries and meta-analyses; 2) studies that did not give genotype frequencies or could not determine the number of genotypes and alleles; 3) overlapping data or incomplete data.

Data extraction
Two investigators extracted and double-checked data. When their conclusions differed, through discussion were reached a consensus. The first author's name, year, country, source of control, match, race and other information were recorded, among which race was divided into "Asian," "African," "European," "Caucasian," "Indian" and "Mixed."

Quality score assessment
The quality assessment scale was independently assessed by two investigators (Supplementary Table S2). The total score was 18 points, studies scoring >12 were high, while scores of ≤12 were regarded as low quality in this meta-analysis. The scale was influenced by selection (two points), control and case source (four points), diagnostic criteria for T2DM (two points), Ascertainment of control (two points), matching (two points), Genotyping examination (two points), Association assessment (two points), and Hardy-Weinberg balance (HWE).

Statistical analysis
Stata 12.0 software (Stata Corporation, College Station, TX) was used to calculate all statistical analyses. The strength of the association of GSTM1, GSTT1, and GSTP1 polymorphisms with T2DM and T2DM complications risk were estimated by calculating the crude OR and their 95% CI. We used the Chisquare-based Q-test and I2 value to evaluate heterogeneity (Higgins et al., 2003). The random-effects model was used when p was less than 0.10 and/or I2 was larger than 50%. if not, the fixed-effects model was applied (DerSimonian and Laird, 2015). Furthermore, a meta-regression analysis was used to investigate sources of heterogeneity (Baker et al., 2009). Subgroups were carried out based on ethnicity, control source, control type, and matching. Sensitivity analysis was performed using the following methods: 1) excluded one study at a time; 2) removed low-quality or Hardy-Weinberg Disequilibrium (HWD) studies; 3) retained only high quality and HWE research. HWE was evaluated by Chi-square goodness-of-fit test, p < 0.05 was defined as HWE, otherwise as HWD. To determine publication bias, Begg's funnel plot (Begg and Mazumdar, 1994) and Egger's test (Egger et al., 1997) were conducted. If publication bias existed, the missing studies would be filled employing a nonparametric "trim and fill" method (Dual and Tweedie, 2000).

FIGURE 1
Flow diagram for searching studies in the current meta-analysis.
Combinations of GSTT1 and GSTP1 polymorphism were not associated with T2DM risk in overall analysis. However, when subgroup analysis for ethnicity was stratified, important association were found in Indians (model 6: OR = 1.75, 95% CI = 1.10-2.78, Table 6).

Heterogeneity and sensitivity analyses
We assessed heterogeneity source by applying a metaregression analysis. The source of heterogeneity for the GSTT1 null genotype and the combined effects of GSTT1 and GSTP1 was revealed to be the source of controls (p = 0.028 and model 1: p = 0.037). HWE (model 7: p = 0.045) was the source of heterogeneity between the combined effects of GSTM1, GSTT1, and GSTP1 gene and T2DM risk.

Publication bias
No publication bias was observed between GSTs polymorphisms with T2DM risk through the Begg's funnel plot and Egger's test (Figures not shown).

Discussion
Oxidative stress is considered to play an important role in the pathogenesis of diabetes mellitus and related complications. GST is a member of phase II metabolic isoenzyme group, which has the effect of regulating various cytotoxicity, genotoxicity and antioxidant. The risk genotypes expressed by each gene (GSTM1 null, GSTT1 null, and GSTP1 Val/Val) may decrease the antioxidant activity of the enzyme, thus increasing the susceptibility of T2DM and other diseases. There was much important evidence showing that GSTM1, GSTT1, and GSTP1 gene polymorphisms were potential genetic factors for T2DM. Yet, previous research remains contradictory. In addition, seven previously published

FIGURE 3
Forest plot of the association between GSTM1 polymorphism and T2DM complications risk.
Frontiers in Genetics frontiersin.org meta-analyses also had certain limitations, therefore the association of the GSTs polymorphisms with T2DM risk needs to be further evaluated. Furthermore, to our knowledge, this research is the first large-scale metaanalysis to investigate the relationship between individual and combined GSTM1, GSTT1, and GSTP1 genes Frontiers in Genetics frontiersin.org 15 polymorphism of T2DM susceptibility, and utilizing FPRP and BFDP to explore the validity of positive results.
Overall, on the individual and combined effects of GSTM1, GSTT1, and GSTP1 polymorphisms, a statistically significantly increased T2DM risk was found. We found that GSTT1 null genotype increased the risk of T2DM complications.
However, GSTM1 and GSTP1 polymorphism were not associated with T2DM complications risk. The GSTM1 null genotype was correlated to a significantly increased risk of T2DM in Caucasians, Asia, hospital-based studies, and non-diabetic controls. The GSTT1 null genotype was connected with an increased risk of T2DM in Asians, Indians, Asia, and so on. GSTP1 IIe105Val polymorphism was correlated with an increased risk of T2DM in Indians, Asia, population-based studies, and healthy controls. The GSTM1-GSTT1 double null genotype significantly increased the risk of T2DM in both subgroup and sensitivity analysis. The effects of GSTM1-GSTP1 gene polymorphisms were related to the risk of T2DM in Caucasians and matching. However, no relationship was identified between GSTT1-GSTP1 polymorphisms and T2DM susceptibility, and, interestingly, Indians had a significantly increased risk. For the triple genotype combinations, we also discovered an association with T2DM. These findings showed that the same gene has variable susceptibility to T2DM in different race groups. It is worth considering that the obvious inconsistencies could be due to sample size, GSTs prevalence and environmental factors (lifestyle). Additionally, because T2DM is a chronic metabolic disease with multiple factors and genes, and different genetic backgrounds may lead to such differences. We found that GSTs gene polymorphism increased the risk of T2DM in Asia, but not in Europe. Most of the original research were conducted in Asia, which may lead to this difference. Therefore, a larger sample size is needed to study the relationship between GSTs gene polymorphism and T2DM in Europe, America, Africa and other geographic regions. The GSTM1-GSTT1, GSTM1-GSTP1 and GSTM1-GSTT1-GSTP1 gene polymorphisms combination are all associated with the increased susceptibility of T2DM. These results indicate that polygenes may have cumulative effects, but the potential gene-gene interactions still need to be studied with a larger sample size.
However, the present study used several subgroups and diverse genetic models, resulting in multiple comparisons, which necessitated adjusting the pooled p value (Attia et al., 2003). FPRP is regarded an acceptable tool for calculating the likelihood of positive findings in multiple hypothesis testing of molecular epidemiology investigations (Wacholder et al., 2004). Wakefield provided a much more reliable Bayesian metric of false positives in genetic epidemiology research (Wakefield, 2008). Potential errors and bias may be caused by genotyping errors or high inter-study heterogeneity (I 2 > 50%) (Page et al., 2003;Clayton et al., 2005;Pompanon et al., 2005;Balding, 2006;Zhang et al., 2013), while large volumes evidence (statistical power >80%) have stricter statistical significance levels or lower false-discovery rate (Benjamini and Hochberg, 1995). Therefore, we consider using FPRP tests, BFDP tests and the Venice criteria to evaluate the false findings in present studies. The current meta-analysis found that all substantial relationship were deemed as less credible during did the reliability analysis. According to the results of the meta-regression analysis, the HWD studies and sources of controls were the main causes of heterogeneity. Bias and errors are prevalent in some low-quality, small-sample HWD studies, making the conclusions of these original studies undependable, molecular biology and disease susceptibility study, in particular. Furthermore, research with a small sample size and substantial results may be likely to swallow than studies with negative outcomes. When they produce positive results, however, their research might not have been rigorous and typically of poor quality.
Therefore, we assessed the sensitivity analysis using matching, high-quality, and HWE.
These three GSTs gene polymorphisms and their interactions were discovered to be strongly related to diabetes susceptibility in this investigation. From 2012 to 2019, a total of five related metastudies were published. Recent meta-studies by Nath et al. (2019) showed that the GSTM1, GSTT1 null genotype and GSTM1-GSTT1 double null genotype increased the risk of T2DM in both Asian and Caucasian populations, the results are consistent with this study. An analysis of eleven studies by Yi et al. (2013) revealed that the GSTM1 gene polymorphisms increases the susceptibility to T2DM in Asians, Caucasians, and Africans. However, Saadat (2017) selected 18 studies (2,595 patients and 2,888 controls) and suggested that GSTP1 IIe105Val polymorphism was not connected with the risk of T2DM, which was contrary to the results of this study. Tang et al. (2013) also analyzed the three genes and suggested that GSTM1 and GSTT1 null genotypes were correlated to an increased risk of T2DM in Asians, while GSTP1 gene polymorphism was not. For the GSTP1 gene, the opposite result was obtained. Three related meta-analysis were published from 2015 to 2019. Nath et al. showed that GSTT1 null genotype increased the risk of T2DM related complications, but GSTM1 null genotype did not, which was consistent with this study. On the contrary, Orlewski and Orlewska (2015) found that GSTM1 gene polymorphisms increased the risk of diabetic nephropathy, while GSTT1 and GSTP1 gene polymorphisms were not statistically significant. Sun et al. (2015) suggested that a strong association GSTM1 and GSTP1 polymorphisms and diabetic retinopathy susceptibility. Further investigate the effect of gene-gene interaction on T2DM, two previous studies reported the correlation between GSTM1-GSTT1 gene polymorphism and T2DM susceptibility, but did not study the combined effect between the remaining genes. In Frontiers in Genetics frontiersin.org this study, all except GSTT1-GSTP1 combination were found to be related. Moreover, these published meta-analyses are obviously inconsistent with current meta-analyses in terms of racial classification. The variability in results may be due to a variety of factors, such as small sample sizes, ethnic differences, failure to assess quality in most eligible studies and failure to develop a more complete genetic model. There are several advantages to the newer meta-analysis. First, The FPRP, BFDP test, and Venice criterion were used to assess the credibility with positive outcomes. Secondly, quality assessment of qualified research is carried out. Third, the sample size was much larger and more detailed data and information were collected. Fourth, quality assessment of qualified research. Fifthly, a relatively complete genetic model is established. Sixth, more subgroup analyses were performed. Seventh, to our knowledge, this is the first metaanalysis to explore the combined effects of GSTM1-GSTP1, GSTT1-GSTP1, and GSTM1-GSTT1-GSTP1 gene polymorphisms on the risk of T2DM. There are some limitations of the present meta-analysis. To begin with, latest meta-analyses only include published articles, and positive results are frequently more likely to be published than negative consequences. The zero effect of GSTM1 may be overestimated if negative results are included. Second, because we could not calculate the HWE of these two genes, we did not evaluate whether the genotype distributions of GSTM1 and GSTT1 polymorphisms in the control group was in HWE. Third, heterogeneity between studies was much large in overall and among several subgroups. According to the results of the meta-regression analysis, the HWD studies and sources of controls were the main causes of heterogeneity, and adequate statistical correction and analysis did not significantly improve the heterogeneity. The source of heterogeneity may also be due to some other factors, such as study area, lifestyle and so on. Fourth, we conducted a subgroup analysis of race and region, and population heterogeneity has been resolved. However, T2DM is closely associated with the environment, the original study did not provide data on lifestyle, so we failed to solve the heterogeneity caused by lifestyle, and did not explore the effect of the geneenvironment combination. Fifth, we failed to further investigate the effect of GSTs gene polymorphisms on the risk of T2DM related complications.
In conclusion, the current study revealed the relationship between GSTM1, GSTT1, and GSTP1 gene polymorphisms with the risk of T2DM, which may be attributed to false positive results rather than actual correlations or biological variables. Larger epidemiological studies should be carried out in the future to confirm or deny our findings.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author contributions L-sL: designed research, performed research, collected data, analyzed data, wrote paper. DW and RT: collected data. QW, LZ and JW: check the data. YL and X-fH: designed research and revised article.