System biology approach to delineate expressional difference in the blood mononuclear cells between healthy and Turner syndrome individuals

Abstract

Background

Turner syndrome (TS) is a rare disorder associated either with complete or partial loss of one X chromosome in women. The information on the genotype–phenotype relationship in TS is inadequate. Comparing the healthy and Turner syndrome patients may help elucidate the mechanisms involved in TS pathophysiology. Gene expression differences between healthy and individuals with Turner syndrome were characterized using the systems-biology approach of weighted gene coexpression network analysis (WGCNA) on 182 microarray peripheral mononuclear blood samples (PBMC).

Results

The coexpression networks of healthy and TS had scale-free topology that ensures network robustness. In the process, five modules were preserved between healthy and TS, which carry several genes common in each module. Two of them, SMCHD1 and PGK1, have already been reported to be involved in TS. Previously reported genes of TS, specifically, PTPN22, RPS4X, CSF2RA, and TIMP1, were missing in their respective modules. Dysfunction, differential expression, or absence of these genes could lead to a progressive disruption of molecular pathways leading to the pathophysiology of TS. Indeed, we observed a significant difference in the functions of these modules when compared within and across the healthy and TS samples. We identified four clusters in the PPI network constructed from the top 15 KME enriched in significant functions.

Conclusion

Overall, our work highlights the potential molecular functions, pathways, and molecular targets of TS that can be exploited therapeutically in the human healthcare system.

Background

This era of genomic technologies provides an immense wealth of molecular information offering a remarkable opportunity for new target discovery [1]. The microarray gene expression data provides an understanding of the underlying system-level view of the disconcerted networks of the disease cells. This, in turn, allows recognition of the potentially responsible key genes or pathways that can be targeted therapeutically [2]. However, the identification of key molecular targets remains a challenge for syndromes whose pathophysiology is mainly due to global genomic imbalance rather than defect in a single gene. Turner syndrome (TS) is a rare disorder associated either with complete or partial loss of one X chromosome in women. Often the women with TS have mosaic karyotypes [3]. Comorbidities associated with TS include short stature, absent or incomplete central puberty, ovarian dysgenesis, infertility, congenital malformations of the heart, endocrine disorders such as type 1 and type 2 diabetes mellitus, osteoporosis, recurrent miscarriages, and autoimmune disorders [4]. A multidisciplinary approach to care is required because Turner syndrome can affect multiple organs at any stage of life [5].

Previously, we demonstrated that the TS networks, constructed from manually curated genes from literature searches and microarray gene expression data, carry scale-free hierarchical fractal attributes [6, 7]. This organization showed that the TS networks were dynamic and heterogeneous, suggesting this to be a compact and self-organized and had efficient information processing potential. We identified few key regulators and essential genes (interologs) from these TS networks that we envisage having potential as therapeutic targets.

Weighted gene coexpression network analysis (WGCNA) is another powerful approach that lets us identify the higher-order relationships among the genes built on their coexpression relationships, demarcating highly similar coexpressed genes into modules, and allowing a robust view of transcriptome organization [8]. Exploring TS through WGCNA may explain a whole new aspect of this syndrome, and this may shed light on the physiology and gene network perturbations that occur in TS.

Our present work highlights the potential molecular functions, pathways, and molecular targets when genomic data from TS and healthy women are compared with each other and analyzed at the level of gene coexpression modules. Integrating this type of data, in which a comparative analysis of modules of healthy and TS samples can be studied in response to relevant molecular perturbations, may further facilitate the identification and validation of novel molecular targets.

Methods

Data set acquisition

The raw gene expression data were downloaded from the GEO data repository (9) (http://www.ncbi.nlm. nih.gov/geo/). As we aimed to compile a wide set of comparable data, we collected as many relevant data sets as we could find. We then filtered out all barring a core collection of data sets that were similar enough to be used for bioinformatics comparison. We considered four inclusion criteria. First, we considered only those data sets that were run on the Affymetrix platform HG-U133 Plus 2. Second, we included only those samples in each data set that were healthy (control) from human peripheral blood mononuclear Cells (PBMC). Third, only female samples were considered. And fourth, we made the correlations between genes more comparable across the studies by omitting all the data sets with fewer than five samples. Qualifying these inclusion criteria, a total of 12 data sets were considered for further analysis (Table 1). These expression data were combined, and prior to the analysis, batch corrected. The genes detected by multiple probes were represented by their mean value. These expression data were then subjected to WGCNA analysis.

Table 1 List of microarray series used in the study

Data preprocessing

The datasets that we collected individually had a small sample size. Performing the statistical analysis individually on each series may not give robust results [19]. Thus, we performed a combined analysis of the collected gene expression data from different series (Table 1). However, combining these datasets can lead to batch effects. So, we removed the batch effects between these datasets through the Empirical Bayes method (ComBat) [20] using the GeneSpring Software 14.9.1 (Agilent Technologies, Inc., Santa Clara, CA, USA). This algorithm adjusts the batch effect in microarray expression data. We assumed that the measured expression value of gene i in sample j of batch g can be expressed in a general form as follows:

$${{\varvec{x}}}_{{\varvec{i}}{\varvec{j}}}= {{\varvec{x}}}_{{\varvec{i}}{\varvec{j}}}^{\mathrm{^{\prime}}}+{{\varvec{b}}}_{{\varvec{i}}{\varvec{j}}}^{{\varvec{x}}}+{{\varvec{\varepsilon}}}_{{\varvec{i}}{\varvec{j}}}^{{\varvec{x}}}$$

(1)

Then, to remove the batch effect, we standardize it:

$$\hat{\user2{x}}_{ij} = \frac{{{\varvec{x}}_{{{\varvec{ij}}}} - \overline{\user2{x}}_{{{\varvec{ij}}}} }}{{{\varvec{\sigma}}_{{{\varvec{x}}_{{\varvec{i}}} }} }}$$

(2)

Here, x_i is the average value of the expression value from gene i in all the samples; σ x_i is the standard deviation for x_i of the expression value from gene i in all the samples.

Coexpression network construction: a system biology approach

Analyzing the coexpression patterns of the genes can provide an understanding of the underlying cellular processes because the coordinated coexpression of the genes encodes interacting proteins [21]. Therefore, we used the collected data sets (Table 1) to construct the coexpression network following the protocols of WGCNA [8, 22] to create a coexpression network of healthy humans (females) from PBMC samples. Separately, we again constructed a coexpression network of Turner syndrome patients from PBMC samples (GSE46687). We did this to get an overview of the comparison between healthy and TS patients from the perspective of a coexpression network.

To construct the gene coexpression network and identify its modules, we used the WGCNA (weighted gene coexpression network analysis) implemented in the R software package (http://www.r-project.org/). We used the Pearson correlation coefficient as a coexpression measure for cluster analysis. In these networks, a node represents the gene expression profile of a given gene. Thus, nodes are connected only if they have a significant pairwise expression profile cutoff across the samples. In our study, for all the pairwise comparisons of gene expression, the absolute value of the Pearson correlation coefficient was calculated across all microarray samples, and the appropriate soft threshold power was chosen. Then, the Pearson correlation matrix was transformed into an adjacency matrix. An adjacency matrix is a matrix of connection strengths through a power function. Thus, the connection strength (adjacency) a_ij between gene expressions x_i and x_j is defined as:

$${{\varvec{a}}}_{{\varvec{i}}{\varvec{j}}}={|{\varvec{c}}{\varvec{o}}{\varvec{r}}({{\varvec{x}}}_{{\varvec{i}}},{{\varvec{x}}}_{{\varvec{j}}})|}^{{\varvec{\beta}}}$$

(3)

As known, in general, all types of biological networks (gene expression networks) have been found to exhibit an approximate scale-free topology (2–6). Scale-free networks follow the power law and decay as P(k) ~ k^−g, where g is some exponent. Depending on the graph, we used the scale-free topology criterion (β). It is expected that the weighted gene coexpression network has an advantage over the unweighted coexpression network as the former one is more robust and the continuous nature of the gene coexpression information is preserved in it. For weighted gene coexpression network, whole network connectivity k_i of the ith gene expression profile x_i is the sum of the connection strengths with all other genes in the network, i.e.,

$${{\varvec{k}}}_{{\varvec{i}}}={\sum }_{\mathbf{j}\in \mathbf{N},\mathbf{j}\ne \mathbf{i}}{|{\varvec{c}}{\varvec{o}}{\varvec{r}}({{\varvec{x}}}_{{\varvec{i}}},{{\varvec{x}}}_{{\varvec{j}}})|}^{{\varvec{\upbeta}}}$$

(4)

where N refers to the set of network genes.

Then, modules were identified following the dynamic tree cut method. The modules are the groups of genes with high topological overlap (1, 7). Through this method, the spurious or isolated connections were excluded during module formation.

The topological overlap for weighted networks (7) is given by,

$${{\varvec{t}}}_{{\varvec{i}}{\varvec{j}}}=\frac{\sum_{{\varvec{u}}}{{\varvec{a}}}_{{\varvec{i}}{\varvec{u}}}{{\varvec{a}}}_{{\varvec{u}}{\varvec{j}}}+{{\varvec{a}}}_{{\varvec{i}}{\varvec{j}}}}{\mathbf{min}\left({\sum_{{\varvec{u}}}{\varvec{a}}}_{{\varvec{i}}{\varvec{u}}},\sum_{{\varvec{u}}}{{\varvec{a}}}_{{\varvec{j}}{\varvec{u}}}\right)+1-{{\varvec{a}}}_{{\varvec{i}}{\varvec{j}}}}$$

(5)

where an adjacency matrix is represented by,

$${\varvec{A}}=\left[{{\varvec{a}}}_{{\varvec{i}}{\varvec{j}}}\right], 0\le {{\varvec{a}}}_{{\varvec{i}}{\varvec{j}}}\le 1$$

(6)

We used the topological overlap dissimilarity measure (1—topological overlap matrix (TOM)) in an average linkage hierarchical clustering to define the modules as branches of the resulting dendrogram. The genes of the modules were then assigned the eigenvalues based on which the significant module eigengenes (ME) were identified in each module.

Functional and pathway enrichment

The functional enrichment analysis was conducted on the genes in the identified modules. GO biological process term and Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses [23] were conducted using DAVID online server [24]. Functional enrichment analysis was based on the cutoff value of Benjamini P-value < 0.05.

PPI network of TS and MODULE identification

To further explore the significant genes of each module in terms of protein–protein interactions, the top 15 K_ME genes in each module were extracted. These genes were then entered into STRING database [25] to form the PPI network. The Cytoscape plug-in molecular complex detection (MCODE) [26, 27] was used to analyze the most notable clustering module. Parameters for MCODE were set with Degree Cutoff as 2, Node Score Cutoff as 0.2, K-Core as 2, and Max. Depth from Seed as 100. Thus, system biology approach may be exploited in human health care system including complex polygenic diseases.

Results

The workflow of the work presented here is illustrated in Fig. 1.

Fig. 1

Schematic illustration of the workflow

Data set acquisition and pre-processing

By filtering the datasets with our inclusion criteria (mentioned in the methodology), we collected a total of 12 microarray series which gave us 182 healthy PBMC female samples (see Table 1).

Module formation and characterization

Coexpression networks were generated using the weighted gene coexpression network (WGCNA) package, in which all the coexpressed genes were connected with varying correlation strengths. This was accomplished using soft thresholding, thereby preserving the continuous nature of the data set and eliminating the need to set an arbitrary correlation score cutoff. The soft threshold power β is a value used to power the correlation of the genes to that threshold. It is assumed that raising the correlation to a power (β) will reduce the noise of the correlations in the adjacency matrix. To pick up one threshold, we used the pickSoftThreshold function, which calculates if the network resembles a scale-free graph, for each power. The lowest power that resulted in a scale-free topology is the one we used. In a scale-free network, small-degree nodes are the most abundant, but the frequency of high-degree nodes decreases relatively slowly. Thus, nodes that have degrees much higher than average, so-called hubs, exist. Because of the heterogeneity of scale-free networks, random node disruptions do not lead to a major loss of connectivity, but the loss of the hubs causes the breakdown of the network into isolated clusters [28]. The biological networks are considered scale-free networks. This is the reason, in case of WGCNA, we selected that β value, which gives scale-free topology. The genes that had highly similar expression levels, clustered together [29]. Such highly correlated genes may shed light on the shared biological processes or shared regulatory mechanisms that could be targeted. Therefore, we performed combined global gene expression profiling on 182 patient samples from PBMCs Set 1, n = 182 (combined 12 datasets in Table 1) and Set 2, n = 26 (only Turner syndrome samples in GSE46687)) and separately constructed a weighted gene coexpression network of Set 1 and Set 2 based on pairwise Pearson correlations between the expression profiles for the identification of gene modules. A gene module is a group of highly coexpressed genes. These modules were detected through unsupervised hierarchical clustering [30]. In a coexpression network, for fulfilling the property of scale-free topology, the parameter β is considered very crucial. Most likely, the biological networks which are based on gene expression data are scale-free [31]. Therefore, β = 6 (for Set 1) and β = 4 (for Set 2) were considered to obtain scale-free topology by the fit index greater than 0.87 and 0.84, respectively. Figures 2 and 3 show the result of several powers for finding a network with scale-free topology properties of Sets 1 and 2, respectively. We found that the coexpression networks had a scale-free topology (Figs. 2C and 3C), which is believed to ensure network robustness and thus resisting random node attacks [32].

Fig. 2

Determination of soft-thresholding powers (β) in the weighted gene coexpression network analysis (WGCNA) for Set 1. A Analysis of the scale-free fit index for various soft-thresholding powers (β). B Analysis of the mean connectivity for various soft-thresholding powers (β). C Log–log plot of whole network connectivity distribution. The x-axis shows the logarithm of whole network connectivity, y-axis the logarithm of the corresponding frequency distribution. On this plot, the distribution approximately follows a straight line, which is referred to as approximately scale-free topology when β = 6

Fig. 3

Determination of soft-thresholding power (β) in the weighted gene coexpression network analysis (WGCNA) for Set 2. A Analysis of the scale-free fit index for various soft-thresholding powers (β). B Analysis of the mean connectivity for various soft-thresholding powers (β). C Log–log plot of whole network connectivity distribution. The x-axis shows the logarithm of whole network connectivity, and y-axis the logarithm of the corresponding frequency distribution. On this plot, the distribution approximately follows a straight line, which is referred to as approximately scale-free topology when β = 4

After identifying modules of coexpressed genes, each module in effect becomes a new network, and a new measure of connectivity (intramodular connectivity, or kin) is defined as the sum of a gene’s connection strengths with all other genes in its module. Intramodular connectivity measures how connected, or coexpressed, a given gene is with respect to the genes of a particular module. The module network dendrogram, constructed by clustering ME distances, showed that modules with high Kin were positioned at the tip of the branches since they exhibit the highest interconnectedness with the rest of the module. After highly similar modules were merged, eleven gene expression modules were determined in Set 1 (Fig. 4A), and five gene expression modules were determined in Set 2 (Fig. 4B). The colors are assigned based on module size. Turquoise (others refer to it as cyan) colors the largest module, next comes blue, next brown, next yellow, then green, and so on. A gray color module is reserved for unassigned genes. Eigengenes, the first principal component of a cluster, is thought to be a representative of a cluster’s expression profile [22]. To get a sense of how related the modules are, one can summarize each module by its first eigengene (referred to as principal components), and then correlate these module eigengenes with each other based on K_ME value (Fig. 4C and D). The module eigengene-based intramodular connectivity measure kME roughly approximates the standard intramodular connectivity kIN. This measure is determined by correlating the expression profile of a gene i with the module eigengene of its resident module: kMEi =|cor(x i, ME)|. This gives top 15 most connected genes of each module.

Fig. 4

WGCNA network and module detection A and B (Upper) Cluster dendrograms in the human Sets 1 (A) and 2 (B). The y-axes correspond to distance (1 − TO). (Lower) Dynamic tree cutting was used to determine modules, generally by dividing the dendrogram at significant branch points. Modules with significant overlap were assigned the same labels. Colors in the horizontal bar represent the modules. C and D. The module network dendrogram was constructed by clustering module eigengene distances in Sets 1 C and 2 D

The multiple dimensional scaling (MDS) plot provides an alternate visualization of the module structure. Figure 5 represents multidimensional scaling plots (MDS) of coexpressed modules in Set 1 and Set 2. Modules tend to form separate ‘fingers’ in this plot. The color denotes genes of that module. Intramodular hub genes are located at the finger tips. The MDS plots are shown in Fig. 5. All the plots are color-coded according to the modules concerning respective Sets 1 or 2. The points with the same color are almost always clustered together which means that the relative position of the points is well-preserved. However, the spatial distributions of the points vary to a large extent when we compare MDS plots of healthy vs TS. This is because the expression pattern of genes is disrupted due to TS.

Fig. 5

A multidimensional scaling plots (MDS) of coexpressed modules in Set 1 and Set 2. Modules tend to form separate ‘fingers’ in this plot. The color denotes genes of that module. Intramodular hub genes are located at the finger tips. A An MDS plot displaying expression data of genes in different modules in Set 1. B An MDS plot displaying expression data of genes in different modules in Set 2

The detailed list of genes of each module in each set is listed in electronic supplementary material, (Additional file 1: Table S1. and S2).

Comparison of healthy and TS patients modules provide insights into TS

A comparison of coexpression networks between healthy and TS human samples could provide valuable insight into the syndrome. We found that five modules were common between Set 1 (Healthy) and Set 2 (Turner Syndrome). The genes that were common between these modules are listed in Table 2. It was identified that the genes highlighted in bold were differentially expressed in TS. Two of them, namely SMCHD1 [33] and PGK1 [34], have already been reported to be involved in TS. Thus, the functions performed by these genes (highlighted in bold) in their respective modules in healthy individuals are perturbed in TS patients due to their differential expression. This may contribute to the pathophysiology of TS. Therefore, as studied previously, it can also be concluded that the TS phenotype is caused due to the additive effect of genes from different loci and due to global genomic imbalance [7].

Table 2 Common genes between the common modules of Set 1 (Healthy) and Set 2 (TS)

Earlier, we reported seed genes, few key regulatory genes, and signature genes involved in Turner syndrome from manually curated genes of TS from the literature [6] and microarray analysis [7] employing the network approach. Table 3 contains the list of those seed genes that are present in the gene expression modules of the healthy and TS coexpression network. We found that many of these genes were present in the different modules of Set 1 and Set 2. Since these reported genes are expected to play an important role in the pathophysiology of TS, we can see that many important genes (present in the modules of healthy samples coexpression network) are missing in their respective modules of TS samples (Table 3). This shows that due to the TS state, the important genes start to function in a different way resulting in disturbed gene functioning. The genes that coexpressed to share some biological process or regulatory mechanisms in healthy females changed their pattern of functioning, thus, causing the perturbation in the normal functioning of an individual with TS.

Table 3 List of previously reported seed genes present in the gene expression modules of the healthy and TS coexpression network

Specifically, four of the genes, namely PTPN22 [35], RPS4X [36, 37], CSF2RA [38, 39], and TIMP1 [40, 41] that we identified in our previous studies [6] as the key regulatory genes of the TS network were found to be missing in their respective modules of TS. Thus, dysfunction, differential expression, or absence of these genes could then lead to a progressive disruption of the molecular pathways in a predictable manner leading to the pathophysiology of TS. This also suggests a high level of genetical, organizational and expressional heterogeneity among the TS patients.

Comparison of functional and pathway enrichment of different modules of healthy vs Turner Syndrome

To compare the potential biological functions of critical modules between healthy and TS samples, we performed the functional enrichment analysis of the GO and KEGG categories. There was a significant difference in the functions that different modules were enriched in when compared within and across the healthy and TS samples (Table 4). However, few biological functions were consistent within and across the healthy and TS samples.

Table 4 Comparative analysis of Functional Enrichment of coexpression modules of Healthy with Turner Syndrome

We identified five modules common between healthy and TS coexpression networks. In the turquoise module of both the healthy and TS coexpression network, the common molecular functions identified were protein binding, RNA binding, and poly(A) RNA binding. The common enriched cellular components were found to be the cytosol, extracellular exosome, and cytoplasm and none of the biological functions were found to be common in turquoise module. In the blue module of both, the healthy and TS coexpression network, the common molecular function identified was protein binding, while none of the biological function and cellular component was found to be common. In the green module, extracellular exosome while in the yellow module, cytosol was found to be the common enriched cellular component of both the healthy and TS coexpression network. None of the functions were common in the brown module of both the healthy and TS coexpression network. None of the KEGG pathways was found to be common in the common modules of healthy and TS coexpression networks.

Though there were few enriched functions and cellular components that were common, there was a significant difference in the enriched functions and pathways among these five common modules. The functional enrichment and KEGG pathway of different modules of healthy and TS coexpression networks are listed in Tables 4 and 5.

Table 5 Top 15 K_ME in each module of TS samples

It can be explained that the similar enriched functions though similar but the set of genes that perform these functions are slightly different (in each module) in TS samples as compared to healthy ones. This means that the functions being similar are somehow being performed by a different set of genes, and due to differential gene expression in the TS state, these functions may not be performed in synchronized and a proper way resulting in the pathophysiology of Turner syndrome. And obviously, there are many enriched molecular functions, biological functions, cellular components, and KEGG pathways that are significantly different among the five common modules between healthy and TS coexpression networks.

While many genes (Table 2) are common between common modules of the healthy and TS patients, most of the genes are different in these modules between healthy and TS. This clearly makes sense that in a healthy state, these common genes function in synchronization and combination of the genes listed in electronic supplementary material, (Additional file 1: Table S1). These genes get replaced by a new set of genes in the TS state (electronic supplementary material, (Additional file 1: Table S2)), and the complete functional and biological process is performed with a whole new set of genes. Many of these genes are differentially expressed in TS. Due to this differential expression, the genes perform the biological and molecular functions in a perturbed way. This explains the mechanism of Turner syndrome and pathophysiology.

PPI network of TS and identification of its clusters

The eigengenes are the principal component of a module, and they can be a representative of a cluster’s expression profile. To understand how these modules interact with one another, we selected the top 15 K_ME genes of each module. These K_ME genes of each module are listed in Table 5. Being the eigengenes of each module, the genes listed in Table 5 may have a potential role in TS. These genes were then entered into the STRING database to construct the protein–protein interaction (PPI) network. The confidence score was set at 0.7, and the genes in the first shell was set at 100. The main network constructed had a clustering coefficient of 0.750 with 129 nodes and 990 edges. In total, we considered the top four sub-networks (modules 1, 2, 3, and 4) with scores 34, 12.5, 7, and 7, respectively, that were detected by MCODE (Fig. 6). Each cluster had the highest scoring node called a seed. It is the node from which the cluster was derived. The genes NDUFV1, EFTUD2, POLR1C, and STT3A are the highest scoring node, i.e., seed genes of clusters 1, 2, 3, and 4, respectively. They act as the eigengene of those cluster. The list of genes of each cluster is listed in electronic supplementary S3. The functional and pathway enrichment of modules of the PPI network of TS was analyzed, and it was found that cluster 1 is mainly enriched in the processes and reactions related to mitochondrial phosphorylation, cluster 2 is mainly enriched in golgi-mediated protein transport and reactions related to it, cluster 3 is mainly enriched in processes related to transcription, and cluster 4 is mainly enriched in post-translational modification processes like glycosylation. The list of different molecular functions, biological functions, cellular components, and KEGG pathways in which these four clusters are enriched is listed in the electronic supplementary material (Additional file 1: Table S4).

Fig. 6

A Protein–protein interaction network of top 15 K_ME genes of each module of TS coexpression network. B-E Modules identified by MCODE where B is cluster 1, C is cluster 2, D is cluster 3, and E is cluster 4. The highest scoring node in the cluster is called a seed. It is the node from which the cluster was derived and is represented by a square shape (see square yellow in B, C, D and E)

Discussion

Turner syndrome (TS) is a rare disorder that is associated with either complete or partial loss of one X chromosome in women, with an approximate occurrence of 1 in 2500 females [42]. Surprisingly, our knowledge of genotype–phenotype relationships in TS is rather limited, with very few specific candidate genes linked to its clinical features. The genotype–phenotype relationship in TS refers to the relationship between karyotype and phenotypic features of TS. TS is defined by the presence of an abnormal X chromosome in the form of monosomy, mosaicism of a 45X cell line with another cell line, which could be 46XX, 46XY, or an abnormal sex chromosome rearrangement [43, 44]. In TS, there is a link between the exact cytogenetic appearance and the phenotype. Pure 45, X monosomy is the most common karyotype and is associated with the most abnormal phenotype. Although the clinical features of TS are well-defined, the severity of the phenotype in TS individuals varies according to the underlying chromosomal constitution [44]. To better characterize gene expression differences between healthy and TS, we took a systems-biology approach by using weighted gene coexpression network analysis on 182 microarray peripheral mononuclear blood samples (PBMC). By comparing the modules in both healthy and TS, we aim at finding more relations between both healthy and TS PBMC samples including both the similarities and the differences.

The coexpression networks of healthy and TS samples had a scale-free topology that ensured network robustness which means they can resist random node attacks. Five modules were preserved between healthy and Turner Syndrome that carry many genes that were common (listed in Table 2) in each module. This means that these five modules are TS specific and are related to TS progression, which implicates many genes that were previously known and many that were unknown that function in TS. We found that many genes were common between the modules of TS and healthy coexpressed networks. Many of these common genes were differentially expressed in TS modules. This signifies that most of these genes being differentially expressed fail to perform the molecular functions of their respective modules. Two of these genes, namely SMCHD1 and PGK1 that were common between turquoise and brown modules, respectively, of healthy and TS have already been reported to be involved in TS. SMCHD1 (OMIM 614982) plays a role in X chromosome inactivation [45]. It was reported earlier in a 13-year-old girl with Turner syndrome-like clinical features in association with intellectual disability, facial dysmorphism, and the psychomotor developmental delay had chromosome deletion syndromes and also had SMCHD1 polymorphism [33]. This SMCHD1 may correlate with TS. In another case study of a 28-year-old Turner syndrome patient, it was found through the DNA studies that the locus controlling X-inactivation is proximal to PGK1 suggesting the role of PGK1, as well, in TS [34]. There were few more common genes between healthy TS modules, namely CDC42SE2, ARPC4, RUNX1-IT1, ZBTB20, USP34, EMILIN2, GLIPR2 that were differentially expressed in TS and may prove to be significant in the progression of TS pathophysiology.

Many genes were different among the common healthy and TS modules. Few more previously reported genes of TS, specifically, PTPN22, RPS4X, CSF2RA, and TIMP1 were found to be missing in their respective modules in TS when compared with healthy ones. Dysfunction or absence of these genes (listed in Tables 1 and 2) could then lead to a progressive disruption of the molecular pathways in a predictable sequence leading to the pathophysiology of TS.

It is unsurprising that the presence of 45,X monosomy results in dysregulation of cellular growth and repair pathways, leading to significant detrimental effects during both embryonic development and later stages. The extended cell cycle observed may mark the initiation of Turner syndrome (TS) pathogenesis, giving rise to a cascade of phenotypic consequences throughout various life stages, including embryonic/fetal, neonatal, pediatric, adolescence, and adulthood phases. The X chromosome harbors transcription factors influencing cell cycle duration and haploinsufficiency may contribute to cell cycle delays. Moreover, the insufficient expression of genes associated with DNA replication or repair on the X chromosome could further prolong the cell cycle in 45,X cells, impacting their ability to respond to S replication signals. TS patients often experience embryonic lethality, growth retardation, short stature, osteopenia/osteoporosis, congenital heart disease, gonadal dysgenesis, impaired pancreatic β-cell function, and neurologic deficits, manifesting as compromised neuronal microstructural integrity (connectivity) in white matter pathways.

While there were few molecular, biological functions and cellular pathways common among these preserved modules of healthy and TS, we also observed a significant difference in the functions of these modules when compared within and across healthy and TS samples. These genes mentioned in Tables 1 and 2 are likely to be involved in molecular and biological functions and pathways that are syndrome-related, thus providing a new window of an investigation into TS pathophysiology. The identified pathways in individuals with Turner syndrome (TS) offer valuable insights into the underlying molecular processes associated with this condition. The dysregulation observed in these pathways has specific implications for various aspects of TS pathology. The spliceosome pathway, crucial for RNA splicing, may influence gene expression patterns, contributing to the diverse clinical manifestations seen in TS. Transcriptional misregulation, involving errors in gene transcription control, further accentuates the complexity of TS molecular dysregulation and its potential impact on aberrant gene activity. Platelet activation and vascular smooth muscle contraction dysregulations hint at possible connections to cardiovascular complications observed in TS. Abnormalities in osteoclast differentiation suggest potential implications for bone health, possibly contributing to conditions like osteoporosis. Lysosome and phagosome dysregulation may influence cellular homeostasis and immune responses in TS. Oxidative phosphorylation, central to energy production, highlights potential metabolic perturbations in TS. The presence of the non-alcoholic fatty liver disease (NAFLD) pathway points to potential liver-related complications in individuals with TS. Collectively, the identification of these pathways not only provides a more elaborated understanding of TS at the molecular level, but also lays the groundwork for future research and therapeutic strategies targeting the specific biological dysregulations associated with this genetic condition.

Further, we identified four clusters in the PPI network constructed from the top 15 K_ME in each module of the TS coexpression network. These modules were enriched in different significant molecular and biological functions, cellular components, and KEGG pathways. Each cluster has the highest scoring node called a seed. It is the node from which the cluster was derived. The genes NDUFV1, EFTUD2, POLR1C, and STT3A are the highest scoring node representing the seed genes of cluster 1, 2, 3, and 4, respectively, as mentioned earlier. They act as the eigengene of that cluster.

This comprehensive network-based meta-analysis methodology has previously been implemented by Jeremy et al. on human and mouse brain transcriptome [19]. This methodology of WGCNA has many advantages over traditional transcriptional analyses, leading to more reliable results compared to previous studies. Although, through this approach we uncover many facets of TS, there are several limitations such as the limited size of TS samples. The currently available datasets of PBMCs samples of TS do not allow a more elaborated study at this moment. Thus, a larger sample size would provide more elaborate results. Overall, our work highlights the potential molecular functions, pathways, and molecular targets when genomic data of TS and healthy women are compared and analyzed at the level of gene coexpression modules that can be targeted therapeutically in the future.

This work also provides an opportunity to ascertain if these genes originating from different modules have lost their functional attributes owing to mutational load or aberrant signal transduction or defunct interactomes. Since the number of interactomes are expected to be more for all the genes, their own detailed characterization would prove to be a rewarding proposition. From these interactomes, regulatory genes may be identified using network-based approach. In any case, these observations are likely to fuel the thought along the line for precise identification of genes actually implicated in giving rise to TS. Another important aspect is that after the implantation of the blastocyst, at what stage TS sets in and how the molecular heterogeneity is accentuated and prevailed. If we know what compels such an events, we may use the information for better management of TS patients including that of prenatal diagnosis.

Conclusion

In this study, gene expression differences between healthy and individuals with Turner syndrome were characterized using the systems-biology approach of weighted gene coexpression network analysis (WGCNA) on 182 microarray peripheral mononuclear blood samples (PBMC). While many genes were found to be common between the modules of TS and healthy coexpressed networks, some of the genes were different among the common healthy and TS modules. These genes may prove to be significant in the progression of TS pathophysiology. Dysfunction or absence of these genes could lead to a progressive disruption of the molecular pathways in a predictable sequence leading to the pathophysiology of TS.

Overall, our scientific work holds potential benefits for various groups within society. Primarily, biotechnologists, geneticists, and doctors stand to gain valuable insights from our analysis, as it may contribute to early diagnosis and proactive intervention in individuals at risk of Turner Syndrome. By identifying specific genes associated with the syndrome, our research could aid geneticists in establishing reference panels for more accurate and efficient genetic testing. This has the potential to facilitate early detection and intervention before symptoms manifest or worsen.

Moreover, our findings could pave the way for advancements in personalized medicine. Pharmaceutical companies might leverage the identified genes to develop Turner Syndrome-specific drugs aimed at reducing the severity of symptoms. This personalized approach has the potential to significantly improve the quality of life for individuals with Turner syndrome by addressing the underlying genetic factors. In essence, our work carries implications for both diagnostic practices and the development of targeted therapeutic strategies, offering a comprehensive contribution to the fields of genetics, medicine, and biotechnology.