Abstract
Background
Both lncRNAs and glycolysis are considered to be key influencing
factors in the progression of bladder cancer (BCa). Studies have shown
that glycolysis-related lncRNAs are an important factor affecting the
overall survival and prognosis of patients with bladder cancer. In this
study, a prognostic model of BCa patients was constructed based on
glycolysis-related lncRNAs to provide a point of reference for clinical
diagnosis and treatment decisions.
Methods
The transcriptome, clinical data, and glycolysis-related pathway gene
sets of BCa patients were obtained from The Cancer Genome Atlas (TCGA)
database and the Gene Set Enrichment Analysis (GSEA) official website.
Next, differentially expressed glycolysis-related lncRNAs were screened
out, glycolysis-related lncRNAs with prognostic significance were
identified through LASSO regression analysis, and a risk scoring model
was constructed through multivariate Cox regression analysis. Then,
based on the median of the risk scores, all BCa patients were divided
into either a high-risk or low-risk group. Kaplan-Meier (KM) survival
analysis and the receiver operating characteristic (ROC) curve were
used to evaluate the predictive power of the model. A nomogram
prognostic model was then constructed based on clinical indicators and
risk scores. A calibration chart, clinical decision curve, and ROC
curve analysis were used to evaluate the predictive performance of the
model, and the risk score of the prognostic model was verified using
the TCGA data set. Finally, Gene Set Enrichment Analysis (GSEA) was
performed on glycolysis-related lncRNAs.
Results
A total of 59 differentially expressed glycolysis-related lncRNAs were
obtained from 411 bladder tumor tissues and 19 pericarcinomatous
tissues, and 9 of those glycolysis-related lncRNAs ([33]AC099850.3,
[34]AL589843.1, MAFG-DT, [35]AC011503.2, NR2F1-AS1, [36]AC078778.1,
ZNF667-AS1, MNX1-AS1, and [37]AC105942.1) were found to have prognostic
significance. A signature was then constructed for predicting survival
in BCa based on those 9 glycolysis-related lncRNAs. ROC curve analysis
and a nomogram verified the accuracy of the signature.
Conclusion
Through this study, a novel prognostic prediction model for BCa was
established based on 9 glycolysis-related lncRNAs that could
effectively distinguish high-risk and low-risk BCa patients, and also
provide a new point of reference for clinicians to make individualized
treatment and review plans for patients with different levels of risk.
Keywords: bladder cancer, glycolysis, lncRNAs, prognosis, signature
Introduction
Bladder cancer (BCa) is one of the most prevalent malignancies of the
urinary system, according to statistics from 185 countries worldwide,
in 2020, it ranked 9th in incidence and 13th in mortality ([38]Sung et
al., 2021). BCa is a heterogeneous disease and can be divided into two
clinical types depending on whether it infiltrates the muscle layer:
25% of cases are muscle invasive BCa (MIBC), and 75% are non-muscle
invasive BCa (NMIBC) ([39]Cumberbatch et al., 2018). NMIBC has a low
mortality rate, while 70% of these tumors will recur and 15% will
progress in stage and grade. In contrast, the long-term survival rate
of MIBC is unsatisfactory because it is prone to develop distant
metastasis ([40]Knowles and Hurst, 2015). Recent studies suggested that
the poor prognosis of MIBC may also be related to the different
molecular mechanisms of its occurrence and progression; some scholars
even believed that MIBC and NMIBC were not the same kind of disease
([41]Gui et al., 2011). With the increasing in-depth study of molecular
biological mechanisms and the rapid development of gene detection
technology, molecular typing and prognostic evaluation of BCa by
genetic testing is becoming a new diagnostic and therapeutic target.
Glycolysis is a metabolic pathway that provides energy to the body.
Changes in the glycolytic state may signify the presence of cancer. One
study by Otto Warburg found that the glycolysis of tumor cells produced
more lactic acid and less ATP, which resulted in tissue hypoxia, even
in an environment with a high oxygen content. This phenomenon is known
as the “Warburg effect” ([42]Warburg, 1956). Enhanced or abnormally
activated glycolysis promoted tumor progression, which has been proven
in breast cancer, colorectal cancer, liver cancer, etc. ([43]Chen et
al., 2019a; [44]Weng et al., 2020; [45]Zhang et al., 2020). In BCa, an
increase of the Warburg effect has been proven to contribute to the
aggressiveness of BCa tumor cells and accelerated the proliferation
rate of BCa cells ([46]Ritterson Lew et al., 2015). Therefore, changes
in glycolytic status may be an emerging marker of malignancy and a
potential prognostic target for patients with BCa.
LncRNAs are a class of nucleotides with a transcriptional length of
more than 200 nucleotides. Although lncRNAs do not possess
protein-coding ability, they can regulate gene expression through
multiple layers and pathways, thus affecting the occurrence and
progression of tumors, and can be used as markers for tumor diagnosis
and prognosis ([47]Evans et al., 2016). There are more than 50,000
lncRNAs in the human genome, and they can play different roles in
different tumors, acting as tumor-suppressor or tumor-promoter genes
([48]Derrien et al., 2012). To effectively predict the prognosis of
patients with BCa, the sensitivity and specificity of a single lncRNA
was not sufficient; it was necessary to combine the high or low
expression of multiple lncRNAs to have stronger diagnostic value.
Hyperglycolysis is a common phenomenon in the rapid proliferation of
tumor cells and is also a marker of tumorigenesis. Recently, increasing
amounts of data have shown that lncRNAs were a crucial regulator of
glycolysis, and some studies have attempted to prevent tumor
progression by inhibiting the activity of key enzymes in the tumor
glycolysis pathway ([49]Yang et al., 2014; [50]Chen et al., 2019a;
[51]Liao et al., 2019; [52]Wang et al., 2019a). The functional study of
glycolysis-related lncRNAs has found that different lncRNAs could
promote glycolysis in various ways, which in turn promoted
tumorigenesis, progression, and metastasis. For example, lncRNA-SNHG7
regulated by c-MYC could increase glycolysis by inhibiting the
expression of miR-34a-5p, thus promoting the progression of breast
cancer ([53]Zhang et al., 2019), while the lncRNA LINRIS stabilized
IGF2BP2 and promoted aerobic glycolysis in colorectal cancer ([54]Wang
et al., 2019b).
The functional study of glycolysis-related lncRNAs in bladder cancer
has found that lncRNA-SLC16A1-AS1 could be used as a target of E2F1 and
a co-activator to induce metabolic reprogramming in the progression of
bladder cancer ([55]Logotheti et al., 2020), and lncRNA UCA1
up-regulated hexokinase 2 in bladder tumor cells to promote glycolysis
through the mTOR–STAT3/microRNA143 pathway ([56]Li et al., 2014).
Prognostic models constructed by differentially expressed
glycolysis-related lncRNAs have been proven to have good predictive
performance in gastrointestinal tumors, breast cancer, gliomas, and
other tumors. In bladder cancer, there is currently no known
glycolysis-related lncRNA prognostic model.
This study aimed to clarify the relationship between glycolysis-related
lncRNAs and the prognosis of BCa. First, a total of 9
glycolysis-related lncRNAs were identified to have prognostic
significance. Then, based on these 9 glycolysis-related lncRNAs, a
novel prognostic prediction model for BCa was constructed and verified,
which proved that the model possessed good predictive power. Finally,
gene function enrichment analysis was performed on the differentially
expressed glycolysis-related lncRNAs identified through this study to
explore their potential functions in the progression of BCa.
Materials and Methods
Sample Sources and Processing
Transcriptome sequencing data and clinical data of BCa patients were
downloaded from The Cancer Genome Atlas (TCGA) database^[57]1, and gene
sets of glycolysis related functions and glycolysis related pathways
were obtained from the Gene Set Enrichment Analysis (GSEA) official
website^[58]2. The “edgeR” package of R software was used to normalize
the entire data set, set | log2FC | > 1 and FDR < 0.05 as the threshold
to construct the volcano map, and obtain the tumor tissue and
pericarcinomatous differentially expressed glycolysis-related genes and
lncRNA. The correlation between the lncRNAs and glycolysis-related
genes was calculated using Pearson correlation. If the correlation
coefficient of lncRNA | R2 | > 0.3 and P < 0.05, it was considered to
be related to glycolysis and used for further analysis.
Construction of the Risk Score Model
Univariate Cox regression and Least Absolute Shrinkage and Selection
Operator (LASSO) regression analysis were used to screen out 9 lncRNAs
with prognostic significance. Genes with P < 0.05 were considered to be
independent glycolysis-related lncRNAs with prognostic significance.
Multivariate Cox regression was used to calculate their respective
coefficients (βI). Then, the formula of risk scoring model composed of
βI and lncRNA expression level (EXPI) was obtained as follows: Risk
score =
[MATH: ∑i=19(βEi*<
mi>xpi) :MATH]
. Then the risk score formula was used to calculate the risk score for
each BCa patient. Subsequently, the median of the risk score was set as
the cutoff value, and all eligible patients with bladder cancer were
divided into high-risk and low-risk groups. Next, Kaplan-Meier survival
analysis was performed to compare the differences in Overall Survival
(OS) between the two groups of BCa patients. Univariate Cox regression
and multivariate Cox regression analysis were used to evaluate whether
risk score, age, gender, and TNM were independent prognostic factors
for BCa. Clinical survival analysis was used to evaluate the predictive
power of the risk scoring model for different clinical subgroups. ROC
curve analysis was used to evaluate the efficiency of the prognostic
model. These results were also verified by the validation data.
Establishment and Evaluation of the Prognostic Model
A nomogram was then constructed to obtain the 3-year and 5-year OS of
BCa patients. This nomogram included the results of multivariate Cox
regression, risk scores, and clinical information such as gender, age,
and TNM stage. The calibration chart generated by the “rms” package of
the R software verified the predictive performance of the nomogram, and
the ROC curve analysis was used to evaluate the accuracy of the
nomogram. Then, a decision curve analysis (DCA) was performed to verify
the predictive effect of the prognostic model.
Gene Set Enrichment Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed
on the differentially expressed genes of BCa patients between the
high-risk and low-risk groups. By analyzing the tumor signaling
pathways and tumor progression processes involved in the differentially
expressed genes, their possible roles and mechanisms in the occurrence
and progression of BCa was further explored.
Results
Differentially Expressed Glycolysis-Related Genes and LncRNAs in Bladder
Cancer
The flow chart in [59]Figure 1 shows the general step-by-step process
of this study. First, the gene expression data, lncRNA sequencing data
and clinical data of the corresponding BCa patients (411 tumor tissues,
19 pericarcinomatous tissues) were downloaded from the TCGA database,
and glycolysis-related gene sets (GLYCOLYTIC_PRCOESS,
HALLMARK_GLYCOLYSIS, REACTOME_GLYCOLYSIS, GLYCOLYSIS_GLUCONEOGENE SIS,
and BIOCARTA_GLYCOLYSIS_PATHWAY) were downloaded from the GSEA
database. The “edgeR” package of R software was used to normalize the
entire data set, Set | log2FC | > 1 and false discovery rate (FDR) <
0.05 as the threshold to obtain glycolysis-related genes and lncRNAs
that were differentially expressed between tumors and pericarcinomatous
tissues ([60]Figure 2).
FIGURE 1.
[61]FIGURE 1
[62]Open in a new tab
Flow chart of key steps during the process.
FIGURE 2.
[63]FIGURE 2
[64]Open in a new tab
Differentially expressed glycolysis-related genes and lncRNAs in
bladder cancer compared to pericarcinomatous tissues. (A,B) Volcano
plot and heat map of glycolysis-related genes differentially expressed
in bladder cancer tissue and pericarcinomatous tissue. (C,D) Volcano
plot and heat map of lncRNAs differentially expressed in bladder cancer
tissues and pericarcinomatous tissue.
Next, Pearson correlation analysis was used to analyze the correlation
between the differentially expressed glycolysis-related genes and
lncRNAs, and to set correlation coefficients | R2 | > 0.3 and P < 0.05
for lncRNAs that were considered to be related to glycolysis. A total
of 59 glycolysis-related lncRNAs were identified as differentially
expressed in BCa tumors and pericarcinomatous tissues
([65]Supplementary Table 1). Subsequently, univariate Cox regression
analysis was performed, through which 10 glycolysis-related lncRNAs
related to the prognosis of BCa were identified ([66]Figure 3A).
Lastly, 9 key glycolysis-related lncRNAs associated with BCa prognosis
were further screened by LASSO regression ([67]Figures 3B,C).
FIGURE 3.
[68]FIGURE 3
[69]Open in a new tab
Screening glycolysis-related lncRNAs with prognostic value. (A) Risk
ratio Forest plot identified 10 glycolysis-related lncRNAs with
prognostic significance (P < 0.05). (B,C) LASSO regression further
screened out 9 glycolysis-related lncRNAs that were significantly
related to prognosis.
Construction of a Risk Scoring Model for 9 Glycolysis-Related LncRNAs
The data of bladder cancer patients with a follow-up time of more than
30 days was selected to construct a risk scoring model, where their
clinical parameters and pathological stages were showed in
[70]Supplementary Table 2. Patients were divided into training set and
validation set at a ratio of 2:1. The risk score formula (Risk score =
[MATH: ∑i=19(βiE*xpi) :MATH]
) was used to calculate the risk scores of all BCa patients. The median
of the risk score was set as the cutoff value, and all eligible
patients with BCa were divided into either high-risk or low-risk
groups. Results of Kaplan-Meier survival analysis showed that the
Overall Survival (OS) of the high-risk group of BCa patients was
significantly lower than that of the low-risk group, both in the
training set and the validation set (P = 1.546e-10, [71]Figures 4A,B).
The glycolysis-related lncRNAs [72]AC011503.2, [73]AC078778.1, and
ZNF667-AS1 were expressed at a lesser extent in tumor tissues of
patients with BCa and were thus considered to be protective factors for
the prognosis of bladder cancer, while glycolysis-related lncRNA
[74]AC099850.3, [75]AL589843.1, MAFG -DT, NR2F1-AS1, MNX1-AS1, and
[76]AC105942.1 were highly expressed in tumor tissues and were
therefore considered to be risk factors for the prognosis of BCa. ROC
curve analysis proved that the model has a good predictive performance
(3-year AUC = 0.753, 5-year AUC = 0.797) ([77]Figure 4C). We evaluated
the survival times of patients in the high- and low-risk groups and
found that mortality rates for patients with high-risk scores were
higher than those with low-risk scores ([78]Figures 4E,G). Heatmap
analysis was performed to reveal the expression profiles of the 9
glycolysis-related lncRNAs in high-risk and low-risk groups ([79]Figure
4I). These results were then validated using validation set data
([80]Figures 4B,D,F,H,J).
FIGURE 4.
[81]FIGURE 4
[82]Open in a new tab
Construction of risk score model. Kaplan-Meier survival analysis showed
that the Overall Survival (OS) of patients in the high-risk training
set (A) and validation set (B) was significantly lower than that in the
low-risk group. The ROC curve of the training set (C) showed that the
AUCs of the 3-year and 5-year OS were 0.753 and 0.797, respectively,
while the validation set (D) were 0.713 and 0.705, respectively.
Survival rate and survival status of bladder cancer patients in the
training set (E) and validation set (F). The distribution of 9
glycolysis-related lncRNAs risk scores for each patient in the training
set (G) and validation set (H). Heatmap of 9 glycolysis-related lncRNAs
in the low-risk group and the high-risk group in the training set (I)
and validation set (J).
Construction and Evaluation of Prognostic Model of 9 Glycolysis-Related
LncRNAs
A nomogram was then constructed to obtain the 3-year and 5-year OS of
BCa patients ([83]Figure 5A). This nomogram included the results of
multivariate Cox regression, risk scores, and clinical information such
as gender, age, and TNM stage. The calibration chart ([84]Figures 5D,F)
generated by using the “rms” package of the R software helped verify
the performance of the nomogram. Next, the ROC curve ([85]Figure 5B)
was used to obtain the area under the curve (3-year AUC = 0.781, 5-year
AUC = 0.821) to evaluate the accuracy of the nomogram. These results
further proved that the prognostic signatures of the 9
glycolysis-related lncRNAs were independent prognostic factors for
bladder cancer. Subsequently, the C index was calculated (training set
0.79, validation set 0.724), and Decision Curve Analysis ([86]Figures
5H,I) was conducted to verify the predictive effect of the prognostic
model. Finally, the validation set data was used to verify these
results ([87]Figures 5C,E,G). Subgroup analysis results ([88]Figure 6)
showed that the model had good predictive efficacy in all age groups,
in the T1-T2 and T3-T4 stages, with or without lymph node metastasis,
and with or without distant metastasis (P < 0.05) ([89]Figures 6A–I).
There was no significant difference observed in the female group (P =
0.107) ([90]Figure 6J). Possibly because of the small number of cases
in the group.
FIGURE 5.
[91]FIGURE 5
[92]Open in a new tab
Establishment and evaluation of the nomogram. (A) Nomogram for
predicting the 3-year and 5-year survival rates of bladder cancer
patients. (B,C) ROC curve analysis results showed that the AUC of the
training set (B) for 3-year OS and 5-year OS were 0.781 and 0.821,
respectively, and those in the validation set (C) were 0.783 and 0.748,
respectively. (D–G) The calibration curves for the 3-year and 5-year OS
of the nomogram in the training set (D,F) and validation set (E,G).
(H,I) The clinical decision curve for the 3-year OS (H) and 5-year OS
(I) of the nomogram.
FIGURE 6.
[93]FIGURE 6
[94]Open in a new tab
The results of subgroup analysis. Analysis showed that the OS of
bladder cancer patients in the high-risk group was significantly lower
than that in the low-risk group based on age and TNM stage (P < 0.05).
(A–I) In male patients, the OS in the high-risk group was also
significantly shorter than that in the low-risk group (P < 0.001),
while no significant difference was shown in female patients (P =
0.107) (J), possibly due to the relatively small number of cases.
Gene Function and Pathway Enrichment Analysis of 9 Glycolysis Related LncRNAs
Signatures
Kyoto Encyclopedia of Genes and Genomes enrichment analysis was
performed to explore the possible roles and mechanisms of the
differentially expressed glycolysis-related lncRNAs in the occurrence
and progression of tumors in the high-risk group and the low-risk group
of BCa patients, Analysis results ([95]Figures 7A–L) showed that these
glycolysis-related lncRNAs were involved in multiple signaling
pathways, such as: B-cell receptor, T-cell receptor, Hedgehog, MAPK,
WNT, Calcium, and Chemokine. They were also involved in oxidative
phosphorylation, ECM receptor interaction, and the development of basal
cell carcinoma. These findings may help researchers determine the
direction of further in-depth research in order to continue studying
the mechanism of how glycolysis-related lncRNAs affects the progression
of BCa.
FIGURE 7.
[96]FIGURE 7
[97]Open in a new tab
Enrichment analysis. (A–L) KEGG pathway analysis indicates that these
glycolysis related lncRNAs were involved in the B cell and T cell
receptor signaling pathway, Hedgehog signaling pathway, MAPK signaling
pathway, TGF-β signaling pathway, chemokine signaling pathway, EMC
receptor interaction, and oxidative phosphorylation processes.
Discussion
Accurately predicting the prognosis of BCa patients is of great concern
for clinicians in order to develop personalized treatment and review
plans. For patients with a high risk of recurrence, we must be more
cautious in the treatment strategy of bladder preservation, and for
low-risk patients, we need to preserve the bladder as much as possible
for the patients to help improve their postoperative quality of life.
To determine whether a BCa patient was of a high-risk type, current
clinical practice was mainly based on TNM staging, postoperative
pathological diagnosis, and molecular classification, but there was no
quantitative risk score or predictive model. Prognostic models
constructed by differentially expressed glycolysis-related lncRNAs have
been proven to have good predictive performance in gastrointestinal
tumors, breast cancer, gliomas, and other tumors. [98]Ho et al. (2021)
constructed a glycolysis-related lncRNA prognosis model that could
identify a subgroup of patients with different prognosis and different
degrees of immune infiltration. In bladder cancer, there is currently
no known glycolysis-related lncRNA prognostic model.
This prediction model was established based on the following
theoretical basis and previous research results: First, an increasing
amount of evidence shows that changes in glucose metabolism are a sign
of tumorigenesis. Some research has pointed out that the aerobic
glycolysis of bladder cancer cells is more active than that of normal
urothelial cells, and the increase in aerobic glycolysis is beneficial
to the proliferation, invasion, and metastasis of bladder tumor cells
([99]Berrondo et al., 2016). Other studies have found that the
activation of the Warburg effect was positively correlated with the
degree of tumor malignancy ([100]Agnihotri and Zadeh, 2016). Secondly,
previous studies have confirmed that lncRNAs were also involved in
glycometabolic reprogramming of tumor cells in BCa. Different lncRNAs
play different roles in the glycolysis process of bladder tumor cells.
Some lncRNAs accelerate tumor progression by promoting glycolysis, and
some of them inhibit glycolysis to suppress cancer. For instance,
[101]Li et al. (2014) found that the lncRNA UCA1 could promote
glycolysis by upregulation of HK2 through the mTOR-STAT3/miRNA and
miR-143 pathways. [102]Hu et al. (2017) found that the lncRNA CASC 8
may act as a tumor suppressor by reducing glycolysis in bladder tumor
cells.
In this study, 9 glycolysis-related lncRNAs were screened in relation
to the prognosis of BCa patients ([103]AC099850.3, [104]AL589843.1,
MAFG-DT, [105]AC011503.2, NR2F1-AS1, [106]AC078778.1, ZNF667-AS1,
MNX1-AS1, and [107]AC105942.1). Among them, glycolysis-related lncRNA
[108]AC011503.2, [109]AC078778.1, and ZNF667-AS1 were protective
factors for the prognosis of BCa, and glycolysis-related lncRNAs
[110]AC099850.3, [111]AL589843.1, MAFG-DT, NR2F1-AS1, MNX1-AS1, and
[112]AC105942.1 were risk factors for the prognosis of bladder cancer.
MNX1-AS1 regulates the expression of RAB1A in bladder tumor cells by
competitively binding with miR-218-5p, thereby promoting the
proliferation, migration, invasion, and epithelial-mesenchymal
transition of bladder tumor cells, which contributes to the growth and
metastasis of bladder cancer cells ([113]Wang et al., 2020).
Multivariate Cox analysis was used to calculate the regression
coefficient and construct the prognostic model. According to the risk
score, the median was used to divide the bladder cancer patients into
high-risk and low-risk groups. Analysis showed that the OS in the
low-risk group was longer than that in the high-risk group. According
to the results of multivariate Cox regression, a histogram was
established, including age, sex, TNM stage, and risk score. The
prediction performance of the model was verified by both a calibration
diagram and decision curve analysis, and the results demonstrated that
the model had good predictive performance.
In terms of GSEA functional enrichment analysis, it was found that
these glycolysis related lncRNAs were significantly enriched for the
regulation of multiple processes, such as the B cell and T cell
receptor signaling pathway, Hedgehog signaling pathway, MAPK signaling
pathway, TGF-β signaling pathway, chemokine signaling pathway, EMC
receptor interaction, and oxidative phosphorylation processes. [114]Li
et al. (2016) found that GALNT1-mediated glycosylation and activation
of Sonic Hedgehog signaling were involved in tumor-initiation and
self-renewal of bladder cancer stem cells ([115]Li et al., 2016). In
another study, Teng Hu et al. found that RSPO3 promoted the
aggressiveness of bladder cancer via the Wnt/β-catenin and Hedgehog
signaling pathways ([116]Chen et al., 2019b). Transcriptional
regulatory factor YAP was overexpressed in bladder cancer tissue and
promoted the spread of bladder cancer by affecting the MAPK pathway
([117]Qiu et al., 2020). All of these prior results further confirmed
that the glycolysis associated lncRNAs screened in the present study
play an important role in the occurrence and progression of BCa.
It should be pointed out that this study has some limitations: First,
this study was a retrospective study of patient data from the TCGA
database, which may lead to selection bias. Secondly, the patient
information obtained for this study lacks longer follow-up time;
therefore, the predictive model still needs to be further verified in
large-scale prospective clinical trials.
Conclusion
This study established a novel prognostic model of BCa based on 9
glycolysis-related lncRNAs, which can effectively distinguish high-risk
and low-risk bladder cancer patients from a new perspective. It can
also provide a new point of reference of ideas for clinicians to
formulate individualized treatment and review plans for patients with
different levels of risk. This study also found that some new
glycolysis-related lncRNAs were related to the prognosis of bladder
cancer, which points out the direction for the next step of
verification experiments and exploration of mechanisms.
Data Availability Statement
The original contributions presented in the study are included in the
article/[118]Supplementary Material, further inquiries can be directed
to the corresponding author.
Author Contributions
YY designed the study and took overall control of the manuscript. ZZ
analyzed and interpreted the data and drafted the manuscript. CL
performed the data analysis and helped write the manuscript. WL
collected relevant references, sorted out the data, and helped draft