Downregulation of Brain Enriched Type 2 MAGEs Is Associated With Immune Infiltration and Poor Prognosis in Glioma

Melanoma associated antigen (MAGE) is an extensively studied family of tumor-associated genes that share a common MAGE homology domain (MHD). Based upon their expression pattern, MAGE genes have been broadly classified into type 1 MAGEs (T1Ms) and type 2 MAGEs (T2Ms) categories. Interestingly, several T2Ms are highly expressed in the brain and involved in the regulation of neuronal development, differentiation, and survival. Available literature suggests possible tumor suppressor functions of a few T2Ms, while information available about their expression, regulation, and clinical significance in glioma is scanty. This prompted us to perform a comprehensive analysis of T2M expression in glioma. Gene expression data from glioma datasets: Oncomine, TCGA, and REMBRANDT study, were used to assess the mRNA expression of T2M genes (MAGED1, MAGED2, MAGED3, MAGED4, MAGED4B, MAGEE1, MAGEE2, MAGEF1, MAGEH1, MAGEL2, NSMCE3, and NDN), and their association with clinical characteristics and composition of the tumor microenvironment. Further, mutation, copy number alteration, and DNA methylation data from TCGA were assessed for determining potential mechanisms of T2Ms expression in glioma. Expression analysis revealed overexpression of MAGED subfamily genes in glioma, while other genes of this family exhibited reduced expression in advanced grades of this malignancy. Further, the expression of T2Ms exhibited varying extent of positive correlations with each other. Amongst downregulated T2Ms, MAGEH1 expression exhibited negative correlations with DNA methylation. Additionally, genes associated with MAGEH1 were enriched in Myc and Hedgehog signaling. Furthermore, T2Ms downregulation was associated with immune infiltration in glioma tissues and poor overall survival of glioma patients. In multivariate Cox regression analysis, MAGEH1 emerged as an independent prognosticator in lower grade glioma. Conclusively, these results suggest that expression of T2Ms is associated with important clinical and molecular features in glioma. Mechanistic studies may further provide novel insights into their role in glioma progression.


INTRODUCTION
Gliomas are a group of heterogeneous primary malignant brain tumors with a dismal outcome (1). Histopathological classification of glioma has been constantly evolving with the most recent WHO classification, 2016, classifying gliomas as astrocytic, oligodendroglial or mixed oligo-astrocytic, and further graded as WHO grade I and II (low grades), III (anaplastic) or IV (glioblastoma multiforme) (2). Molecular alterations in glioma including genetic as well as epigenetic changes have been constantly determined and integrated into common clinical practice (2,3). Isocitrate dehydrogenase gene (IDH) mutation and co-deletion of chromosome 1p and 19q (1p/ 19q codeletion) are frequently observed in low grade gliomas and associated with better patient prognosis (4). IDH mutation confers global epigenomic changes, commonly called the CpG island methylation phenotype (G-CIMP) leading to suppression of tumor suppressor genes (5). Also, methylation O-6methylguanine-DNA methyltransferase (MGMT) gene promoter is associated with better response to temozolomide (TMZ) chemotherapy (6). While these biomarkers help in the management of glioma with the current treatment regimen, identification of novel molecular features may provide better therapeutic opportunities in glioma.
Melanoma associated antigen (MAGE) gene family consists of more than 50 genes that share a common MAGE homology domain (MHD) (7,8). MAGE genes evolved through transposition and segmental duplications in the genome (9). Based on chromosomal location and pattern of gene expression, members of this family are further subdivided into two types. Type 1 MAGEs (T1Ms) include MAGE-A, B, and C subfamilies. The human T1Ms are present as clusters on the X chromosome. T1Ms display germ cell and cancer specific expression, and thereby fulfill the criteria of cancer testis antigens. T1M peptides are recognized by cytotoxic T lymphocytes (CTLs) in a variety of cancers and serve as ideal candidate antigens for tumor vaccines (10)(11)(12). Human Genome Organisation (HUGO) has approved the inclusion of twelve Type 2 MAGEs (T2Ms), including MAGED1, MAGED2, MAGED3 (TRO), MAGED4, MAGED4B, MAGEE1, MAGEE2, MAGEF1, MAGEH1, MAGEL2, NSMCE3, and NDN (Necdin). T2Ms are not restricted to X chromosome and are expressed in a variety of tissues.
The current knowledge of the physiological functions of the MAGE genes is mostly related to T1Ms. Several T1Ms form complex with E3 RING ubiquitin ligases to form MAGE-RING ligases (MRLs) that are involved in ubiquitination mediated protein degradation and other cellular processes such as regulation of transcription and cell cycle (8). T1Ms have also been implicated in apoptosis, cancer cell invasion (13), stem cell maintenance, and DNA repair (14). However, information about the functions of T2Ms is limited. MAGED1, MAGEH1, and NDN are highly expressed in the brain and involved in neuronal development (8,15,16). MAGEE1 and MAGED4 were initially identified as brain specific members of the MAGE gene family and were significantly enriched in glial tissues (17,18). However, Necdin is a neuron specific growth suppressor, is downregulated in tumors, including glioma and tumor cell lines, thereby confirming its tumor suppressor functions (19,20). Low frequency mutations of MAGEH1 in glioma have been shown to affect its nuclear localization (21). Furthermore, aberrant MAGEH1 expression has been linked to dementia (22). These studies suggest crucial roles of T2Ms in neuronal growth, survival and possibly in the pathogenesis of central nervous system diseases, including glioma.
Although wide variations in the nucleotide sequences of regulatory regions of T1Ms enables their differential expression (9), epigenetic alterations including DNA methylation and histone modification underlie the co-expression of several T1Ms and other tumor associated antigens in cancers (23). On the contrary, all T2Ms display striking sequence homology in their regulatory regions which possibly permit a common mechanism to control their expression (9). Although contradictory reports also exist regarding maternal imprinting of the Necdin gene (24,25), epigenetic regulation of MAGED4 in glioma is documented (18). While current evidence suggests crucial roles of T2Ms in normal brain functions and neurological disorders, no systematic study has been done to investigate the expression, regulation, cellular function, and clinical significance of T2Ms in glioma. Therefore, the current study was undertaken to investigate the regulatory mechanism responsible for altered expression of T2Ms, their cellular functions, and clinical significance in glioma using large scale multi-omics datasets.

5'Azacytidine Treatment and qRT-PCR
Human glioblastoma cells (2X 10 4 ) were plated in each well of 6 well dish. The next day the cells were treated with 1 µM Azacytidine (#A2385, Sigma-Aldrich). After 48 h the cells were washed twice with ice-cold PBS and processed for total cellular RNA isolation using Trizol reagent (Thermo Fisher Scientific Inc.) as per manufacturer's protocol. One µg of total RNA was reverse transcribed using random hexamers, dNTPs, and M-MuLV reverse transcriptase enzyme (Fermentas, USA). The MAGEH1 expression was determined by quantitative Real-time PCR using primers Forward, 5'-CGGAGCAATTTTCAGG GCAC-3'; Reverse, 5'-AGCACTTTCCAGACCAGAGC-3'. The mRNA expression of MAGEH1 was normalized to 18S ribosomal RNA using primers: Forward, 5'-GTAACCCGTTG AACCCCATT-3'; Reverse, 5'-CCATCCAATCGGTAGTAGCG-3'. Ct values for each PCR were analyzed by the 2 -DCt method. Total cellular RNA isolated from vehicle-treated glioblastoma cells was processed identically and served as control. Depicted results were drawn based on three biological replicates.

Correlation and Pathway Enrichment Analysis
For correlation analysis in TCGA-LGG and GBM data, the correlation module of the GlioVis tool was used. Gene expression data from three randomly chosen TCGA datasets including TCGA skin cutaneous melanoma, TCGA breast cancer and TCGA lung squamous cell carcinoma, along with CCLE pan-cancer cell line data of 1156 cell lines, were downloaded from cBioPortal website (https://www.cbioportal.org/) (29,30). Spearman's correlation analysis was performed followed by heatmap generation and hierarchical clustering using HemI software (31). The default parameters of hierarchical clustering using the average linkage method and Pearson distance were used. Similarly, whole transcriptome correlations of MAGEH1 were extracted for TCGA-LGG dataset using cBioPortal. After applying a filter for a cutoff of FDR corrected p-value of 0.05 for Spearman's r value, a total of 13,536 genes with Spearman's r value ranging from 0.689 to -0.660 were available for gene set enrichment analysis in GSEA software (Broad Institute, http:// www.broad.mit.edu/gsea/) (32). Predefined molecular signature database hallmark gene set (version 7.1) was used as a reference gene set for pathway enrichment (33).

DNA Methylation Analysis
DNA methylation of T2Ms in TCGA cancer datasets was estimated and visualized using MEXPRESS web server (https:// mexpress.be) (34,35). The MEXPRESS web server uses DNA methylation data of cancer and normal tissues from TCGA datasets, which were originally developed on the "Illumina 450k Beadchip" platform. The predesignated methylation probes for each gene were taken into consideration.

Survival Analysis
Kaplan-Meier survival analysis was performed using a survival tool available at the GlioVis web server, which utilized the "survival" package in R to generate Kaplan-Meier plots. Hazard ratios are determined to utilize the "coxph" function from the "survival" package. For Kaplan-Meier analysis, Patients were distributed in high and low expression groups based on optimal cutoff determined using maximally selected rank statistics (maxstat) function for continuous variables, as provided in the "survminer" package. For survival analysis using univariate and multivariate Cox regression, MAGEH1 and MAGED1 expression was taken as a continuous variable. For TCGA-LGG and TCGA-GBM datasets, information on overall survival, disease-specific survival, progression-free interval, and the disease-free interval was available. For, TCGA-GBM data, disease-free interval events were excluded as suggested previously (36).

TIMER Analysis
Tumor immune estimation score (TIMER) database (https:// cistrome.shinyapps.io/timer/), which utilized the RNA sequencing data from TCGA for estimating the correlation between gene expression and level of tumor-infiltrating immune cells (37). We utilized TIMER to calculate the association between gene expression of T2Ms with tumor purity and infiltration of immune cells including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells in LGG and GBM datasets.

Statistical Analysis
Data analysis was performed using Graphpad (version 6) and Stata software (version 11). Chi-square test was used to calculate the expression fold change with threshold p-value <0.001) between normal and glioma in Oncomine datasets analysis. Mann-Whitney U-test was used for comparison among histological subtypes, molecular subtype and grades (**p<0.001; **p<0.01; *p<0.05; ns, p>0.05). Pearson correlation was used to calculate the correlation of DNA methylation of T2M genes to its expression in LGG and GBM dataset. Kaplan-Meier survival analysis was performed using Wilcoxon and logrank test, p-value <0.05 was considered as statistically significant.

Expression Pattern of T2Ms in Glioma
To compare the expression of T2Ms in glioma with normal brain tissue, we utilized publicly available gene expression datasets: Oncomine, TCGA, and REMBRANDT study. The outline of the study has been depicted in Figure 1A. Oncomine analysis was performed to compare gene expression in multiple datasets in parallel to get reliable information regarding the change in T2Ms gene expression ( Figure 1B). Interestingly, Oncomine analysis revealed that in glioma, several MAGE genes, including MAGED1, MAGED4, and MAGED4B exhibited elevated expression ( Figure 1B). In contrast, an opposite pattern was observed for TRO (MAGED3), MAGEE1, MAGEE2, MAGEH1, MAGEL2, and NDN, which exhibited reduced expression in tumors compared to normal brain tissues. NSMCE3, MAGED2, and MAGEF1 did not exhibit a clear pattern of overexpression/ downregulation in Oncomine analysis.  We further analyzed the coexpression pattern of T2Ms in TCGA datasets of glioma and the results are presented in Figure 1C. This analysis revealed a strong positive correlation among MAGED1, MAGED4, and MAGED4B (r>0.4) while they exhibited a comparatively weak correlation (r <0.3) to other MAGEs. Similarly, TRO, MAGEE1, MAGEH1, MAGEL2, and NDN exhibited positive correlations (r >0.3), while others exhibited varying levels of weaker correlations (r <0.3) ( Figure 1C). We further performed similar co-expression analysis on gene expression datasets for lower grade glioma (LGG) and glioblastoma multiforme (GBM), separately, and observed broadly similar coexpression pattern of T2Ms in LGG and GBM (Supplementary Figures S1A, B). Furthermore, we also assessed the correlation of T2Ms expression in other TCGA cancer datasets to determine whether these coexpression patterns are unique to glioma or are common to most cancers. Interestingly, we observed that in TCGA datasets of cutaneous melanoma (TCGA-SKCM), breast cancer (TCGA-BRCA), and lung squamous cell carcinoma (TCGA-LUSC), T1Ms and T2Ms form distinct clusters (Supplementary Figures S2A-C). However, within the cluster, they display varying degrees of positive correlations, thereby supporting the possibility of common regulatory mechanisms for their expression. These co-expression patterns also reflected well in a pan-cancer cell line gene expression dataset from Cancer Cell Line Encyclopedia (CCLE, Supplementary Figure S2D). Therefore, expression of T2Ms is highly coordinated in two distinct subgroups of T2Ms in glioma, one of which is overexpressed while the other exhibits downregulation in glioma. In view of the reported overexpression of MAGED4 in glioma (18), we focused our further analysis on downregulated T2Ms using MAGEH1 as a representative member of this highly co-expressed T2M subgroup. Interestingly, all T2Ms which were found to be downregulated in the expression analysis exhibited preferential enrichment in the brain, as determined by the Genotype-Tissue Expression (GTEx) project whole-body gene expression data (Supplementary Figures S3A-E).

Expression of T2Ms in Different Histological Subtypes
REMBRANDT dataset contains normal brain tissues, therefore it was used to validate the expression of different T2Ms in glioma compared to the normal glial tissue. For all subsequent analyses, TCGA along with the REMBRANDT dataset was utilized for clinical correlations. REMBRANDT dataset revealed reduced MAGEH1 expression in glioma compared to its normal counterpart. Interestingly, analysis of both REMBRANDT and TCGA datasets brought out a steady decrease in MAGEH1 expression with advancing grades of glioma (Figures 2A, B). Furthermore, the analysis of histopathological data revealed a strong correlation of MAGEH1 downregulation with astrocytoma histology (Figures 2C, D). Among TCGA defined molecular subtypes of GBM described by Wang et al. 2017 (38), MAGEH1 exhibited the highest expression in the proneural subtype while classical and mesenchymal subtypes exhibited its comparable expression ( Figures 2E, F). Regarding epigenetic alterations, higher MAGEH1 expression was observed in G-CIMP GBMs compared to non-G-CIMP GBMs in the REMBRANDT dataset (p<0.01, Figure 2G). In TCGA dataset, IDH mutant tumors, which are closely associated with the G-CIMP phenotype, exhibited higher MAGEH1 expression than non-IDH mutant gliomas (p<0.001, Figure 2H). Furthermore, IDH mutant tumors that harbor 1p/19q codeletion expressed higher levels of MAGEH1 compared to the non co-deletion tumors with IDH mutation (p<0.001), which together indicate epigenetic regulation as possible mechanisms of T2Ms expression.

Alterations in T2M Genes in Glioma
To determine the possible contribution of genetic alterations in altered expression of T2Ms in glioma, copy number variation (CNV), and mutation data of TCGA-LGG (lower grade glioma consisting grade II +III) and TCGA-GBM (glioblastoma) dataset were assessed. Our results showed heterozygous deletion and copy number gain to be common events in T2M genes. However, they rarely undergo mutations, or display amplification and/or deep deletion events ( Figure 3A). Interestingly, it was also observed that some glioma patients harbor heterozygous deletion of multiple T2Ms. Further, to determine the effect of copy number variation (CNV) on T2Ms expression, we compared their expression in glioma patients based on CNV data and observed no difference in expression of MAGEH1 ( Figure 3B) and MAGEE1 ( Figure 3C) between tumors with diploid and shallow deletion status, and between diploid and copy number gain status for the respective genes. In TCGA-LGG dataset, shallow copy number deletion was indeed associated with higher expression of MAGEH1 compared to tumors with diploid status for this gene (p=0.0009, Figure 3D), while no such difference was observed for MAGEE1 ( Figure 3E).
Further, we separately analyzed the prognostic significance of MAGEH1 in GBM using REMBRANDT and TCGA datasets. This analysis also brought out the association of reduced MAGEH1 expression with poor OS in both REMBRANDT (HR=1.55, 95% CI=1-2.39, p<0.05, Figure 7A), and TCGA (HR=2.09, 95% CI=1.18-3.68, p<0.01, Figure 7B) dataset. Since IDH mutation are associated with poor patient survival   Figure 7C) or mutant (HR= 1.94, 95% CI=1.2-3.16, p<0.01, Figure 7D) genotypes. Additionally, we performed univariate and multivariate Cox regression analysis for prognostic significance of relevant clinicopathological features along with MAGED1 and MAGEH1 expression in TCGA-LGG (a combined group of grade II and grade III) and TCGA-GBM datasets, separately. Overall survival (OS), disease specific survival (DSS), disease free interval (DFI) and progression-free interval (PFI) were taken as outcome measures. In the univariate analysis for LGG, we used MAGEH1 expression, MAGED1 expression, age, gender, tumor histology, WHO grade, IDH mutation status, 1p/19q codeletion status, and MGMT promoter methylation status to analyze the association of these variables to patient outcome. Among these, age, tumor grade and astrocytoma histology were associated with poor OS, DSS and PFI, while IDH mutation status, 1p/19q codeletion status, and MGMT promoter methylation status were associated with better OS, DSS, and PFI (Table 1). Additionally, IDH mutation status was also associated with DFI. Interestingly, MAGEH1 expression, but not MAGED1 expression was associated with better OS, DSS, and PFI in LGG. Furthermore, multivariate analysis revealed that higher MAGEH1 expression is also an independent prognostic indicator of the better OS, DSS, and PFI in LGG ( Table 2). On the other hand, multivariate analysis using LGG dataset for MAGED1 suggested no association of MAGED1 with patient survival ( Table 3).
Next, we performed survival analysis for TCGA-GBM data with MAGEH1 and MAGED1 expression, along with age, gender, molecular subtypes (38), IDH mutation status, and MGMT promoter methylation status ( Table 4). In univariate analysis, we observed an association of higher MAGEH1 expression with favorable PFI, while no association of MAGED1 was observed with the disease outcome. Additionally, multivariate analysis of the association of MAGEH1 with PFI revealed only marginal significance (p=0.07), while no association of MAGEH1 was observed with OS, DSS, and DFI ( Table 5). Similar to univariate analysis, multivariate analysis for prognostic significance of MAGED1 in TCGA-GBM dataset revealed no association of MAGED1 expression with patient survival ( Table 6).

MAGEH1 Associated Pathways in Glioma
Considering its prognostic significance in LGG, we determined the potential functional associations of MAGEH1 in LGG by performing gene set enrichment analysis over whole transcriptome correlations of MAGEH1 expression in TCGA-LGG dataset. Interestingly, positively correlated genes were    highly enriched in Myc target genes ( Figures 8A, B) and Hedgehog signaling ( Figure 8C). Further, genes negatively correlated to MAGEH1 expression were enriched in different immune associated pathways ( Figures 8D-G, I-K) and cancer associated pathways including epithelial to mesenchymal transition ( Figure 8H) and apoptosis ( Figure 8L).

DISCUSSION
Despite of the recent advances in the understanding of glioma biology, limited improvements have been observed in patient     survival (1). The current study was undertaken to analyze the clinical utility of T2Ms expression in glioma. We further focused our analysis to determine the coexpression pattern, underlying regulatory mechanism, and association of T2Ms with tumor immunity. Oncomine analysis provided the benefit of comparing gene expression between tumor and normal brain tissues from multiple datasets, in parallel. This would compensate for the dataset dependent variations and only highly significant observations were considered using a stringent cutoff value of p<0.001. This analysis revealed that some MAGEs including MAGED subfamily are overexpressed, while most others including MAGEE1, MAGEE2, MAGEL2, MAGEH1, NDN, and TRO are downregulated in glioma. The downregulation of these T2Ms suggests their potential tumor suppressor functions in glioma. Downregulation and tumor suppressor functions of Necdin in glioma have been previously reported in detail (19). MAGEH1 expression has been reported to be downregulated in hepatocellular carcinoma (39) and cholangiocarcinoma (40). Further, MAGEH1 overexpression inhibits proliferation, cell migration, and invasion in hepatocellular carcinoma (39), induces apoptosis in melanoma cell lines (41), and cell cycle arrest in multiple cancer cell lines (40,42). In agreement with these studies, we also observed that MAGEH1 is significantly downregulated with advanced grade of glioma. Further, higher MAGEH1 was also associated with age, IDH mutation, and 1p/19q codeletion status. Similar to T1Ms, which are coordinately expressed (23), our correlation analysis revealed that T2Ms are also co-expressed in glioma tissues. Further, we observed highly coordinated expression within two major groups of the MAGE gene family in multiple cancers, including glioma. Interestingly, while MAGED subfamily members showed upregulation in glioma, they also exhibited strong positive correlation among themselves, therefore suggesting common regulation and functions. Theoretically, co-expression of structurally similar MAGE proteins can regulate a large number of oncogenic functions. Nevertheless, functional interaction of two MAGE proteins to perform a single oncogenic function has also been reported (43). Furthermore, both MAGEH1 and Necdin have been shown to interact with p75 neurotrophin receptors, a class of receptors that regulate cell growth and neuronal survival (44).
We hypothesized that the coexpression pattern observed for the downregulated T2Ms is guided by a common regulatory mechanism, including DNA methylation. DNA hypomethylation associated overexpression of MAGED4 has previously been reported in glioma (18). Therefore, it was to our interest whether the downregulated T2Ms that emerged in the expression analysis are also regulated by DNA methylation. Interestingly, we identified several intragenic CpG sites in MAGEH1 gene, where DNA methylation was closely associated with its gene expression. Although the IDH mutation in glioma leads to a hypermethylation phenotype and global gene repression, interestingly, MAGEH1 exhibited higher expression in IDH mutant tumors. This is in agreement with previous reports which demonstrated that DNA hypermethylation may also upregulate some genes by promoting transcription accessibility at some genomic regions (45).
The prognostic utility of T2Ms was revealed by survival analysis in TCGA datasets. Interestingly, the association of MAGEH1 with better survival of patients was separately validated in major histological and molecular subtypes of glioma. Further, this association was independent of IDH status and other clinically relevant variables. The association of MAGED4 overexpression with poor prognosis of glioma patients has been previously reported (46). In agreement to this, our analysis revealed that while MAGED1, which is co-expressed with MAGED4, was associated with poor prognosis in glioma, it did not emerge as an independent prognostic factor in multivariate analysis. Using a similar analysis of cancer-testis antigens in TCGA datasets, association of Necdin expression with favorable overall survival in glioma has also been reported (20). Our study revealed that NDN expression is downregulated in a coordinated manner with other T2Ms and higher expression of these T2Ms is associated with better patient prognosis in lower grade glioma.
While the current pieces of evidence strongly suggest tumor suppressor functions of some T2Ms, limited information is available for detailed functions of these proteins in glioma. We, therefore, performed pathway analysis for MAGEH1 co-expressed genes, which revealed strong positive enrichment of Myc and Hedgehog signaling genes. While both of these signaling pathways promote glioma progression (47,48), their functional association with MAGEH1 remains unclear. Nevertheless, a higher c-Myc expression is indeed associated with better prognosis in glioma patients (49). This might be partially due to less characterized proapoptotic properties of c-Myc, which may sensitize Myc overexpressing glioma cells to DNA damaging agents, such as chemotherapy and radiotherapy (49).
Recent evidence suggests that immune response plays a crucial role in glioma and immunity associated pathways are associated with therapeutic response and clinical outcome. Further, analysis of immunogenomic profiles has led to the identification of several immunomodulatory pathways, such as DNA damage repair (50) and tumor-associated immunomodulatory proteins, such as Galectin-1 and IGFBP2 (51,52). Pathway analysis revealed that the expression of MAGEH1 was negatively correlated with several immunity associated pathways along with level of immune infiltration in LGG. Contrary to this, a positive correlation of MAGEH1 with different immune cell population was observed in GBM. These results suggest that T2Ms might be involved in regulating tumor immune infiltration distinctly in LGG compared to GBM. LGG have been shown to exhibit distinct immune composition compared to GBMs (53,54). In addition to expected intratumor heterogeneity, the microenvironment composition variably influences gene expression profiles (55,56).
In conclusion, the current study provides intrinsic details of coordinated expression, epigenetic regulation, and prognostic significance of T2Ms in glioma. These findings will provide a strong basis and resource for designing future research on Type 2 MAGE proteins. However, this study is limited by the utilization of gene expression as the primary measure for alterations of T2Ms. Therefore, a protein level study with mechanistic exploration may provide better insights into the role of T2M proteins in glioma biology.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
MA conceptualized the study. SC supervised the study and provided essential infrastructure. MA, SK, and JS and AC performed data curation, interpretation, and statistical analysis. MA and SK wrote the original manuscript.
All authors contributed to the article and approved the submitted version.