*Article* **Identification of Antitumor** *miR-30e-5p* **Controlled Genes; Diagnostic and Prognostic Biomarkers for Head and Neck Squamous Cell Carcinoma**

**Chikashi Minemura 1,† , Shunichi Asai 2,3,†, Ayaka Koma <sup>1</sup> , Naoko Kikkawa 2,3, Mayuko Kato <sup>2</sup> , Atsushi Kasamatsu <sup>1</sup> , Katsuhiro Uzawa <sup>1</sup> , Toyoyuki Hanazawa <sup>3</sup> and Naohiko Seki 2,\***


**Abstract:** Analysis of microRNA (miRNA) expression signatures in head and neck squamous cell carcinoma (HNSCC) has revealed that the *miR-30* family is frequently downregulated in cancer tissues. The Cancer Genome Atlas (TCGA) database confirms that all members of the *miR-30* family (except *miR-30c-5p*) are downregulated in HNSCC tissues. Moreover, low expression of *miR-30e-5p* and *miR-30c-1-3p* significantly predicts shorter survival of HNSCC patients (*p* = 0.0081 and *p* = 0.0224, respectively). In this study, we focused on *miR-30e-5p* to investigate its tumor-suppressive roles and its control of oncogenic genes in HNSCC cells. Transient expression of *miR-30e-5p* significantly attenuated cancer cell migration and invasive abilities in HNSCC cells. Nine genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2*, and *RTN4R*) were identified as putative targets of *miR-30e-5p* control. Their expression levels significantly predicted shorter survival of HNSCC patients (*p* < 0.05). Among those targets, *FOXD1* expression appeared to be an independent factor predicting patient survival according to multivariate Cox regression analysis (*p* = 0.049). Knockdown assays using siRNAs corresponding to *FOXD1* showed that malignant phenotypes (e.g., cell proliferation, migration, and invasive abilities) of HNSCC cells were significantly suppressed. Overexpression of *FOXD1* was confirmed by immunostaining of HNSCC clinical specimens. Our miRNA-based approach is an effective strategy for the identification of prognostic markers and therapeutic target molecules in HNSCC. Moreover, these findings led to insights into the molecular pathogenesis of HNSCC.

**Keywords:** microRNA; HNSCC; tumor-suppressor; *miR-30e-5p*; *FOXD1*; TCGA

### **1. Introduction**

A head and neck squamous cell carcinoma (HNSCC) may arise in the oral cavity, hypopharynx, nasopharynx, and larynx. Based on Global Cancer Statistics 2018, HNSCC is the eighth most common cancer in the world [1]. Approximately 890,000 people are diagnosed with HNSCC annually, and 450,000 die from the disease [2]. A notable characteristic of HNSCC is that many patients are already in the advanced stage of disease at the time of diagnosis. Treatment strategies are limited for patients at such an advanced stage [3]. Cisplatin-based treatment is selected for patients for whom surgery is not indicated. However, cancer cells develop drug-resistance in the course of treatment [4]. No effective treatments have been reported for patients with HNSCC who have acquired drug

**Citation:** Minemura, C.; Asai, S.; Koma, A.; Kikkawa, N.; Kato, M.; Kasamatsu, A.; Uzawa, K.; Hanazawa, T.; Seki, N. Identification of Antitumor *miR-30e-5p* Controlled Genes; Diagnostic and Prognostic Biomarkers for Head and Neck Squamous Cell Carcinoma. *Genes* **2022**, *13*, 1225. https://doi.org/ 10.3390/genes13071225

Academic Editors: Giuseppe Iacomino and Fabio Lauria

Received: 12 June 2022 Accepted: 7 July 2022 Published: 9 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

resistance [5]. Thus, there is a pressing need for improved understanding of the molecular pathogenesis of HNSCC using the latest genomic science.

The human genome project showed that a vast number of non-coding RNA molecules (ncRNAs) are transcribed and function in normal and diseased cells [6,7]. MicroRNAs (miR-NAs), which constitute a class of ncRNA, are single strands of RNA, 19 to 23 bases in length. The unique nature of miRNA is that a single miRNA negatively controls a large number of RNA transcripts (both protein coding RNAs and ncRNAs) in a sequence-dependent manner [8,9]. Therefore, dysregulated miRNA expression can disrupt tightly controlled RNA networks in normal cells. The involvement of aberrantly expressed miRNAs has been reported in a variety of human diseases, including cancer [10].

Over the past decade, a large number of studies have shown that aberrantly expressed miRNAs are closely connected to human oncogenesis [11]. Identifying miRNAs with altered expression in cancer cells is a first step in the characterization of their role. The latest RNA-sequence technology is accelerating the creation of miRNA expression signatures in a variety of cancer types [12]. We have established miRNA expression signatures in several types of cancers, including HNSCC [13,14].

The analysis of miRNA signatures has revealed that members of the *miR-30* family are frequently downregulated in several types of cancers, including HNSCC [15,16]. In the human genome, the *miR-30* family is composed of five species, i.e., *miR-30a* to *miR-30e*. *miR-30c* is further subdivided into *miR-30c-1* and *miR-30c-2*. The *miR-30* family is encoded by six genes located on human chromosomes 1, 6, and 8 [17]. The guide strands of the *miR-30* family share the same seed sequence (GUAAACA). The passenger strands of the *miR-30* family are divided into two groups according to their seed sequences. The first group (UUUCAGU) includes *miR-30a-3p*, *miR-30d-3p*, and *miR-30e-3p*. The second group (UGGGAG) encompasses *miR-30b-3p*, *miR-30c-1-3p*, and *miR-30c-2-3p*. We theorized that clarifying the targets/pathways controlled by miRNAs for each cancer type could greatly enhance our understanding of the molecular pathogenesis of cancer cells.

The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga), a landmark cancer genomics program, molecularly characterizes over 20,000 primary cancers (and matched normal samples) spanning 33 cancer types. Based on a TCGA analysis, all members of the *miR-30* family (except *miR-30c-5p*) are downregulated in HNSCC tissues. Moreover, the low expression of two miRNAs, *miR-30e-5p* and *miR-30c-1-3p*, significantly predicts shorter survival of patients with HNSCC (*p* = 0.0081 and *p* = 0.0224, respectively).

In this study, we focused on *miR-30e-5p* and investigated its control of genes involved in the molecular pathogenesis of HNSCC. A total of nine genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2*, and *RTN4R*) were identified as putative targets of *miR-30e-5p*. Their expression level predicted shorter survival of the patients (*p* < 0.05).

Among these genes, we focused on *FOXD1* (Forkhead Box D1), which belongs to the forkhead family of transcription factors. Here, we investigated its functional significance in HNSCC cells.

RNA immunoprecipitation (RIP) assays revealed that *FOXD1* mRNA was incorporated into the RNA-induced silencing complex (RISC). Additionally, dual luciferase reporter assays showed that *miR-30e-5p* was directly bound to the 30 -UTR of *FOXD1*. Functional assays using siRNAs targeting *FOXD1* showed that the knockdown of *FOXD1* expression attenuated HNSCC cell malignant behaviors (migration and invasion). Our approach will reveal the new insights into the involvement of miRNA and the molecular pathogenesis of HNSCC.

### **2. Materials and Methods**

### *2.1. Analysis of miRNAs and miRNA Target Genes in HNSCC Patients*

The sequences of the *miR-30* family were confirmed using miRbase ver. 22.1 (https: //www.mirbase.org, accessed on 10 July 2020) [18]. TCGA-HNSC (TCGA, Firehose Legacy) was used to investigate the clinical significance of miRNAs and their target genes. Clinical parameters and gene expression data were obtained from cBioPortal

(http://www.cbioportal.org/, data downloaded on 20 August 2020) and OncoLnc (http: //www.oncolnc.org/, data downloaded on 20 August 2020) [19–21].

The putative target genes with the *miR-30-5p* binding sites were selected by using TargetScanHuman ver. 7.2 (http://www.targetscan.org/vert\_72/, accessed on 22 November 2020) [22].

To identify differentially expressed genes in HNSCC tissues, a newly created microarray data (accession number: GSE180077) was used in this study. Three hypopharyngeal squamous cell carcinoma (HSCC) tissues, three normal hypopharyngeal tissues, and two cervical lymph nodes harvested from one HSCC patient who underwent surgical resection at Chiba University Hospital were subjected to Agilent whole genome microarrays. In this study, we compared gene expressions in cancer tissues with those in normal tissues. The clinical information of this patient was summarized in Table S1.

### *2.2. HNSCC Cell Lines*

Two human HNSCC cell lines were obtained from the RIKEN BioResource Center (Tsukuba, Ibaraki, Japan). Sa3 was harvested from upper gingiva cancer of 63y male. SAS was harvested from oral tongue cancer of 69y female.

### *2.3. RNA Extraction and Quantitative Real-Time Reverse-Transcription PCR (qRT-PCR)*

Total RNA was isolated using TRIzol reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer's protocol. qRT-PCR was performed using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Waltham, MA, USA) and the StepOnePlus™ Real-Time PCR System (Applied Biosystems). Gene expressions were quantified relatively by the delta-delta Ct method (used *GAPDH* as internal control). Taqman assays used in this study are summarized in Table S2.

### *2.4. Transfection of miRNAs and siRNAs into HNSCC Cells*

The protocols used for transient transfection of miRNAs and siRNAs were described in our previous studies [23–26]. The miRNA precursors and siRNAs used in this report were detailed in Table S2. The miRNAs and siRNAs were transfected into HNSCC cell lines using Opti-MEM (Gibco, Carlsbad, CA, USA) and LipofectamineTM RNAiMax Transfection Reagent (Invitrogen, Waltham, MA, USA). All miRNA precursors were transfected into HNSCC cell lines at 10 nM, and siRNAs were transfected at 5 nM. Mock transfection consisted of cells without precursors or siRNAs. Control groups were transfected with the negative control precursor.

### *2.5. RIP Assay*

The procedures for the RIP assay have been described previously [23]. Briefly, SAS cells were cultured in 10-cm plates at 80% confluency. Negative control miRNA precursors and *miR-30e-5p* precursors were transfected. After 6 h, immunoprecipitation was performed using the MagCaptureTM microRNA Isolation Kit, Human Ago2, obtained from FUJIFILM Wako Pure Chemical Corporation (Wako, Osaka, Japan) according to the manufacturer's protocol. Expression levels of *FOXD1* bound to Ago2 were measured by qRT-PCR. Taqman primers used in the assay are summarized in Table S2.

### *2.6. Functional Assays of HNSCC Cells (Cell Proliferation, Migration, and Invasion Assays)*

The procedures for conducting functional assays (cell proliferation, migration, and invasion assays) with HNSCC cells were described in earlier publications [23–26]. In brief, for proliferation assays, Sa3 or SAS cells were transferred to 96-well plates at 3.0 <sup>×</sup> <sup>10</sup><sup>3</sup> cells per well. Cell proliferation was evaluated using the XTT assay kit II (Sigma–Aldrich, St. Louis, MO, USA) 72 h after the transfection procedure. For migration and invasion assays, Sa3 and SAS cells were transfected in 6-well plates at 2.0 <sup>×</sup> <sup>10</sup><sup>5</sup> cells per well. After 48 h, transfected Sa3 and SAS cells were added into each chamber at 1.0 <sup>×</sup> <sup>10</sup><sup>5</sup> per well. Corning BioCoatTM cell culture chambers (Corning, Corning, NY, USA) were used for

migration assays whereas Corning BioCoat Matrigel Invasion Chambers were used for invasion assays. After 48 h, the cells on the lower surface of chamber membranes were stained and counted for analysis. All experiments were performed in triplicate.

### *2.7. Plasmid Construction and Dual-Luciferase Reporter Assays*

The procedures for plasmid construction and the dual-luciferase reporter assays were outlined previously [23–26]. Briefly, a partial wild-type sequence, including the seed sequence of *FOXD1* within the 30 -untranslated region (30 -UTR), was inserted into the psiCHECK-2 vector (C8021; Promega, Madison, WI, USA). Alternatively, a deletion type that was missing the *miR-30e-5p* target site was synthesized. These synthesized vectors (50 ng) were transfected into 1.0 <sup>×</sup> <sup>10</sup><sup>5</sup> cells in each well using Lipofectamine 2000 (Invitrogen, Waltham, MA, USA) and Opti-MEM. After 48 h of transfection, dual luciferase reporter assays using the Dual Luciferase Reporter Assay System (Promega, Madison, WI, USA) were conducted. Luminescence data were presented as the *Renilla/Firefly* luciferase activity ratio.

### *2.8. Immunohistochemistry*

The procedures for immunohistochemistry were described in our previous studies [23–26]. The clinical samples were obtained from HNSCC cases who received surgical treatment at Chiba University Hospital. The slides were incubated with primary antibody FOXD1 (1: 50, PA5-27142, Thermo Fisher Scientific, Waltham, MA, USA). The clinical features are shown in Table S3. The reagents used in the analysis are listed in Table S2.

### *2.9. Gene Set Enrichment Analysis (GSEA)*

To analyze the molecular pathways related to *FOXD1* (regulated by *miR-30e-5p*), GSEA was performed. Using TCGA-HNSC data, HNSCC patients were divided into high and low expression groups according to the Z-score of the *FOXD1* expression level. A ranked list of genes was generated by the log<sup>2</sup> ratio comparing the expression levels of each gene between the two groups. The ranked gene list was uploaded into GSEA software [27,28] and we applied the Hallmark gene set in The Molecular Signatures Database [27,29].

### *2.10. Statistical Analysis*

Statistical analyses were determined using JMP Pro 15 (SAS Institute Inc., Cary, NC, USA). Welch's *t*-tests were performed to determine the significance of differences between two groups. Dunnett's tests were applied for comparisons among multiple groups. For correlation analyses, Spearman's test was applied. Survival analyses were analyzed by log-rank test. Each expression levels of target genes, age, disease stage, and pathological grade in TCGA-HNSC were used as variables for Cox's proportional hazards model. A *p*-value less than 0.05 was considered statistically significant.

### **3. Results**

### *3.1. Expression Levels and the Clinical Significance of the miR-30 Family in HNSCC Clinical Specimens Assessed by TCGA Analysis*

The human *miR-30* family consists of 12 members. They include *miR-30a-5p*, *miR-30a-3p*, *miR-30b-5p*, *miR-30b-3p*, *miR-30c-1-5p*, *miR-30c-1-3p*, *miR-30c-2-5p*, *miR-30c-2-3p*, *miR-30d-5p*, *miR-30d-3p*, *miR-30e-5p*, and *miR-30e-3p*. Seed sequences of *miR-30-5p* are identical. In contrast, two types of seed sequences are present in *miR-30-3p* (Figure 1A). Human *miR-30e* and *miR-30c-1* are located on chromosome 1p34.2, whereas *miR-30c-2* and *miR-30a* are on chromosome 6q13, and *miR-30b* and *miR-30d* are found on chromosome 8q24.22 (Figure 1B).

*Genes* **2022**, *13*, 1225 5 of 22

(Figure 1B).

(Figure 1B).

**Figure 1.** Twelve miRNAs are present in the human genome as members of the *miR-30* family. (**A**) Seed sequences of *miR-30-5p* are identical. In contrast, two types of seed sequences are present in *miR-30-3p*. (**B**) The locations of each microRNA family on chromosome. **Figure 1.** Twelve miRNAs are present in the human genome as members of the *miR-30* family. (**A**) Seed sequences of *miR-30-5p* are identical. In contrast, two types of seed sequences are present in *miR-30-3p*. (**B**) The locations of each microRNA family on chromosome. **Figure 1.** Twelve miRNAs are present in the human genome as members of the *miR-30* family. (**A**) Seed sequences of *miR-30-5p* are identical. In contrast, two types of seed sequences are present in *miR-30-3p*. (**B**) The locations of each microRNA family on chromosome.

*miR-30e* and *miR-30c-1* are located on chromosome 1p34.2, whereas *miR-30c-2* and *miR-30a* are on chromosome 6q13, and *miR-30b* and *miR-30d* are found on chromosome 8q24.22

*miR-30e* and *miR-30c-1* are located on chromosome 1p34.2, whereas *miR-30c-2* and *miR-30a* are on chromosome 6q13, and *miR-30b* and *miR-30d* are found on chromosome 8q24.22

Expression levels of all members of the *miR-30* family were validated by TCGA database analyses. All members of the *miR-30* family (except for *miR-30c-5p*) were significantly downregulated in HNSCC tissues (*n* = 484) compared to normal tissues (*n* = 44) (Figure 2). Expression levels of all members of the *miR-30* family were validated by TCGA database analyses. All members of the *miR-30* family (except for *miR-30c-5p*) were significantly downregulated in HNSCC tissues (*n* = 484) compared to normal tissues (*n* = 44) (Figure 2). Expression levels of all members of the *miR-30* family were validated by TCGA database analyses. All members of the *miR-30* family (except for *miR-30c-5p*) were significantly downregulated in HNSCC tissues (*n* = 484) compared to normal tissues (*n* = 44) (Figure 2).

**Figure 2.** The expression level of the *miR-30* family was analyzed using TCGA-HNSC database. A total of 484 HNSCC tissues and 44 normal epithelial tissues were evaluated (N.S.: not significant).

*Genes* **2022**, *13*, 1225 6 of 22

Moreover, low expression of *miR-30e-5p* and *miR-30c-1-3p* significantly predicted poor prognosis of patients with HNSCC (Figure 3). Expression levels of other members of the *miR-30* family were not related to patient prognosis (Figure 3). poor prognosis of patients with HNSCC (Figure 3). Expression levels of other members of the *miR-30* family were not related to patient prognosis (Figure 3).

**Figure 2.** The expression level of the *miR-30* family was analyzed using TCGA-HNSC database. A total of 484 HNSCC tissues and 44 normal epithelial tissues were evaluated (N.S.: not significant).

Moreover, low expression of *miR-30e-5p* and *miR-30c-1-3p* significantly predicted

**Figure 3.** Kaplan–Meier survival analyses of HNSCC patients using data from TCGA-HNSC. Patients were divided into two groups according to the median miRNA expression level: high and low expression groups. The red and blue lines represent the high and low expression groups, **Figure 3.** Kaplan–Meier survival analyses of HNSCC patients using data from TCGA-HNSC. Patients were divided into two groups according to the median miRNA expression level: high and low expression groups. The red and blue lines represent the high and low expression groups, respectively.

#### respectively. *3.2. Effect of Transient Transfection of miR-30e-5p on HNSCC Cell Proliferation, Migration and Invasion*

*3.2. Effect of Transient Transfection of miR-30e-5p on HNSCC Cell Proliferation*, *Migration and Invasion*  In this study, we focused on *miR-30e-5p* because its expression was significantly downregulated in HNSCC tissues, and it was closely associated with poor prognosis of the patients, suggesting that *miR-30e-5p* acts as a tumor-suppressive miRNA in HNSCC cells.

In this study, we focused on *miR-30e-5p* because its expression was significantly downregulated in HNSCC tissues, and it was closely associated with poor prognosis of the patients, suggesting that *miR-30e-5p* acts as a tumor-suppressive miRNA in HNSCC cells. The tumor-suppressive roles of *miR-30e-5p* were assessed by transient transfection of The tumor-suppressive roles of *miR-30e-5p* were assessed by transient transfection of *miR-30e-5p* in two cell lines, Sa3 and SAS. Transient transfection of *miR-30e-5p* inhibited HNSCC cell proliferation (Figure 4A). Cancer cell migration and invasive abilities were markedly blocked by *miR-30e-5p* expression in Sa3 and SAS cells (Figure 4B,C). Photographs of typical results from the migration and invasion assays are shown in Figure S1.

*miR-30e-5p* in two cell lines, Sa3 and SAS. Transient transfection of *miR-30e-5p* inhibited HNSCC cell proliferation (Figure 4A). Cancer cell migration and invasive abilities were markedly blocked by *miR-30e-5p* expression in Sa3 and SAS cells (Figure 4B,C). Photographs of typical results from the migration and invasion assays are shown in Figure

S1.

**Figure 4.** Functional assays of *miR-30e-5p* in HNSCC cell lines (Sa3 and SAS)**.** (**A**) Cell proliferation was assessed using XTT assays 72 h after miRNA transfection. (**B**) Cell migration was assessed using a membrane culture system 48 h after seeding miRNA-transfected cells into the chambers. (**C**) Cell invasion was determined using Matrigel invasion assays 48 h after seeding miRNA-transfected cells into the chambers. **Figure 4.** Functional assays of *miR-30e-5p* in HNSCC cell lines (Sa3 and SAS)**.** (**A**) Cell proliferation was assessed using XTT assays 72 h after miRNA transfection. (**B**) Cell migration was assessed using a membrane culture system 48 h after seeding miRNA-transfected cells into the chambers. (**C**) Cell invasion was determined using Matrigel invasion assays 48 h after seeding miRNA-transfected cells into the chambers.

#### *3.3. Screening for Oncogenic Targets of miR-30e-5p in HNSCC 3.3. Screening for Oncogenic Targets of miR-30e-5p in HNSCC*

To identify genes modulated by *miR-30e-5p* that were closely involved in HNSCC molecular pathogenesis, we performed in silico database analyses combined with our gene expression data. We created new gene expression data using clinical specimens of hypopharyngeal squamous cell carcinoma (HSCC). Three HSCC tissues, three normal hypopharyngeal tissues, and two cervical lymph nodes harvested from one HSCC patient who underwent surgical resection at Chiba University Hospital were subjected to Agilent whole genome microarrays. In this study, we compared gene expressions in cancer tissues with those in normal tissues. The clinical information of this patient was summarized in Table S1. Expression data were deposited in the GEO database (accession number: GSE180077). Our selection strategy is shown in Figure 5. To identify genes modulated by *miR-30e-5p* that were closely involved in HNSCC molecular pathogenesis, we performed in silico database analyses combined with our gene expression data. We created new gene expression data using clinical specimens of hypopharyngeal squamous cell carcinoma (HSCC). Three HSCC tissues, three normal hypopharyngeal tissues, and two cervical lymph nodes harvested from one HSCC patient who underwent surgical resection at Chiba University Hospital were subjected to Agilent whole genome microarrays. In this study, we compared gene expressions in cancer tissues with those in normal tissues. The clinical information of this patient was summarized in Table S1. Expression data were deposited in the GEO database (accession number:GSE180077). Our selection strategy is shown in Figure 5.

**Figure 5.** Flow chart of the strategy used to identify putative tumor suppressor genes regulated by *miR-30e-5p*.

**Figure 5.** Flow chart of the strategy used to identify putative tumor suppressor genes regulated by

*miR-30e-5p*. TargetScan Human database (release 7.2) provides data on the putative targets of *miR-30e-5p*. A total of 1576 genes are listed. Among these targets, 97 genes were TargetScan Human database (release 7.2) provides data on the putative targets of *miR-30e-5p*. A total of 1576 genes are listed. Among these targets, 97 genes were upregulated in HSCC patients. Furthermore, the increased expression of 54 genes was confirmed in TCGA database analyses (Table 1).


upregulated in HSCC patients. Furthermore, the increased expression of 54 genes was confirmed in TCGA database analyses (Table 1). **Table 1.** Candidate target genes regulated by *miR-30-5p*.

84733 *CBX2* Chromobox homolog 2 1 0.005 3.00 0.133 27 *ABL2* ABL proto-oncogene 2, non-receptor tyrosine kinase 1 0.008 2.19 0.155


**Table 1.** *Cont.*

<sup>1</sup> Gene Expression Omnibus, <sup>2</sup> Fold Change, <sup>3</sup> 5-year Overall Survival.

### *3.4. Clinical Significance of miR-30e-5p Targets in Patients with HNSCC Determined by TCGA Analysis*

Among the 54 putative target genes, high expression of nine genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2*, and *RTN4R*) showed statistically significant correlations with the 5-year overall survival frequencies of patients with HNSCC (*p* < 0.05; Figures 6 and 7).

Therefore, univariate analysis for 5-year overall survival was conducted first, and then multivariate analysis was performed for the statistically significant variables (*p* < 0.05). Each expression levels of target genes, age, disease stage, and pathological grade in TCGA-HNSC were used as variables for Cox's proportional hazards model. As a result, the high expression level of *FOXD1* (HR: 1.374, 95% CI: 1.002–1.890, *p* = 0.049), age (≥70) (HR: 1.922, 95% CI: 1.369–2.698, *p* < 0.001) and disease stage (III and IV) (HR: 1.774, 95% CI: 1.159–2.716, *p* = 0.008) were independent prognostic factors (Table 2).

The correlations of the expression levels between these nine genes and *miR-30e-5p* were evaluated by TCGA-HNSC. A Spearman's rank test confirmed weak negative correlations

in seven genes (*DDIT4*, *FOXD1*, *FZD2*, *MINPP1*, *PAWR*, *PFN2*, and *RTN4R*) (*p* < 0.001, *r* < −0.2, Figure 8), whereas the correlations were not confirmed in *FXR1* and *HMGB3*.

*Genes* **2022**, *13*, 1225 10 of 22

**Figure 6.** Expression levels of nine target genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2,* and *RTN4R*) in HNSCC clinical specimens from TCGA-HNSC. All genes were found to be upregulated in HNSCC tissues (*n* = 518) compared with normal tissues (*n* = 44). **Figure 6.** Expression levels of nine target genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2*, and *RTN4R*) in HNSCC clinical specimens from TCGA-HNSC. All genes were found to be upregulated in HNSCC tissues (*n* = 518) compared with normal tissues (*n* = 44).

**Figure 7.** Clinical significance of nine target genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2,* and *RTN4R*) according to TCGA-HNSC data analysis. Kaplan–Meier curves of the 5 year overall survival rates according to the expression of each gene are presented. High expression levels of all nine genes were significantly predictive of a poorer prognosis in patients with HNSCC. Patients were divided into two groups according to the median gene expression level: high and low expression groups. The red and blue lines represent the high and low expression groups, **Figure 7.** Clinical significance of nine target genes (*DDIT4*, *FOXD1*, *FXR1*, *FZD2*, *HMGB3*, *MINPP1*, *PAWR*, *PFN2*, and *RTN4R*) according to TCGA-HNSC data analysis. Kaplan–Meier curves of the 5-year overall survival rates according to the expression of each gene are presented. High expression levels of all nine genes were significantly predictive of a poorer prognosis in patients with HNSCC. Patients were divided into two groups according to the median gene expression level: high and low expression groups. The red and blue lines represent the high and low expression groups, respectively.



*FXR1* (High vs. Low expression) 1.651 1.242–2.193 0.001 1.303 0.930–1.825 0.124 HR: hazard ratio, CI: confidence interval.

respectively.

*FZD2* (High vs. Low expression) 1.337 1.011–1.768 0.042 1.069 0.776–1.473 0.684 *HMGB3* (High vs. Low expression) 1.413 1.069–1.867 0.015 1.186 0.861–1.633 0.297 *MINPP1* (High vs. Low expression) 1.451 1.096–1.921 0.009 1.308 0.940–1.818 0.111

*PAWR* (High vs. Low expression) 1.438 1.084–1.908 0.012 1.255 0.913–1.727 0.162 *PFN2* (High vs. Low expression) 1.642 1.238–2.178 0.001 1.238 0.885–1.731 0.213 *RTN4R* (High vs. Low expression) 1.326 1.003–1.752 0.048 0.962 0.699–1.323 0.810

Disease Stage (III, IV vs. I, II) 1.746 1.158–2.633 0.008 1.774 1.159–2.716 0.008 Pathological Grade (3, 4 vs. 1, 2) 0.901 0.653–1.245 0.529 - - -

HR: hazard ratio, CI: confidence interval.

*HMGB3*.

Age (≥70 vs. <70) 1.628 1.197–2.213 0.002 1.922 1.369–2.698 < 0.001

The correlations of the expression levels between these nine genes and *miR-30e-5p*  were evaluated by TCGA-HNSC. A Spearman's rank test confirmed weak negative correlations in seven genes (*DDIT4*, *FOXD1*, *FZD2*, *MINPP1*, *PAWR*, *PFN2,* and *RTN4R*) (*p* < 0.001, *r* < −0.2, Figure 8), whereas the correlations were not confirmed in *FXR1* and

**Figure 8.** Correlation analysis by TCGA-HNSC for nine genes. Negative correlation of expression levels between *miR-30e-5p* and nine genes in HNSCC clinical specimens. **Figure 8.** Correlation analysis by TCGA-HNSC for nine genes. Negative correlation of expression levels between *miR-30e-5p* and nine genes in HNSCC clinical specimens.

### *3.5. Regulated Expression of the Nine Identified Genes by miR-30e-5p in HNSCC Cells*

qRT-PCR revealed that the mRNA expression levels of *DDIT4*, *FOXD1*, *MINPP1*, *PAWR*, and *PFN2* were significantly suppressed in *miR-30e-5p*-transfected HNSCC cells. The other genes were not significantly suppressed by *miR-30e-5p* (Figure 9).

Considering the multivariate analysis, the correlation analysis, and the results of qRT-PCR, *FOXD1* was of the greatest interest among the nine targets. Thus, we focused on *FOXD1* as the target gene of *miR-30e-5p* in this study.

*3.5. Regulated Expression of the Nine Identified Genes by miR-30e-5p in HNSCC Cells* 

The other genes were not significantly suppressed by *miR-30e-5p* (Figure 9).

qRT-PCR revealed that the mRNA expression levels of *DDIT4*, *FOXD1*, *MINPP1*, *PAWR,* and *PFN2* were significantly suppressed in *miR-30e-5p*-transfected HNSCC cells.

**Figure 9.** Regulation of the expression of nine genes by *miR-30e-5p* in HNSCC cells. qRT-PCR showing significantly reduced expression of *FOXD1* mRNA (top, middle) 72 h after *miR-30e-5p* **Figure 9.** Regulation of the expression of nine genes by *miR-30e-5p* in HNSCC cells. qRT-PCR showing significantly reduced expression of *FOXD1* mRNA (top, middle) 72 h after *miR-30e-5p* transfection in Sa3 and SAS cells (N.S.: not significant compared to control group).

#### transfection in Sa3 and SAS cells (N.S.: not significant compared to control group). *3.6. Incorporation of FOXD1 mRNA into the RNA-Induced Silencing Complex (RISC) and Direct Control of FOXD1 Expression by miR-30e-5p in HNSCC Cells*

Considering the multivariate analysis, the correlation analysis, and the results of qRT-PCR, *FOXD1* was of the greatest interest among the nine targets. Thus, we focused on *FOXD1* as the target gene of *miR-30e-5p* in this study. *3.6. Incorporation of FOXD1 mRNA into the RNA-Induced Silencing Complex (RISC) and Direct Control of FOXD1 Expression by miR-30e-5p in HNSCC Cells*  To confirm incorporation of *FOXD1* mRNA into RISC, RIP assays were performed (Figure 10A,B). The schematic illustration displays the concept of RIP assays. Ago2-bound miRNA and mRNA were isolated by the immunoprecipitation of Ago2, the protein that plays a central role in RISC (Figure 10A). Using samples isolated by immunoprecipitation, qRT-PCR showed that the *FOXD1* mRNA level was significantly higher than those of mock or miR control transfected cells (*p* < 0.01; Figure 10B), suggesting a significant incorporation into RISC.

To confirm incorporation of *FOXD1* mRNA into RISC, RIP assays were performed (Figure 10A,B). The schematic illustration displays the concept of RIP assays. Ago2-bound miRNA and mRNA were isolated by the immunoprecipitation of Ago2, the protein that plays a central role in RISC (Figure 10A). Using samples isolated by immunoprecipitation, qRT-PCR showed that the *FOXD1* mRNA level was significantly higher than those of mock or miR control transfected cells (*p* < 0.01; Figure 10B), suggesting a significant To confirm that *miR-30e-5p* bound directly to the 30 -UTR of *FOXD1*, a dual-luciferase reporter assay was performed. Luciferase activity was significantly reduced following co-transfection with *miR-30e-5p* and a vector containing the *miR-30e-5p*-binding site in the 3 0 -UTR of *FOXD1* (Figure 10C). In contrast, co-transfection with a vector containing the *FOXD1* 3 0 -UTR, in which the *miR-30e-5p*-binding site was deleted, resulted in no change in luciferase activity (Figure 10C).

To confirm that *miR-30e-5p* bound directly to the 3′-UTR of *FOXD1*, a dual-luciferase reporter assay was performed. Luciferase activity was significantly reduced following cotransfection with *miR-30e-5p* and a vector containing the *miR-30e-5p*-binding site in the 3′- UTR of *FOXD1* (Figure 10C). In contrast, co-transfection with a vector containing the

incorporation into RISC.

in luciferase activity (Figure 10C).

**Figure 10.** Isolation of RISC-incorporated *FOXD1* mRNA by Ago2 immunoprecipitation. Direct regulation of *FOXD1* expression by *miR-30e-5p* in HNSCC cells. (**A**) Schematic illustration of RIP assay. (**B**) qRT-PCR suggested *FOXD1* mRNA was significantly incorporated into RISC. (**C**) TargetScan database analysis predicting putative *miR-30e-5p*-binding sites in the 3′-UTR of *FOXD1* (upper panel). Dual-luciferase reporter assays showed reduced luminescence activity after cotransfection of the wild-type vector and *miR-30e-5p* in Sa3 cells (lower panel). Normalized data were calculated as the *Renilla/Firefly* luciferase activity ratio (N.S.: not significant compared to control group). **Figure 10.** Isolation of RISC-incorporated *FOXD1* mRNA by Ago2 immunoprecipitation. Direct regulation of *FOXD1* expression by *miR-30e-5p* in HNSCC cells. (**A**) Schematic illustration of RIP assay. (**B**) qRT-PCR suggested *FOXD1* mRNA was significantly incorporated into RISC. (**C**) TargetScan database analysis predicting putative *miR-30e-5p*-binding sites in the 30 -UTR of *FOXD1* (upper panel). Dual-luciferase reporter assays showed reduced luminescence activity after co-transfection of the wild-type vector and *miR-30e-5p* in Sa3 cells (lower panel). Normalized data were calculated as the *Renilla/Firefly* luciferase activity ratio (N.S.: not significant compared to control group).

*FOXD1* 3′-UTR, in which the *miR-30e-5p*-binding site was deleted, resulted in no change

### *3.7. Expression of FOXD1 in HNSCC Clinical Specimens*

*3.7. Expression of FOXD1 in HNSCC Clinical Specimens*  HNSCC clinical specimens displayed moderate immunoreactivity in the cytoplasm (Figure 11B,D), whereas normal epithelium showed no expression of *FOXD1* (Figure HNSCC clinical specimens displayed moderate immunoreactivity in the cytoplasm (Figure 11B,D), whereas normal epithelium showed no expression of *FOXD1* (Figure 11A,C). The clinical features of HNSCC specimens are summarized in Table S3.

### *3.8. Effects of FOXD1 Knockdown on the Proliferation, Migration, and Invasion of HNSCC Cells*

11A,C). The clinical features of HNSCC specimens are summarized in Table S3.

To assess the tumor-promoting effect of *FOXD1* in HNSCC cells, we performed knockdown assays using siRNAs.

First, the inhibitory effects of two different siRNAs targeting *FOXD1* (si*FOXD1*-1 and si*FOXD1*-2) expression were examined. The *FOXD1* mRNA levels were effectively inhibited by each siRNA (Figure S2).

The knockdown of *FOXD1* slightly inhibited Sa3 and SAS cell proliferation (Figure 12A). In contrast, cell migration and invasion were significantly inhibited after si*FOXD1*-1 and si*FOXD1*-2 transfection into Sa3 and SAS cells (Figure 12B,C). Photographs of typical results from the migration and invasion assays are shown in Figure S3.

**Figure 11.** Immunohistochemical staining of *FOXD1* in HNSCC clinical specimens. Weak to moderate immunoreactivity of *FOXD1* was observed in the cancer lesions (**B**,**D**) whereas negative **Figure 11.** Immunohistochemical staining of *FOXD1* in HNSCC clinical specimens. Weak to moderate immunoreactivity of *FOXD1* was observed in the cancer lesions (**B**,**D**) whereas negative immunoreactivity was shown in normal mucosa (**A**,**C**). *Genes* **2022**, *13*, 1225 16 of 22

a membrane culture system 48 h after seeding miRNA-transfected cells into the chambers. (**C**) Cell invasion assessed by Matrigel invasion assays 48 h after seeding miRNA-transfected cells into

We performed gene set enrichment analysis (GSEA) to identify genes that were differentially expressed in the high and low expression groups using TCGA-HNSCC data. A GSEA analysis of *FOXD1* showed that the most enriched molecular pathway in the high expression groups of *FOXD1* was "epithelial-mesenchymal transition" (Table 3, Figure 13). Additional pathways (MYC targets, TNFα signaling, Hypoxia, E2F targets, Glycolysis, and DNA repair) were also associated with groups expressing high levels of *FOXD1* (Table 3). The aberrant expression and activation of these pathways were closely associated with the downregulation of *miR-30e-5p* in HNSCC cells and contributed to

**Enriched Gene Sets in High** *FOXD1* **Expression Group** 

Epithelial mesenchymal transition 3.674 *q* < 0.001

**Name Normalized Enrichment Score FDR** *q***-Value** 

*MYC* targets V1 3.412 *q* < 0.001 *TNFα* signaling via NFκB 3.262 *q* < 0.001

Hypoxia 3.044 *q* < 0.001 E2F targets 2.988 *q* < 0.001 Glycolysis 2.797 *q* < 0.001 DNA repair 2.602 *q* < 0.001

*3.9. FOXD1-Mediated Molecular Pathways in HNSCC Cells* 

chambers.

HNSCC oncogenesis.

**Table 3.** Results of gene set enrichment analysis.

assessed by XTT assay 72 h after siRNA transfection. (**B**) Cell migration assessed using a membrane culture system 48 h after seeding miRNA-transfected cells into the chambers. (**C**) Cell invasion assessed by Matrigel invasion assays 48 h after seeding miRNA-transfected cells into chambers.

### *3.9. FOXD1-Mediated Molecular Pathways in HNSCC Cells*

We performed gene set enrichment analysis (GSEA) to identify genes that were differentially expressed in the high and low expression groups using TCGA-HNSCC data. A GSEA analysis of *FOXD1* showed that the most enriched molecular pathway in the high expression groups of *FOXD1* was "epithelial-mesenchymal transition" (Table 3, Figure 13). Additional pathways (MYC targets, TNFα signaling, Hypoxia, E2F targets, Glycolysis, and DNA repair) were also associated with groups expressing high levels of *FOXD1* (Table 3). The aberrant expression and activation of these pathways were closely associated with the downregulation of *miR-30e-5p* in HNSCC cells and contributed to HNSCC oncogenesis.


**Table 3.** Results of gene set enrichment analysis.

**Figure 13.** *Cont*.

**4. Discussion** 

the *miR-199* family and the *miR-216* family [14,23,30].

**Figure 13.** Pathways enriched among the differentially expressed genes in the high *FOXD1* expression group according to gene set enrichment analysis. The nine significantly enriched pathways (top 9) are shown. Most enriched pathway was "Epithelial Mesenchymal Transition".

In the human genome, numerous miRNAs have very similar sequences and therefore constitute a "family". Because the seed sequences of the miRNA family are the same, each miRNA family controls the expression of the same set of genes. Therefore, aberrant expression of the miRNA family members will cause the disruption of intracellular molecular networks, and these events can initiate the transformation of normal cells to cancer cells. Based on our miRNA signatures, we previously focused on several miRNA families and identified the molecular pathways that were controlled by the *miR-29* family,

Unfolded protein response 2.275 0.001

Inflammatory response 1.753 0.014 Apical junction 1.739 0.015

Oxidative phosphorylation 1.713 0.015

*IL6/JAK/STAT3* signaling 1.617 0.027

*P53* pathway 1.702 0.016

Apoptosis 1.601 0.029

G2M checkpoint 2.248 0.001 *TGFβ* signaling 2.025 0.001 KRAS signaling up 1.946 0.005 Coagulation 1.848 0.010

**Figure 13.** Pathways enriched among the differentially expressed genes in the high *FOXD1* expression group according to gene set enrichment analysis. The nine significantly enriched **Figure 13.** Pathways enriched among the differentially expressed genes in the high *FOXD1* expression group according to gene set enrichment analysis. The nine significantly enriched pathways (top 9) are shown. Most enriched pathway was "Epithelial Mesenchymal Transition".

#### pathways (top 9) are shown. Most enriched pathway was "Epithelial Mesenchymal Transition". **4. Discussion**

**4. Discussion**  In the human genome, numerous miRNAs have very similar sequences and therefore constitute a "family". Because the seed sequences of the miRNA family are the same, each miRNA family controls the expression of the same set of genes. Therefore, aberrant expression of the miRNA family members will cause the disruption of intracellular molecular networks, and these events can initiate the transformation of normal cells to In the human genome, numerous miRNAs have very similar sequences and therefore constitute a "family". Because the seed sequences of the miRNA family are the same, each miRNA family controls the expression of the same set of genes. Therefore, aberrant expression of the miRNA family members will cause the disruption of intracellular molecular networks, and these events can initiate the transformation of normal cells to cancer cells. Based on our miRNA signatures, we previously focused on several miRNA families and identified the molecular pathways that were controlled by the *miR-29* family, the *miR-199* family and the *miR-216* family [14,23,30].

cancer cells. Based on our miRNA signatures, we previously focused on several miRNA families and identified the molecular pathways that were controlled by the *miR-29* family, the *miR-199* family and the *miR-216* family [14,23,30]. Our RNA sequence-based miRNA signatures, including those in HNSCC, showed that some *miR-30* family members were frequently downregulated in cancer tissues. Those observations suggest that these miRNAs control targets with pivotal roles in cancer progression, metastasis, and drug-resistance [31]. Here, we first investigated the expression levels and clinical significance of all members of the *miR-30* family using TCGA database. The family consists of 12 miRNAs: guide and passenger strands of *miR-30a*, *miR-30b*, *miR-30c-1*, *miR-30c-2*, *miR-30d*, and *miR-30e*). Among them, the expression levels of *miR-30e-5p* and *miR-30c-1-3p* were closely associated with HNSCC molecular pathogenesis. Identifying molecular networks controlled by these miRNAs is important for understanding HNSCC oncogenesis.

In this study, we focused on *miR-30e-5p* to investigate the targets that it controls in HNSCC. The seed sequences of the guide strands of the *miR-30a* family are identical, suggesting that these miRNAs may control common targets in a sequence-dependent manner. Previous studies reported that *miR-30e-5p* had tumor-suppressive functions in several types of cancers [32–34]. A previous study showed that *miR-30e-5p* suppressed astrocyte elevated gene-1 (*AEG-1*), an oncogene that contributes to angiogenesis, metastasis, and EMT processes [35]. Another study showed that *miR-30e-5p* was a direct transcriptional target of P53 in colorectal cancer. Expression of *miR-30e-5p* blocked tumor cell migration, invasion, and in vivo metastasis by directly controlling integrin molecules [36].

A single miRNA controls a large number of genes, and targeted genes can differ depending on the type of cancer. Defining miRNA-targeted genes is an important focus in miRNA research. Our target identification strategy successfully identified nine genes that significantly predicted 5-year survival frequencies of HNSCC patients. Unfortunately, none of those genes were independent prognostic factors contributing to 5-year overall survival rates. According to multivariate analysis, a high expression level of *FOXD1* was an independent prognostic factor (HR: 1.374, 95% CI: 1.002–1.890, *p* = 0.049). Furthermore, the negative correlation between the expression of *FOXD1* and *miR-30e-5p* was confirmed. For those reasons, we focused on *FOXD1* as the target gene in this study. Future studies may examine the other genes in the belief that resultant data will improve our understandings of the molecular pathogenesis of HNSCC.

For example, *FXR1* is an RNA-binding protein that regulates co-transcriptional and post-transcriptional gene expression. It is a member of the Fragile X-mental retardation (*FXR*) family of proteins, which includes *FMR1* and *FXR2* [37]. Recent studies showed that the overexpression of *FRX1* was observed in several types of cancers, including HNSCC, and its expression correlates with poor prognosis of the patients [38]. In oral cancer cells, *FXR1* stabilized *miR-301a-3p* and *miR-301a-3p*, both of which target *p21* [39]. Another study showed that the overexpression of *FXR1* facilitates the bypass of senescence and tumor progression [40]. Aberrant expression of *FXR1*-enhanced HNSCC progression and its expression is closely associated with the molecular pathogenesis of HNSCC.

*PFN2* was identified as a *miR-30e-5p* target in this study. *PFN2* is a member of the profilin family, namely, profilin 1, 2, 3, and 4 [41]. Profilins are actin-binding proteins involved in the regulation of cytoskeletal dynamics. Overexpression of *PFN2* was reported in several types of cancers, including HNSCC [42]. Recently, we revealed that *PFN2* was directly regulated by tumor-suppressive *miR-1*/*miR-133* clustered miRNAs in HNSCC cells, and its overexpression promoted cancer cell malignant transformation [25]. Another study showed that overexpression of *PFN2* enhanced the aggressiveness of small cell lung cancer (SCLC), including cell proliferation, migration, and invasion. Knockdown of *PFN2* decreased the cells' aggressive nature. Moreover, a mouse xenograft model demonstrated that the overexpression of *PFN2* dramatically elevated SCLC growth and vasculature formation as well as lung metastasis [43]. Notably, a previous study showed that *miR-30a-5p* suppressed the expression of *PFN2* in lung cancer cell lines [44].

In this study, we confirmed that the overexpression of *FOXD1* facilitated cancer cell malignant transformation, e.g., enhanced migratory and invasive abilities in HNSCC cells. *FOXD1* is a member of the forkhead family of transcription factors [45]. *FOXD1* expression is upregulated in several types of cancers [46]. Moreover, the knockdown of *FOXD1* impairs the colony-forming abilities of oral cancer cells after radiation treatment [47]. Another study showed that upregulated *FOXD1* contributed to melanoma cells' resistance to vemurafenib via the recruited expression of connective tissue growth factor [48]. In gastric cancer (GC) cells, *FOXD1*-AS1 expression induced *FOXD1* translation, and these events enhanced GC cell progression and cisplatin resistance [49]. Importantly, *miR-30a-5p* directly binds to the 3 0 -UTR of *FOXD1* and suppresses its expression in several types of cancer cells, e.g., lung squamous cell carcinoma, pancreatic ductal adenocarcinoma, osteosarcoma, and ovarian cancer [50–53]. Our present data and previous studies indicate that the overexpression of *FOXD1* is closely involved in aggressive cancer cell transformation in a wide range of cancers, including HNSCC.

These studies show that searching for target genes of tumor-suppressive miRNAs is an attractive strategy for exploring the molecular mechanisms of HNSCC.

### **5. Conclusions**

We showed that the *FOXD1* gene is directly controlled by tumor-suppressive *miR-30e-5p* in HNSCC cells. The aberrant expression of *FOXD1* facilitated cancer cell migration and invasion, and its expression was closely associated with the prognosis of HNSCC patients. Our strategy (analysis of tumor-suppressive miRNAs and their controlled genes) can identify genes that are deeply involved in the molecular pathogenesis of HNSCC.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/genes13071225/s1, Figure S1: Photomicrographs of cells subjected to migration and invasion assays. Figure S2: Efficiency of siRNA mediated *FOXD1* knockdown in HNSCC cell lines (Sa3 and SAS cells). Figure S3: Photomicrographs of cells subjected to migration and invasion assays. Table S1: Clinical features of the patient in GSE180077. Table S2: Used reagents in this study. Table S3: Clinical features of HNSCC patients for IHC.

**Author Contributions:** Conceptualization, N.S., C.M and S.A.; methodology, N.S., C.M. and S.A.; software, S.A.; validation, C.M. and S.A.; formal analysis, C.M. and S.A.; investigation, C.M., S.A. and A.K. (Ayaka Koma); resources, M.K. and N.K.; data curation, C.M. and S.A.; writing—original draft preparation, N.S., C.M. and S.A.; writing—review and editing, A.K. (Atsushi Kasamatsu); visualization, C.M. and S.A.; supervision, T.H. and K.U.; project administration, N.S.; funding acquisition, N.S., A.K. (Atsushi Kasamatsu), M.K., N.K. and K.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by JSPS KAKENHI (grant numbers: 21K09577, 22K09679, 21K09367, 22H03286, 20H03883, 21KK0159).

**Institutional Review Board Statement:** The study was conducted according to the Declaration of Helsinki and approved by the Bioethics Committee of Chiba University (approval number: 915, 11 July 2021).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Our expression data were deposited in the GEO database (accession number: GSE180077).

**Acknowledgments:** The results shown here are in part based upon data generated by TCGA Research Network: https://www.cancer.gov/tcga (accessed on 20 August 2020).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


### *Article*
