Transcriptomic analysis of elderly women with low muscle mass: association with immune system pathway

Sample characteristics

Using the definition of Newman et al. [1] and the residual previously defined for our population, we included forty elderly women in the study, of which 20 had LMM (residual <−1.32) and 20 age and race-matched controls without LMM. Lower weight, total and appendicular lean mass and lower BMI were observed in the LMM group (Table 1). No difference in abdominal circumference was found between the groups. A higher percentage of fat was observed in the group with LMM (median 41.63 (3.94) vs. 38.61 (4.18); p = 0.02 for LMM and controls respectively). Regarding associated chronic diseases and medications for chronic use, no differences were found between groups. Comparing the biochemistry between the groups, a significantly lower serum level of vitamin D (25OHD) (18.75 (7.06) vs. 24.15 (6.74) ng/ml; p = 0.018) and serum insulin (11.6 (5.7) vs. 17 (14.57) UI/mL; p = 0.02) was observed in the LMM group. The participants included in this study are homogeneous in terms of lifestyle and habits. They reside in the same geographic region, and have the same socioeconomic, educational level and access to the same health system. As part of the SPAH cohort (São Paulo Ageing and Health Study) they were followed up in the city of São Paulo, and are routinely guided by medical doctors, medical students, and other health professionals from the University of São Paulo (USP), School of Medicine. Serum albumin levels were similar in both groups (4.65 ± 0.2 vs. 4.65 ± 0.22 g/dL; p = 0.30 for LMM and control respectively) and were in the normal range (3.5–5.2 g/dL) suggesting a reasonable nutritional level. Regarding consumption of milk and derivative products, no difference was observed between the groups with LMM and control (median 527.1 vs. 553.2 mg/day; p = 0.75 for LMM and control respectively). All other characteristics are presented in Table 1.

Discriminating transcript clusters and differential expression in LMM

In this analysis, for simplicity purposes, transcript clusters are sometimes identified by their corresponding genes symbols. The logistic regression model adjustment identified 821 unique transcript clusters of which 565 with positive coefficients (PPR) and 255 with negative coefficients (PNR) (Supplementary Table 1 and Figure 1A). Differential expression analysis identified 56 and 79 up- and down-regulated transcript clusters respectively (Supplementary Table 2). Among these transcript clusters, fifteen are also identified in the logistic regression model (Table 2). We then used unsupervised clustering of the normalized expression of the top 40 predictors to show their expression pattern between samples. This analysis stratifies the control and LMM groups in two different clusters (Figure 1B). The heat map of the expression profile shows a high homogeneity in each group when we consider this subset of 40 discriminating transcript clusters (Figure 1B).

Figure 1. Gene expression pattern from the logistical regression analysis. (A) The volcano plot was constructed using the full list of 67528 transcript clusters analyzed. The top 821 transcripts were highlighted in blue (556 PPR) and green (255 PNR). A p-value of <0.01 and |beta| >1 was considered statistically significant. (B) The heat map shows the unsupervised clustering of the normalized expression pattern. The dendrogram indicates two clusters that stratified LMM and control group in distinct clusters. The right panel of the heat map shows the bar graph of beta coefficient values. Abbreviations: Beta: Estimated logistic regression coefficient; PPR: Predictor with Positive Relationship also called Up in panel A and B; PNR: Predictor with Negative Relationship also called Down in panel A and B.

Support vector machine-based LMM classification model

RFE algorithm discriminates a list of the 45 most informative variables from the discriminating transcripts clusters identified in the logistic regression adjustment. This subset of transcript clusters was used to train a SVM model with a polynomial kernel of degree 1. Model performance was evaluated in the test cohort and shows an accuracy of 0.70. This list of 45 classifiers is enriched in differentially expressed transcript clusters based on their fold changes (Supplementary Table 3).

Functional gene set enrichment highlights an alteration of muscle development and metabolism pathways in LMM

Genes associated with the subset of transcript clusters with positive and negative regression coefficients were then used to perform to different gene set analyses with ClueGO in Cytoscape software. Using the lists of annotated genes, ClueGO identifies and clusters biological pathways and gene ontology terms that participate in the same biological function, thereby removing redundancy. The top significant and non-redundant biological pathways and gene ontology terms associated with PPR genes were shown in Figure 2A. Muscle development pathways and metabolism ontology terms were overrepresented in this analysis. ClueGO shows few gene interactions between the different enriched gene set groups, suggesting an absence of hub genes involved in the simultaneous regulation of different pathways (Figure 2A). Genes with negative regression coefficients were associated with metabolism deregulation (Supplementary Figure 1). Even though these analyses provide us interesting information about the mechanisms underlying the LMM, they did not reveal the human phenotypes that were associated with the observed alterations.

Figure 2. Functional analysis of genes and human phenotypes associated with LMM. (A) Pathways that are associated with the gene discriminated in the logistic regression (with positive coefficient) are indicated by colored nodes. GSA terms are interconnected with their associated genes. Related GSA terms are indicated by the same color. (B) The human phenotype ontologies are presented in green and shared genes are represented as red (positive coefficient) and blue (negative coefficient) dots. Abbreviation: GSA: Gene Set Analysis.

Human phenotype ontology (HPO) analysis reveals the involvement of immune and inflammatory alterations

To gain additional insights into the phenotypes previously associated with the LMM discriminating genes (identified in the logistic regression analysis), a HPO analysis was performed with the full list of predictors using GenePattern [12]. Figure 2B shows several relevant phenotypes up-regulated in LMM. Phenotypes associated with immune system deregulation, inflammation and pathologies involving muscle disfunction were overrepresented in the HPO network. Interestingly, these phenotypes are highly connected, suggesting a tight relationship between them and highlighting a probable association of their common mediators with the physiological mechanisms that underlie LMM development. In addition, based on network connectivity analysis, the following genes were predicted as hub genes of the HPO network: (i) CLCNKB and KDSR (abnormality of endocrine system, glucose intolerance, abnormality of immune system physiology, abnormal inflammatory response, decreased body weight); (ii) KIAA0586 and OCRL (abnormality of endocrine system, abnormality of immune system physiology, abnormal inflammatory response, decreased body weight, abnormal phalanx morphology); (iii) FANCM (abnormality of endocrine system, abnormality of immune system physiology, decreased body weight, abnormal phalanx morphology); (iv) TGDS (abnormality of immune system physiology, abnormal inflammatory response, decreased body weight, abnormal phalanx morphology); (v) NDN (abnormality of endocrine system, glucose intolerance, abnormality of immune system physiology, decreased body weight). These genes are highly connected with most of the phenotypes and might play prominent role in the regulation of LMM.

Prediction of the regulatory networks

The regulatory network system upstream to the expression of the discriminating genes (with positive regression coefficients) was reconstructed using Expression2Kinase tool. This prediction tool allows us to identify in Figure 3A the top 20 predicted transcription factors (TFs) that can regulate the expression of the genes detected in this study (Figure 3A). The circos diagram in Figure 3B represents the relationship between the top 2 TFs that target the largest number of genes identified in our study (49 and 42 genes for FOXL1 and TCF4 respectively). We then analyzed the subnetworks of regulation and observed a common regulatory network pattern associated with TCF4 and HNF1A that was identified as a PPR (estimated coefficient β = 22.87; p-value = 0.007). The regulatory network in Figure 3C presents these TFs (red nodes); (ii) the intermediate proteins likely to be involved in the formation of a regulatory complex with these TFs (green nodes) and (iii) the kinases that are expected to be the activators of the regulatory complexes (dark blue nodes). This predicted regulatory network connects HNF1A and TCF4 and suggests a possibility of a co-regulatory system involving both TFs (Figure 3C).

Figure 3. Reconstruction of the predicted transcriptional regulatory network. TF regulatory network was reconstructed using Expression2Kinase tool. (A) Top 20 predicted TFs were presented and ranked based on their combined scores. The intensity of the red coloration is proportional to the combined scores. The ratio of the target genes (x axis) indicates the proportion of genes targeted by a determined TF. (B) The circos diagram shows the interaction between the top two predicted TFs (TCF4 in red and FOXL1 in blue), ranked by the ratio of target genes) and their targeted genes (gray). Both TFs target 16 common genes. (C) Co-regulatory network of TCF4 and HNF1A. The transcriptional regulatory network is presented with the transcription factors (TFs) in red, the intermediate protein that are predicted to interact with these TFs in green and the kinases in dark blue.

Gene modules detected by weighted gene co-expression network analysis are correlated with patient's characteristics

Fourteen co-expression modules were identified in the weighted co-expression analysis performed exclusively with the patient group. Correlation of clinical and laboratorial data with gene modules shows tight relationship with relevant parameters (Figure 4A). Pink module is moderately correlated with both serum levels of vitamin D (r = 0.63, p-value = 0.004) and Agatston score (r = 0.51, p-value = 0.02). This module regroups genes involved in immune system regulation including: (I) innate immunity, phagocytosis, and neutrophil degranulation, (ii) B cell receptor and IL1 signaling pathways (Figure 4B). In contrast, green-yellow and light-green modules are negatively correlated with plasma level of vitamin D and body fat percentage (r = −0.58, p-value = 0.009 and −0.51, p-value = 0.05 respectively). Both green-yellow and light-green modules are associated with ubiquitination mechanisms (Figure 4C). In addition, genes of green-yellow are involved in erythrocyte specific metabolism pathways and oxidoreduction processes (Supplementary Figure 1).

Figure 4. Weighted gene co-expression network analysis reveals modules correlated with patient´s characteristics. (A) Correlation matrix of detected gene modules with the characteristics of patients. Red color indicates positive correlation and blue color indicates negative correlation. Pink module is positively correlated with plasma levels of vitamin D and Agatston score. (B) Pink module is enriched with immune system pathways. (C) Green-yellow module is enriched with pathways associated with ubiquitination processes.

Subjects and microarray dataset

The subjects of this study were selected from the population-based survey (São Paulo Ageing and Health [SPAH] study) followed by the bone metabolism outpatient clinic of the Faculty of Medicine of the University of Sao Paulo (FMUSP). Transcriptome analysis of a total of 90 elderly women from this population were performed and described in our previous studies [33, 34] using microarray technology. From this microarray dataset publicly deposited in the GEO NCBI database under the accession number GSE152073, we selected 40 elderly women, of which 20 were with LMM (defined by a Newman's residual <−1.32) and 20 age and race-matched controls (residual >−1.32). Clinical and laboratory characteristics of this cohort were extracted from our database. Race was defined based on the self-reported race of second-generation ancestors, an approach previously used for the Brazilian population [35]. Physical activity was classified as (a) low, when not even housework is performed; (b) moderate, for regular housework, non-regular walking, gardening; and (c) high, for regular housework and daily regular physical activity at least twice a week for 30 min [36]. This research was approved by the local Ethics in Research Committee of the Medicine Faculty of São Paulo University/Brazil.

Transcriptomic analysis

Functional gene set analysis

To facilitate the interpretation of the biological relevance of discriminating genes selected from the logistic regression adjustment, a functional gene set analysis (GSA) was performed using ClueGO [39]. ClueGO, a plug-in of Cytoscape [40], predicts the functional gene ontology terms and biological pathways associated with our genes list. Then, it organizes the detected pathways in functionally grouped networks to highlight the relevant biological relationship between them. We applied a fusion criterion to reduce the redundancy of the terms that have similar associated gene sets. To further filter the most important terms and pathways, two-sided hyper-geometric distribution tests was used and terms at a significant level Bonferroni adjusted p-value of ≤0.05 were selected. To gain more insights into the human phenotype previously associated with the gene sets, a single-sample GSEA (ssGSEA) was also performed with GenePattern [12, 41] using the Human Phenotype Ontology (HPO) database [42]. Differentially activated HPO were then identified by comparing ssGSA score between LMM and control using RankProd package [43].

Regulatory gene network prediction

To refine the exploratory analysis, the prediction of the regulatory networks of transcription factors (TFs) and kinases underlying gene expression of the PPR gene list was performed using Expression2Kinases (X2K) tool [44]. Only the top 20 list of human TFs and kinases ranked based on p-value were selected. Transcriptional regulatory network was visualized using Cytoscape software [40].

Predictive model of LMM using support vector machine (SVM)

We submitted the discriminating transcript clusters identified by logistic regression to Recursive Feature Elimination (RFE) algorithm using e1071 and caret package. RFE was applied to the same subset of microarray dataset included in the logistic regression analysis (20 LMM and 20 controls) to select the most informative features (transcript clusters) for the SVM model training. The Training cohort was consisted of 70% of the SPAH dataset randomly selected and the remaining sample (test cohort) was used to validate the model. SVM model was trained with a polynomial kernel of degree 1. Model performance was evaluated based on its accuracy.

Weighted gene co-expression network analysis (WGCNA)

Gene co-expression analysis was performed using WGCNA method [45] as previously described elsewhere [26]. We included in this analysis the 20 LMM patients previously analyzed in the logistic regression analysis. To construct the weighted gene co-expression networks, WGCNA determines the connection strengths between nodes (protein coding genes) by computing the pairwise correlations between their expression profiles and using a soft thresholding of the Pearson correlation. Modules of the weighted gene co-expression network constituted of groups of highly correlated genes were identified and correlated with relevant laboratorial and clinical data of the patients. Module highly correlated with patient's data were selected for further gene set enrichment analysis using ClueGO [39] plugin in Cytoscape [40].