Comprehensive Analysis of the Immune and Prognostic Implication of COL6A6 in Lung Adenocarcinoma

Collagen type VI alpha 6 chain (COL6A6), a novel collagen, has been considered as a tumor suppressor and therapeutic target in several tumors. However, the functional role of COL6A6 in immune cell infiltration and prognostic value in lung adenocarcinoma (LUAD) remains unknown. Here, we evaluated COL6A6 expression and its impact on survival among LUAD patients from The Cancer Genome Atlas (TCGA) and several other databases. COL6A6 was downregulated in LUAD tissues compared to normal tissues at both mRNA and protein levels. COL6A6 expression was negatively associated with pathological stage, tumor stage, and lymph node metastasis. High COL6A6 expression was a favorable prognostic factor in LUAD. Next, we explored the associations between COL6A6 expression and immune cell infiltration. COL6A6 expression was positively associated with the infiltration of B cells, T cells, neutrophils and dendritic cells. Additionally, the immune cell infiltration levels were associated with COL6A6 gene copy number in LUAD. Consistently, gene set enrichment analysis showed that various immune pathways were enriched in the LUAD samples with high COL6A6 expression, including pathways related to T cell activation and T cell receptor signaling. The impacts of COL6A6 on immune activity were further assessed by enrichment analysis of 50 COL6A6-associated immunomodulators. Thereafter, using Cox regression, we identified a seven-gene risk prediction signature based on the COL6A6-associated immunomodulators. The resulting risk score was an independent prognostic predictor in LUAD. Receiver operating characteristic curve analysis confirmed that the seven-gene signature had good prognostic accuracy in the TCGA-LUAD cohort and a Gene Expression Omnibus dataset. Finally, we constructed a clinical nomogram to predict long-term survival probabilities, and calibration curves verified its accuracy. Our findings highlight that COL6A6 is involved in tumor immunity, suggesting COL6A6 may be a potential immunotherapeutic target in LUAD. The proposed seven-gene signature is a promising prognostic biomarker in LUAD.


INTRODUCTION
Lung cancer is the leading cause of cancer-related death globally and its incidence continues to increase (1). Non-small cell lung cancer (NSCLC) accounts for 85% of lung cancer, and lung adenocarcinoma (LUAD) is the major NSCLC subtype (2). The 5-year overall survival (OS) rate of lung cancer is approximately 15%, mainly as a result of late diagnosis and the limited treatment options (3). Therefore, novel treatment strategies and biomarkers are urgently needed.
In recent years, immunotherapy has emerged as a promising strategy for various cancers, especially lung cancer (4). The success of immunotherapy highlights the importance of the tumor microenvironment (TME). Immune cell infiltration of the TME has been shown to be associated with immunotherapy outcomes and survival among patients with solid tumors (5). For instance, tumors with T cell infiltration have a favorable response to anti-PD-1/PD-L1 therapy (6). However, only about 20% of NSCLC cases respond to immunotherapy. Further exploration of the interplay between immune cell infiltration and tumor cells will provide insights into the treatment of LUAD.
Collagen type VI alpha 6 chain (COL6A6) is a protein-coding gene. COL6A6 belongs to the collagen VI (COL6) family which acts a vital role in the extracellular matrix. COL6A6 encodes a 2,262-amino acid protein that contains multiple von Willebrand factor domains and forms a component of the basal lamina of epithelial cells (7). This protein may regulate epithelial cell-fibronectin interactions and participate in cell adhesion (8,9). Cell adhesion has a strong association with cancer metastasis (10). COL6A6 overexpression significantly suppressed tumor growth and metastasis capacity in pituitary adenoma (11). The existence of COL6 family members in the extracellular matrix indicates the potential function of these proteins in the TME. In addition, COL6A6 was found to be highly expressed in various human organs, especially in lung. However, no previous studies have reported on COL6A6 in lung cancer.
In this study, we explored the prognostic and immune implication of COL6A6 in LUAD. We downloaded data from The Cancer Genome Atlas (TCGA) database for detailed analysis. The underlying roles and prognostic value of COL6A6 in LUAD were confirmed through multiple bioinformatics analyses. Further, we systematically evaluated the association between COL6A6 and immune cell infiltration, as well as the signaling pathways regulating the COL6A6-mediated immune response. Moreover, we identified an immune prognostic signature using the COL6A6-associated immunomodulators, and we then validated its prognostic accuracy in a LUAD dataset from the Gene Expression Omnibus (GEO) database. Finally, a nomogram was established integrating the immune signature and clinical pathological features.

Acquisition of Gene Expression and Clinical Data
The mRNA sequencing data (HTSeq-FPKM) and associated clinical information for LUAD were downloaded from TCGA database (https://portal.gdc.cancer.gov/, up to February 23, 2020) using TCGAbiolinks R package (12), which involved 522 cases and 535 tumor tissues. Nine cases without corresponding tumor samples were eliminated. Subsequent processing excluded cases with insufficient or missing data on age, sex, race, overall survival time, local invasion, lymph node metastasis, distant metastasis, and TNM stage. Normalized gene expression data in the form of transcripts per million (TPM) were log2 (TPM+1) transformed for survival analysis. Finally, 485 cases with eligible clinical information were included in the Cox regression analysis. Moreover, 535 tumor tissues were retained to perform Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORTx) analysis to figure out the influence of COL6A6 expression on the immune microenvironment. Furthermore, after creating a prognostic signature using the TCGA data, it was validated using mRNA expression data and related clinical data in the GSE26939 dataset from GEO database (http://www.ncbi.nlm.nih.gov/geo/). This dataset included 116 LUAD samples and the detailed clinical data are shown in Supplementary Table 1. Ethical approval and informed consent were not required in this study due to the data being publicly available in the TCGA and GEO databases.

Validation of COL6A6 Expression in Tissues and Cell Lines
Protein expression of COL6A6 in LUAD and normal lung tissues were evaluated based on immunohistochemistry data from the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) (16). COL6A6 expression in various LUAD cell lines were also analyzed using GENEVESTIGATOR 7.6.0, which enables analysis of deeply curated transcriptomic data from public repositories (17).

Construction of Protein-Protein Interaction Networks
PPI networks were constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (version 11.0, https://string-db.org/) (18).The minimum required interaction score was set at 0.900 (highest confidence). The results were visualized in Cytoscape software (version 3.7.2) (19). Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of genes in the PPI networks were performed using the clusterProfiler R package (20). The GO annotation involved biological process (BP), molecular function (MF) and cellular component (CC) annotation. P < 0.05 was set as the cut-off criterion.

Gene Set Enrichment Analysis
GSEA was performed to identify the biological processes influenced by the genes in the prognostic signature (21). TCGA-LUAD samples with the top 25% and lowest 25% COL6A6 expression were used as the high and low expression groups, respectively. GSEA was accomplished in GSEA software version 4.0.3 based on the Molecular Signatures Database version 7.1. C2 (curated gene sets) and C5 (GO gene sets) were searched to identify enriched KEGG pathways and GO terms associated with COL6A6 expression. P<0.05 and false discovery rate (FDR) <0.05 were considered to indicate statistical significance.

Associations Between COL6A6 and Tumor-Infiltrating Immune Cells
Tumor Immune Estimation Resource (TIMER) is a comprehensive resource for systematically analysis of immune infiltrates across diverse cancer types (https://cistrome.shinyapps.io/timer/) (22). The correlations between COL6A6 expression and the abundances of six immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were estimated using the "Gene" module. Additionally, the "SCNA" module was used to explore the correlations between somatic copy number alteration and the abundances of immune infiltrates.
We further determined the associations between COL6A6 expression and the infiltration of 22 types of tumor-infiltrating lymphocytes (TILs) in LUAD using CIBERSORTx. CIBERSORTx (https://cibersortx.stanford.edu/) is the nextgeneration version of CIBERSORT, which supports deconvolution of bulk RNA-Seq data (23). Among the 535 LUAD tumor tissues with complete gene expression data in the TCGA database, samples with the top 25% and lowest 25% COL6A6 expression were included as the high and low expression groups, respectively. One thousand permutations were used. P <0.05 was set as the criterion to select the immune infiltrates that may be affected by COL6A6 expression. In addition, the correlations between COL6A6 expression and gene markers of TILs were also analyzed using the "correlation analysis" module in GEPIA2.

Analysis of COL6A6-Associated Immunomodulators
Immunomodulators associated with COL6A6 were extracted from the integrated repository portal TISIDB (http://cis.hku. hk/TISIDB/) (24). This database aims to elucidate tumor and immune system interactions by integrating multiple heterogeneous data types. Immunostimulators and immunoinhibitors that were significantly correlated with COL6A6 expression (Spearman correlation test, P<0.05) were subjected to further analysis. More specifically, the COL6A6associated immunomodulators were used to build a PPI network and were subjected to KEGG and GO analyses, as mentioned above.

Establishment of a Prognostic Signature and Survival Analysis
A prognostic immune gene signature was identified based on the COL6A6-associated immunomodulators. Forward stepwise variable selection was performed using the Akaike Information Criterion in Cox models (25). Based on the immune gene signature, the risk score was generated as follows: is the coefficient of each gene derived from the Cox regression, and x i is the expression level of each gene. Kaplan-Meier survival curve and log-rank test were used to estimate the associations of the signature and clinical features with OS. The time-dependent receiver operating characteristic (ROC) curve analysis was performed to appraise the prognostic accuracy of the risk score using the timeROC R package (26). The TCGA-LUAD dataset was used as the training dataset and the GSE26939 dataset was used as a validation dataset.

Construction and Evaluation of a Nomogram
As a widely used predictive model, a nomogram can provide an individualized prognostic risk assessment intuitively and visually (27). Based on clinical characteristics and the risk score, we constructed a nomogram to predict the probability of 1-, 3-, and 5-year OS using the rms R package. The concordance index (Cindex) was calculated to assess the predictive accuracy of the nomogram based on a bootstrap method with 1,000 replicates. Calibration curves were plotted to compare the predicted OS with actual OS rates.

Statistical Analysis
Statistical analysis was performed in R software version 3.6.1 and IBM SPSS Statistics version 22.0. Logistic regression analysis was used to identify the associations between clinical characteristics and COL6A6 expression. Univariate and multivariate Cox regression analyses were conducted to assess clinical factors associated with OS. P<0.05 was considered statistically significant unless otherwise indicated.

Expression and Function of COL6A6
This study was carried out according to the flow chart shown in Figure 1. According to the GEPIA2 results (Figure 2A), the COL6A6 mRNA expression in LUAD tumor tissues was significantly lower than that in normal tissues (|log2(fold change) |>2, P<0.01). The differential expression was also confirmed by a meta-analysis using the LCE web portal ( Figure 2B). Furthermore, the COL6A6 protein expression was explored using HPA database. The typical immunohistochemistry result revealed downregulated COL6A6 expression in tumor tissues ( Figure 2C). The protein concentration in the plasma was 490ng/L, which was detected by mass spectrometry and estimated from spectral counts in a publicly available dataset from the PeptideAtlas (Supplementary Figure 1). In addition, the GENEVESTIGATOR analysis showed that 92.7% (51 of 55) LUAD cell lines had low COL6A6 expression ( Figure  2D). According to cBioPortal, 13% of TCGA-LUAD patients had genomic COL6A6 alterations, and missense mutation was the most common type ( Figure 2E). As shown in Figure 3A, a PPI network of COL6A6 was constructed using STRING, with the highest confidence threshold (0.900). GO annotation showed that the 11 genes in this PPI network were enriched in 29 GO terms (13 BP terms, 10 CC terms, and 6 MF terms). These terms indicated that COL6A6 mainly plays a role in the extracellular matrix ( Figure 3B). According to the KEGG analysis, COL6A6 is involved in five pathways, including "Focal adhesion" and "PI3K-Akt signaling pathway" ( Figure 3C).

Survival Analysis and Prognostic Factors
GEPIA2 survival analysis indicated that lower COL6A6 expression was significantly associated with poor OS (P=0.014, group cutoff=quartile; Figure 4A) and advanced pathological stage (P=0.015; Figure 4B). Moreover, the prognostic value of COL6A6 was also substantiated by the survival meta-analysis using LCE web portal [hazard ratio (HR)=0.83, Figure 4C]. As shown in Table 1, we performed Cox regression analysis to   Figure 4D) were significantly associated with OS. Multivariate analysis identified up-regulated COL6A6 expression and lower pathological stage as independent prognostic factors that predicted favorable OS ( Figure 4E).

Association Between COL6A6 Expression and Clinical Characteristics
We investigated the mechanism underlying the effect of COL6A6 expression in LUAD by analyzing its relationship with clinical factors. The median expression value was used to create a categorical dependent variable based on COL6A6 expression. As shown in  lymphocyte differentiation, regulation of T cell activation, and chemokine signaling pathway ( Figure 5). These results indicate that COL6A6 might play important roles in the tumor immune microenvironment.

Relationship Between COL6A6 Expression and Tumor Immune Infiltrates
Previous research has demonstrated that TILs act as independent predictors of survival and immunotherapy efficacy in lung cancer (28). Thus, we sought to determine whether COL6A6 expression was related to immune cell infiltration in LUAD. The "Gene" module of TIMER was used to roughly analyze the correlations. Upregulated COL6A6 was positively associated with the infiltration levels of all six immune cells assessed (B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages and dendritic cells) ( Figure 6A, P<0.0001). As shown in Figure  6B, with chromosome arm-level deletion of COL6A6, infiltration levels of CD8 + T cells (P= 0.027), macrophages (P=0.001), neutrophils (P<0.001), and dendritic cells (P<0.001) were significantly lower. In addition, we also found COL6A6 expression varied among different immune subtypes ( Figure 6C).
COL6A6 expression was highest in the inflammatory subtype, which indicates a better prognosis in LUAD (29).
To obtain a deeper understanding of the relationship between COL6A6 expression and tumor immune infiltrates, we further determined the proportions of 22 types of TILs in LUAD using CIBERSORTx. Among the 134 TCGA-LUAD samples in each group, 125 samples in the high expression group and 107 samples in the low expression group met the screening criterion. The differences in the proportions of the 22 subpopulations of immune cells in these two groups are shown in Figure 7.  Table 3). Taking these results together, it can be concluded that COL6A6 may play a potential role in modulating the abundance of B cells, T cells, neutrophils, and dendritic cells.

Establishment and Validation of a Seven-Gene Signature
To determine the prognostic value of the COL6A6-associated immunomodulators in LUAD, we performed a stepwise Cox multivariate regression analysis. A prognostic gene signature consisting of seven genes (CD160, CD40LG, CD86, IL10, PVR, TNFSF13B, and TNFSF14) was identified based on the TCGA-LUAD dataset. The biological functions of the seven genes and the risk coefficients are shown in Table 4. Univariate Cox regression analysis was performed to evaluate the associations between these genes and OS (Supplementary Figure 2A). The risk scores were then calculated based on the proposed formula. Patients in the training dataset were then divided into high-and low-risk groups according to the optimal risk score cut-off. The Kaplan-Meier survival curve revealed that the survival time of the high-risk group was significantly lower (P < 0.001) ( Figure  9A). The area under the curve (AUC) values of the risk score and pathological stage were 0.721 and 0.683, respectively. Combining the risk score and pathological stage, an AUC of 0.77 was achieved ( Figure 9C). The risk scores, survival status distribution of the patients, and gene expression profiles related to the seven-gene signature in the training dataset are shown in Figure 9E. In addition, univariate Cox regression analysis showed that the risk score was significantly associated with OS (HR= 2.718, 95% CI=2.143-3.448, P<0.001). As shown in Supplementary Figure 2B, the risk score was an independent predictor of prognosis in the multivariate Cox regression model (HR=2.462, 95% CI=1.905-3.183, P<0.001) after adjusting for age, gender, pathological stage, T, N, and M stage. The GSE26939 dataset was used to validate the performance of the seven-gene signature. A risk score for each patient was calculated using the same method. The optimal cut-off was also determined and patients in the validation dataset were then divided into highand low-risk groups accordingly. Consistently, the patients in the high-risk group had notably poorer prognosis (P =0.014) ( Figure  9B). In the validation dataset, the AUCs for 1-, 3-, and 5-year OS prediction based on the gene signature were 0.65, 0.625, and 0.666, respectively ( Figure 9D). The distribution of the risk scores and the gene expression profiles are shown in Figure  9F. Collectively, these results indicated that the seven-gene signature had good performance in predicting the OS of LUAD patients.

Construction and Evaluation of a Prognostic Nomogram
A prognostic nomogram to predict the 1-, 3-, and 5-year OS of LUAD patients was established based on the TCGA-LUAD dataset using Cox regression ( Figure 10A). Risk score, age, gender and stage were features that were included in the nomogram. By calculating the score of each feature for each patient, we can predict the 1-, 3-, and 5-year OS probability, contributing to personalized precision treatment. The C-index of the prognostic nomogram was 0.723. The calibration curves further revealed that the nomogram had good performance in predicting the OS of LUAD patients ( Figures 10B-D).

DISCUSSION
Lung cancer is an increasingly prevalent disease that threatens global health. Late-stage diagnosis and the limitations of classical therapies lead to a poor OS. Recent evidence suggests that immunotherapy is a powerful tool for lung cancer treatment (4)(5)(6). Therefore, the identification of novel biomarkers and promising immune-related therapeutic targets for lung cancer has become imperative in clinical practice. In this study, we found that the expression of COL6A6 was downregulated at both mRNA and protein levels in LUAD tissues. High COL6A6 expression was an independent predictor of a favorable prognosis. Additionally, COL6A6 expression was significantly associated with several clinical characteristics, including pathological stage, T stage, and N stage. Moreover, COL6A6 was found to be associated with immune cell infiltration levels and various immune pathways in LUAD. Furthermore, we identified a seven-gene risk signature based on COL6A6-associated immunomodulators using stepwise Cox regression model. This signature showed good predictive performance and a prognostic nomogram was constructed based on the signature for clinical application.
Previous studies revealed that COL6A6 expression was distinctly higher in lung tissue than in other tissues (Supplementary Figure  3) (30, 31). However, as a new member of the COL6 family, the role of COL6A6 in lung cancer is unclear. We found that COL6A6 expression in LUAD tissues was significantly lower than in normal tissues. Additionally, we found that upregulated expression of COL6A6 was associated with a favorable prognosis. Logistic regression analysis revealed that COL6A6 expression was related to pathological stage, T stage, and N stage. Meta-analysis effectively combines statistical strength from multiple data sets which allows for greater precision compared to when using any of the single studies. The survival meta-analysis in this study suggested that high COL6A6 expression was a protective factor for LUAD patients. Therefore, we speculated that COL6A6 may act as a tumor suppressor in LUAD.
Another important finding of this study was the effects of COL6A6 on immune cell infiltration in LUAD. The TIMER analysis showed close correlations between COL6A6 expression and abundances of B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. Decreased copy number may lead to reduced infiltration of CD8+ T cells, macrophages, neutrophils and dendritic cells. The abovementioned evidence indicated the potential role of COL6A6 in the immune microenvironment. As some immune infiltrates, such as macrophages, have been found to exist in both antitumorigenic and protumorigenic forms, we used CIBERSORTx for further analysis. Consistent with the TIMER results, we found that the proportions of several types of B cells, T cells, neutrophils, and dendritic cells were apparently increased in the high expression group while the proportions of macrophages M0 and Tregs were lower. It has been reported FIGURE 7 | Effects on immune cell infiltration of COL6A6 in lung adenocarcinoma by CIBERSORTx. The proportions of B cells naive (P = 0.041), dendritic cells resting (P = 0.038), mast cells resting (P < 0.0001), neutrophils (P = 0.043), and T cells CD4 memory resting (P < 0.0001) were higher in the high expression group than the low expression group. In contrast, the proportion of macrophages M0 (P < 0.0001), plasma cells (P < 0.001), and Tregs (P < 0.001) were lower. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. that the immune microenvironment is frequently associated with lung cancer outcome. T cells, especially cytotoxic and memory T cells, have been consistently shown to have positive effects on prognosis (28,(32)(33)(34). B cells, dendritic cells, and eosinophils have also been reported to be associated with a favorable prognosis (32,33,35). In contrast, high proportions of Tregs and macrophages M0 have been associated with a poor prognosis among LUAD patients (28,36). Thorsson et al. identified six immune subtypes that span cancer tissue types and molecular subtypes, and the inflammatory subtype had the best prognosis (29). In this study, COL6A6 expression was highest in the inflammatory immune subtype, indicating a better prognosis in LUAD; the positive impact of COL6A6 expression prognosis among LUAD patients was consistent with the potential roles of diverse immune cell types. This further confirms the association between COL6A6 expression and survival of LUAD patients.
Intriguingly, the GSEA results suggested that COL6A6 was associated with differentiation and activation of several lymphocytes including B and T cells. Consistently, functional analysis of the COL6A6-associated immunomodulators indicated that regulation of T cell activation, leukocyte cell−cell adhesion, and T cell receptor signaling pathway were related to COL6A6-mediated immune events. COL6A6 participates in the cell adhesion process, and many other cell adhesion molecules, especially integrins, are known to play key roles in the regulation of immune cell infiltration into the TME (37). Integrins have various effects on T cell infiltration, including T cell recruitment into tumors, priming, and effector function (37,38). Multiple integrin pairs, including    a1b1, a2b1, a3b1, a10b1, and avb3 integrin, have been identified as potential binding partners for COL6 (39,40). The PI3K-Akt pathway is another crucial pathway involved in tumor immunity. As a major pathway in tumorigenesis, its role in tumor immunity has also been reported. Okkenhaug demonstrated that the pathway contributes to the development and differentiation of B cells and T cell subsets (41). PI3K inhibition in dendritic cells enhances type I immune responses. Moreover, inhibition of the PI3K pathway enhances CD8+ T cell infiltration into tumor tissue (42). PI3K-AKT-mTOR inhibition can regulate immunosuppressive cytokine secretions, effect Treg infiltrations into tumor tissues and promote the development of memory T cells (43). The PI3K pathway may also regulate the expression of PD-L1 (44). Furthermore, COL6A6 was reported to inhibit the PI3K-Akt pathway; Long et al. found that COL6A6 suppressed the growth and metastasis of pituitary adenocarcinoma by blocking the PI3K-Akt pathway (11). Together, these findings show that COL6A6 may play a pivotal role in the tumor immune microenvironment of LUAD. Gene signatures have been widely used as prognostic predictors in various cancers (45,46). Zhuang et al. analyzed M1 macrophage- related genes and identified a four-gene signature for predicting thyroid cancer prognosis (45). Li et al. analyzed the gene expression profiles of frozen tumor tissue samples from 19 public NSCLC cohorts and identified a prognostic immune-related gene signature involving 25 gene pairs (46). Their immune signature achieved moderate prognostic accuracy (C-index=0.64). In this study, we established a seven-gene signature for LUAD based on COL6A6associated immunomodulators. The risk score derived from the immune gene signature was an independent predictor of LUAD prognosis. The signature had high accuracy in both the training dataset and validation dataset. Combined with clinical features, we further built a prognostic nomogram for personalized prediction, which was found to have a C-index of 0.723. Our findings may provide clinicians with a convenient and accurate way to evaluate the prognosis of patients with LUAD after surgery. However, there are several limitations in this study. First, our study was mainly based on online databases. Nonetheless, integrating bioinformatic analyses, such as meta-analysis add strength to this study. We will further verify the expression and function of COL6A6 as well as its correlation with immune cell infiltration in the laboratory. Second, we were unable to analyze the relationship between COL6A6 expression and immunotherapy response, as these data in the TCGA and other databases were insufficient. Lastly, the prognostic value of the seven-gene signature should be further validated in an in-house patient population.
In summary, COL6A6 is downregulated in LUAD tissues, while increased COL6A6 expression predicts a favorable prognosis. COL6A6 expression is associated with the infiltration of various immune cells, such as B cells, T cells, and dendritic cells, in LUAD. COL6A6 may play a role in the tumor immune microenvironment of LUAD. The seven-gene signature derived from COL6A6associated immunomodulators is an independent predictor of OS in LUAD. The prognostic nomogram showed good performance in predicting the OS of LUAD and might thus be beneficial for individualized treatment and medical decision making.

AUTHOR CONTRIBUTIONS
FY, MQ, XL, and YM: study concept and design. YM, HG, HC, and JL: acquisition and analysis of the data. YM and MQ: drafting and revising of the manuscript. All authors contributed to the article and approved the submitted version.