Identification of novel prostate cancer genes in patients stratified by Gleason classification: Role of antitumoral genes.

Prostate cancer (PCa) is a tumor with a great heterogeneity, both at a molecular and clinical level. Despite its global good prognosis, cases can vary from indolent to lethal metastatic and scientific efforts are aimed to discern those with worse outcomes. Current prognostic markers, as Gleason score, fall short when it comes to distinguishing these cases. Identification of new early biomarkers to enable a better PCa distinction and classification remains a challenge. In order to identify new genes implicated in PCa progression we conducted several differential gene expression analyses over paired samples comparing primary PCa tissue against healthy prostatic tissue of PCa patients. The results obtained show that this approach is a serious alternative to overcome patient heterogeneity. We were able to identify 250 genes whose expression varies along with tissue differentiation-healthy to tumor tissue, 161 of these genes are described here for the first time to be related to PCa. The further manual curation of these genes allowed to annotate 39 genes with antitumoral activity, 22 of them described for the first time to be related to PCa proliferation and metastasis. These findings could be replicated in different cohorts for most genes. Results obtained considering paired differential expression, functional annotation and replication results point to: CGREF1, UNC5A, C16orf74, LGR6, IGSF1, QPRT and CA14 as possible new early markers in PCa. These genes may prevent the progression of the disease and their expression should be studied in patients with different outcomes.

differential expression analysis

principal component analysis

prostate cancer

prostate antigen levels

The Cancer Genome Atlas

trimmed mean of m method

INTRODUCTION

After lung cancer, prostate cancer (PCa) is the second most commonly diagnosed tumor in males worldwide, with 1 276 106 new cases in 2018. In recent years, an increase in cancer cases has been observed in Europe, with PCa the most frequent cancer in males and the third most common cause of cancer mortality.

PCa is the neoplasia with the highest hereditary component. Up to 15% of all PCa is attributable to high‐risk hereditary factors. Moreover, numerous genetic susceptibility markers for PCa have been identified using different approaches such as family‐based studies, candidate gene association studies, genome‐wide association studies and emerging RNA‐Seq technology. Some gene examples of the most widely identified in the expression profile of PCa include BRCA1, BRCA2, PTEN and HOXB13. , , However, only a few have been associated with cancer outcomes and often with inconsistent results. Specifically, antitumor genes could reveal new insights into molecular mechanisms of PCa pathogenesis and be a promising approach for gene therapy of PCa. , However, compared to other cancers such as breast or colorectal cancers, , little has been described regarding antitumor genes in PCa.

These knowledge gaps could be due to the high intra‐ and intertumor heterogeneity of PCa, both at a molecular and clinical level. , In this sense, the clinical behavior of PCa ranges from indolent to lethal metastatic disease, with the Gleason score as the strongest prognostic marker identified to date. , However, Gleason grading presents moderate reproducibility when used for the pathological classification of PCa in routine clinical practice. Thus, the development of new biomarkers to enable better PCa distinction and classification is a sizable challenge. , Identifying gene expression profiles could be useful in this context to overcome individual tumor variability by comparing healthy prostate tissue as control. ,

In this article, we analyze the genomic expression of PCa to identify genes implicated in the proliferation and metastasis progression of PCa by comparing gene expression of paired samples between prostate tumoral and healthy tissues. The final objective was to find new candidates for early PCa markers.

MATERIAL AND METHODS

Dataset description

The TCGA‐PRAD dataset was downloaded from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) with TCGAbiolinks package version 2.14.1. This dataset included RNA‐Seq expression data from 56 499 genes for 551 samples of primary PCa adenocarcinoma, from 495 patients, in addition to social patient data, Gleason related measures, prostate antigen levels (PSA), mutational landscape information and sample collection data. Procedures of samples are specified in the original paper. Methodological, biological gene expression and vial (aliquots) variability were assessed by principal component analysis (PCA) using the NOISeq package, version 2.30.0. Sample variability assessment by PCA dismissed batch effect and biological variability but showed an aliquot variability, thus using only vial A samples for further analysis (Figure S1).

Differential expression analysis

Gene expression analysis was performed using R language version 3.6.2 (https://cran.r-project.org/bin/windows/) and R Studio 1.2.5019; ggplot2 package version 3.3.3 (https://rstudio.com/) and pheatmap version 1.0.12 was used for graphical representation. For the discovery process, only paired samples (51) were selected from the TCGA‐PRAD dataset. Additionally, this set was divided into three further subsets according to the Gleason Score (“Gleason 6,” “Gleason 7” and “Gleason ≥8”) (Figure 1).

FIGURE 1

Methodology flowchart [Color figure can be viewed at wileyonlinelibrary.com]

The four subsets, the whole of the 51 paired samples and the three subsets after stratifying by Gleason (Figure 1) were independently analyzed using gene filtering with filterByExpr from the edgeR package with default parameters, version 3.28.1. Normalization was done using the trimmed mean of m method (TMM) with the NOISeq package, and values were later transformed in the logarithmic scale. Differentially expressed genes between tumor and healthy adjacent samples for each subset were obtained using a two groups single‐channel experimental design from the limma package, version 3.42.2. Inclusion criteria included a false discovery rate‐adjusted P‐value <.01 and at least a 2‐fold change.

Functional annotation

Annotation was carried out over the total pool of differential expressed genes obtained from the four subsets described in Section 2.2 (Figure 1).

Automatic annotation and enrichment analysis was done using ConsensusPathDB , across all available databases and selecting Gene Ontology categories at level 4.

Manually curated annotation was done using own‐made scripts for mining the Pubmed database (https://www.ncbi.nlm.nih.gov/pubmed/). The search terms for each gene differentially expressed were: “[‘prostate cancer’ AND ‘gen symbol’],” “[‘prostate adenocarcinoma’ AND ‘gen symbol’].” If no significant results were found, a secondary search was performed using the terms: [“cancer” AND “gen symbol.”] The articles retrieved were manually inspected, and gene functions were categorized according to known Prostate Cancer Hallmarks and stored in a database.

Replication datasets

No available datasets with paired data were found for PCa primary tumor. Consequently, to replicate our findings, three independent cohorts from the GEO database (http://www.ncbi.nlm.nih.gov/geo): GDS2545, GDS2546, GDS2547; were used. They included data from healthy tissue, healthy adjacent tissue, primary PCa tissue and metastatic tissue with 171, 167 and 164 samples, respectively (Table 1). Only normalized gene expression data was available, but data could undergo logarithmic transformation. Moreover, TCGA‐PRAD remaining samples were normalized and scaled, using the same parameters as described in Section 2.2 (Figure 1). Replication analysis was done with all samples, but the 51 paired samples were used for discovering the antitumoral differentially expressed gene set.

TABLE 1

Study and validation datasets number of samples

	Source	Healthy prostate tissue	Healthy prostate adjacent tissue	Primary prostate tissue	Metastatic prostate tissue
Discovery dataset	TCGA‐PRAD	–	51	51	–
Replication datasets	TCGA‐PRAD	–	–	444	–
	GDS2545	18	63	65	25
	GDS2546	17	59	66	25
	GDS2547	17	58	64	25

Replication samples' antitumoral expression tendencies were analyzed for different tumor stages. Statistical tests were performed using the statistical rstatix package version 0.6.0, Welch's ANOVA with Bonferroni correction was used for groups' comparison. Differences were considered statistically significant at the P < 0.01 level. Joint plots from TCGA‐PRAD and GEO datasets showing expression changes and statistical analysis are in Supporting Information 1.

RESULTS

As previously described in Section 2.2, in order to identify clear expression patterns related to PCa progression only paired samples were used for the differential expression analysis (DEA) (Figure 1). Paired samples analysis involves comparing healthy adjacent and tumor samples from the same patient, avoiding interpersonal variability interferences in detecting differential gene expressions due to disease progress. DEA analysis between healthy and tumoral samples were performed overall for 51 patients, subset named “All,” and for three additional subsets resulting from the stratification of the 51 patients according to their Gleason score and named: “Gleason 6” subset, “Gleason 7” and “Gleason ≥8” (Figure 1).

DEA over the subset “All” showed 195 differentially expressed genes, from which 107 were overexpressed and 88 underexpressed in tumoral tissue with regard to healthy adjacent tissue (Table S1). The subset “Gleason 6” showed 54 differentially expressed genes, 45 of them overexpressed. The number of differentially expressed genes was much higher in the subset “Gleason 7” reaching 192 genes, 115 overexpressed. However, the subset “Gleason ≥8” did not present any significant differentially expressed genes between tumoral tissue and adjacent tissue. The differentially expressed genes found in subsets “All,” “Gleason 6” and “Gleason 7” were represented in heatmaps and clustered samples by healthy and tumor tissue (Figures S2‐S4). In all cases, differential expression between healthy tissue and the corresponding subgroup tumoral tissue is visible. The interpersonal variability in expression independent of the tissue, is also noteworthy.

Functional annotation was done in a two‐step process. First using automatic annotation and then a manually curated process. The annotation was performed over the differentially expressed genes found in the subsets “All”, “Gleason 6” and “Gleason 7” (Figure 2), which made a total of 250 genes, 140 overexpressed and 110 underexpressed in tumor tissue vs paired healthy samples. Most of the 250 genes were protein‐coding genes (82.8%), followed by long noncoding RNA, lncRNA, (14%) and pseudogenes (3.2%).

FIGURE 2

Venn diagram of the 250 significant (P‐value <.01) differentially expressed genes between paired samples—healthy and tumor. The figure shows results found in the subsets “All,” “Gleason 6” and “Gleason 7.” (A) Results across the three datasets and their overlap for overexpressed genes; (B) Results across the three datasets and their overlap for underexpressed genes [Color figure can be viewed at wileyonlinelibrary.com]

Automatic functional analysis evidenced that, in general, PCa is related to the endocrine process and that over and underexpressed genes are associated with different biological processes. Moreover, most of the differentially expressed genes detected are specific to a particular Gleason score and are related to different biological processes. For example, concerning the six overexpressed genes from the subset “Gleason 6,” they were linked to cell cycle (G2/M transition, checkpoints and Rho GTPases); the 26 specific genes to subset “Gleason 7” were related to lipase activity and rhodopsin‐like receptors.

In the case of underexpressed genes, the differences were significant in the number of specific genes related to the Gleason score. While “Gleason 6” presented two exclusive underexpressed genes, “Gleason 7” presented 20. The two genes specific for “Gleason 6” could not be related to any biological process so far (C16orf74 and PHYHIPL), while the 20 genes of the “Gleason 7” were involved in androgen receptor, skeletal system development and neuron fate commitment pathways. Additional information to Figure 2 is available in Table S2.

Results obtained with the automatic approach highlight that PCa, regardless of the Gleason score, is associated with changes in the endocrine system. Moreover, bring the insight that genes differentially expressed in patients with Gleason score 6 are related to proliferation while in the case of Gleason score 7 they are mainly related to the progression of the disease.

Taking into account the delay in information flow from published resources to the specialized databases, we performed a manually curated functional annotation (Supporting Information 2). This approach allowed us to have a detailed understanding of the DEA genes' functions, and classify them according to cancer hallmarks for PCa and their antitumoral activity. It is noteworthy that using manually curated functional annotation there is still little information or a complete lack of information for 45 of these genes in relation to their function or specific role in cancer development (Figure 3). Moreover, we found that 43 genes had been previously described to be related to cancer (eg, omics experiments), but there was not enough information related to specific cancer hallmarks. However, there were 162 genes for which the manually curated annotation encountered new information and allowed their classification according to those hallmarks. From these 162 genes, 89 were directly linked to Prostate Cancer Hallmarks, whereas 73 were found to be related to hallmarks described in other cancers (eg, breast cancer, colorectal cancer or cancer in general). Most of 162 genes were related to metastasis (96) and proliferation (75), followed by resisting cell death (Figure 3). There were 110 genes strictly related to one hallmark and 52 linked to more than one.

FIGURE 3

Results of the manually curated functional annotation over the 250 significant (P‐value <.01) differentially expressed genes between paired samples—healthy and tumor. (A) Percentage of genes merged by the available information linked to their activity related to cancer and prostate cancer. (B) Number of genes whose activity could be pooled into cancer hallmarks in prostate cancer (inner circle) and cancer in general (outer circle), which includes prostate cancer. Genes could be related to more than one hallmark [Color figure can be viewed at wileyonlinelibrary.com]

Manually curated functional annotation revealed that, among the 162 genes with information, there were 39 antitumoral genes (24%). From these genes, 17 were identified to have known antitumoral activity in PCa by manual functional annotation (Table S5). Remarkably, the remaining 22 antitumoral genes found in our study have not been previously related to PCa, though all of them have been previously described to have antitumoral activity in other cancers (Table 2).

TABLE 2

Antitumoral genes described for the first time to be involved in PCa

Gene	Function	Cancer hallmark	PCa exp.
CGREF1	Inhibits AP‐1, c‐Jun, c‐Fos, p42/44 and p38 suppressing proliferation	Prol.	OE
FFAR2	Its loss promotes colon cancer by an epigenetic dysregulation of inflammation and its OE induced apoptosis in leukemia	A.I.D.	OE
SRARP	Its expression is inversely correlated with genes that promote cell proliferation	Prol.	OE
UNC5A	UE in most cancers. Reduces apoptosis when unbound to its ligand	R.C.D.	OE
TGM3	Its OE acts repressing EMT in colorectal cancer but promoting it in hepatocellular carcinoma	Met.	OE
TOX3	Inhibits cancer migration and invasion via transcriptional regulation of SNAI1 and SNAI2 but seems to play a critical role in progression of breast cancer	Met.	OE
PGM5‐AS1	Its downregulation inhibits the proliferation and metastasis via increasing miRNAs expression in colorectal cancer and esophageal squamous cell carcinoma	Met. Prol.	UE
CRABP2	Its loss reduces viability and proliferation and induces apoptosis, cytotoxicity and interferon‐signaling in malignant peripheral nerve sheath tumors	Prol., R.C.D.	UE
C16orf74	OE in pancreatic cancer, plays a crucial role in proliferation and invasion	Met., Prol.	UE
P2RX6	Promotes renal cancer cells migration and invasion modulating ERK1/2 phosphorylation and MMP9 signaling pathway	Met.	UE
LGR6	Downregulation inhibits proliferation and invasion and increases apoptosis by increased expression of Bcl‐2 and caspase‐3 and inhibition of PI3K/AKT	Met., Prol.	UE
MSLN	OE in multiple cancers, activates the NFB, MAPK and PI3K pathways and subsequently induce resistance to apoptosis or promote cell proliferation, migration and metastasis	Met., Prol., R.C.D.	UE
PDE1C	Silencing PDE1C significantly mitigates proliferation and EMT in glioblastoma	Met., Prol.	UE
ACTC1	ACTC1‐negative glioma had better prognosis than ACTC1 positive, and it was related to epithelial mesenchymal transition (EMT)	Met.	UE
EMX2OS	Induces proliferation, invasion and sphere formation in ovarian cancer via regulating the mir‐654‐3P/AKT3/PD‐L1 axis	Met., Prol.	UE
LINC00958	Its OE facilitates cell proliferation, migration and invasion	Met., Prol.	UE
IGSF1	Its knockdown could inhibit cell proliferation and significantly impair the migration and invasion in vitro in thyroid cancer	Met., Prol.	UE
SYT8	Inhibition was correlated with decreased invasion, migration and fluorouracil resistance in gastric cancer	Met.	UE
CHP2	Its OE promotes proliferation by activating AKT and suppression of FOXO3 transcription factor in breast cancer	Prol.	UE
CRABP2	Its loss reduces viability and proliferation and induces apoptosis, cytotoxicity and interferon‐signaling in malignant peripheral nerve sheath tumors	Prol., R.C.D.	UE
QPRT	It suppresses spontaneous cell death by inhibiting overproduction of active‐caspase‐3	R.C.D.	UE
PON3	Enhances cell death resistance	R.C.D.	UE
CA14	It is usually upregulated in cancer and linked with deadification	M.R.	UE

The antitumoral activity was scored according to genetic manipulation experiments and pathways activity annotations (Tables S3 and S4). From the 39 antitumoral genes, 34 are protein‐coding genes and 5 are regulatory lncRNAs, 29 presented differential underexpression.

The expression pattern found in the paired sample analysis for most of the 39 antitumoral genes agrees with the expression pattern (over‐ or underexpressed) previously described in the literature for other cancers. However, there were eight genes (TRPM8, CCN5, TGM3, TOX3, ACTC1, EMX2OS, PGM5‐AS1 and CRABP2) whose activity can be antitumoral or carcinogenic depending on the cancer type (Tables S3 and S4). For example, the overexpression of TGM3 (Table 2) in colorectal cancer presents antitumoral activity, while, in the case of hepatocellular carcinoma, it has a carcinogenic effect. We also found an increase in the overall number of antitumoral genes with cancer progression, Gleason 6 vs Gleason 7, as shown in Figure 4, and most were related to proliferation and metastasis inhibition.

FIGURE 4

Results summary of antitumoral differentially expressed genes found in “Gleason 6” and “Gleason 7” subsets merged by cancer hallmarks manually curated annotation [Color figure can be viewed at wileyonlinelibrary.com]

To corroborate our findings, we studied the expression levels of these antitumoral genes in the remaining tumoral samples from TCGA‐PRAD dataset not used for discovery and in the GEO datasets GDS2545, GDS2546 and GDS2547 (Figure 1). We could replicate the results in the remaining TCGA‐PRAD nonpaired samples in all cases but one, SYT7.

Regarding replicate analysis in the GEO samples, 31 of the 39 antitumoral genes were present in at least one of the three GEO replicate datasets (Table S6). Unfortunately, none of these experiments were paired samples. The replication percentage was lower in the case of the GEO expression datasets but still reached more than 74% replication in all samples where the gene was present. As an example of variable replication, we have PDE1C and ACTC1, where the expression profile of one replication sample agreed with TCGA_PRAD paired analysis results, but the other did not. In other cases, differences were so small between healthy adjacent tissue and primary PCa tissue that no conclusions could be made. However, in most cases, the tendency was to agree with the TCGA_PRAD paired sample analysis results. Results of all antitumoral genes are provided in Supporting Information 3.

Considering paired differential expression, manually curated functional annotation and replication results, we propose seven new potential PCa early‐stage antitumoral genes: CGREF1, UNC5A, C16orf74, LGR6, IGSF1, QPRT and CA14. Taking into account that some of the remaining genes could not be assessed in replica datasets, thus not being discarded for future analysis. Additionally, we found three genes—TGM3, TOX3 and CRABP—that present significant differential expressions between healthy and tumor samples and have been proved to play a role in other cancer processes (eg, colorectal, breast cancer). However, this evidence shows contradictory results (antitumoral or carcinogenic activity) depending on the type of cancer and thus the role of these three genes in PCa should be further investigated.

DISCUSSION

The use of differential gene expression between tumoral and healthy tissue over paired samples points to this approach as a serious alternative to overcome patient heterogeneity. It allowed the identification of 250 genes whose expression varies along with tissue differentiation—healthy to tumor tissue—, 161 of these genes are described here for the first time to be related to PCa. All these genes were found to be implied in diverse molecular mechanisms related to cancer. It is noteworthy that 15% of these 250 genes could be related to antitumoral activity, allowing the identification of new potential early PCa markers.

Paired DEA strategy highlighted clear differences between healthy and tumor tissue when considering all available paired samples. These differences persisted when samples were stratified by Gleason score. There were no genes differentially expressed between the “healthy tissue” and “Gleason ≥8” tumor when using this analysis. This result is not surprising due to the “field effect” that occurs within the normal epithelium adjacent to cancer. The absence of differentially expressed genes in this comparison could be explained by the histological changes and altered gene expression in the healthy tissue, an effect that has been already described to be more pronounced in aggressive PCa disease. , , Thus, we focused on the paired comparisons of healthy tissue with earlier differentiation stage tumors (“Gleason 6” and “Gleason 7”) in order to identify potential early PCa markers. This stratification revealed that the number of differentially expressed genes depended on the tumor differentiation stage, with a higher number in “Gleason 7” compared to “Gleason 6.” This result is in line with previous findings for PCa and breast cancer, relating gene expression changes to different tumor profiles and behavior. ,

Automatic functional analysis over the identified genes did not bring other insights than those previously known to be related to PCa in the literature. For example, the well‐known role of the endocrine process in the development and progression of PCa. The “Gleason 6” subset was related to proliferation mechanisms, and the “Gleason 7” subset with progression indicating a loss of local cell control in the latter. These functions agree with described changes in behavior between Gleason 6 stage and Gleason 7, specifically in Gleason 7 (4 + 3) in PCa.

The use of automatic functional annotation implies some drawbacks due to the delay of information flow between published results and databases. Thus, manually curated functional analysis allowed us, on the one hand, to have a deeper and actualized insight about the newly identified genes, and on the other hand, to link their activity, when data was available, to cancer hallmarks. This strategy identified the antitumoral activity of many newly discovered differentially expressed genes in PCa, such as the underexpressed C16orf74. There was no available information for this gene when using automatic functional annotation, but manual annotation showed that it plays a crucial role in proliferation and invasion when overexpressed in pancreatic cancer. This strategy also revealed that almost 20% of the genes were not protein‐coding but regulatory genes, lncRNA and pseudogenes. This finding is in line with the increasingly recognized role of noncoding genes in the course of PCa, and it supports the integration of these genes into precision medicine. However, to date, their role has been reserved only for research purposes and not for translational medicine.

More interesting is that this exhaustive approach revealed a high number of genes with potential antitumoral activity (n = 39) based on findings in the literature using manually curated annotation. Among them, we found 22 genes that have been identified previously as antitumoral in other cancers (eg, breast and colorectal cancer). However, they are described here for the first time to be differentially expressed between healthy and tumor tissue and in different grades of PCa differentiation. For example, CGREF1, whose overexpression in cancer suppresses proliferation. This does not necessarily mean that the 22 new antitumoral genes will act in PCa in the same way as described in the literature in other cancers. In fact, among these 22 genes, there are six (TGM3, TOX3, ACTC1, EMX2OS, PGM5‐AS1 and CRABP2) whose activity varies in the literature between carcinogenic or antitumoral, depending on the cancer type. However, this opens a new line of study for potential antitumoral genes in PCa in addition to the 17 already reported for this cancer.

In global terms, PCa has a good prognosis, with 5‐year survival data reaching more than 90%, but some cases have the potential to progress to a lethal outcome. , In fact, among the 495 patients in the TCGA‐PRAD subset with primary PCa, only 2% (10 cases) died: Gleason 6 (1 case) or Gleason 7 (4 cases 3 + 4, 1 case 4 + 3), Gleason ≥8 (2 cases), and unclassified (2 cases). Until now, it remains unclear how to identify between “low” and “high” risk subgroups through traditional clinical and pathologic criteria. , The search for new biomarkers has been suggested as a possible new alternative for classifying PCa risk. , We hypothesize that the antitumoral differentially expressed genes identified in our study could be playing a relevant role in the body defense against the tumoral process. We found an increase of antitumoral gene's expression with higher Gleason scores, mainly related to cancer hallmarks such as proliferation, metastasis, resisting cell death and castration resistance (Figure 4). Thus, these genes may be potential early biomarkers of changes in the clinical course of the disease.

A limitation of the present study is the low number of paired samples we used to identify new genes involved in PCa. However, this approach has been shown to overcome interpersonal PCa heterogeneity, and results could help find other ways of classifying the risk of PCa patients. A major challenge was to find replication datasets containing information about the identified genes; however, the results could be replicated in GEO databases up to 74%, and in the remaining TCGA‐PRAD samples could be achieved for all but one gene (SYT8). Another strength of our study is the use of manually curated functional annotation that allowed the discovery of antitumoral activity of genes differentially expressed in PCa with no annotation in current databases (eg, C16orf74).

In conclusion, taking into account paired differential expression, functional annotation and replication results, we propose the following genes as possible early markers in PCa: CGREF1, UNC5A, C16orf74, LGR6, IGSF1, QPRT and CA14. These genes may prevent the progression of the disease, and their expression should be studied in patients with different outcomes. Moreover, currently Gleason 6 and Gleason 7 treatment consists mainly in active surveillance, but these results already show molecular changes related to cancer proliferation and progression in these stages. Further studies with suitable datasets from large consortia are needed to confirm the role of these genes as early biomarkers.

CONFLICT OF INTEREST

The authors declare no potential conflict of interests.

AUTHOR CONTRIBUTIONS

Elisa Díaz de la Guardia‐Bolívar contributed to the conceptualization, collection of data, programming, preparation of the figures and tables, draft writing and editing of the manuscript. Rocío Barrios‐Rodríguez contributed writing, review and editing of the manuscript. Igor Zwir contributed reviewing and editing the manuscript and resources. José Juan Jiménez‐Moleón contributed drafting and reviewing the manuscript. Coral del Val acted as main supervisor, conceptualization, visualization and preparation of figures, draft, writing, review and edition of the manuscript. All authors read and approved the final manuscript. The work reported in the article has been performed by the authors, unless clearly specified in the text.

Appendix S1. Supporting Information.

Click here for additional data file.