Gene expression profiling in the Cynomolgus macaque Macaca fascicularis shows variation within the normal birth range

Abstract

Background

Although an adverse early-life environment has been linked to an increased risk of developing the metabolic syndrome, the molecular mechanisms underlying altered disease susceptibility as well as their relevance to humans are largely unknown. Importantly, emerging evidence suggests that these effects operate within the normal range of birth weights and involve mechanisms of developmental palsticity rather than pathology.

Method

To explore this further, we utilised a non-human primate model Macaca fascicularis (Cynomolgus macaque) which shares with humans the same progressive history of the metabolic syndrome. Using microarray we compared tissues from neonates in the average birth weight (50-75^th centile) to those of lower birth weight (5-25^th centile) and studied the effect of different growth trajectories within the normal range on gene expression levels in the umbilical cord, neonatal liver and skeletal muscle.

Results

We identified 1973 genes which were differentially expressed in the three tissue types between average and low birth weight animals (P < 0.05). Gene ontology analysis identified that these genes were involved in metabolic processes including cellular lipid metabolism, cellular biosynthesis, cellular macromolecule synthesis, cellular nitrogen metabolism, cellular carbohydrate metabolism, cellular catabolism, nucleotide and nucleic acid metabolism, regulation of molecular functions, biological adhesion and development.

Conclusion

These differences in gene expression levels between animals in the upper and lower percentiles of the normal birth weight range may point towards early life metabolic adaptations that in later life result in differences in disease risk.

Background

Clinical, experimental and epidemiological studies have highlighted a link between the early-life environment and the health and well-being of offspring in later life. An adverse maternal environment has been linked to an increased risk of developing metabolic and cardiovascular disorders including type 2 diabetes, obesity, hyperlipidemia, insulin resistance and hypertension [1–7]. An important feature of these studies is that these relationships exist within the normative birth range and do not depend on extremes of birth weight. This has led to the proposal that later life disease risk is the result of maladaptive consequences of plastic mechanisms which would normally be adaptive [8, 9].

It is proposed that developmental plasticity determines the trajectory of development through epigenetic processes such that the fetus attempts to match its later phenotype to the environment [10]. It has been proposed that low birth weight is a marker of a poor early life nutritional environment [11] and thus a smaller fetus is more likely to develop a metabolic capacity appropriate for a low nutrient postnatal environment. But, if faced with a high nutrient environment it is more likely to become obese and insulin resistant [12]. Although, epigenetic processes have been increasingly implicated largely from rodent studies involving nutritional manipulation of the dam [13, 14] the molecular mechanisms underlying altered disease susceptibility are largely unknown. There is also some evidence that these developmental trajectories, and associated long-term gene expression and epigenetic changes can be reversed by the administration of the adipokine leptin to the neonatal rat although the concentrations used were higher than physiological levels [10, 12, 15, 16]. These data suggest that a better understanding of the molecular events associated with impaired early life development may help in designing future intervention strategies.

To identify the possible molecular pathways associated with variations in the fetal environment, we have utilised a non-human primate (NHP) model, the Macaca fascicularis (Cynomolgus macaque) to elucidate whether variations within the normal birth weight range are associated with differential gene expressions patterns. Cynomolgus macaques share with humans the same progressive history of the metabolic syndrome [17] which makes this model directly relevant to humans and importantly, Cynomolgus macaque is a monotocous species in which spontaneous variation in fetal growth rather than experimental manipulation can be investigated. This study therefore we have investigated the effect of spontaneous lower birth weight on gene expression in key tissues (umbilical cord, hepatic tissue and skeletal muscle) from female Cynomolgus macaque neonates.

Methods

Collection of Umblical cords

Sixty-five pregnant Cynomolgus macaque dams, sired naturally by one male, were monitored prior to delivery at the Vietnam Primate Breeding and Development Corporation. After birth, dams were sedated (ketamine-HCl; 7 mg/kg) to facilitate collection of the umbilical cord. The cords were collected and immediately snap-frozen in liquid nitrogen and stored at -80°C for later analyses. Neonates were weighed at birth and promptly returned to the dams. All animal procedures were approved by Nafovanny, subsidiary of the Ministry of Forestry, Vietnam, and performed in accordance with the guidelines set by the national advisory committee for laboratory animal research (NACLAR) of Singapore.

Collection of hepatic and skeletal muscle samples

The normative birth range was assessed from these 65 pregnancies and 8 neonates were selected based on their birth weights to comprise 2 groups: 1) lower birth weight group (LBW); n = 4 classified as those that were within the 5^th to 25^th birth weight percentile, birth weight range 299-317 g and 2) average birth weight group (ABW); n = 4 classified as those that were within the 50^th to 75^th birth weight percentile, birth weight range 358-398 g. The normal gestation of Cynomolgus macaque is approximately 155-170 days [18]. We have estimated the gestational age based on early ultrasound measurements (greatest length of the embryo at the time of pregnancy detection) and used those pregnancies where fetuses were within normal distribution for full term Cynomolgus macaques [18]. On postnatal day 5, neonates were sedated with an intramuscular injection of ketamine-HCl (15 mg/kg), and exsanguinated under anesthesia. Liver and skeletal muscle (biceps femoris) were collected and immediately snap frozen in liquid nitrogen and stored at -80°C for later analyses

Preparation of Total RNA

Total RNA was isolated from umbilical cords and neonatal tissues using TRIzol reagent according to the manufacturer's instructions (Invitrogen). RNA integrity was confirmed by bio analyzer 2100 (Agilent Technologies, Santa Clara, USA). An RNA Integrity Number (RIN) value of 7.5 above was considered acceptable and used in further experiments.

Microarray analysis

For microarray analysis, RNA from 6 groups of samples (Cord: ABW and LBW; Liver: ABW and LBW; skeletal muscle: ABW and LBW) were labeled using QuickAmp 1-color labeling kit (Agilent Technologies) according to manufacturer's protocol. The Cy3 labeled cRNA were subsequently hybridized to Agilent Rhesus Macaque (G2519F-015421) Gene Expression microarray. The Rhesus macaque gene expression microarray used in this study represented 43,803 Rhesus monkey probes synthesized as 60-mers spotted using the Agilent SurePrint technology (Agilent Technologies). The microarrays were scanned with Agilent High resolution Scanner and the images were feature extracted using FE software version 10.5 (Agilent Technologies).

Data analysis was performed using Genespring GX ver10 (Agilent Technologies). The raw signal intensity from each samples is global normalized to 75th percentile and base-line transformed. Probes flagged with present call in at least 75% of the samples in any of the 6 groups were used for subsequent data analysis.

Two-way ANOVA was performed with p value cut-off at 0.05 to identify genes that are differentially expressed in the tissue type and birth weight. Due to limited annotation of Rhesus Monkey genome, the microarray probes are annotated against human genes and ontology. For mapping against the human genome, the probes from the monkey microarray were re-annotated using Agilent eARRAY. The probe sequences were aligned to human genome (hg18) using BLAST based algorithm and the associated human annotation was extracted from the database.

To study the gene expression profile of the birth weight in each of the tissue type, Welch T-Test with p-value cut off of 0.05 and fold change of 1.5X was performed between the ABW and LBW samples in each of the tissue groups.

Complete microarray data is available at the Gene Expression Omnibus (GEO) database under the accession number GSE32069.

Network Analysis

The microarray data was imported into Pathway Studio version 7 (Ariadne Genomics, Rockville, USA) for network analysis. Gene Set Enrichment Analysis (GSEA) was performed on ABW vs LBW in the respective tissue using Kolmogorov-Smirnov algorithm with p-value cut-off at 0.05. In addition, Network Enrichment Analysis (NEA) was performed to identify expression hub of the treatment. Sub network was generated by connecting entities to their expression target network using the Resnet 7 Mammalian database.

Quantitative RT-PCR

Five up-regulated and two down regulated genes were selected for verification using qRT-PCR. The primers were designed using the Cynomolgus if available or Rhesus macaque sequences using the primer 3 software if not available [19]. The gene symbols and the primers are listed in Table 1.

Table 1 Sequences of primers used for qRT-PCR:

We used skeletal muscle RNA to verify the array. Briefly, total RNA was extracted as mentioned above from four ABW and four LBW neonates and were reverse transcribed using Applied Biosystem's high-capacity cDNA reverse transcription kit using 1 μg of total RNA in a reaction volume of 20 μl. The PCR reactions were carried out in 25 μl of SYBR Green Master Mix with 200 ng of cDNA using 7500 real time PCR system (Applied Biosystems, CA, USA). The comparative Ct method was used to calculate the relative gene expression [20]. 18S RNA was used as the internal control which was validated using the method described in Schmittgen and Livak [20] and found to be stable and consistent.

Results

Microarray analysis

Two-way factorial ANOVA identified 1973 genes which were differentially expressed in the three tissue types between ABW and LBW neonates (P < 0.05). Of these, 1141 genes were up regulated in the LBW samples while 832 genes were down-regulated in the LBW samples compared to ABW (Figure 1).

Figure 1

Hierarchical clustergram of 1973 genes (1141 up regulated in LBW and 832 down regulated in LBW) identified by ANOVA (P < 0.05) in all the three tissues (umbilical cord (C), liver (L) and skeletal muscle (B) analyzed. The relative expression is reflected by the intensity of the color (Green = down regulated, red = up regulated)

Gene expression profiles of umbilical cord, liver and skeletal muscle were also compared using Welch-t-test. There were 250 genes significantly and differently expressed in the liver, 850 genes significantly and differently expressed in the skeletal muscle and 891 genes significantly and differently expressed in cord samples (P < 0.05, >1.5 fold difference) (Figure 2, 3). The top 50 genes whose gene expression levels changed in each tissue based on their fold difference based on Welch T-test are given in Tables 2, 3, 4.

Figure 2

Heat map depicting the differentially expressed genes in skeletal muscle (850 genes), liver (210 genes) and cord (891 genes) (P < 0.05, >= 1.5 fold difference). A. 733 genes (584 genes up regulated in LBW skeletal muscle, 149 genes down regulated in LBW skeletal muscle), B. 19 genes (16 genes up regulated in LBW skeletal muscle and 15 genes up regulated in LBW liver, 3 genes down regulated in LBW skeletal muscle and 4 genes down regulated in liver). C 182 genes (115 genes up regulated in LBW liver and 67 genes down regulated in LBW liver). D 5 genes (2 genes up regulated in LBW cords and 3 genes up regulated LBW liver, 3 genes down regulated in LBW cords and 2 genes down regulated in LBW liver). E. 788 genes (338 genes up regulated in LBW cord and 450 genes down regulated LBW cord). F. 94 genes (60 genes up regulated in LBW skeletal muscle and 22 up regulated in LBW cord, 34 genes down regulated in LBW skeletal muscle and 72 genes down regulated in LBW cord. G. 4 genes (4 genes up regulated in LBW skeletal muscle and liver, 4 genes down regulated in LBW cord).

Figure 3

Venn diagram depicting the differentially expressed genes in skeletal muscle (850 genes), liver (210 genes) and cord (891 genes) t-test unpaired unequal variance, LBW vs ABW, p < 0.05, >= 1.5 fold difference.

Table 2 List of top 50 genes significantly differentially regulated in the skeletal muscle based on fold difference (t-test).

Table 3 List of top 50 genes significantly differentially regulated in liver based on fold difference (t-test).

Table 4 Top 50 genes significantly differentially regulated in cord based on fold difference (t-test).

Of the 250 genes which were differently expressed in the liver, 182 genes were unique to the liver (115 genes up regulated in the LBW group and 67 genes down-regulated in the LBW group). There were 19 genes which were significantly and differently expressed between liver and skeletal muscle (16 up regulated in LBW skeletal muscle and 15 genes up regulated in liver and 3 down-regulated in the LBW skeletal muscle and 4 down regulated in LBW liver). There were 5 genes which are significantly and differently expressed between liver and cord (3 up regulated in LBW liver and 2 genes up regulated in LBW cord and 2 down regulated in LBW liver and 3 down regulated in the cord).

Of the 850 genes significantly and differently expressed in skeletal muscle, 733 genes were specific to the skeletal muscle; i.e. showed altered regulation only in the skeletal muscle (584 genes up regulated in the LBW samples and 149 genes down regulated in the LBW group). There were 94 genes which were significantly and differently expressed between skeletal muscle and cord (60 up regulated in LBW skeletal muscle and 22 genes up regulated in LBW cord and 34 genes down regulated in the LBW skeletal muscle and 72 genes down regulated in the LBW cord).

Of the 891 genes significantly and differently expressed in umbilical cord; 788 genes showed altered regulation only in umbilical cord (338 genes up regulated and 450 genes down regulated in LBW group). There were 4 genes which are significantly and differently expressed between liver, skeletal muscle and umbilical cord (4 genes up regulated in the LBW skeletal muscle and liver while the same four genes were down regulated in LBW cord) Figure 2, 3.

Functional classification of genes

Gene ontology was used to classify genes based on functional significance. The main component of the Gene Ontology (GO) annotation taken into consideration was metabolic function. Genes were classified into fifteen functional categories: Cellular lipid metabolic process (43 genes), Cellular biosynthetic process (95 genes), Cellular macromolecule synthesis (222 genes), Cellular nitrogen compound metabolic process (10 genes), Cellular carbohydrate metabolic process (6 genes), Cellular catabolic process (9 genes), Nucelobase, Nucleoside, nucleotide and nucleic acid metabolic process (216 genes), Other cellular metabolic process (29 genes), Other metabolic process (36 genes), Transport (141 genes), Regulation of molecular functions (28 genes), Biological adhesion (27 genes), Developmental process (74 genes) Other biological processes (252 genes) and Non classified (795) with p-value of (p > 0.05). The genes enriched for each GO term were further classified into the number of genes up regulated or down regulated in each tissue with a fold difference of 1.5 or greater (Table 5) (Additional file 1).

Table 5 Gene ontology classification to group genes using Genespring ver10 (Agilent Tech, Santa Clara) of similar functional families.

Quantitative RT-PCR

To validate the micoarray results we carried out quantitative RT-PCR (qRT-PCR) using the same RNA samples in the microarray analysis. We selected seven genes (a novel gene XM_116936, PAS domain containing serine/threonine kinase (PASK), Adenisine kinase (ADK) transcript variant ADK-short, ELMO/CED-12 domain containing 1 (ELMOD1), Sine oculins homeobox homolog 1 (SIX1) Retinoblastoma like-1 (RBL1) and Solute carrier family 12 (potassium/chloride transporters) member 9 (SLC12A9) for this validation of the array using skeletal muscle cDNA. Of these 5 genes were significantly up regulated in the array and 2 genes were significantly down regulated in the array. Our results from the qPCR complement our results from the microarray (Table 6). The fold differences along with the values which derived from the microarray are presented in Table 6.

Table 6 Verification of seven genes from the microarray using Real-time RT-PCR analysis in skeletal muscle.

Discussion

In the present study, we have identified genes involved in key metabolic signaling pathways in three tissue types in a non-human primate model, that were differentially expressed according to the birth weight of the animal. Importantly, this differential expression was across the normal birth weight spectrum and therefore likely to represent adaptive pathways that the fetuses uses to predict its postnatal environment. The identification of novel signaling pathways that appear to be regulated by the early life environment is a key step in designing future experimental paradigms to understand the association between birth weight and disease risk. Metabolic disease particularly, has been strongly with early life adversity [21, 22]. Our data begin to shed light on the key signaling pathways that are vulnerable to subtle changes in the early life environment.

The strength of our study, despite its small size, is that we have focused on infants whose growth was not experimentally manipulated but lay within normal birth weight range. Many experimental models have manipulated pregnancy in an effort to produce fetal growth restriction. Such studies have shown that offspring which are born growth restricted catch-up in growth with their normally nourished counterparts and in adulthood are obese, hypertensive, hyperinsulinemic, leptin resistant and display sedentary behavior [23–25]. Investigations into underlying mechanisms and the determination of gene expression levels that may explain these altered phenotypes have produced conflicting results [26–28] which may reflect variations in the model systems used and the gender of the animals used [28]. Taken together, although these studies have established the link between early life nutritional adversity to later pathophysiology, there are limitations in the interpretation of rodent studies in development as applied to humans.

In the present study we aimed to study the molecular associates of growth variation within the normative range and exclude pathology. This is because the growing literature on developmental outcomes highlights that the importance of variation in risk is associated with non-pathological developmental environments. Accordingly we studied relatively small infants born between the 5^th and 25^th centile but excluded the smallest neonates, which may reflect obstetrical abnormalities. These infants were compared to infants in the middle of the normative range (50^th-75^th centile) and accordingly we excluded infants who may have had macrosomia as a result of the mother's being over-nourished by being maintained in captivity. Thus we are confident that we have excluded pathological influences and demonstrated that within the normative range patterns of gene expression may vary considerably with variable birth weights. Indeed we found a number of genes with more than a 10 fold shift in expression levels. There are important implications to this observation. Historically, experimental and epidemiological focus has been on the extremes of birth weight (either small for gestational age or large for gestational age) and there has been little focus within the normal birth weight range continuum. What is evident from the present study is that relatively small changes in the birth phenotype may be associated with profound changes in molecular physiology. In turn this also suggests the presence of highly evolved processes by which the fetus can adjust its development in response to subtle cues from mother [29].

The Cynomolgus, as in the human, has monotocous pregnancies with haemochorial placentae; they have omnivorous diets and monogastric digestion. They also share with humans the same progressive history of the metabolic syndrome [17]. We have identified alterations in the levels and expression patterns of a number of genes involved in different metabolic processes including cellular lipid metabolism, cellular biosynthesis, cellular macromolecule synthesis, cellular nitrogen metabolism, cellular carbohydrate metabolism, cellular catabolism, nucleotide and nucleic acid metabolism, biological adhesion and development.Recently, transcriptional profiling in rats subjected to gestational under nutrition was performed in young adult male rats where 249 genes were shown to be differentially expressed in the liver [28]. We compared the genes which are significantly altered in our array with those identified in the rat array study and have identified twelve similar or closely related genes from those identified in the rat: Tribbles homolog 2 (Trib2); 3-hydroxyanthanillate dioxygenase (Haao), transmembrane serine protease 6 (Tmmprss6); Thioredoxin domain containing10 (Txndc10); tralation initiation factor 4A3 (Eifa3); Ribosomal protein L31 (Rpl31); Danse 1-like 2 (Dnase1l2); Quinolate phosphoribosyl transferase (Qprt); general transcription factor IIa 2 (Gtf2a3); General transcription factor II H 3 [28]. Only four of these genes (Trib2, Trip10, EIF4A3 and Dnase1l2) were altered in the same direction in the livers of LBW macaques as in the rat array. We have also compared our observations with those identified from LBW term placentas by McCarthy and colleagues [30] and found similarities in expression changes in the genes Procollagen-lysine, 2-oxoglutarate 5-dioygense (PLOD2); Soluble interleukin-1 receptor accessory protein (IL1RAP); Solute carrier family 2 (facilitated glucose transporter) member3 (SLC2A3); Myosin V1(MYH6); Ribosomal protein S6(RPS6) and Latexin (LXN). The Tribbles homolog 2 (Trib2) is also increased in the LBW infants and suggests another possible way in way insulin/IGF-1 signalling might be impaired during development. Tribbles belongs to a family of kinase-like proteins and are reported to be negative regulators of Akt, the principle target of insulin signaling [31, 32].

Studies using animal models such as rodents to understand the developmental origins of metabolic diseases have shown that epigenetic changes in genes correlates with expression changes including metabolic enzymes such as PEPCK, transcriptional factors such as PPARα which regulate fat metabolism and factors associated with insulin action (e.g. PI3 kinase, PKC-ζ), the key regulatory genes which are responsible for bringing these changes are not known [15, 33]. One aim in our study was to identify whether there were shifts within the expression of key early regulators and from the array we have identified one such key regulatory gene the PAS Kinase (PASK), an evolutionarily conserved PAS domain containing serine/threonine kinase whose expression is altered as a result of adverse early developmental conditions. PAS domains are evolutionarily conserved and appear from archaea, bacteria to eukaryotes and are present in many signaling proteins where they act as signal sensor domain [34]. The PASK gene, whereby expression was up regulated in the LBW animals, has been shown to be a metabolic sensor based on mice knock-out studies; mice lacking PASK are resistant to high-fat induced obesity, hepatic steatosis and are resistant to insulin [35].

Conclusions

In summary, this paper has identified significant variation in gene expression in multiple tissues in primate newborns of different growth trajectories but within the normative range of birth size. Further detailed analyses may improve our understanding of how alterations in such genes due to adverse early life environment predisposes towards metabolic syndrome. These data give strength to the hypothesis that developmental plasticity operating within the normative range of birth weights can influence metabolic and other physiological systems in a way that might have later health consequences. It emphasizes that the concept of developmental programming need not involve pathological changes in growth trajectories to have molecular and presumably functional consequences.