Identification of biomarkers for childhood obesity based on expressional correlation and functional similarity.

The aim of the current study was to identify potential biomarkers of childhood obesity, and investigate molecular mechanisms and candidate agents in order to improve therapeutic strategies for childhood obesity. The GSE9624 gene expression profile was downloaded from the Gene Expression Omnibus database. The differentially expressed genes (DEGs) in omental adipose tissues were analyzed with limma package by comparing samples from obese and normal control children. Two‑way hierarchical clustering was applied using the pheatmap package. The co‑expression (CE) analysis was performed using online CoExpress software. Subsequent to functional classification via the GOSim package, the gene network enriched by DEGs was visualized using the Cytoscape package. The codon usage bias of the DEGs was then examined using the CAI program from the European Molecular Biology Open Software Suite. In total, 583 DEGs (273 upregulated genes and 310 downregulated genes) were observed in the omental adipose tissues between samples from obese and normal control children. Hierarchical clustering identified a significant difference between samples from obese and normal control children. Subsequent to CE analysis, 130 DEGs, which were classified into 4 clusters, were selected. The following 3 upregulated and 2 downregulated genes were identified to be significant: Upregulated genes, microtubule‑associated protein tau (MAPT), destrin (actin depolymerizing factor) (DSTN) and spectrin, β, non‑erythrocytic 1 (SPTBN1); downregulated genes, Rho/Rac guanine nucleotide exchange factor 2 (ARHGEF2) and spindle and kinetochore associated complex subunit 1 (SKA1). The top 3 amino acids were identified to be glycine, leucine and serine with a high bias. The DEGs MAPT, DSTN, SPTBN1, ARHGEF2 and SKA1 are suggested to be candidate biomarkers for childhood obesity.

Introduction

Childhood obesity is a prevalent disease worldwide (1)and has been recognized as a serious public health concern due to its increasing morbidity in addition to adverse health effects (2). Compared with normal-weight children, obese children have an increased probability of obesity in adulthood (3), which may lead to hypertension, diabetes, dyslipidemia, chronic heart disease and other disorders (4). In addition, the prevalence of depression in obese children is greater than those of normal-weight children (5).

Obesity in children is a multifactorial disease, to which numerous genetic factors contribute (6). Previous genetic studies have identified that the mutation in melanocortin-4 receptor (MC4R) is associated with childhood obesity through binding with α-melanocyte-stimulating hormone (α-MSH). It has been identified that α-MSH could inhibit feeding, thus leading to an increasing risk of obesity (7–10). A missense amino acid substitution (R236G) has been reported to contribute to an inherited susceptibility to obesity. The R236G mutation contributes to an aberrant fusion protein, which is capable of interfering with central melanocortin signaling and appears to increase the risk of early-onset obesity (11–13). The fat mass and obesity associated gene has been recognized as another gene relevant to childhood obesity by influencing appetite and body composition (14,15). However, the effective approaches of preventing obesity development in children remain limited and further studies of the molecular mechanisms are required.

The aim of the present study was to identify potential biomarkers of childhood obesity. Furthermore, the current study aimed to elucidate the molecular mechanism for the therapy of childhood obesity. The significant differentially expressed genes (DEGs) in omental adipose tissues were screened by comparison of samples between obese and normal-weight children with two-way hierarchical clustering, co-expression (CE) analysis and gene network construction. The notable amino acids were identified by measuring of codon usage bias.

Materials and methods

Source of data

The microarray expression profile of GSE9624 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE9624) was extracted from the Gene Expression Omnibus database, based on the Affymetrix Human Genome U133 Plus 2.0 Array GPL570 (HG-U133_Plus_2). The data used was originally from 11 samples of omental adipose tissues, with 5 from obese children and 6 from normal-weight children.

Data preprocessing and DEGs analysis

The probe-level data in the CEL files were converted into expression measures using the Affy package (version 1.48.0) (16) with R software (https://www.r-project.org/). Missing data were inputted and data normalization was performed using quartile normalization (17). A boxplot graph was produced to present the chip data distribution and the median values. When the median values were consistent subsequent to normalization, the chip data were used for DEGs analysis.

The omental adipose tissues of normal-weight children were considered as the controls. DEGs in omental adipose tissues between obese children samples and normal controls were identified by the limma package (version 3.26.8) (https://bioconductor.org/packages/release/bioc/html/limma.html) (18). The P-values and log fold change (logFC) were calculated. P<0.01 and |logFC| >1 were considered as the cut-off criteria.

Hierarchical clustering of DEGs

Two-way hierarchical clustering was applied using the pheatmap package in R language (version 1.0.8) (https://cran.r-project.org/web/packages/pheatmap/index.html) (19). Subsequently, the heat map of the DEGs was generated. The clustering was performed using the Euclidean distance and Ward's method (20).

CE analysis of DEGs

The CE analysis was performed using online CoExpress software (http://www.bioinformatics.lu/CoExpress/). To identify CE DEGs, the correlation coefficient of each gene pair was calculated using the Pearson's correlation coefficient (21). The correlation coefficient ranged from −1 to +1, ‘−’ represented a negative correlation and ‘+ represented a positive correlation. When CE<0.95, the gene pairs were cut off. Subsequently, the CE gene network was constructed by connecting the CE genes with straight lines.

Gene network construction

The functional similarity of DEGs was calculated using GOSim package (version 1.8.0) (https://www.bioconductor.org/packages/release/bioc/html/GOSim.html) (22), which functions based on various information theoretic similarity concepts for gene ontology (GO) terms. P<0.05 was used as the cut-off criterion.

Following functional classification, the gene network enriched by the screened DEGs was constructed. The Cytoscape package (version 3.2.1) (http://www.cytoscape.org/) was applied for the network visualization (23).

Measures of codon usage bias

It has been previously identified that there are large differences in codon usage bias among genes with different functions (24). The program of codon usage analysis has been previously used for analyzing codon and amino acid usage patterns (25). In the current study, the codon usage bias of the DEGs in each cluster was examined using the codon adaptation index (CAI) value, which was calculated with the CAI program from the European Molecular Biology Open Software Suite (version 6.5.0) (http://emboss.sourceforge.net/index.html) (26). CAI values between 0 and 1 indicated a positive correlation between the codon usage bias and the CAI value.

Results

DEGs analysis

The boxplot graph demonstrated that the chip data had been normalized and were available for DEGs selection (Fig. 1).

Figure 1.

Boxplot graph of chip data distribution. The data were normalized and available for differentially expressed gene selection. The pink and white boxes represent samples of omental adipose tissues in obese and normal children, respectively. The x-axis represents samples; the y-axis represents gene expression values. The black horizontal line in each box represents the median value of the chip data, and the median indicates the gene expression value subsequent to normalization.

Based on the cut-off criteria of P<0.01 and |log FC| >1, 583 DEGs were identified in omental adipose tissues between the samples from obese children and those of normal-weight controls, including 310 upregulated genes and 273 downregulated genes.

The heat map of DEGs (Fig. 2) indicated a marked difference in the DEGs of omental adipose tissues between the samples from obese children and those of normal-weight controls, observed by the clear color difference.

Figure 2.

Heat map of differentially expressed genes. Red indicates high expression and blue indicates low expression. The color change from blue to red indicates movement from low to high levels of expression.

A total of 130 DEGs were selected using the cut-off criterion of CE>0.95. The CE network of those DEGs is presented in Fig. 3.

Figure 3.

Co-expression network of the differentially expressed genes. The downward and upward triangles refer to down- and upregulated genes, respectively. The edges represent the interactions.

Functional classification of the 130 DEGs was performed. According to the GO functional nodes of these DEGs, they were classified into 4 clusters (summarized in Table I).

Table I.

Functional classification of correlated differentially expressed genes.

Cluster name	Genes	Enrichment P-value	Biological function
C1 (BP)	ARHGEF2, MAPT, SPTBN1, SKA1, DSTN	0.0154	GO:0051493~regulation of cytoskeleton organization
	ARHGEF2, MAPT, SKA1	0.0199	GO:0031110~regulation of microtubule polymerization or depolymerization
	ARHGEF2, MAPT, SKA1	0.0350	GO:0070507~regulation of microtubule cytoskeleton organization
	ARHGEF2, MAPT, SKA1	0.0464	GO:0032886~regulation of microtubule-based process
C2 (BP)	ARHGEF2, MAPT, EIF5A, SPTBN1, DSTN	0.0004	GO:0043244~regulation of protein complex disassembly
	ARHGEF2, MAPT, SPTBN1, SKA1, DSTN	0.0154	GO:0051493~regulation of cytoskeleton organization
	ARHGEF2, MAPT, SPTBN1	0.0335	GO:0043242~negative regulation of protein complex disassembly
C3 (BP)	NRP2, LEP, SLITRK3, IRX5, EFNA2, EPHB3, SNAP25 NRP2, SLITRK3, EFNA2	0.0320	GO:0048666~neuron development
	EPHB3, SNAP25	0.0469	GO:0007409~axonogenesis
C4 (CC)	CENPN, CENPM, SKA1	0.0064	GO:0000777~condensed chromosome kinetochore
	CENPN, CENPM, SKA1	0.0080	GO:0000779~condensed chromosome, centromeric region
	CENPN, CENPM, SKA1	0.0104	GO:0000776~kinetochore
	CENPN, CENPM, SKA1	0.0221	GO:0000775~chromosome, centromeric region
	CENPN, CENPM, SKA1	0.0234	GO:0000793~condensed chromosome

BP, biological process; GO, gene ontology; C1-4, cluster 1–4; CC, cellular component; ARHGEF, 2 Rho/Rac guanine nucleotide exchange factor 2; MAPT, microtubule-associated protein tau; SPTBN1, spectrin, β, non-erythrocytic 1; SKA1, spindle and kinetochore associated complex subunit 1; DSTN, destrin (actin depolymerizing factor); EIF5A, eukaryotic translation initiation factor 5A-1; NRP2, neuropilin 2; LEP, leptin; SLITRK3, SLIT and NTRK-like family member 3; IRX5, iroquois homeobox 5; EFNA2, ephrin-A2; EPHB3, ephrin type-B receptor 3; SNAP25, synaptosomal-associated protein 25; CENPN, centromere protein N; CENPM, centromere protein M.

The gene network of the DEGs with similarities of expression and function is presented in Fig. 4. The gene network included 38 nodes, 34 DEGs and 35 edges. There were 3 significantly upregulated genes [microtubule-associated protein tau (MAPT), destrin (actin depolymerizing factor) (DSTN) and spectrin, β, non-erythrocytic 1 (SPTBN1)] and 2 significantly downregulated genes [Rho/Rac guanine nucleotide exchange factor 2 (ARHGEF2) and spindle and kinetochore associated complex subunit 1 (SKA1)].

Figure 4.

Gene network of differentially expressed genes with similarity of expression and function. The downward and upward triangles refer to down- and upregulated genes, respectively. The red circular nodes represent functional clusters. The edges represent the interactions.

Codon usage bias

The CAI values of the DEGs in the 4 clusters were 0.724 (cluster 1), 0.687 (cluster 2), 0.712 (cluster 3) and 0.705 (cluster 4), indicating a significant codon usage bias. To confirm the codon usage bias, amino acids translated by the codon in each sequence were summarized in Table II. The top 3 amino acids were glycine, leucine and serine, which exhibited a high bias.

Table II.

Codon usage of differential expression genes sequence in the 4 clusters.

						AA Sum
Codon	AA	Cluster 1	Cluster 2	Cluster 3	Cluster 4	Cluster 1	Cluster 2	Cluster 3	Cluster 4
GCA	A	194	26	695	190	712	148	2147	666
GCC	A	253	40	641	210
GCG	A	45	48	163	58
GCT	A	220	34	648	208
TGC	C	212	24	771	198	376	51	1579	427
TGT	C	164	27	808	229
GAC	D	128	20	469	117	270	40	1058	255
GAT	D	142	20	589	138
GAA	E	212	18	880	187	573	60	1731	436
GAG	E	361	42	851	249
TTC	F	203	46	857	216	531	108	2242	616
TTT	F	328	62	1385	400
GGA	G	281	46	842	198	1146	185	2907	765
GGC	G	255	43	658	203
GGG	G	387	62	825	180
GGT	G	223	34	582	184
CAC	H	179	22	693	200	341	43	1519	355
CAT	H	162	21	826	155
ATA	I	114	11	752	171	406	52	2302	528
ATC	I	127	19	605	153
ATT	I	165	22	945	204
AAA	K	294	18	1478	406	491	44	2328	609
AAG	K	197	26	850	203
CTA	L	134	18	470	136	1322	212	4625	1282
CTC	L	275	49	804	249
CTG	L	372	42	958	288
CTT	L	220	45	851	230
TTA	L	135	21	778	173
TTG	L	186	37	764	206
ATG	M	165	22	776	183	165	22	776	183
AAC	N	109	10	617	163	244	22	1550	384
AAT	N	135	12	933	221
CCA	P	273	35	880	214	1029	157	2750	682
CCC	P	339	57	831	170
CCG	P	86	17	213	69
CCT	P	331	48	826	229
CAA	Q	178	24	832	177	495	60	1815	448
CAG	Q	317	36	983	271
AGA	R	289	30	1020	228	846	137	2503	663
AGG	R	337	41	843	233
CGA	R	35	9	128	48
CGC	R	64	16	143	58
CGG	R	80	28	233	51
CGT	R	41	13	136	45
AGC	S	264	30	662	194	1235	182	3930	1118
AGT	S	207	23	666	188
TCA	S	191	31	795	227
TCC	S	252	44	771	202
TCG	S	47	12	126	40
TCT	S	274	42	910	267
ACA	T	164	12	845	206	538	72	2187	587
ACC	T	197	25	564	134
ACG	T	33	10	131	54
ACT	T	144	25	647	193
GTA	V	124	13	435	136	642	112	2144	634
GTC	V	136	24	433	126
GTG	V	246	39	680	212
GTT	V	136	36	596	160
TGG	W	331	42	945	233	331	42	945	233
TAC	Y	95	11	397	132	179	24	1090	305
TAT	Y	84	13	693	173
TAA	*	108	9	819	214
TAG	*	134	17	526	149
TGA	*	208	22	827	206

Columns 3–6 indicate the number of codons of the differentially expressed genes in each cluster. The last column shows the number of AA encoded by different codons. *No corresponding amino acid. AA, amino acid; A, alanine; C, cysteine; D, aspartate; E, glutamate; F, phenylalanine; G, glycine; H, histidine; I, isoleucine; K, lysine; L, leucine; M, methionine; N, asparagine; P, proline; Q, glutamine; R, arginine; S, serine; T, threonine; V, valine; W, tryptophan; Y, tyrosine.

Discussion

Childhood obesity is a multisystem disease with potentially life-threatening consequences (27). In spite of numerous genetic studies aiming to elucidate the pathogenesis of childhood obesity, the molecular mechanisms in the development and progression of this disease remain unclear. In the present study, 3 notably upregulated genes (MAPT, DSTN and SPTBN1) and 2 notable downregulated genes (ARHGEF2 and SKA1) were identified. Several functions, including microtubule polymerization or depolymerization, condensed chromosome kinetochore, regulation of cytoskeleton organization and regulation of microtubule cytoskeleton organization, were observed to be significantly enriched by these DEGs.

SKA1 is a protein-coding gene. The SKA1 complex is part of the conserved kinetochore-microtubule interface and directly associates with microtubules as part of oligomeric assemblies. It has been previously demonstrated to serve a critical role in interacting with dynamic microtubules at the outer kinetochore by depolymerizing microtubule ends (28). As presented in Table I and Fig. 4, SKA1 was significantly enriched in regulation of microtubule polymerization or depolymerization. The depolymerizing microtubules attached by macromolecular kinetochores are necessary for the chromosome movements to ensure regular chromosome segregation (29,30). Notably, SKA1 was observed in the current study to also be enriched in cluster 4, which consisted of biological functions such as condensed chromosome kinetochore. Previous studies have identified that chromosomal (such as those on chromosome 11p14-p12 and chromosome 16p11.2) deletions may result in obesity (31–34). Therefore, it is suggested that SKA1 may be associated with childhood obesity by depolymerizing microtubules and disturbing chromosome segregation.

An additional significantly downregulated gene identified in the present study is ARHGEF2. ARHGEF2 encodes guanine nucleotide exchange factor H1, which is a microtubule-associated exchange factor (35,36). ARHGEF2 activity is downregulated by interaction of its C-terminus with microtubules (37). Rho GEF is associated with microtubules and becomes active upon microtubule depolymerization (38). Through functional analysis, ARHGEF2 was observed to be significantly enriched in the regulation of microtubule polymerization or depolymerization. Hence, similar to that of SKA1, it is suggested that ARHGEF2 may be a candidate gene in childhood obesity, acting via targeting the regulation of microtubule polymerization or depolymerization.

In the present study, MAPT was observed as a significantly upregulated gene. MAPT encodes the microtubule-associated protein (MAP) tau. MAP has been observed to serve an important role in the promotion of microtubule assembly in vitro (39). The results of the current study indicated that MAPT targeted regulation of microtubule cytoskeleton organization. Previous studies have demonstrated that disruption of microtubule assembly inhibited the translocation of the insulin responsive glucose transporter isoform (GLUT4) (40,41). GLUT4 translocation promotes insulin to stimulate glucose uptake in adipose tissue, ultimately resulting in obesity (42). Thus, MAPT may serve a stimulative role in the development and progression of childhood obesity.

SPTBN1 is a member of the β-spectrin gene family (43). Spectrin is a protein functioning in actin cross-linking and the molecular scaffold. It connects the plasma membrane to the actin cytoskeleton and determines cytoskeleton organization including cell shape, arrangement of transmembrane proteins and organization of organelles (44). SPTBN1 was identified in the current study as a significantly upregulated gene, enriched in pathway of regulation of cytoskeleton organization. Cytoskeleton organization is a biological process involving various cellular components in adipose tissues. Fat deposition in mature fat cells and cell proliferation have been suggested to accelerate pre-adipocytes to differentiate into adipocytes, serving a role in the regulation of body weight (45). Thus, it is suggested that SPTBN1 may be a candidate gene for childhood obesity, acting via the regulation of cytoskeleton organization function.

DSTN encodes for the protein destrin. Destrin is a pH-independent protein, with the capacity to sever actin filaments (F-actin) (46). Destrin has been previously reported to serve a role in the reorganization of the actin cytoskeleton in response to stress and cell stimuli (47). In addition, destrin belongs to the actin-depolymerizing factor family, which mediates a pH-sensitive destruction of actin filaments (48–51). Filament networks, such as peripheral actin filaments, have been identified to have a mechanical connection to various cytoskeletal structures (52). By reviewing the functional classifications and the gene network in the current study, it is suggested that DSTN participates in the regulation of cytoskeletal organization. Previously, organization of the actin cytoskeleton was suggested to be necessary for early adipocyte differentiation, resulting in hyperplasia of adipose tissues which is a critical event for the development of obesity (53,54). This led to it being hypothesized that DSTN may be a biomarkers of childhood obesity.

However, there were limitations in the present study. The results obtained were web-based and were not verified by any biological experiments. In addition, the data downloaded were from a European database, it remains unclear whether these results would be consistent for children from other continents. Thus, further experimental studies based on the observations of the current study are required.

In conclusion, MAPT, DSTN, SPTBN1, ARHGEF2 and SKA1 may be candidate biomarkers of childhood obesity. Microtubule polymerization or depolymerization, condensed chromosome kinetochore, regulation of cytoskeleton organization and regulation of microtubule cytoskeleton organization may be important biological pathways in the progression of childhood obesity. The observations of the current study may have implications for the understanding of the molecular mechanisms, and identification of candidate therapeutic agents, for childhood obesity. However, further studies are required to confirm the preliminary observations outlined in the current study, with the aim of this information being used in a clinical setting.