Topological analysis of metabolic networks integrating co-segregating transcriptomes and metabolomes in type 2 diabetic rat congenic series
Animals
A colony of GK/Ox rats bred locally and derived in 1995 from the GK/Par colony was
used to produce the congenic strains. BN rats were obtained from a commercial supplier
(Charles River Laboratories, Margate, UK). All congenics were derived from these strains
using a genetic marker-assisted breeding strategy (“speed congenics”) as previously
described 11], 44] and maintained by brother–sister mating. They were specifically designed to contain
GK alleles over genomic regions of various lengths (1–176 Mb) of rat chromosome 1,
introgressed onto the genetic background of the BN strain (Table 1, Fig. 1). The targeted GK genomic intervals between markers D1Rat27 (90.3 Mb) and D1Got254
(264.37 Mb) cover several cardiometabolic disease (CMD) -relevant QTLs originally
mapped in F2 (GKxBN) genetic crosses 9], 45], 46]. GK alleles on the X chromosome were lost early in the breeding program by two consecutive
breedings of male backcross animals to BN females. All animals used in this study
were systematically genotyped with markers chosen to accurately monitor retention
of GK alleles across the congenic intervals and absence of GK allele contaminants
from the genetic background 47], 48]. All strains were co-housed in order to avoid cage-specific microbiome selection
14], 49].
Table 1. Genomic details of the BN.GK congenic strains
Fig. 1. Phenotypic features of the congenic strains. Association for body weight (a), adiposity index (b), QTL for adiposity index (c), association for cumulative plasma glucose (d), and insulin (e) during the intravenous glucose tolerance and insulin secretion tests were measured
in male rats. Solid bars represent the GK genomic segments of chromosome 1 of each congenic strain introgressed
onto the genetic background of the BN strain. The y-axis shows genomic length (Mb) and boundaries of the genomic region of GK origin. The
location of the adipose tissue QTL mapped to chromosome 1 in the GK?×?BN F2 cross
9] is reported with significance threshold shown with a dashed line (c). Details of GK chromosomal regions introgressed in each congenic strain are given
in Table 1. Significantly (P??0.05) increased and decreased values of the phenotypes between congenic strains
and the BN control are indicated in red and green, respectively. Phenotype data are available in Additional file 1: Table S1
Rats were allowed free access to tap water and standard laboratory chow pellets (BK
Universal Ltd, Grimston, Aldbrough, Hull, UK) and were maintained on a 12-h light–dark
cycle. All rats were identified using a microchip (identity chip, Animal Care Ltd,
York, UK) linked to a database specifically designed to administer the project (husbandry,
phenotype scheduling, and data storage) and store genetic information and phentoypic
data 48], 50].
Phenotype analysis
Three-month-old male congenic rats and BN controls were used for all experiments.
Intravenous glucose tolerance and insulin secretion tests (IVGTTs) were performed
following procedures strictly identical to those consistently applied in both F2 (GK?×?BN)
hybrids 9], 20], which we used to map diabetes QTLs in the GK, and BN.GK congenic strains derived
for several GK QTLs 11], 22], 51]–54]. Briefly, rats in the post-absorptive state at the end of the post-prandial glycemic
response (4.5 to 5 h fasting from 9–9:30 am when food was removed until 2 pm when
the IVGTT procedures were initiated) were anesthetized by injection of 95 mg/kg body
weight intraperitoneal ketamine hydrochloride (Ketalar, Parke-Davies, UK). Rats were
injected with a solution of 0.8 g glucose/kg body weight via the saphenous vein. Blood
samples were collected into heparinized tubes before glucose injection and 5, 10,
15, 20, and 30 min afterward. Plasma was separated by centrifugation prior to glucose
assays using a diagnostic kit (ABX, Shefford, UK) on a Cobas Mira Plus automatic analyzer
(ABX, Shefford, UK) and assay of immunoreactive insulin (IRI) using an ELISA kit (Mercodia,
Uppsala, Sweden). Cumulative values of plasma glucose and plasma insulin during the
IVGTTs were calculated to evaluate overall glucose tolerance and insulin secretion
capacity in response to glucose, respectively. At 6 months, rats were killed by terminal
anesthesia following an overnight fast (16–18 h). Retroperitoneal fat pads (RFPs)
were collected, weighed, snap frozen in liquid nitrogen, and stored at ?80 °C until
preparation of tissue extracts and RNA for analysis of the metabolome and the transcriptome,
respectively. The adiposity index was calculated as the ratio between RFP weight and
body weight.
Metabotyping of white adipose tissue by
1
H NMR spectroscopy
Tissues samples (30–50 mg) were weighed into 2-mL Eppendorf tubes and were each homogenized
in 1.5 ml 50 % methanol using TissueLyser (5 min at 25 Hz; QIAGEN, Germany). The mixtures
were each transferred into 3-mL glass tubes and 0.7 mL chloroform was added into each
sample. The mixtures were vortexed for 1 min followed by centrifugation at ~3500 g
for 25 min at 10 °C. The aqueous phase was decanted and the methanol was removed under
a fume cupboard before freeze drying. The lipid phase was pipetted out and chloroform
was removed under a fume cupboard. Dried extracts were reconstituted using 500 ?L
of 0.1 M phosphate buffer solution (10 %
2
H
2
O/H
2
O v/v, with 0.05 % sodium 3-trimethylsilyl-(2,2,3,3-
2
H4)-1-propionate for chemical shift reference at ?0.0) in 5-mm tubes for NMR acquisition.
Standard
1
H NMR spectra were measured on a Bruker spectrometer (Rheinstetten, Germany) operating
at 600.22 MHz
1
H frequency, as described previously 18], 55]. The
1
H NMR spectra were imported into Matlab and phase- and baseline-corrected at high
resolution. The region ?5.0–4.5 was removed to eliminate baseline effects of imperfect
water signal pre-saturation. Each spectrum was normalized to a constant intensity
sum and each variable was mean centered. Analyses were carried out using R and Matlab.
Orthogonal partial least squares discriminant analysis
The method allows enhanced focus on strain and diet intervention while minimizing
other biological/analytical variation. Sample classes were modeled using the orthogonal
partial least squares (O-PLS) algorithm. This algorithm derives from the partial least
squares (PLS) regression method. In linkage analysis version, the model explains the
maximum separation between genotypes Y (coded as 0, 1, 2 for GK allele numbers) using the NMR data X. Further details on standard O-PLS implementation in metabonomics
have been given previously 56], 57]. The model coefficients locate the NMR signals significantly associated with genotypic
variation in a specific genomic region Y.
RNA preparation and Illumina bead array hybridization
Total RNA was individually isolated from 100 mg of RFP (four biological replicates
per strain) using the RNeasy® 96 Universal Tissue kit (Qiagen, Crawley, UK): frozen
tissue samples were transferred into cooled RNeasy® 96 Universal Tissue plates and
homogenized in QIAzol Lysis Reagent using a Qiagen Tissue Lyser. Following phase separation
after addition of chloroform, total RNA was purified with RNeasy columns using a spin
technology according to the manufacturer’s guidelines and eluted in RNase-free water.
RNA concentrations were determined using a NanoDrop spectrophotometer and RNA quality,
purity, and integrity were assessed using an Agilent 2100 Bioanalyser (Agilent Technologies,
Waldbronn, Germany).
Samples were independently used to hybridize Gene Expression Sentrix® BeadChip RatRef-12
v1 arrays (Illumina Inc., San Diego, CA, USA) containing 22,523 oligonucleotide probes
(replicated 30 times on average). They allowed interrogation of transcript levels
for 21,910 genes (6274 RefSeq NM transcripts, 15,983 Refseq XM transcripts, 12 Refseq
XR transcripts, 250 Unigene clusters). Double-stranded cDNA and purified biotin-labeled
cRNA were synthesized from 300 ng high quality total RNA using the Illumina® TotalPrep
RNA Amplification Kit (Ambion Inc., Austin, TX, USA). cRNA concentrations were determined
using a NanoDrop spectrophotometer whilst cRNA quality and integrity were assessed
using an Agilent 2100 Bioanalyser (Agilent Technologies, Waldbronn, Germany). Hybridizations
onto the arrays were carried out using 750 ng of each (132) biotinylated cRNA in a
58 °C hybridization oven for 18 h. Following washing and staining with Streptavidin-Cy3,
the BeadChip Arrays were scanned on the Illumia® BeadArray Reader (Illumina Inc.,
San Diego, CA, USA). Resulting data were then preliminarily analyzed using the Illumina®
BeadStudio Application software before undergoing comprehensive statistical analysis.
Particular attention was given to the following quality control parameters: 0???G
sat???1; Green 95 percentile (GP95) for consistency between arrays (around 2000);
GP5 background level in the range of the low 100 s or below.
Microarray experiments were compliant with MIAME (Minimum Information About a Microarray
Experiment) and both protocol details and raw data have been deposited in ArrayExpress
(http://www.ebi.ac.uk/arrayexpress/) under accession number E-MTAB-1048.
Network representation of genome–metabolome associations
In order to explore genome–metabolome associations, a functional association network
was derived from metabotype and genotype correlation coefficients using the bipartite
graph Rgraphsviz package from R to represent the O-PLS correlation matrix derived
from the linkage analysis between NMR variables and genotypes. A cutoff was then applied
to the P value of the Pearson’s correlation coefficient adjusted using Benjamini and Hochberg’s
multiple testing correction (P BHadj
??0.05); significant correlations were set to 1 and non-significant correlations
were set to 0, defining the adjacency matrix for the graph. Hence, G?=?(N, E) specifies
a graph G with N denoting the two node sets (two types of nodes: genomic regions and
metabolites) and E the edge set (link between nodes, here a correlation between metabotypes
and genotypes above the cutoff, i.e., P BHadj
??0.05).
Integration of eQTL-responsive genes and mQTL-responsive metabolites on KEGG metabolic
networks to identify functional candidates
Differentially regulated genes and metabolites that could be associated with GK haplotypes
were mapped to KEGG pathways. Metabolic pathways were imported in R using KEGGgraph
55], 58]. An in-house python script was written to generate an adjacency table. The resulting
region-specific networks were mathematically formalized as a directed multilabeled
graph Gm,n?=?(Vm,E) composed of two types of nodes Vm, where m?=?(“metabolic reactions”,
“metabolites”), and functional equally-weighted edges corresponding to network connectivity
(i.e., metabolic reactions) between enzymes and metabolites. To identify pairs of
eQTL-associated genes coding for enzymes that are metabolically connected with mQTL-associated
metabolites, we computed shortest paths from the eQTL genes mapped on the KEGG network
to the target mQTL metabolites using the igraph R package. Defining shortest path
lengths (spl) corresponds to counting the minimal number of additional reactions required
to connect a given gene and a given metabolite across the metabolic network. For example,
Galm and ?-D-glucose have a spl of 0 as ? -D-glucose is the product of the reaction catalyzed
by Galm (rn:R01602). The gene Asns is annotated as AsnsA and AsnsB since it involves two catalytic sites for two different
reactions (rn:R00256 for AsnsB and rn:R00578 for AsnsA).
shRNA-based inhibition of Galm and Asns expression in vitro in 3T3-L1 cells
3T3-L1 fibroblasts (ATCC, Molsheim, France) were cultured in 10 % calf serum (PAA,
Velizy, France) containing DMEM high glucose (Life Technologies, Saint Aubin, France).
Cells were plated at 10
5
cells/well density until confluence and differentiated into adipocytes in 10 % fetal
bovine serum (FBS; Life Technologies, Saint Aubin, France) containing DMEM high glucose,
IBMX, dexamethasone (Sigma-Aldrich, Saint-Quentin, France), and insulin. Differentiated
3T3-L1 cells were maintained in 5 % FBS and DMEM high glucose.
We used pGFP-V-RS-shRNA plasmids (Origene, Rockville, MD, USA) containing short hairpin
RNA (shRNA) sequences specifically designed to target Asns or Galm and a puromycin resistance gene cloned between integrative long terminal repeat sequences.
The Platinum-Ecotropic Retroviral Packaging Cell Line (Cell Biolabs, San Diego, CA,
USA) producing host range recombinant ?-retroviruses was used for shRNA-containing
viral production. Platinum E cells were maintained in DMEM supplemented with glucose
and 10 % FBS, puromycin, and blasticidin (Sigma-Aldrich, Saint-Quentin, France). Transfection
was performed after adapting culture medium for 3T3-L1 cells (DMEM high glucose and
10 % calf serum without antibiotics). Fugene 6 HD® (Roche, Boulogne, France) was used
according to the manufacturer’s recommendation in 3T3-L1 medium. Plasmid transfection
efficiency was determined by GFP fluorescence. Supernatants were collected 48 h after
transfection to transduce 3T3-L1 in the presence of polybrene (Sigma-Aldrich, Saint-Quentin,
France). Puromycin selection started 24 h post-transduction.
Differentiation and glucose transport in shRNA-transfected 3T3-L1 adipocytes
For differentiation analysis, cells were first incubated in 10 % formaldehyde (Sigma-Aldrich,
Saint-Quentin, France), washed with 60 % isopropanol, and dried. A solution Oil Red
O (Sigma-Aldrich, Saint-Quentin, France) was added and dishes were washed with distilled
water. Quantification of coloration was performed by spectrophotometry at 490 nm with
a plate reader (Perkin Elmer, Villebon, France). A separate batch of cells was used
for glucose transport analysis, which was determined by incubation with a solution
containing 0.5 ?Ci tritium labeled 2-deoxy-D-glucose. Briefly, adipocytes were cultured
in DMEM high glucose without FBS for 4 h and washed with a buffer containing CaCl
2
, MgCl
2
, fatty acid free BSA (PAA, Velizy, France) in PBS (Life Technologies, Saint Aubin,
France), followed by 20 min of insulin stimulation (100 nM). After 10 min of incubation
with labeled 2-deoxy-D-glucose, cells were washed in ice-cold PBS and collected in
a solution of NaOH for radioactivity recording and protein content quantification.