Identification of therapeutic targets and prognostic biomarkers from the hnRNP family in invasive breast carcinoma

Expression levels of the hnRNP family members in patients with BRCA

The mRNA expression patterns of the hnRNPs were determined using the Oncomine and The Cancer Genome Atlas (TCGA) databases. Oncomine analysis revealed that there were significant changes in hnRNP expression at the transcriptional level between different types of cancer and normal tissues following the threshold search criteria (P value < 0.05; fold change > 2, and top 10% gene rank). The significant unique analyses between cancer and normal tissues that met our selection criteria in pan-cancer analysis are shown in Figure 1A and Supplementary Table 1. HNRNPA0, HNRNPA1, HNRNPA2B1, PCBP1, HNRNPL, HNRNPM, SYNCRIP, and HNRNPU were significantly upregulated in the two main subtypes of BRCA (IDC and ILC) compared with in normal tissues in the datasets by Gluck et al. [20], Radvanyi et al. [21], Ma et al. [22] and Turashvili et al. [23]. As shown in Figure 1B, in the TCGA cohort, HNRNPA0, HNRNPA2B1, HNRNPC, HNRNPD, PCBP1, HNRNPF, HNRNPK, HNRNPL, HNRNPM, PTBP1, SYNCRIP, HNRNPR, and HNRNPU were found to be upregulated in 1,093 invasive BRCA tissue specimens compared with in 112 normal tissue specimens from the TCGA cohort (P < 0.05). The results of Radvanyi et al. showed that the mRNA expression levels of HNRNPA1, HNRNPL, and HNRNPR were notably increased in invasive mixed BRCA tissues than in normal tissues. Radvanyi et al. and Ma et al. revealed that HNRNPA, PCBP1, HNRNPL, and SYNCRIP were upregulated in ductal BRCA tissues in situ than in normal tissues. However, Finak et al. [24] reported significantly lower expression levels of HNRNPA0, HNRNPA1, HNRNPA2B1, HNRNPC, HNRNPD, PCBP1, HNRNPF, RBMX, HNRNPK, and HNRNPU in BRCA-associated stroma than in normal tissues. The results of Curtis from the European Genome-phenome Archive demonstrated that HNRNPU was downregulated in invasive BRCA tissues than in normal tissues. These findings suggest that the expression of hnRNP family members is associated with patients in different areas or different histological types of BRCA.

Figure 1. Aberrant mRNA expression of hnRNPs in BRCA. (A) mRNA expression of hnRNPs in various cancer types based on the Oncomine database. The Figure shows the number of datasets with statistically significant mRNA overexpression (red) or downregulated expression (blue) of hnRNPs in different types of cancer versus normal tissues. The threshold was designed with a P value of 0.05, fold change of 2, and gene ranking of 10%. The gene-rank percentile was analyzed for the top target gene from all genes measured in each research. Cell color was determined by the best gene-rank percentile for analysis within the cell. (B) Box plots showing the distribution of the expression of the hnRNPs across BRCA and normal tissue samples in the TCGA-BRCA dataset.*P < 0.05; **P < 0.01; ***P < 0.001.

Correlation of hnRNP mRNA expression with individual stages of cancer, lymph node metastases, and pathological subtypes in TCGA-BRCA

We evaluated the relationship between the mRNA expression levels of hnRNPs and the different clinicopathological subgroups of patients in the TCGA-BRCA cohort. As shown in Figure 2, significant differences between the groups were assessed to generate P values and then illustrated using heat maps. These data indicated that the upregulated mRNA expression levels of HNRNPA0, HNRNPA2B1, HNRNPC, HNRNPD, PCBP1, HNRNPF, PTBP1, HNRNPK, HNRNPL, HNRNPM, SYNCRIP, HNRNPR, and HNRNPU were significantly associated with high cancer stages and lymph node metastases (LNM) with few exceptions. However, the mRNA expression level of HNRNPA1 was not associated with cancer stages and LNM in patients with BRCA (P > 0.05; Figure 2A, 2B). Except for HNRNPA1 and RBMX, the expression levels of the remaining members were higher in IDC tissues than in normal tissues (P < 0.0001; Figure 2C). In addition, except for HNRNPA1, HNRNPD, and RBMX, the expression levels of the remaining 12 members were higher in ILC tissues than in normal tissues (P < 0.05; Figure 2C). These results suggest that the mRNA expression levels of HNRNPA0, HNRNPA1, HNRNPA2B1, HNRNPC, HNRNPD, HNRNPF, HNRNPK, HNRNPL, SYNCRIP, HNRNPR, and HNRNPU are markedly higher in IDC tissues than in ILC tissues (P< 0.001; Figure 2C), indicating that hnRNPs play significant roles in the tumorigenesis and progression of BRCA, particularly IDC and ILC. However, owing to the limited sample size of some histological groups (mixed, other, mucinous, medullary, or metaplastic), we did not present an in-depth analysis of their subgroup comparison data.

Figure 2. The relationship between the mRNA expression of hnRNPs and tumor stage, nodal metastasis status, and pathological subtypes of patients with TCGA-BRCA. (A) Heat map displaying the correlation between the expression of hnRNPs and tumor stage. The boxplot shows the correlation between HNRNPU expression and tumor stages. (B) Heat map displaying the association between the expression of hnRNPs and nodal metastasis status. The boxplot shows the association between HNRNPU expression and the nodal metastasis status. (C) Heat map displaying the correlation between the expression of hnRNPs and tumor pathological subtypes. The boxplot shows the association between hnRNPU expression and tumor pathological subtypes. Gray represents *P < 0.05, blue represents **P < 0.01, green represents ***P < 0.001, red represents ****P < 0.0001, and white represents no significant difference.

Prognostic value of hnRNPs in BRCA

To date, the potential prognostic value of only few members of the hnRNP family have been investigated in certain cohorts. Therefore, to analyze the independent prognostic value of the hnRNP family members, we analyzed the prognostic role of 15 hnRNP family members in patients with BRCA using the Kaplan–Meier (KM) method. As shown in Figure 3A, 3B, 10 of the 15 tested genes significantly correlated with OS (P < 0.05). Among these genes, HNRNPD, HNRNPL, and SYNCRIP were risky genes, with a hazard ratio (HR) of >1 and log-rank P value of <0.05. In contrast, HNRNPA0, HNRNPA2B1, HNRNPC, PCBP1, HNRNPF, RBMX, and HNRNPR were protective genes, with an HR of <1 and log-rank P value of <0.05. In addition, the mRNA expression levels of HNRNPA0, HNRNPA2B1, HNRNPC, PCBP1, RBMX, PTBP1, HNRNPK, HNRNPL, SYNCRIP, HNRNPR, and HNRNPU were associated with recurrence-free survival (RFS) in patients with BRCA (P < 0.05; Figure 3C). In particular, among these 11 genes, increased mRNA expression levels of HNRNPA0, PTBP1, HNRNPK, and SYNCRIP were associated with poor RFS in patients with BRCA under the selected best cutoff value (P < 0.05; Figure 3D).

Figure 3. Correlation between hnRNP mRNA expression and clinical prognosis in patients with BRCA. (A) Forest plots showing the correlation between the expression of hnRNPs and patient OS, as examined via Kaplan–Meier analysis. (B) High hnRNPD, hnRNPL, or SYNCRIP expression was significantly associated with poor OS. (C) Forest plots showing the association between the mRNA expression levels of hnRNPs and RFS. (D) High expression of HNRNPA0, PTBP1, HNRNPK, and SYNCRIP was significantly associated with poor RFS.

Despite recent advances in BRCA therapy, metastasis remains the leading cause of mortality in patients with BRCA. Subgroup analysis revealed that the increased expression of HNRNPU significantly indicated a poor OS in patients with BRCA with LNM (HR, 2.11; P = 0.029), However, the increased expression of HNRNPA0, HNRNPC, HNRNPF, and HNRNPK predicted improved survival in patients with BRCA with LNM (HR, <1; P < 0.05). RFS data analysis revealed that the increased expression levels of HNRNPA0, HNRNPA2B1, HNRNPC, HNRNPD, HNRNPF, HNRNPK, HNRNPR, and HNRNPU were associated with improved RFS in patients with BRCA with LNM (HR, <1; P < 0.05). In contrast, the increased expression of PTBP1 (HR, 1.25; P = 0.032), HNRNPL (HR, 1.73; P <.0001), HNRNPM (HR, 1.42; P = 0.01), and SYNCRIP (HR, 1.84; P <.0001) were associated with poor RFS in patients with BRCA with LNM (Supplementary Table 2).

Protein expression levels of hnRNPs in BRCA and their association with the clinicopathological features of patients with BRCA

Because HNRNPA1 and RBMX were not altered in 1,093 patients with invasive BRCA compared with in 112 normal tissues from the TCGA cohort, the mRNA expression level of PCBP1 is not associated with LNM in BRCA. Therefore, we further identified the protein expression levels of the other 12 genes in patients with BRCA using the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database. Our results indicated that HNRNPA2B1, HNRNPC, HNRNPD, HNRNPF, PTBP1, HNRNPK, HNRNPL, HNRNPM, SYNCRIP, and HNRNPU were upregulated in BRCA tissues compared with in normal tissues (Supplementary Figure 1), which corroborated the mRNA expression results detailed in the previous section of this manuscript. We analyzed the association between the protein expression profile and clinicopathological data and found that increased protein expression levels of HNRNPA2B1, HNRNPC, HNRNPD, HNRNPF, PTBP1, HNRNPK, HNRNPL, HNRNPM, and HNRNPU were associated with high cancer stages (Figure 4A). Further, we found that the protein expression levels of HNRNPA2B1, HNRNPC, HNRNPD, HNRNPF, PTBP1, HNRNPK, HNRNPL, HNRNPM, SYNCRIP, and HNRNPU were higher in IDC tissues than in normal tissues and that the protein expression levels of HNRNPC, HNRNPF, PTBP1, HNRNPL, and HNRNPM were higher in ILC tissues than in normal tissues. In addition, the expression levels of HNRNPC, HNRNPF, HNRNPL, HNRNPM, SYNCRIP, HNRNPR, and HNRNPU were significantly higher in IDC tissues than in ILC tissues (Figure 4B; P < 0.05).

Figure 4. Correlation between hnRNP protein expression, tumor stage, and pathological subtypes. (A) Heat map displaying the correlation between hnRNP protein expression and tumor stage. The boxplot shows the association between HNRNPU protein expression and tumor stage. (B) Heat map displaying the correlation between hnRNP protein expression and tumor pathological subtypes. The boxplot shows the association between HNRNPU protein expression and tumor pathological subtypes. Gray represents *P < 0.05, blue represents **P < 0.01, green represents ***P < 0.001, red represents ****P < 0.0001, and white represents no significant difference.

Genetic mutation, interaction, and coexpression analyses of hnRNPs

The different mutation spectra seen in the geographical and ethnic populations can be used to identify the environmental exposure contributing toward BRCA development. We considered their possible roles in BRCA risk and their association with the type of cancer developed. Therefore, we calculated the alteration frequency of hnRNPs across the TCGA-BRCA cohort using the cBioPortal online tool. hnRNPs were altered in 595 samples from 1065 patients (60%) in the TCGA-BRCA cohort (Figure 5A). Notably, HNRNPU, HNRNPL, and HNRNPA2B1 were the three most frequently altered genes (14%, 2.9%, and 2.5%, respectively), and amplification, mutation, and deep deletion were the three most common types of mutations. We also confirmed the alteration frequency of hnRNPs in the two other datasets (Metastatic Breast Cancer Project Provisional and Breast Cancer; METABRIC) (Supplementary Figure 2A, 2B). Similar to the results of the TCGA-BRCA cohort, HNRNPU was the most frequently altered gene in these two datasets (13% and 23%, respectively). Next, we explored the potential coexpression genes of hnRNPs. Pearson’s correlation analysis revealed the following significant correlations: HNRNPA0 with PTBP1 and HNRNPM; HNRNPA2B1 with HNRNPC, HNRNPD, HNRNPL, and HNRNPU; HNRNPC with HNRNPF, HNRNPK, HNRNPL, and HNRNPU; HNRNPD with HNRNPL; HNRNPF with HNRNPK; PTBP1 with HNRNPL and HNRNPM; and HNRNPK with HNRNPR and HNRNPU (R > 0.4; N = 1093; Figure 5B). These results indicate that several hnRNPs interact with each other during the pathogenesis of BRCA. We conducted a protein–protein interaction (PPI) network analysis of hnRNPs to present the interaction among the 12 genes using STRING. As shown in Figure 5C, the purple connecting line shows the correlation among the genes when experimentally determined. We then constructed another network for hnRNPs and the most frequently altered coexpressed genes. The results showed that 43 genes, including BRCA2, TP53, ZNF532, CYLC1, COL12A1, SEMA5A, PLXNB3, and PCM1, were closely related with alterations in hnRNPs. Using Cytoscape, we found that of the 12 hnRNPs, HNRNPU, SYNCRIP, HNRNPL, and PTBP1 were the hub genes in the network (Figure 5D). Next, we studied the potential biological processes of the top 50 genes that positively and negatively correlated with HNRNPU mutation in BRCA using Linkomics (Supplementary Figure 2C–2E). The Gene Ontology (GO) term functional enrichment analyses of these genes are summarized (P < 0.01) in Supplementary Figure 2F–2H. Among the GO terms for biological processes, the top five terms were “protein binding,” “ion binding,” “molecular transducer activity,” “hydrolase activity,” and “transporter activity.”

Figure 5. Genetic mutation, interaction, and coexpression analyses of hnRNPs in patients with BRCA. (A) Summary of the alteration frequencies of nodal metastasis-related hnRNPs in patients with TCGA-BRCA (Breast Invasive Carcinoma, TCGA). The number of patients for each BRCA subtype is indicated at the top of each column and each type of alteration, including deletions (blue), amplification (red), multiple alterations (gray), or mutations (green), is shown using the cBioPortal. (B) Correlation plot based on Spearman’s correlation test results to show the correlation of gene expression among the 12 hnRNP family members in BRCA. (C) PPI among hnRNPs was derived from STRING. The purple connection indicates the protein interaction determined by the experiment and the blue connection indicates that the proteins are homologous. (D) The 12 hnRNP family members and coexpressed gene networks were screened using Cytoscape. Red and yellow represent the top 10 hub genes.

Immunohistochemistry (IHC)-predicted patterns for potential clinical applications

Because HNRNPU and SYNCRIP were identified as the core molecular genes and SYNCRIP positively correlated with OS and HNRNPU positively correlated with tumor stage and metastasis, we further validated the IHC pattern from the Human Protein Atlas database to simulate clinical applications on the basis of pathology. Several IHC images of hnRNP proteins displayed a trend toward differential expression in BRCA subtype samples (IDC or ILC) compared with in normal tissue samples. These results were in accordance with the transcriptional level of hnRNPs. SYNCRIP staining showed negative or low staining in the cytoplasm and cell membranes of normal tissues, and HNRNPU staining showed negative staining in the cytoplasm and cell membrane. In BRCA (IDC and ILC) tissues, strong and moderate staining patterns of SYNCRIP were clearly located in the nucleus, cytoplasm, and cell membrane. For HNRNPU, moderate and strong patterns were mainly located in the nucleus, cytoplasm, and cell membrane. This unique IHC staining of BRCA could help distinguish cancer tissues from normal tissues as well as predict clinical outcomes (Figure 6A, 6B). The relative staining intensity of SYNCRIP was negative (two cases) for normal tissues; weak (three cases), moderate (four cases), and strong (one case) for IDC tissues; and weak (one case) and moderate (two cases) for ILC tissues (Figure 6C). For HNRNPU staining, both cases were weak for normal tissue; moderate (two cases) and strong (six cases) for IDC tissues; and moderate (one case) and strong (three cases) for ILC tissues (Figure 6D).

Figure 6. The IHC expression pattern of SYNCRIP and HNRNPU in normal breast tissues and BRCA. (A, B) Representative images are shown for strong, moderate, weak and negative expression of SYNCRIP and HNRNPU in normal breast tissues and BRCA (IDC and ILC). The black rectangle indicates a higher magnification of the indicated area in the image. (C) Bar chart of the IHC staining intensity of SYNCRIP for BRCA (a total of 13 cases). (D) Bar chart of the IHC staining intensity of hnRNPU for BRCA (a total of 14 cases).

HNRNPU promoted breast cancer cell proliferation, migration, and invasion

Because HNRNPU was discovered as a key gene in BRCA among the hnRNP family members and was associated with nodal metastasis in BRCA, we speculated that HNRNPU is involved in the regulation of BRCA metastasis. To test this hypothesis, we first investigated the expression of HNRNPU in 32 breast IDC and adjacent normal tissues. The expression level of HNRNPU was higher in IDC tissues than in normal tissues (Figure 7A). We then examined the expression level of HNRNPU in a series of BRCA cell lines (MCF-7, ZR-75-1, MDA-MB-330, and MDA231-LM2) and a normal breast epithelial cell line (MCF-10A). The expression level of HNRNPU was higher in the BRCA cell lines than in MCF-10A cells (Figure 7B). To investigate the functional roles of HNRNPU in invasive BRCA progression, HNRNPU expression was knocked down using shRNAs in ZR-75-1 and MDA-MB-330 cells (Figure 7C, 7D). MTT assay revealed that HNRNPU knockdown substantially reduced the proliferation ability of ZR-75-1 and MDA-MB-330 cells (Figure 8A, 8B). Moreover, the wound-healing and Matrigel invasion assays revealed that HNRNPU knockdown significantly inhibited the migration and invasion abilities of ZR-75-1 and MDA-MB-330 cells (P < 0.05, Figure 8C–8H). Further analysis of HNRNPU knockdown in ZR-75-1 and MDA-MB-330 cells revealed that HNRNPU depletion can increase the expression level of the epithelial marker E-cadherin and decrease the expression levels of the mesenchymal markers vimentin, N-cadherin, and matrix metallopeptidase 9, suggesting that HNRNPM is important for the EMT phenotype (Figure 8I–8L). Taken together, our results suggest that HNRNPU knockdown significantly inhibits the proliferation, migration, and invasion of BRCA cells in vitro. Therefore, we reached the tentative conclusion that HNRNPU is the key molecule that plays a fundamental role in BRCA metastasis.

Figure 7. HNRNPU expression is upregulated in BRCA cell lines and clinical tumor samples. Altering expression of HNRNPU proteins in clinical tumor samples (A) and four human BRCA cell lines (B), as determined via western blot analysis. Glyceraldehyde 3-phosphate dehydrogenase was used as the control. ZR-75-1 and MDA-MB-330 cells were transfected with shRNA and the mRNA (C), and protein (D) expression levels of HNRNPU were measured via reverse transcription quantitative real-time polymerase chain reaction and western blotting. *P < 0.05; **P < 0.01; ***P < 0.001.

Figure 8. HNRNPU improved metastatic ability of cells by activating EMT phenotypes in BRCA cells. (A, B) Effects of HNRNPU knockdown on the proliferation of ZR-75-1 and MDA-MB-330 cells. (C–E) Representative images from the Transwell migration assays for ZR-75-1 and MDA-MB-330 cells transfected with sh-HNRNPU and sh-NT. (F–H) Representative images from the wound-healing assays for ZR-75-1 and MDA-MB-330 cells transfected with sh-HNRNPU and sh-NT. (I–L) Effects of HNRNPU knockdown on the expression of EMT signal protein detection in ZR-75-1 (I, J) and MDA-MB-330 (K, L) cells using western blotting. *P < 0.05; **P < 0.01; ***P < 0.001 compared with sh-NT.

Ethics statement

The study has been conducted in accordance with the ethical standards and the Declaration of Helsinki as well as according to national and international guidelines. All procedures were approved by the Ethics Committee of the Institute of Clinical Pharmacology. Informed consent has been obtained.

Bioinformatics analysis

During the validation process, the transcription levels of hnRNP in different types of cancer were analyzed using the Oncomine 4.5 (https://www.oncomine.org) gene expression array dataset [46] and TIMER (https://cistrome.shinyapps.io/timer/) [47]. The mutation frequency of hnRNPs in a variety of cancers was analyzed based on TCGA, and a detailed diagram was constructed for hnRNPs in BRCA samples based on mutations. CNA data were analyzed using the online database cBioPortal (https://www.cbioportal.org/) [48]. To intuitively analyze the relationship between the mRNA and protein expression levels of hnRNPs and clinicopathology, we compared the patterns of the expression levels of hnRNPs in major BRCA types and downloaded a boxplot based on these patterns from the UALCAN web portal (http://ualcan.path.uab.edu/), which is an interactive Web portal that provides gene expression analyses based on TCGA and MET500 cohort data as well as protein expression analysis based on CPTAC data [49]. The comparison between hnRNPs expression and TCGA-BRCA tumor pathological subtypes, lymph metastasis status, and major subclasses were selected to create heat maps. Additionally, LinkedOmics (http://www.linkes.org/login.php) was used to analyze the mRNA expression levels of hnRNP-coexpression genes in the TCGA-BRCA cohort based on Pearson’s correlation coefficient [50]. All results were graphically presented in a heat map. To analyze OS and RFS, we divided the patients into two groups based on the best-performing threshold using KM plotter (https://www.kmplot.com). Next, the distribution of HR with 95% confidence intervals of OS and RFS across the hnRNPs was selected to create forest plots. Further, to examine the differences in the protein expression levels of hnRNPs between IDC, ILC, and normal tissues, we obtained the immunohistochemistry data of SYNCRIP and HNRNPU from the HPA database (https://www.proteinatlas.org/) [51].

Patient samples

All samples (IDC and adjacent normal tissues) were obtained from Nanyang Second People's Hospital (Nanyang, Henan, China) between October 2019 and June 2020 (n = 32). All pathological features were confirmed by experienced pathologists, and none of the patients received preoperative anticancer treatment.

Cell culture

The human breast cell lines MCF-10A, ZR-75-1, MCF-7, MDA231-LM2, and MDA-MB-330 were obtained from the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in Dulbecco’s modified Eagle medium (Invitrogen, Shanghai, China) supplemented with 10% fetal bovine serum (Gibco, Logan, UT) and 100 U of penicillin and streptomycin at 37° C in a humidified atmosphere containing 5% CO2.

Cell transfection

The shRNA for the specific inhibition of HNRNPU expression and a negative control shRNA were synthesized by GenePharma Co., Ltd. (Shanghai, China). Exponentially growing untreated cells were plated 24 h before transfection. The plated cells were transfected with HNRNPU shRNA using 5 μL of Lipo-RNAiMAX following the manufacturer’s instructions (Invitrogen, Carlsbad, CA, USA). After treatment, the cells were processed for further analysis.

RNA extraction and quantitative real-time polymerase chain reaction

Total RNA was extracted from the cultured cells or tissue samples using the TRIzol reagent (Invitrogen) following the manufacturer’s instructions. First-strand cDNA synthesis was performed using the Prime Script RT Master Mix (TaKaRa, Dalian, China) and reverse transcription quantitative real-time polymerase chain reaction was performed using the SYBR Green PCR kit (TaKaRa) according to the manufacturer’s instructions. Glyceraldehyde 3-phosphate dehydrogenase was used as an endogenous control. The gene expression was calculated using the 2−ΔΔCt method. All data represent the average of three replicates.

Cell viability assays

One day after transfection, cells were re-seeded into 96-well plates. Cell viability was measured using a Cell Counting Kit 8 assay (Dojindo Laboratories, Kumamoto, Japan), according to the manufacturer's instructions. The absorbance value of each well was measured at 450 nm. Five replicates were set for each group and the whole experiment was repeated three times.

Wound-healing assay

Cells were seeded in 6-well plates and were transfected with the shRNA after 24 h. After culturing the cells for an appropriate amount of time, artificial wounds were made using a micropipette tip. The wounded monolayers were washed with phosphate-buffered saline (PBS) to remove cell debris. The distance between the two edges of the wound was calculated and imaged at 0 and 24 h under a microscope.

Transwell invasion assays

For the invasion assays, cells were seeded in a 24-well Corning FluoroBlok chamber precoated with Matrigel (BD Biosciences, Billerica, MA, USA). A medium containing 20% FBS was added to the lower chamber in a 5% CO₂ incubator at 37° C for 24 h. Then, the cells remaining on the lower side of the membrane were stained with 4% paraformaldehyde for 15 min, stained with 0.1% cresyl violet, washed three times with PBS, and air-dried. Five fields were randomly selected for counting the number of migrated cells and images (200× magnification) were captured under a phase-contrast microscope.

Western blot analysis

Standard Western blotting protocols were followed [18]. Whole-cell lysates were probed with antibodies specific to GAPDH (Santa Cruz, CA, USA), HNRNPU, MMP9, E-Cadherin, Vimentin, and N- Cadherin (Cell Signaling, UK).

Statistical analysis

Data are expressed as the mean ± standard error of the mean for at least three independent experiments. Analysis of variance was used to distinguish the differences in each group. P < 0.05 was considered statistically significant. All analyses were performed using GraphPad Prism version 8.0 (Graphpad Software, CA, USA).