Evaluation of reference genes for gene expression studies in human brown adipose tissue

Abstract

Human brown adipose tissue (BAT) has during the last 5 year been subjected to an increasing research interest, due to its putative function as a target for future obesity treatments. The most commonly used method for molecular studies of human BAT is the quantitative polymerase chain reaction (qPCR). This method requires normalization to a reference gene (genes with uniform expression under different experimental conditions, e.g. similar expression levels between human BAT and WAT), but so far no evaluation of reference genes for human BAT has been performed. Two different microarray datasets with samples containing human BAT were used to search for genes with low variability in expression levels. Seven genes (FAM96B, GNB1, GNB2, HUWE1, PSMB2, RING1 and TPT1) identified by microarray analysis, and 8 commonly used reference genes (18S, B2M, GAPDH, LRP10, PPIA, RPLP0, UBC, and YWHAZ) were selected and further analyzed by quantitative PCR in both BAT containing perirenal adipose tissue and subcutaneous adipose tissue. Results were analyzed using 2 different algorithms (Normfinder and geNorm). Most of the commonly used reference genes displayed acceptably low variability (geNorm M-values <0.5) in the samples analyzed, but the novel reference genes identified by microarray displayed an even lower variability (M-values <0.25). Our data suggests that PSMB2, GNB2 and GNB1 are suitable novel reference genes for qPCR analysis of human BAT and we recommend that they are included in future gene expression studies of human BAT.

Keywords:

• brown adipose tissue
• human
• quantitative PCR
• reference genes
• white adipose tissue

Abbreviations

BAT	=	Brown adipose tissue
CV	=	Coefficient of variation
PCR	=	Polymerase chain reaction
PR-BAT	=	Perirenal BAT
PT-BAT	=	Perithyroid BAT
qPCR	=	Quantitative PCR
RMA	=	Robust multi-array average
SOS	=	Swedish Obese Subjects
WAT	=	White adipose tissue

Introduction

It is generally considered that 2 different types of adipose tissues exists; white adipose tissue (WAT) and brown adipose tissue (BAT). They are also considered to have distinct and reciprocal biological roles. WAT is the main site for energy storage, whereas BAT has the unique ability to convert chemical energy to heat through uncoupling of the oxidative phosphorylation.¹

BAT has been extensively studied in rodents. In humans, it has been known since the seventies that brown adipocytes exist in adults,² but it is only in the recent years that the amount of BAT found in human adults have been considered large enough to play an important metabolic role.³ The conversion of excess calories into heat via BAT has made manipulations of this tissue a very attractive treatment strategy for obesity and obesity related-co morbidities.⁴ However, there is still relatively little knowledge of this “new human organ”, especially at the molecular level. The expression levels of relatively few genes have been investigated in human BAT and quantitative polymerase chain reaction (qPCR) has been the method of choice in most of these cases.

Quantitative PCR is a well-established technique for quantification of messenger ribonucleic acids (mRNA). The sensitivity and accuracy of the method makes it suitable for mRNA analysis in various biological samples. However, many factors such as the reverse transcriptase step, PCR efficiency, sample loading and tissue type, can affect the qPCR result. Therefore, the method depends on reference genes (genes with uniform expression under different experimental conditions) to correct for these variations. Some specific genes, like 18S, PPIA and GAPDH, are frequently used as reference genes and even though they are considered to be stably expressed, high variability in expression of these genes has been reported.^5-7 In order to select the most stable reference gene in a specific experimental setting, statistical algorithms like geNorm and Normfinder have been developed. With these methods and other similar techniques, researchers have published suitable reference genes in different human tissues.^5,8-12 Our research group has previously evaluated reference genes in WAT under different experimental conditions.¹¹ Chechi and coworkers recently published a validation study on reference genes in human epicardial adipose tissue,¹² an adipose tissue depot shown to contain brown adipocytes.¹³ Apart from this, there are no available data on the variability of reference genes in the study of human BAT.

The aim of this study was to identify and evaluate reference genes that can be used in real-time PCR analysis of human BAT. To this end, we used microarray data to find new reference gene candidates. These novel gene candidates and commonly used reference genes were evaluated with qPCR in 2 human adipose tissue settings using 2 different algorithms.

Results

Microarray-based selection of novel candidate reference genes

Initially, “genes with intermediate high expression levels” were selected separately from the 2 different microarray datasets (3304 genes from the PT-BAT dataset and 3346 genes from the PR-BAT data set) to avoid genes with very high or low expression levels. These genes were ranked based on their coefficient of variation (CV) values and the 100 genes with the lowest CV values from each dataset were cross-referenced in the 2 datasets. A CV-score was calculated by multiplying the 2 CV values. The 10 genes with the lowest CV-score are listed in Table 1. The genes with the lowest CV-scores were HUWE1 and RING1 for the PR-BAT and the PT-BAT datasets, respectively. Seven genes with the lowest CV-score were selected for qPCR analysis.

Table 1. Microarray based evaluation and selection of putative reference genes

Gene Title	Gene Symbol*	PR-BAT CV	PT-BAT CV	CV score**
Commonly used reference genes
Beta-2-microglobulin	B2M	0.075	0.111	0.0084
Ubiquitin C	UBC	0.065	0.140	0.0090
Ribosomal protein, large, P0	RPLP0	0.101	0.117	0.0118
Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide	YWHAZ	0.095	0.132	0.0125
Peptidylprolyl isomerase A (cyclophilin A)	PPIA	0.127	0.110	0.0140
Low density lipoprotein receptor-related protein 10	LRP10	0.100	0.176	0.0176
Polymerase (RNA) II (DNA directed) polypeptide A, 220kDa	POLR2A	0.083	0.215	0.0179
Phosphoglycerate kinase 1	PGK1	0.163	0.136	0.0222
Actin, β	ACTB	0.101	0.239	0.0242
TATA box binding protein	TBP	0.145	0.179	0.0259
Importin 8	IPO8	0.168	0.171	0.0287
Glucuronidase, β	GUSB	0.138	0.249	0.0344
Glyceraldehyde-3-phosphate dehydrogenase	GAPDH	0.152	0.283	0.0429
Hydroxymethylbilane synthase	HMBS	0.267	0.240	0.0642
Hypoxanthine phosphoribosyltransferase 1	HPRT1	0.274	0.417	0.1144
Transferrin receptor (p90, CD71)	TFRC	0.382	0.530	0.2022
18S rRNA	18S	n.a.	0.245	n.a.
PR-BAT dataset selected putative reference genes
HECT, UBA and WWE domain containing 1	HUWE1	0.044	0.082	0.0036
Guanine nucleotide binding protein (G protein), β polypeptide 2	GNB2	0.046	0.105	0.0048
Family with sequence similarity 96, member B	FAM96B	0.063	0.090	0.0057
Proteasome (prosome, macropain) subunit, β type, 2	PSMB2	0.061	0.097	0.0059
Tumor protein, translationally-controlled 1	TPT1	0.064	0.098	0.0063
Guanine nucleotide binding protein (G protein), β polypeptide 1	GNB1	0.060	0.105	0.0063
Transmembrane protein 248	TMEM248	0.046	0.151	0.0069
Nuclear receptor binding protein 1	NRBP1	0.053	0.136	0.0072
SCY1-like 1 (S. cerevisiae)	SCYL1	0.056	0.131	0.0074
ADP-ribosylation factor 1	ARF1	0.062	0.123	0.0076
PT-BAT data set selected putative reference genes
Ring finger protein 1	RING1	0.065	0.069	0.0045
Family with sequence similarity 96, member B	FAM96B	0.063	0.090	0.0057
Proteasome (prosome, macropain) subunit, β type, 2	PSMB2	0.061	0.097	0.0059
Guanine nucleotide binding protein (G protein), β polypeptide 1	GNB1	0.060	0.105	0.0063
GNAS complex locus	GNAS	0.081	0.091	0.0074
ADP-ribosylation factor 1	ARF1	0.062	0.123	0.0076
Ring finger protein 1	RING1	0.065	0.127	0.0082
RNA binding protein S1, serine-rich domain	RNPS1	0.079	0.110	0.0087
Transcription factor 25 (basic helix-loop-helix)	TCF25	0.077	0.115	0.0089
ELAV (embryonic lethal, abnormal vision, Drosophila)-like 1 (Hu antigen R)	ELAVL1	0.073	0.122	0.0089

*Gene symbols marked in bold text were selected for qPCR analysis.

**CV-score is the product of the CV values from the PT-BAT and PR-BAT datasets.

n.a., not available.

Microarray based analysis of commonly used reference genes

In addition to the unbiased candidate reference gene search, CV-scores for 17 commonly used reference genes were also calculated using the PR-BAT and PT-BAT datasets (Table 1). B2M, UBC and RPLP0 displayed the lowest CV-scores of the commonly used reference genes. However, the CV scores for these reference genes were in general higher compared to the top 10 candidate reference genes from the unbiased candidate gene search (Table 1). Probe sets for the commonly used reference gene 18S was not included in the PR-BAT dataset and the 18S could therefore not be included in the comparison. Eight of the commonly used references genes were selected for qPCR analysis based on a low CV-score or that they have previously been used in analysis of human BAT.

qPCR evaluation of selected reference genes

To evaluate the stability of the candidate reference genes in perirenal adipose tissue (n = 20), 2 different algorithms were used. The M-values generated by geNorm and S-values generated by Normfinder software are illustrated in Figure 1A-B. The genes are ranked from the most stable, i.e. lowest M or S value, to the least stable reference genes. The results from both analyses were quite similar; GNB1, GNB2, PSMB2, RING1 and TPT1 were the most stable genes in both analyses. GeNorm ranked the combination of GNB1 and PSMB2 (M value = 0.10) as the most stable, whereas Normfinder favored the combination of GNB2 and PSMB2 (S value = 0.01). The commonly used reference genes displayed less stable M and S values (Fig. 1).

Figure 1. Evaluation of stability- and M-values for reference genes in adipose tissue (AT) using qPCR. The evaluation was performed in 2 sets of adipose tissue samples. The first set (A and B) included 20 samples of perirenal AT with half of them containing variable levels of brown adipocytes. The second set (C and D) included the previously analyzed preirenal AT samples as well as 20 s.c. AT samples from both lean and obese subjects. qPCR results from both sets of samples were analyzed using 2 different algorithms; the normfinder algorithm (A and C) and the geNorm algorithm (B and D). Gene symbols in bold text are classified as commonly used reference genes.

When subcutaneous adipose tissue samples from lean and obese individuals were included in the analysis together with the PR-BAT samples (n = 40), the increased heterogeneity of the tissue samples were reflected in increased M and S values (Fig. 1C-D). Like the geNorm analysis of PR-BAT the combination of GNB1 and PSMB2 (M = 0.12) was the most stable. However, in the NormFinder algorithm PSMB2 and RPLP0 (S = 0.13) was suggested as the best combination of reference genes.

Among the commonly used reference genes, RPLP0, LRP10, YWHAZ and GAPDH generally displayed the lowest variability (S or M values) in the 2 groups of samples analyzed (Fig. 1). In all analysis, 18S displayed distinctly higher variability than the other reference genes analyzed (Fig. 1).

Discussion

In this study, we have identified reference genes that may be suitable for qPCR analysis in human BAT. Most of the investigated commonly used reference genes display adequately low variability but novel reference genes identified in this study, such as GNB1, GNB2 and PSMB2, display even lower variability.

There are several challenges in the molecular and functional characterization of human BAT and brown adipocytes. The first challenge is that none of the sites where BAT is commonly found in humans are easily assessable. In general, some type of surgical procedure is needed to obtain the BAT-containing samples. Secondly, the samples obtained are usually very heterogeneous in the amount of BAT that they contain, analyzed as UCP1 expression levels.^14-17 If the samples used are very heterogeneous in their amount of BAT, this will increase the variability of any parameter analyzed. Hence, discrete differences in brown adipocytes compared to white adipocytes may be missed due to sample heterogeneity. If the parameter analyzed is gene expression by qPCR, variability within the reference genes may increase the total variability of the results. Consequently, finding stable reference genes is important for the future molecular evaluation of human BAT.

In the current literature regarding qPCR analysis of human BAT, a large number of different reference genes have been used. Some examples of reference genes that has been used is TATA-binding protein (TBP),^14,18 cyclophilin A (PPIA),^13,19 β actin (ACTB),^15,20 glyceraldehyde-3-phosphate dehydrogenase (GAPDH),¹³ large ribosomal protein P0 (RPLP0)^16,17 and 18S rRNA.²¹ In our study, all these reference genes, with the exception of 18S, seemed reasonable robust. A relatively high degree of variability for 18S expression has also been demonstrated in brown adipocyte containing human epicardial adipose tissue.¹² For the geNorm algorithm, a M-value below 0.5 has been described as typical for stable genes in a homogeneous set of samples.²² In our analysis, the majority of genes meet this criterion (Fig. 1B and D). However, a stronger consensus of what references genes to use in the analysis of human BAT would increase comparability of different gene expression studies as well as the fidelity of the results. It would also be interesting to analyses if the newly identified reference genes also are stable in rodent models of brow/beige adipose tissue.

In addition, the stability of the selected reference genes has also been analyzed separately in the 3 different groups of WAT samples (white UCP1 negative perirenal adipose tissue, and s.c WAT from lean and obese individuals) included in the qPCR study. These genes show similar stability also in this setting (data not shown) and could therefore also be considered as stable reference genes in different WAT depots and in s.c. WAT from individuals with different degree of obesity. However, specifically designed studies to verify and extend these findings are desirable.

One of the limitations of this study is that BAT from only one human depot (perirenal adipose tissue) was used in the qPCR evaluation of the reference genes (we have no perithyroid BAT samples left for additional analysis). Another limitation is that none of these samples are entirely composed of BAT but rather represents islets of brown adipocytes within a majority of white adipocytes. One additional limitation is that the reference genes have only been evaluated in a BAT/WAT context and not compared to any other human tissues. The strength of the study lies within the possibility to identify novel reference genes using human BAT expression profiles, since the expression profiling procedure is not dependent on reference genes for normalization.

We conclude that most of the commonly used reference genes displayed acceptably low variability, whereas 18S, in this study, displayed much higher variability. In addition, we have identified novel reference genes, such as GNB1, GNB2 and PSMB2, which are even more stable than the commonly used reference genes. The use of reference genes with low variability will increase the fidelity in qPCR results from human BAT.

Subjects, Materials and Methods

Subjects

All study participants received written and oral information before giving written informed consent. The Regional Ethics Committee in Gothenburg approved these studies and the studies followed the Helsinki declaration.

Perithyroid BAT (PT-BAT) Study

Perithyroid adipose tissue biopsies were obtained during surgery in the thyroid area at the Sahlgrenska University Hospital (Gothenburg, Sweden) as previously described.¹⁶ Paired biopsies were also taken from the subcutaneous depot in the surgical incision area. Nine paired BAT containing perithyroid adipose tissue and subcutaneous adipose tissue biopsies were analyzed by microarray (U133Plus2.0 microarray, Affymetrix, #900470).

Perirenal BAT (PR-BAT) Study

Perirenal adipose tissue biopsies were obtained by laparoscopic technique from kidney donors at the department of transplantation at the Sahlgrenska University Hospital (Gothenburg, Sweden) as previously described.¹⁷ In this study, samples from 10 study participants with biopsies containing BAT (UCP1+) and 10 sex, age, and BMI matched controls with biopsy samples not containing BAT (UCP1-) were used. The samples were subjected to both microarray (Expression assay Gene 1.0 ST, Affymetrix, #901085) and qPCR analysis. Brief descriptions on the study participants can be found in Table 2.

Table 2. Characteristics of the study participants

Samples	BMI (Mean± SD)	Age (Mean± SD)	Gender (Males/Females)	Biopsy procedure
Perirenal (UCP1+)	25.7 ± 2.2	42.4 ± 13.8	6M/4F	Surgery
Perirenal (UCP1−)	25.8 ± 3.3	44.0 ± 9.2	6M/4F	Surgery
Subcutaneous (Lean)	22.5 ± 1.0	35.2 ± 6.8	5M/5F	Needle
Subcutaneous (Obese)	35.1 ± 3.3	40.2 ± 5.5	5M/5F	Needle

The Swedish Obese Subjects (SOS) SibPair Study

The SOS Sib Pair study consists of 732 study participants from 154 nuclear families.^23,24 Each family contained BMI discordant sibling pairs (BMI difference >10 kg/m2). Subcutaneous adipose tissue biopsies were obtained by needle aspirations. In this study, biopsies from 10 lean and 10 obese study participants were used for qPCR analysis. Brief descriptions on the study participants can be found in Table 2.

Microarray based selection of candidate reference genes

Total adipose tissue RNA was prepared with the RNeasy lipid tissue kit (QIAGEN, # 74804), or using the phenol-chloroform extraction method of Chomczynski and Sacchi.²⁵ Preparation of cDNA and hybridization to DNA microarray was performed according to standard Affymetrix protocols (Affymetrix). Microarray raw data was analyzed using the RMA (robust multi-array average) method.

Initially, signal values (RMA values converted to linear scale) and CV values for all genes were calculated in both the PT-BAT and the PR-BAT datasets. Selection of genes with intermediate high expression levels was performed after graphic illustration of the distribution of signal values. Different signal cut off values were used for the 2 different datasets due to differences in overall signal intensities (signal values between 50-500 for the PT-BAT data set and 500-5000 for the PR-BAT dataset). Genes with signal values within these ranges were selected for further investigation. The selected genes were then ranked based on their CV values and the 100 genes with the lowest CV values from each data set were cross-referenced to the corresponding other dataset. A CV-score (the product of the CV values from the 2 datasets) was calculated and genes were selected for qPCR analysis bases on a low CV-score.

Selection of commonly used reference genes

Sixteen genes included in the TaqMan Human Endogenous Control Arrays (http://tools.lifetechnologies.com/content/sfs/manuals/cms_042704.pdf) were selected as commonly used reference genes. In addition, LRP10 was also included in this analysis since our research group has previously shown that this is a suitable reference gene in human WAT.¹¹ CV-scores for these 17 genes were calculated as described above using both the PT-BAT and the PR-BAT datasets. If multiple probe sets for a selected gene were present in the microarray data set, the probe set with the highest signal value was used in the CV-score calculation.

qPCR analysis

Total RNA was reversed transcribed using High capacity cDNA Reverse Transcription kit (Applied Biosystems, #4368814) according to the manufacturer's instruction. Recombinant RNaseout® Ribonuclease Inhibitor (Invitrogen, #10777-019) was added to prevent RNase-mediated degradation. All the cDNA-reactions were run in duplicate and thereafter pooled for the qPCR analysis.

The qPCR was performed using TaqMan® Custom Array cards (Applied Biosystems). The arrays were designed as CustomFormat 32 with TaqMan probe and primer sets for reference candidate genes in triplicates (Table 3). Each port on the TaqMan® Arrays was loaded with cDNA corresponding to 100 ng total RNA, combined with nuclease free water and 50 μl TaqMan® Gene Expression Master Mix (Applied Biosystems, #4369016) to a final volume of 100 μl. The TaqMan® Arrays were analyzed using the 7900HT system with a TaqMan Array Upgrade (Applied Biosystems). Thermal cycling conditions were: 50°C for 2 min, 94.5°C for 10 min, followed by 40 cycles of 97°C for 30 s and 59.7°C for 1 min.

Table 3. List of candidate reference genes with gene symbol, Assay ID and mean Cq value

Gene Symbol	Assay ID	Cq Mean (all groups)
B2M	Hs00984230_m1	21.7
UBC	Hs01871556_s1	27.0
RPLP0	Hs99999902_m1	22.9
YWHAZ	Hs03044281_g1	23.0
PPIA	Hs04194521_s1	28.5
LRP10	Hs00204094_m1	25.9
GAPDH	Hs02758991_g1	23.4
18S	Hs99999901_s1	13.3
HUWE1	Hs00948075_m1	26.2
GNB2	Hs00929275_g1	25.6
TPT1	Hs02621289_g1	21.2
RING1	Hs00968517_m1	27.2
FAM96B	Hs00211406_m1	27.5
PSMB2	Hs01002946_m1	26.2
GNB1	Hs00929799_m1	25.5

The TaqMan® Custom Array cards were analyzed using Expression Suite Software V.1.0.3 (Applied Biosystems). To evaluate the candidate reference genes we used 2 well established algorithms, geNorm and Normfinder.^26,27 geNorm was developed in 2002 by Vandesompele and coworkers. It determines pairwise variation among all candidate reference genes and generates an M value, called gene-stability measurement, a dimensionless parameter that results from calculating the average of all logarithms of expression ratios between a putative gene and all other reference genes being evaluated. Genes with the lowest M values have the most stable expression. The NormFinder algorithm ranks the set of candidate normalization genes according to their expression stability S for each reference. Importantly, if the samples belong to different treatment groups, NormFinder also takes into account the intragroup and intergroup variations.

Disclosure of Potential Conflicts of Interest

MKS is currently employed by Amgen AB, Solna, Sweden. None of the other authors declare any conflict of interest.

Acknowledgment

We gratefully acknowledge the valuable contributions of Jenny Palming.

Funding

This study was supported by Grants from The Sahlgrenska Academy, the Swedish federal government under the LUA/ALF agreement concerning research and education of doctors, research grants from AstraZeneca R&D, the Research Foundation of the Swedish Diabetes Association and 7TH Framework program International Research Staff Exchange Scheme (IRSES), FUEGO project.