Identification of immune-related lncRNA panel for predicting immune checkpoint blockade and prognosis in head and neck squamous cell carcinoma.

PURPOSE: Immunotherapy is changing head and neck squamous cell carcinoma (HNSCC) treatment pattern. According to the Chinese Society of Clinical Oncology (CSCO) guidelines, immunotherapy has been deemed as first-line recommendation for recurrent/metastatic HNSCC, marking that advanced HNSCC has officially entered the era of immunotherapy. Long non-coding RNAs (lncRNAs) impact every step of cancer immunity. Therefore, reliable immune-lncRNAs able to accurately predict the immune landscape and survival of HNSCC are crucial to clinical management.
METHODS: In the current study, we downloaded the transcriptomic and clinical data of HNSCC from The Cancer Genome Atlas and identified differentially expressed immune-related lncRNAs (DEir-lncRNAs). Further then, Cox and least absolute shrinkage and selection operator (LASSO) regression analyses were performed to identify proper DEir-lncRNAs to construct optimal risk model. Low-risk and high-risk groups were classified based on the optimal cut-off value generated by the areas under curve for receiver operating characteristic curves (AUC), and Kaplan-Meier survival curves were utilized to validate the prediction model. We then evaluated the model based on the clinical factors, immune cell infiltration, and chemotherapeutic and immunotherapeutic efficacy between two groups.
RESULTS: In our study, we identified 256 Deir-lncRNAs in HNSCC. A total of 18 Deir-lncRNA pairs (consisting of 35 Deir-lncRNAs) were used to construct a risk model significantly associated with survival of HNSCC. Cox proportional hazard regression analysis confirmed that our risk model could be served as an independent prognostic indicator. Besides, HNSCC patients with low-risk score significantly enriched of CD8+ T cell, and corelated with high chemosensitivity and immunotherapeutic sensitivity.
CONCLUSION: Our risk model could be served as a promising clinical prediction indicator, effective discoursing of the immune cell infiltration of HNSCC patients, and distinguishing patients who could benefit from chemotherapy and immunotherapy.

INTRODUCTION

Head and neck squamous cell carcinoma (HNSCC) mainly arises in the oral cavity, oropharynx, hypopharynx, and larynx. HNSCC was the seventh most common cancer worldwide.  Multiple pathogenic factors (including heavy use of tobacco, alcohol, and chronic infection of human papilloma virus) have been proven in association with the initial of HNSCC. , Smoking could increase the risk of HNSCC from 5‐ to 25‐fold, and the risk increases with the quantity of cigarette. ,  The consumption of alcohol also independently doubles the risk of HNSCC.

Cancer is characterized by the accumulation of various of genetic alterations, which results in the acquisition of ten biological capabilities during the multistep development of human tumors. Subsequent molecular studies have revealed that these genetic or epigenetic alterations often result in mutated cellular antibodies or biomarkers and aberrantly expressed genes. , ,  These abnormal expressed molecules on the surface of cancer cells could be recognized by CD8+ T cell, distinguishing them from their normal counterparts and even inhibiting the outgrowth of cancer cells.  Therefore, in the past decades, researchers were dedicated to efficiently activate cancer immunity system to fight cancer cells, which is called “immunotherapy.” However, the cancer immunity cycle does not perform optimally in cancer patient; sometimes, rarely protective efforts, even promoting assists by T cell were found during the progression of human cancers. , , Either tumor antibodies may not be detectable, or dendritic cells and T cells may classify these antigens as self thereby creating low ratio of cytotoxic T lymphocytes. , , Besides, negative regulators to T cell responses are also responsible for the failure of immune protection.

The goal of cancer immunotherapy is to initiate or reinstate the cycle of cancer immunity to exert anti‐tumor effects. The cancer immunotherapies must be carefully configured to overcome these negative feedback mechanisms. Indeed, immunotherapy exhibits marvelous efficacy in clinical trials ; however, only a minority of patients have incredible response to immunotherapy in clinical practice. Therefore, a deep understanding of immune system and identification of these sorts of HNSCC patients will be helpful to develop effective therapeutic strategies for cancer immunotherapy and to improve the overall survival of HNSCC. In the present study, we aim to comprehensively explore immune‐related lncRNA expression profile and construct a model for predicting outcome, drug sensitivity, and adapting implement clinical strategies in HNSCC.

METHODS

Data download and pretreatment

The online database TCGA program was aimed to molecularly characterize human cancers. In our study, a series of transcriptome RNA‐sequencing data of HNSCC samples were downloaded through GDC portal of TCGA (https://portal.gdc.cancer.gov/), which contained data from 501 HNSCC cases and 44 normal controls. The clinical data were also extracted (including TNM stage, gender, grade, age, survival status, and survival data). The long non‐coding RNA (lncRNA) data and messenger RNA (mRNA) were extracted from transcriptome RNA‐sequencing data, and the annotation file of Ensembl IDs (from http://asia.ensembl.org/index.html) was applied to convert RNA‐seq results into gene symbols. The study was performed based on the TCGA public database. Since there was no personal identifying information was used in the current study, it was granted an exemption from ethics approval from the Institutional Review Board of the Ningbo Medical Center Lihuili Hospital.

Immune‑related differential expressed lncRNAs acquisition

Firstly, immune‐related genes were obtained from the ImmPort database (http://www.immport.org) and was used to screen immune‐related lncRNAs (ir‐lncRNAs) by a co‐expression strategy. Then, correlation analysis was performed between all lncRNAs of HNSCC samples and known immune‐related genes by Pearson's correlation. The cut‐off criteria of adjusted p‐value (adj. p‐value) was set as 0.001 and the criterion of correlation coefficients was set as more than 0.4. Subsequently, we applied limma package of R software to identify the differentially expressed ir‐lncRNAs. The thresholds were set as log fold change (FC) >1 along with false discovery rate (FDR) <0.05. Besides, to analyze the function of ir‐lncRNAs, gene ontology (GO) was performed using the online database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.ncifcrf.gov/, version 6.8).

Pairing DEir‐lncRNAs

The Cancer Genome Atlas (TCGA) is a comprehensive database that includes multi‐layered cancer genome profiles. Large‐scale collection of data inevitably generates batch effects introduced by differences in processing at various stages from sample collection to data generation. To eliminate these factors, we cyclically paired these differentially expressed ir‐lncRNAs (DEir‐lncRNAs). We constructed a new 0‐or‐1 matrix by comparison of all the DEir‐lncRNAs expression in each patient of TCGA. We defined the DEir‐lncRNA pair as 1 if the expression level of lncRNA A is higher than lncRNA B; otherwise, this pair was defined as 0. This new matrix was used as further analyses.

Survival‐related DEir‐lncRNAs pair

Paired DEir‐lncRNAs associated with clinical outcomes in HNSCC patients were identified as survival‐related DEir‐lncRNA pairs. Survival‐related DEir‐lncRNA pairs were selected by univariate Cox regression applying survival package of R software. These survival‐related Deir‐lncRNA pairs were specified for subsequent research.

Establishment of a risk model

The least absolute shrinkage and selection operator (LASSO) is a regression‐based methodology permitting for a large number of covariates in the model, and then regulating the impact of a coefficient may have on the overall regression. The greater the penalization, the greater the shrinkage of coefficients, with some reaching 0, thus automatically removing unnecessary/uninfluential covariates. The survival‐related DEir‐lncRNA pairs were selected for Lasso regression and Cox proportional hazard regression analysis, as well as for construction of the model. Risk score is a simplified version of a prognostic model, in which scores are assigned to each risk factor (on the basis of rounded regression coefficients). The risk score of each sample was calculated by the following formula:

The value of coefficient was generated from Cox proportional hazard regression analysis. The 1‐, 2‐, 3‐, 4‐, 5‐, and 6‐year receiver operating characteristic (ROC) curves of the model were cyclically plotted till the maximum areas under curve for each ROC (AUC) were obtained. Then, the highest point of risk score of maximum AUC were set to the cut‐off criteria to classify HNSCC patients into high or low risk of group.

Validation of the constructed risk model

To validate the constructed risk model, we performed Kaplan–Meier analysis to show the survival difference of high and low risk group. Subsequently, we also re‐assessed the relationship between the model and clinicopathological characteristics. Then, to confirm whether the model can be used as an independent clinical prognostic biomarker, we performed the univariate and multivariate Cox regression analyses, and the results were represented by forest map.

Immune cell abundance analysis

We collected the currently acknowledged methods to calculate the immune infiltration statues from CIBERSORT.  The results were shown in a lollipop diagram, which was finished by violet package of R software.

Prediction of clinical application of the risk model

Recent studies have predicted the response to immune checkpoint inhibitors (ICIs) of cancers based on their genomic features, and the immunophenoscore (IPS) was correlated with responses to ICIs immunotherapy. ,  The score of IPS‐CTLA‐4 and IPS‐PD1 stands for the potential for response rates of ICIs in HNSCC. Besides, we calculated a half inhibitory centration (IC50) of common chemotherapeutic drugs, such as cisplatin, docetaxel, gemcitabine, and paclitaxel. The difference in the IC50 between the high‐ and low‐risk groups were assessed, and the results were shown as box drawings by pRRophetic package of R.  Then, we also analyzed the relationship between the model and the expression of ICI‐related genes.

RESULTS

Acquisition of differential immune‐related lncRNAs

We firstly abstracted the lncRNA and mRNA matrix from transcriptome RNA‐sequencing data of TCGA (https://www.cancer.gov/about‐nci/organization/ccg/research/structural‐genomics/tcga). Subsequently, the immune‐related mRNAs (ir‐mRNAs) of TCGA were achieved by the intersection of mRNA matrix of TCGA and a list of well‐known immune‐related genes from ImmPort database. Then, the immune‐related lncRNAs (ir‐lncRNAs) were obtained by Pearson's analysis between lncRNA matrix and ir‐mRNAs of HNSCC from TCGA. Eventually, we obtained a total of 804 ir‐lncRNAs. Subsequently, we identified differentially expressed ir‐lncRNAs (DEir‐lncRNAs) between 501 HNSCC tumor tissues and 44 normal tissues, the results showing 256 DEir‐lncRNAs with FDR‐adjusted p < 0.05 and absolute log fold change (FC) >1 (Table S1). These DEir‐lncRNAs were represented by heatmap and volcano plots in Figure 1. The functions of lncRNAs are thought to reflected by their associated mRNAs. Therefore, to exploit the underlying functions of immune‐related lncRNA signature, gene ontology (GO) analyses was carried out for differentially co‐expressed mRNAs. The results of the GO analysis are shown in Figure S1; many immunity‐associated processes were identified, including immune response, innate immune response, adaptive immune response, and T and B cell receptor signaling pathway.

FIGURE 1

Volcano (1A) and heatmap (1B) plots of DEir‐lncRNAs of HNSCC

Identification of survival‐related DEir‐lncRNAs pair

We cyclically paired DEir‐lncRNAs pairs through comparison of all the DEir‐lncRNAs expression in each HNSCC patient of TCGA. Then, univariate Cox regression was performed to identify survival‐related DEir‐lncRNA pairs, and 253 DEir‐lncRNA pairs were stepped into LASSO regression analysis. Taking into account the large number of confounding factors caused by multiple variables, we applied LASSO regression to obtain the best DEir‐lncRNAs pairs to improve the effectiveness of subsequent analysis. Eventually, a total of 18 DEir‐lncRNAs pairs were identified and applied for establishment of a risk model of HNSCC, which is composed of 35 lncRNAs (LINC01063, AL139288.1, MIR9‐3HG, AC091563.1, MIR924HG, KDM4A‐AS1, AC106820.3, AC022031.2, AC108488.1, AC004687.1, RUSC1‐AS1, ZNF687‐AS1, AL354733.3, SCAT1, LINC00460, AC104825.1, KTN1‐AS1, AL133243.2, AC022144.1, HS1BP3‐IT1, LINC02454, AC112484.3, C5orf66‐AS1, AC132192.2, AC004148.1, MIR1‐1HG, AC007991.2, AL133353.1, AC083964.1, AP003555.1, SNHG25, AC012236.1, P3H2‐AS1, AC008735.2, and AC060766.4).

Based on the results of LASSO regression, a total of 61 DEir‐lncRNA pairs was stepped into multivariate Cox proportional hazard regression analysis to specify the variables in the model formula, and 18 DEir‐lncRNA pairs were determined. The formula for calculation the risk score is: 0.324 × expression of (LINC01063|AL139288.1)−0.250 × expression of (MIR9‐3HG|AC091563.1)−0.295 × expression of (MIR9‐3HG|MIR924HG)+0.554 × expression of (KDM4A‐AS1|AC106820.3)+0.298 × expression of (AC022031.2|AC108488.1)−0.573 × expression of (AC004687.1|RUSC1‐AS1)+0.403 × expression of (ZNF687‐AS1|AL354733.3)−0.343 × expression of (SCAT1|LINC00460)−0.248 × expression of (AC104825.1|KTN1‐AS1)+0.378 × expression of (AL133243.2|AC022144.1)−0.371 × expression of (HS1BP3‐IT1|LINC02454)+0.4102 × expression of (AC112484.3|C5orf66‐AS1)+0.420 × expression of (AC132192.2|AC004148.1)+0.345 × expression of (MIR1‐1HG|AC007991.2)+0.340 × expression of (AL133353.1|AC083964.1)+0.380 × expression of (AP003555.1|P3H2‐AS1)−0.281 × expression of (AC008735.2|SNHG25)−0.357 × expression of (AC012236.1|AC060766.4). The risk score of each sample of TCGA was computed and then calculated the areas under curve for each receiver operating characteristic (ROC) curve of 18 pairs. The distribution of risk score of each sample was exhibited in Figure 2. The maximum area under curve of ROC (AUC) was 0.810 (Figure 3A), and the optimal cut‐off point to differentiate the high‐ or low risk‐group among patients with HNSCC was 1.429 (Figure 3B).

FIGURE 2

Risk scores of each HNSCC patient based on risk model

FIGURE 3

Time‐dependent receiver operating characteristic (ROC) analysis at 1–6 years to identify the maximum of the AUC (A) and to identify the optimal cut‐off point to group HNSCC case (B)

Validation of the risk model

According to the cut‐off point confirmed above, 195 patients were grouped into high‐risk group and 304 into low‐risk group. Kaplan–Meier analysis showed that patients in the low‐risk group exhibited a longer survival time than those in the high‐risk group (p < 0.0001) (Figure 4). Subsequently, we performed stratified analysis between high‐ and low‐risk group by clinical characteristics (such as gender, age, TNM stage, and grade), and subgroup analysis results robustly supported that better overall survival of patients in low‐risk group than in high‐risk group (Figure S2). In order to avoid basis from clinical features, univariate and multivariate Cox proportional hazard regression analysis were performed to evaluate the independent prognostic value of risk model. The results showed that our risk model, clinical stage, and age could be served as independent prognostic indicators (Figure 5).

FIGURE 4

The Kaplan–Meier plot of HNSCC grouped by risk model

FIGURE 5

Univariate and multivariate Cox proportional hazard regression of risk model and other clinical features

Then, we assessed the relationship between the risk score of HNSCC and clinicopathological characteristics (TNM stage, age, grade). The results showed that patients with advanced T stage more likely gathered in the high‐risk group, indicating that higher scores of the 18 DEir‐lncRNA pairs might be significantly associated with the progression of HNSCC (Figure S3).

Immune cell abundance analysis based on risk model

We analyzed the infiltration difference of immune cells between the low‐ and the high‐risk group classified by the risk model. The immune infiltration estimations for TCGA expression profiles was calculated by CIBERSORT (https://cibersort.stanford.edu/). Then, the intersection of the infiltration estimation file and our risk score file were acquired for further Pearson's correlation analysis. The association between immune cell infiltration difference and different risk group was presented by vioplot (Figure 6). Patients in low‐risk group showed higher infiltration levels of B cell, CD8 T cells, CD4 memory‐activated T cell, regulatory T cells (Tregs), and mast cells, whereas patients in high‐risk group was more correlated with NK cells, macrophages M1, macrophages M2, activated mast cells infiltration.

FIGURE 6

Estimation of tumor‐infiltrating cells by the risk model

Subsequently, we also assessed the association between several ICI biomarkers and our risk model, and the results showed that low risk group were more likely with high expression of cytotoxic T lymphocyte‐associated antigen‐4 (CTLA4) and programmed cell death protein 1 (PD‐1), whereas the expression of programmed cell death ligand 1(PD‐L1) showed no statistical differences (Figure 7).

FIGURE 7

The gene expression of PD1, CTLA4, and PD‐L1 in high‐risk and low‐risk groups

Prediction drug sensitivity based on risk model

Considering the positive association between two important ICI biomarkers, we assess the innate sensitivity or resistance to anti‐PD‐1 therapy. Our results showed that HNSCC patients with low‐risk score represented significantly higher IPS‐CTLA4 and IPS‐PD1 (immunophenoscore, IPS), indicating the potential of ICI application for HNSCC (Figure 8). Besides, we also identified associations between risk and the efficacy of common chemotherapeutics (including cisplatin, docetaxel, gemcitabine, and paclitaxel) in HNSCC (Figure 9). The results indicated that HNSCC patients with low risk represented higher sensitivity to docetaxel and gemcitabine.

FIGURE 8

The association between the response to ICIs and risk model in HNSCC

FIGURE 9

The associations between risk model and the efficacy of common chemotherapeutics (including cisplatin, docetaxel, gemcitabine, and paclitaxel) in HNSCC

DISCUSSION

HBSCCs are a heterogeneous group of cancers originated from oral cavity, pharynx, and larynx, which are characterized by different molecular subgroups and clinical features. The primary treatment of early stage of HNSCC is radical radiotherapy or surgery, and the treatment option of advanced stage of HNSCC is radiotherapy combination with chemotherapy. However, a substantial proportion of HNSCC patients die from their disease after such treatment. Besides, these extensive surgery and high‐dose chemotherapies are also mutilating for HNSCC patients.

The increasing understanding of tumor immunology has led to intense preclinical and clinical research to new therapy, called “immunotherapy.” Advances in immunotherapy have brought hope, especially the favor outcome of clinical trials, such as KEYNOTE‐012, KEYNOTE‐141, KEYNOTE‐055, and KEYNOTE‐040. , , , ,  Two anti‐PD‐1 antibodies (belongs to ICIs), pembrolizumab and nivolumab, are approved for the treatment option of recurrent and metastatic HNSCC. Subgroup analysis showed that HNSCC patients with high expression of PD‐L1 are deemed more sensitivity to above two antibodies. However, only a minority of patients have benefit from these two drugs in the clinical practice. Therefore, to identify a proportion of HNSCC patients susceptible of immunotherapy, researchers took advantage of advanced technology, such as targeted next‐generation sequencing, and had constructed multiple relevant models or biomarkers to prediction. The most representative biomarker is tumor mutation burden (TMB) in lung cancer.  TMB is significantly associated with improved benefit among patients with non‐small‐cell lung cancer (NSCLC) treated with ICIs and is an independent variable.  Therefore, to identify such reliable biomarkers or models for distinguishing this sort of patients or different risk patient subgroup, we performed a comprehensively analysis to the expression of multi‐type of RNA profile in HNSCC.

We analyzed the expression profile of lncRNA in HNSCC for the following reasons: (1) the amount of ncRNAs in the whole genome is 98%, which were once thought to be “garbage sequence”; (2) increasing evidence demonstrated the critical role of lncRNA in biological process and the association between lncRNA and human cancers; (3) recent studies promulgated the eye‐catching roles of lncRNA in anti‐tumor immunity. Firstly, we synthesized the mRNA profile of TCGA HNSCC and well‐known ir‐mRNA from public database. Subsequently, Pearson correlation analysis was performed to identify ir‐lncRNAs. Finally, differentially expressed ir‐lncRNAs (DEir‐lncRNAs) were filtrated through “limma” package of R software, and a total of 256 DEir‐lncRNAs were identified. These 256 DEir‐lncRNAs (including LINC01614, LINC01322, and LINC00460, detail information presented in Table S1) were consist of an amount of well‐studied lncRNAs associated with human cancers, such as LINC01614 in breast cancer, lung cancer, and gastric cancer. , Changes in the composition, physical properties, and spatial conformation of the extracellular matrix (ECM) play an important role in the immune microenvironment of tumor, and these changes are proven having a significant impact on drug sensitivity and patient survival. , , LINC01322 was significantly enriched in pathways related to the ECM, such as in ECM structural constituents, celladhesion–molecule binding, actin binding, and glycosaminoglycan binding.

It was impossible to assess all the expression level of 256 DEir‐lncRNAs in clinical practice. Therefore, we firstly converted these DEir‐lncRNAs into relative expression, namely lncRNA pairs, to avoid examining specific expression values of every lncRNA. And then, we performed a new method for variable selection and shrinkage in Cox's proportional hazards model – the least absolute shrinkage and selection operator (LASSO), which was deemed more accurate than stepwise selection. There were 18 DEir‐lncRNA pairs (consisting of 35 DEir‐lncRNAs, detail information presented in Table S1) included in our risk model for prediction of outcome and drug sensitivity in HNSCC patients. Most of 35 DEir‐lncRNAs were associated with various human cancers, including HNSCC. Besides, multiple lncRNAs were deemed as an important immune‐related biomarker in auxiliary prognosis prediction for cancer. In our study, the areas under curve for each receiver operating characteristic (ROC) curve of 18 pairs was drawn to estimate the maximum AUC for survival outcome and to identify optimal cut‐off point to differentiate the high‐ or low risk‐group among patients with HNSCC. Our results showed the AUC for 3‐year survival of HNSCC is 0.810, indicating the robust prognosis potential of our risk model. Besides, Kaplan–Meier analysis and univariate and multivariate Cox proportional hazard regression analysis authenticated the poor prognosis prediction value of our risk model in HNSCC.

Increasing attention of lncRNAs has been attracted not only for the value of early screening, diagnosis, and prognosis, but also for the prediction of response to drug treatment. , , Long non‐coding RNA SCAT1 could be served as a pretreatment biomarker in predicting pathological response and outcome in esophageal squamous cell carcinoma with neoadjuvant chemoradiotherapy. A stroma‐related lncRNA signature (SLS) has been reported to predict adjuvant chemotherapy benefit in patients with colon cancer. LncRNAs are also participated in the acquired resistance to chemotherapy, , and targeting to these lncRNAs can reverse drug resistance and enhance the sensitivity of cancer cells to chemotherapy. Several studies have indicated that lncRNAs not only are involved in the typical hallmarks of cancer but also are closely correlated with the regulation of immune response. , Given the crucial role of lncRNAs in cancer, we also assessed the drug sensitive of our risk model, which may provide new insights for precise treatment and individualized management of patients with HNSCC. Firstly, we estimated the immune cell infiltration profile, and results showed significant enrichment of B cells naive, plasma cells, CD8+ T cells, T follicular helper cells (Tfh cells), regulatory T cells (Treg), mast cells resting in low‐risk group, while significant enrichment of CD4+ T cell memory resting, NK cells resting, macrophages M0 and M2 in high‐risk group. Tfh cells have been established as a CD4+ subset and expressed a range of cell surface molecules that serve important functions not only in their interactions with B cells but also for their identification like CXC‐chemokine receptor 5 (CXCR5) and the coinhibitory receptor programmed cell death 1 (PD‐1). , ,

The pro‐tumorigenic or anti‐tumorigenic effects of cancer immune system exerted rely on specific microenvironment. ,  Numerous factors in the tumor microenvironment can act to modulate the existing activated anti‐tumor T cell immune response, such as expression of PD‐L1 or PD‐1. , Besides, in tumor immunity, Treg cells are involved in tumor development and progression by inhibiting anti‐tumor immunity. Recent studies have explored that Treg cells is linked to the tumor immune surveillance and carcinogenesis. ,  There is an interesting phenomenon in our study, namely enrichment of Treg cells in low‐risk group of patients with HNSCC. Numerous factors in the tumor microenvironment can modulate the anti‐tumor T cell immune response, such as the expression of IFN‐γ. Patients with high expression of PD‐L1 showed more sensitivity to ICIs, such as pembrolizumab and nivolumab; these patients represented favor outcome, which rationalize the above phenomenon. And our prediction of drug sensitivity also showed that HNSCC patients with low‐risk score represented significantly higher sensitivity to PD‐1 and CTLA4.

In summary, we identified an immune‐related lncRNAs in HNSCC, and constructed a risk model for prediction prognosis for patients with HNSCC and might help in distinguishing patients who could benefit from anti‐tumor immunotherapy.

AUTHOR CONTRIBUTIONS

All authors contributed to the planning and design of the study. CZ and QL were involved in review of the raw data and directly involved in the analysis. ZS, YS, and YMS provided analytical feedback based on aggregated results. CZ, QL, and GZ drafted the manuscript, with input from all authors. JW, HD, and YMS were responsible for the chart making. All authors provided substantive review and commentary on multiple drafts and approved the final version. QL and GZ had full access to all of the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.

CONFLICT OF INTEREST

The authors declare that there are no conflicts of interest.

FigS1

Click here for additional data file.

FigS2

Click here for additional data file.

FigS3

Click here for additional data file.

TableS1

Click here for additional data file.