Identification of novel targets for multiple myeloma through integrative approach with Monte Carlo cross-validation analysis.

More than one pathway is involved in disease development and progression, and two or more pathways may be interconnected to further affect the disease onset, as functional proteins participate in multiple pathways. Thus, identifying cross-talk among pathways is necessary to understand the molecular mechanisms of multiple myeloma (MM). Based on this, this paper looked at extracting potential pathway cross-talk in MM through an integrative approach using Monte Carlo cross-validation analysis. The gene expression library of MM (accession number: GSE6477) was downloaded from the Gene Expression Omnibus (GEO) database. The integrative approach was then used to identify potential pathway cross-talk, and included four steps: Firstly, differential expression analysis was conducted to identify differentially expressed genes (DEGs). Secondly, the DEGs obtained were mapped to the pathways downloaded from an ingenuity pathways analysis (IPA), to reveal the underlying relationship between the DEGs and pathways enriched by these DEGs. A subset of pathways enriched by the DEGs was then obtained. Thirdly, a discriminating score (DS) value for each paired pathway was computed. Lastly, random forest (RF) classification was used to identify the paired pathways based on area under the curve (AUC) and Monte Carlo cross-validation, which was repeated 50 times to explore the best paired pathways. These paired pathways were tested with another independently published MM microarray data (GSE85837), using in silico validation. Overall, 60 DEGs and 19 differential pathways enriched by DEGs were extracted. Each pathway was sorted based on their AUC values. The paired pathways, inhibition of matrix metalloproteases and EIF2 signaling pathway, indicated the best AUC value of 1.000. Paired pathways consisting of IL-8 and EIF2 signaling pathways with higher AUC of 0.975, were involved in 7 runs. Furthermore, it was validated consistently in separate microarray data sets (GSE85837). Paired pathways (inhibition of matrix metalloproteases and EIF2 signaling, IL-8 signaling and EIF2 signaling) exhibited the best AUC values and higher frequency of validation. Two paired pathways (inhibition of matrix metalloproteases and EIF2 signaling, IL-8 signaling and EIF2 signaling) were used to accurately classify MM and control samples. These paired pathways may be potential bio-signatures for diagnosis and management of MM.

Introduction

Multiple myeloma (MM) is an incurable cancer of plasma cells, caused by aberrant expansion of monoclonal plasma B cells in the bone marrow [1]. MM accounts for 10% of all hematological cancers, and is characterized by wide clinical and pathophysiologic heterogeneities, with lethal outcomes. The median survival time of MM patients is 7–8 years [2]. Thus, a better understanding of MM biology will aid in developing new therapeutic modalities that could potentially cure MM.

Many advances in the understanding of MM pathogenesis have been the result of major developments in genomic technologies [3], [4]. In recent years, genomic technologies identified certain disease-related biomarkers [5], [6]. For example, Leone et al. [7] have demonstrated that CDKN2C plays an important role in the progression and clinical outcome of MM. Jagani et al. [8] have indicated that Bmi-1 is crucial for MM growth. However, the reproducibility and overlap of the extracted genes are poor. Generally, gene biomarkers obtained from gene-based classification methods are often produced independently. Due to this, the gene signatures might not synergistically improve the overall classification ability.

In an attempt to overcome these shortcomings, it is important to understand the complicated interactions between genes, to help elucidate essential principles of cellular systems and the disease machinery [9]. To obtain a clear interpretation of genomic results, pathway analysis is the first criteria to identify abnormal pathways, to shed light on the potential biology of genes, thereby decreasing complexity and promoting explanatory power [10]. Moreover, several studies have reported that pathways-based classifiers are more reproducible and usually achieve better results, as compared to the single gene biomarkers-based classifier [11].

Remarkably, more than one pathway is involved in the development of cancer and its progression, due to the complex characterization of biological systems. Two or more pathways may be interconnected to further affect the disease onset, as functional proteins might participate in multiple pathways [12]. Therefore, identifying cross-talk between pathways is important to understand the molecular mechanisms of MM. Intuitively, different pathways influence each other, but at present, there is no reliable method to quantify the amount of cross-talk between paired pathways [13]. An integrative approach using Monte Carlo cross-validation has been created, to quantify the cross-talk between paired pathways.

Therefore, this study considers gene expression profile and biological pathway data as study objects, and utilizes Monte Carlo cross-validation analysis to detect pathway cross-talk in MM. The pathway cross-talk may be potential signatures for early detection and treatment of MM.

Materials and methods

The integrative approach using Monte Carlo cross-validation method comprised of four steps: Firstly, identification of differentially expressed genes (DEGs) between MM and healthy control samples was conducted. Secondly, the DEGs obtained were mapped to the pathways downloaded from the database of ingenuity pathways analysis (IPA), to reveal the relationship between DEGs and pathways enriched by these DEGs. A subset of pathways enriched by DEGs were also required. Thirdly, a discriminating score (DS) value for each pair of pathways was computed. Lastly, random forest (RF) classification was used to identify paired pathways with high cross-talk, and Monte Carlo cross-validation analysis was then repeated 50 times, to find the best paired pathways. During the analysis, all steps were repeated 50 times. After 50 runs, the top 10 paired pathways with the best AUC were extracted and were considered as significant paired pathways. These best paired pathways were tested with another independently published MM microarray data (GSE85837) using in silico validation.

Acquisition of gene expression profile

The gene expression library of MM (accession number: GSE6477) [14] was retrieved from Gene Expression Omnibus (GEO) database, based on the GPL96 platform of [HG-U133A]Affymetrix Human Genome U133A Array. There were 150 MM samples and 12 healthy donor plasma cell samples, in the GSE6477 profile. The repeated probes were first eliminated, after obtaining the microarray profile. The remaining probes were then mapped to the genomics, to further obtain human gene symbols. Finally, 12,437 genes were identified for subsequent analysis.

DEGs identification

During the research, “normalized quantile” was used to conduct normalization on mRNAs, with an average value determined across 162 samples [15]. Genes with values greater than 0.25-fold quantile average across all samples were extracted. Compared to several other estimators, the quantile-adjusted conditional maximum likelihood (qCML) was the most reliable in terms of bias, on a wide range of conditions, and performed best in smaller samples with a common dispersion quotient [16]. The qCML of edgeR package from Bioconductor was used to verify if these genes were expressed in a differential manner. Next, raw P values were corrected using a Benjamini-Hochberg method, based on a false discovery rate (FDR) [17]. Genes were considered differentially expressed when FDR was less than 0.001 and |log fold change (FC)|was more than 2.

Pathway enrichment analysis

Ingenuity pathways analysis (IPA), is widely utilized as a pathway database to analyze gene expression profile in the context of known biological responses and higher-order response pathways. In the present study, pathway enrichment analysis for DEGs was implemented using Fisher's exact test based on IPA tool, with the goal to extract significant pathways enriched by DEGs between MM and control samples. Hence, 589 biological pathways deposited in the IPA database were first downloaded. After the Fisher's exact test was applied to the genes in the IPA pathways and DEGs, the pathways enriched with P-value less than 0.01 were extracted. Later, the Benjamini-Hochberg procedure was used to correct the P values. Differential pathways were extracted based on the significance of FDR < 0.05.

DS calculation for pathway cross-talk

DS is an index used to compare the expression levels in the subgroup of samples showing amplification and in samples without amplification [18]. Thus, DS was employed to analyze the pathway cross-talk, in this paper. The DS was counted by comparing the gene expression levels of each paired pathway enriched by DEG in each sample, based on the description in the study by Cava et al. [19]. Generally, DS score implicates the relationships between paired pathways, and a larger DS suggests higher difference of activity between pathways.

Extracting the best paired pathways

RF created by Breiman [20], is a statistical method used to handle two issues of variable selection. To classify this methodology, an RF classification model was applied on the paired pathways based on the DS values of each sample. This helped to classify both MM and the control samples. AUC was calculated by a 10-fold cross-validation method, based on the following indexes: mtry and ntree. The mtry (the number of variables randomly sampled as candidates at each split) was equivalent to sqrt (p); p was the number of variables in the data matrix; and ntree (the number of trees grown) was equal to 500. All AUC values were then ranked in descending order, and the top 10 paired pathways were selected.

As documented in the validation analysis, the sample size in the training set was generally larger than that in the testing set. Of note, is the ratio of 6 to 4, which is a common distribution proportion. For example, Zhang et al. [21] randomly selected 60% as the training set and the remaining 40% as the testing data. Thus, in this study, the Monte Carlo cross-validation method was employed to randomly select 60% of the original microarray data comprising of training data, and the remaining 40% was assigned to the testing data. This step was repeated 50 times, randomly forming new training and test datasets each time. For each bootstrap, DEGs and pathway enrichment analysis for DEGs, and a DS for the paired pathways was conducted. For each bootstrap, a training set was used to detect a DS for the top 10 paired pathways with the best AUC value between the two groups. For each bootstrap, a testing partition was employed to confirm the top 10 paired pathways. At the end of 50 runs, the list of the top 10 paired pathways sorted by descending rank were selected, such that each pathway pair was extracted in 50 bootstraps. Ultimately, the top 10 paired pathways ranked for all 50 runs, were regarded as significant.

In silico validation with independent MM microarray data

To predict these best paired pathways, other MM data of 15 patients with MM and 13 control patients, was obtained from the publicly available microarray dataset GSE85837. For validation, all steps and selection criteria were the same as the above analysis.

Results

Identification of DEGs and pathway enrichment analysis

After normalization analysis, 1051 genes were obtained that had greater than 0.25-fold quantile average across all samples. From the differential expression analysis of MM and control samples, 60 DEGs with significance set as FDR < 0.001 and |log FC| > 2 (as shown in Table 1) were identified. Table 2 displays 19 significant pathways enriched by DEGs, their FDR values, number of genes for each pathway, and the count of common genes between DEGs and genes in the pathways.

Table 1

List of differentially expressed genes (DEGs).

Genes	LogFC	FDR	Genes	LogFC	FDR
RNASE2	−4.2668	1.25E-24	LMO2	−2.02467	1.26E-04
CLC	−4.94929	1.05E-13	CD24	−2.24512	1.65E-04
PRG3	−3.90837	2.93E-12	MRC1	−2.06384	2.36E-04
RNASE3	−4.31871	5.08E-12	CRISP3	−2.73151	2.55E-04
PRG2	−5.93443	1.01E-09	MS4A3	−2.17214	2.61E-04
MPO	−2.52454	8.15E-09	ZNF358	2.042213	2.66E-04
DEFA4	−2.22389	1.70E-08	MS4A6A	−2.21587	3.54E-04
ELANE	−4.39842	1.94E-08	IL4R	−2.1392	3.74E-04
EPX	−3.18113	9.02E-08	PXDC1	−2.08509	4.67E-04
ARMC7	2.314169	2.21E-07	CD320	2.241042	4.83E-04
CAMP	−3.0549	9.82E-07	CD163	−2.52643	5.17E-04
IGHG1	−3.7307	2.23E-06	LUM	−3.58203	5.20E-04
CTSH	−2.87845	3.55E-06	CXCL8	−2.40243	5.26E-04
AIF1	−2.16656	4.43E-06	NOD2	−2.95453	5.26E-04
LST1	−2.19262	8.38E-06	S100A8	−2.70373	5.61E-04
IGLJ3	−4.49176	1.13E-05	CXCL2	−2.71229	5.61E-04
CLPP	2.32109	1.65E-05	IGHD	−4.07084	5.84E-04
LTF	−3.89554	1.65E-05	HMOX1	−2.36868	5.99E-04
FABP4	−2.79511	3.18E-05	MS4A4A	−3.38126	6.11E-04
LYVE1	−2.9819	3.18E-05	HIF1A	−2.13129	6.18E-04
NRP1	−2.08155	3.34E-05	VNN2	−2.4071	6.18E-04
S100A12	−2.16333	3.34E-05	MAFB	−3.4547	6.41E-04
CEACAM8	−3.52708	4.20E-05	ALDH1A3	−2.01185	6.80E-04
IGLV1-44	−4.8402	6.64E-05	LHFP	−2.18934	6.98E-04
KCTD12	−2.56891	6.99E-05	S100A9	−2.3273	8.20E-04
CD14	−3.15551	7.82E-05	MNDA	−3.14832	9.14E-04
RPS11	2.079562	8.44E-05	HLA-DPA1	−2.64283	9.63E-04
IGK	−3.99784	9.70E-05	TGFBI	−2.93660	9.89E-04
P2RY13	−2.46684	1.05E-04	MAGEA4	2.361871	9.92E-04
DTX2	2.071491	1.05E-04	MROH7	−2.03615	9.95E-04

FC, fold change; FDR, false discovery rate.

Table 2

Pathways enriched by differentially expressed genes (DEGs).

Pathway	FDR	Genes in pathway	Number of common genes
Granulocyte Adhesion and Diapedesis	6.09E-11	163	4
EIF2 Signaling	2.88E-10	173	6
Atherosclerosis Signaling	1.32E-07	119	1
Hepatic Fibrosis / Hepatic Stellate Cell Activation	1.66E-07	137	2
LXR/RXR Activation	1.78E-07	121	1
T Helper Cell Differentiation	2.91E-06	62	3
Bladder Cancer Signaling	5.83E-06	86	5
Role of Macrophages	2.23E-05	282	4
Complement System	4.96E-05	32	3
Altered T Cell and B Cell Signaling in Rheumatoid Arthritis	6.80E-05	76	2
IL-8 Signaling	6.96E-05	183	1
PI3K Signaling in B Lymphocytes	1.01E-04	122	3
Regulation of the Epithelial-Mesenchymal Transition Pathway	1.17E-04	175	1
Neuroprotective Role of THOP1 in Alzheimer's Disease	3.38E-04	40	1
Role of Osteoblasts	5.18E-04	214	2
Acute Myeloid Leukemia Signaling	5.48E-04	76	3
ILK Signaling	7.71E-04	181	2
Macropinocytosis Signaling	8.53E-04	68	2
LPS/IL-1 Mediated Inhibition of RXR Function	9.14E-04	210	2

FDR, false discovery rate; common gene, the overlap between DEGs and genes in the pathway.

Selection of best paired pathways

To evaluate the classification ability of this methodology, an RF model was used on the paired pathways to compute the AUC values based on the DS values, using a 10-fold cross-validation method. Each pathway was then sorted out based on its corresponding AUC values. A total of 32 paired pathways with AUC not less than 0.850 were identified. In literature, AUC greater than 0.7 is regarded as good, and an AUC of 1.0 denotes a perfect classification [22]. Higher AUC values indicate good classification of disease, that is, higher AUC indicates a stronger pathway correlation with the disease. Thus, the focus was only on pathways with the 10 best AUC values, selected by Colaprico et al. [23]. The top 10 paired pathways that had the best classification ability for MM and control samples for all 50 runs, are displayed in Table 3. The paired pathways, inhibition of matrix metalloproteases and EIF2 signaling, obtained the best AUC value of 1.000. Moreover, the paired pathways of IL-8 signaling and EIF2 signaling also revealed good performance, with 0.975 AUC. Similar performance was found in the paired pathway of IL-8 signaling pathway and regulation of eIF4 and p70S6K signaling pathway with AUC of 0.939.

Table 3

Top 10 pairs of pathways with AUC value.

Pairs of pathways	AUC
(1a) Inhibition of matrix metalloproteases	1.000
(1b) EIF2 signaling
(2a) IL-8 signaling	0.975
(2b) EIF2 signaling
(3a) IL-8 signaling	0.939
(3b) Regulation of eIF4 and p70S6K signaling
(4a) IL-8 signaling	0.927
(4b) FGF signaling
(5a) PI3K signaling in B lymphocytes	0.910
(5b) EIF2 signaling
(6a) PI3K signaling in B lymphocytes	0.900
(6b) FGF signaling
(7a) Colorectal cancer metastasis signaling	0.898
(7b) EIF2 signaling
(8a) Altered T cell and B cell signaling in rheumatoid arthritis	0.897
(8b) EIF2 signaling
(9a) Inhibition of matrix metalloproteases	0.892
(9b) Regulation of eIF4 and p70S6K signaling
(10a) EIF2 signaling	0.890
(10b) TREM1 signaling

AUC, area under the curve.

Following this, the top 10 paired pathways with occurrence frequency of not less than 5 in the 50 runs, were identified. Based on this result, it was found that the PI3K signaling pathway of B lymphocytes and EIF2 signaling pathway were involved in 27 bootstraps; the paired pathways (inhibition of matrix metalloproteases and EIF2 signaling) appeared in 10 runs, and the paired pathways (PI3K signaling in B lymphocytes; regulation of eIF4 and p70S6K signaling) existed in 10 bootstraps; and IL-8 signaling and EIF2 signaling pathways were involved in 7 runs. Specific information is shown in Table 4.

Table 4

Top 10 pairs of pathways based on occurrence number not less than 5.

Pathway pairs	Total occurrence number
PI3K Signaling in B Lymphocytes; EIF2 Signaling	27
Inhibition of Matrix Metalloproteases; EIF2 Signaling	10
PI3K Signaling in B Lymphocytes; Regulation of eIF4 and p70S6K Signaling	10
Colorectal Cancer Metastasis Signaling; EIF2 Signaling	8
IL-8 Signaling; EIF2 Signaling	7
PI3K Signaling in B Lymphocytes; GABA Receptor Signaling	7
Role of NFAT in Regulation of the Immune Response; EIF2 Signaling	6
EIF2 Signaling; IL-17A Signaling in Fibroblasts	6
Granulocyte Adhesion and Diapedesis; EIF2 Signaling	6
Leukocyte Extravasation Signaling; EIF2 Signaling	5
EIF2 Signaling; IL-10 Signaling	5
MSP-RON Signaling Pathway; EIF2 Signaling	5

To validate the best paired pathways, independent MM data from the publicly available microarray dataset of GSE85837, was identified.

Table 5 shows the top 10 paired pathways with the best classification for MM and control samples for all 50 runs. The paired pathways, inhibition of matrix metalloproteases and EIF2 signaling pathways had the best AUC value of 0.971. Moreover, the paired TREM1 signaling and EIF2 signaling pathways also indicated high performance, with an AUC of 0.968. Similar performance levels were found in the paired IL-8 signaling and EIF2 signaling pathways, with an AUC of 0.961. These results further demonstrate that paired pathways are useful in diagnosing MM.

Table 5

Top 10 pairs of pathways with AUC value, validated using the other MM microarray data.

Pairs of pathways	AUC
(1a) Inhibition of matrix metalloproteases	0.971
(1b) EIF2 signaling
(2a) TREM1 signaling	0.968
(2b) EIF2 signaling
(3a) EIF2 signaling	0.961
(3b) IL-8 signaling
(4a) IL-10 signaling	0.956
(4b) EIF2 signaling
(5a) Antigen Presentation Pathway	0.955
(5b) LPS/IL-1 Mediated Inhibition of RXR Function
(6a) Caveolar-mediated Endocytosis Signaling	0.955
(6b) Agranulocyte Adhesion and Diapedesis
(7a) Regulation of eIF4 and p70S6K Signaling	0.953
(7b) Actin Nucleation by ARP-WASP Complex
(8a) Antigen Presentation Pathway	0.947
(8b) IL-6 Signaling
(9a) PI3K Signaling in B Lymphocytes	0.945
(9b) GABA Receptor Signaling
(10a) EIF2 Signaling	0.944
(10b) Actin Cytoskeleton Signaling

AUC, area under the curve.

Table 6 shows the top 10 paired pathways with occurrence frequency ≥ 5 in the 50 runs. From this table, it can be seen that the paired pathways of caveolar-mediated endocytosis signaling and agranulocyte adhesion, as well as diapedesis were involved in 45 bootstraps. The paired pathways (inhibition of matrix metalloproteases and EIF2 signaling) appeared in 43 runs, and the paired pathways involved in 37 bootstraps were antigen presentation and IL-6 signaling pathways. Also, IL-8 signaling and EIF2 signaling pathways were involved in 31 runs. Significantly, based on the above results, it was observed that a total of 5 paired pathways were the common ones, which included EIF2 signaling/inhibition of matrix metalloproteases, EIF2 signaling/IL-8 signaling, EIF2 signaling/IL-17A signaling in fibroblasts, EIF2 signaling/IL-10 signaling, and PI3K signaling in B lymphocytes/GABA receptor signaling.

Table 6

Top 10 pairs of pathways based on occurrence number not less than 5, selected from the validated microarray data.

Pathway pairs	Total occurrence number
Caveolar-mediated Endocytosis Signaling; Agranulocyte Adhesion and Diapedesis	45
EIF2 Signaling; Inhibition of Matrix Metalloproteases	43
Antigen Presentation Pathway; IL-6 Signaling	37
EIF2 Signaling; IL-8 Signaling	31
EIF2 Signaling; Regulation of eIF4 and p70S6K Signaling	31
Antigen Presentation Pathway; LPS/IL-1 Mediated Inhibition of RXR Function	21
EIF2 Signaling; IL-17A Signaling in Fibroblasts	11
EIF2 Signaling; IL-10 Signaling	9
Regulation of eIF4 and p70S6K Signaling; Actin Nucleation by ARP-WASP Complex	9
PI3K Signaling in B Lymphocytes; GABA Receptor Signaling	5

Discussion

Pathway analysis has become the preferred approach to understand biological processes of genes, as this kind of analysis promotes scientific validity and explanation [24]. Generally, different pathways involved in several biological processes act in a synergistic way. Extraction of pathway cross-talk makes data available to study pathway functions in greater detail, relative to single pathways [25]. However, traditional pathway analysis methods mainly focus on individual dysregulated pathways, while pathway interactions are frequently not considered [26]. Thus, in this paper, an integrative approach with Monte Carlo cross-validation analysis based on the DS values, was used to identify the best paired pathways that could distinguish MM from control samples. Ultimately, a total of 32 paired pathways with AUC not less than 0.850, among MM and control samples were identified. Nevertheless, only the top 10 paired pathways with higher AUC were focused on.

In this study, inhibition of matrix metalloproteases and EIF2 signaling pathways obtained the best AUC value of 1.000, with an occurrence frequency of 10, in analyzing GSE6477. In silico validation, inhibition of matrix metalloproteases and EIF2 signaling obtained the best AUC value of 0.971, with an occurrence frequency of 43. Matrix metalloproteinases, as a group of zinc-dependent endopeptidases involved in the degradation of extracellular matrix (ECM), exert important functions in tissue remodeling, and are connected with various physiological processes including migration, invasion, angiogenesis, and regulation of inflammatory processes [27], [28], [29]. Angiogenesis and inflammation are the hallmarks of cancer [30], [31]. Moreover, suppressing matrix metalloproteases initially appeared to be a promising therapeutic approach for cancer [32]. Significantly, the functions of matrix metalloproteases are influenced by reactive oxygen species (ROS). The inflammatory response at the tumor site creates large amounts of ROS that are produced by activated neutrophils and macrophages. Moreover, EIF2 alpha kinase has been indicated to control ROS levels [33]. Hence, the cross-talk between inhibition of matrix metalloproteases and EIF2 signaling pathways may play a crucial role in the etiology of MM.

Another set of paired pathways, notably IL-8 signaling and EIF2 signaling pathways, were involved in this study for GSE6477 in 7 runs, with an AUC of 0.975. In silico validation using GSE85837, IL-8 signaling and EIF2 signaling pathways were obtained, with an AUC value of 0.961 and an occurrence frequency of 31. As reported, multiple cytokines or chemokines, play an important role as mediators of paracrine signals between the tumor and multifarious components of the tumor microenvironment, which eventually results in tumor growth and progression [34], [35].

As a proinflammatory cytokine, IL-8 is a soluble mediator released by tumor cells that simultaneously exerts key functions in autocrine and paracrine modes in the tumor microenvironment. Significantly, in cancer tissues, tumor-derived IL-8 has been indicated to enhance tumor cell survival, proliferation, as well as migration based on autocrine activity, while inducing angiogenesis in endothelial cells due to paracrine activity [36], [37]. Moreover, there is significant epithelial–mesenchymal transition (EMT) during carcinoma progression [38]. Of note, Romaine et al. [39] have implicated that IL-8 signaling pathway is critical in the EMT process of human carcinoma cells. In addition, in MM patient plasma cells, IL-8 receptors (CXCR1 and CXCR2) were observed, and IL-8 parallel MM disease activity was related to bone marrow angiogenesis [40]. Furthermore, EIF2 signaling pathway has been reported to mediate proinflammatory cytokine expression [41]. Thus, it is speculated that the regulation of IL-8 and EIF2 signaling pathways may be important in the onset and progression of MM.

However, several limitations should be noted. To begin with, this study was a preliminary study of molecular mechanisms of MM. This work was conducted based on the bioinformatics approach, while the conclusions have not been verified using experiments. Further investigations are required to reveal changes in MM pathways, using animal experiments or patients’ tissues.

In conclusion, this analysis sheds new light on the involvement of pathway cross-talk in the pathology of MM. However, the identified pathway cross-talk require more functional studies in later work.

Conflicts of interest

The authors have declared that they have no conflict of interest.