Gene expression profiles in testis of pigs with extreme high and low levels of androstenone

Abstract

Background:

Boar taint is a major obstacle when using uncastrated male pigs for swine production. One of the main compounds causing this taint is androstenone, a pheromone produced in porcine testis. Here we use microarrays to study the expression of thousands of genes simultaneously in testis of high and low androstenone boars. The study allows identification of genes and pathways associated with elevated androstenone levels, which is essential for recognising potential molecular markers for breeding purposes.

Results:

Testicular tissue was collected from 60 boars, 30 with extreme high and 30 with extreme low levels of androstenone, from each of the two breeds Duroc and Norwegian Landrace. The samples were hybridised to porcine arrays containing 26,877 cDNA clones, detecting 563 and 160 genes that were differentially expressed (p < 0.01) in Duroc and Norwegian Landrace, respectively. Of these significantly up- and down-regulated clones, 72 were found to be common for the two breeds, suggesting the possibility of both general and breed specific mechanisms in regulation of, or response to androstenone levels in boars. Ten genes were chosen for verification of expression patterns by quantitative real competitive PCR and real-time PCR. As expected, our results point towards steroid hormone metabolism and biosynthesis as important biological processes for the androstenone levels, but other potential pathways were identified as well. Among these were oxidoreductase activity, ferric iron binding, iron ion binding and electron transport activities. Genes belonging to the cytochrome P450 and hydroxysteroid dehydrogenase families were highly up-regulated, in addition to several genes encoding different families of conjugation enzymes. Furthermore, a number of genes encoding transcription factors were found both up- and down-regulated. The high number of clones belonging to ferric iron and iron ion binding suggests an importance of these genes, and the association between these pathways and androstenone levels is not previously described.

Conclusion:

This study contributes to the understanding of the complex genetic system controlling and responding to androstenone levels in pig testis. The identification of new pathways and genes involved in the biosynthesis and metabolism of androstenone is an important first step towards finding molecular markers to reduce boar taint.

Background

The production of uncastrated male pigs is profitable to the swine production, because it leads to improved feed efficiency and greater lean yield of the carcass [1]. However, production of entire males is a challenge due to boar taint, an unpleasant odour and flavour often present in the meat of un-castrated male pigs [2]. Castration is undesirable both for ethical as well as for economical reasons and thus alternative methods are needed to prevent tainted meat. Boar taint is primarily caused by high levels of the two components androstenone and/or skatole in the pig carcasses [3, 4]. Moderate to high heritabilities have been shown for both androstenone and skatole levels in fat [5–7]. For example, a recent study of Duroc and Norwegian Landrace showed heritabilities ranging from 0.5–0.6 for androstenone and 0.23–0.56 for skatole [8].

Androstenone (5α-androst-16-en-3-one) is a 16-androstene steroid produced from pregnenolone in the steroid hormone pathway in boar testis near sexual maturity. It is released into the blood and transported to salivary glands [9], where it regulates reproductive functions in female pigs [10]. Due to its lipophilic properties, it is also easily stored in the adipose tissue [11], causing boar taint when the fat is heated. Another steroid produced from pregnenolone in the same pathway is testosterone. Testosterone stimulates growth and fertility and the challenge is to reduce the level of androstenone without affecting the level of testosterone. Genes specifically affecting the production of androstenone in the testis are therefore of interest, as well as genes involved in degradation and elimination of this steroid.

Several candidate genes have been suggested to affect levels of androstenone. The biosynthesis of 16-androstenes from pregnenolone is catalysed by the andien-β synthase enzyme system [12], which consists of cytochrome P450-C17 (CYP17) and cytochrome b5 (CYB5) [13]. Although CYP17 has been proposed as a very potent candidate gene affecting levels of androstenone, studies so far have not found significant effects of SNPs within the CYP17 gene and androstenone levels [14, 15]. Similarly, it has not been possible to find a significant correlation between levels of CYP17 and levels of 16-androstene steroids in fat [14]. CYB5 has been found positively correlated to the production of androstenone [14]. In addition, a SNP in the 5' un-translated region of porcine CYB5 is found to be associated with decreased fat androstenone levels [16].

Another potent class of candidate genes is the sulphotransferase (SULT) enzymes whose main function is to transfer a sulfo group to a range of molecules including for example steroid hormones [17]. The SULTs have been proposed to regulate the amount of unconjugated 5α-androstenone available for accumulation in fat [18], and high proportions of sulphoconjugated to unconjugated 16-androstene steroids are observed in porcine plasma and Leydig cells [18, 19].

The levels of androstenone vary between breeds, with e.g. higher androstenone levels in Duroc and Hampshire compared to Landrace and Yorkshire [20]. It is therefore of interest to investigate gene expression patterns in individuals with high or low levels of androstenone in different breeds. In this study boars with extreme levels of androstenone were selected from the two commercial breeding populations in Norway, Duroc and Norwegian Landrace, for gene expression analysis in testis using microarray technology. The objective of the study was to identify differentially expressed genes that could point towards pathways associated with extreme levels of androstenone in pigs. To our knowledge, this is the first microarray experiment performed in this context.

Results

Porcine cDNA microarrays containing 26,877 clones were used to examine the transcript profile of Duroc (D) and Norwegian Landrace (NL) pigs with high (H) and low (L) levels of androstenone. A total of 1533 NL boars and 1027 D boars were analysed and average androstenone levels were 1.17 ± 1.15 μg/g and 3.22 ± 2.68 μg/g for NL and D, respectively. 30 high and 30 low androstenone boars were selected from each breed and average values were 5.91 ± 2.41 μg/g for NL high androstenone (NLH), 0.15 ± 0.04 μg/g for NL low androstenone (NLL), 10.27 ± 2.68 μg/g for D high androstenone (DH) and 0.42 ± 0.13 μg/g for D low androstenone (DL). Testicle derived samples were hybridised using a balanced block design and significantly affected genes were identified using linear models (Limma). Box plots (See Additional file 1: Boxplot of the arrays) show that the deviation of the raw log ratios from 0 was successfully normalised. The statistical test detected 563 and 160 clones as differentially expressed in D and NL, respectively (p < 0.01) (See Additional file 2: Microarray results analysed using Limma, Duroc and Additional file 3: Microarray results analysed using Limma, Norwegian Landrace). The top 100 genes are presented for D (Table 1) and NL (Table 2). Venn diagrams were created to explore any overlap between breeds and 72 genes were found in common for D and NL (See Additional file 4: Genes common for Duroc and Norwegian Landrace at p < 0.01). Another more robust, but less powerful non-parametric test for statistical analysis, namely Fisher's Sign Test (FST) was performed as well. (See Additional file 5: Microarray results analysed using Fisher's sign test, Duroc and additional file 6: Microarray results analysed using Fisher's sign test, Norwegian Landrace).

Table 1 Top 100 genes identified in Duroc

Table 2 Top 100 genes identified in Norwegian Landrace

The cDNA clones found to be differentially expressed were used to search for statistically overrepresented gene ontology (GO) terms compared with the GOs represented by all the genes on the array. The most significant terms analysed for molecular function in D were ferric iron binding, iron ion binding, oxidoreductase activity and steroid binding (Figure 1). For NL we also found terms related to functions like electron transport and steroid dehydrogenase activity (Figure 1). Furthermore, the genes were classified according to their biological processes (Figure 2) and cellular components (See Additional file 7: Gene ontology (GO) results for the cellular component ontology, Duroc and Additional file 8: Gene ontology (GO) results for the cellular component ontology, Norwegian Landrace). The biological processes most significant in both breeds were steroid biosynthesis and steroid metabolism.

Figure 1

Gene ontology (GO) results for the biological process ontology. The genes differentially expressed at p < 0.01 were analysed for over-represented GO terms in the biological process ontology. Some terms were common for Duroc and Norwegian Landrace (red) while some were only significant in Duroc (yellow) or Norwegian Landrace (blue). Number of genes is the number of genes that was found significantly differentially expressed (p < 0.01) for a term.

Figure 2

Gene ontology (GO) results for the molecular function ontology. The genes differentially expressed at p < 0.01 were analysed for over-represented GO terms in the molecular function ontology. Some terms were common for Duroc and Norwegian Landrace (red) while some were only significant in Norwegian Landrace (blue). Number of genes is the number of genes that was found significantly differentially expressed (p < 0.01) for a term.

Both real competitive PCR (rcPCR) and real-time PCR were used to validate expression profiles for ten of the interesting differentially expressed genes from the microarray experiment. Real-time PCR was used to confirm the expression of CYP17, CYB5, ferritin light polypeptide (FTL) and translation elongation factor 1 alpha 1 (EEF1A1) in the D samples. The NL pigs were not included in this study. The expression levels of all genes were normalized to β-actin (ACTB) and 18S rRNA. Three of the four genes, FTL, CYB5 and CYP17, were confirmed by this method, while for EEF1A1 no significant difference was found (Table 3). rcPCR was used to analyse gene expression profiles of sphingomyelin phosphodiesterase 1 (SMPD1), steroidogenic acute regulatory protein (StAR), sulfotransferase family 2A, dehydroepiandrosterone-preferring, member 1 (SULT2A1), 17-beta-hydroxysteroid dehydrogenase IV (HSD17B4), aldo-keto reductase family 1, member C4 (AKR1C4), and cytochrome P450, subfamily XIA, polypeptide 1 (CYP11A1). The expression levels of all genes were normalised using the housekeeping gene hypoxanthine guanine phosphoribosyltransferase 1 (HPRT). StAR, SULT2A1, HSD17B4, AKR1C4 and CYP11A1 were all differentially expressed in both D and NL (p < 0.01), while SMPD1 was not significantly differentially expressed in either of the two breeds (Table 4 and 5, respectively).

Table 3 Results from Real-time PCR

Table 4 Results from the rcPCR bootstrap statistics (x4000), Duroc

Table 5 Results from the rcPCR bootstrap statistics (x4000), Landrace

Discussion

The present study describes the transcriptional profiles of testicle from animals with extreme high and low levels of androstenone in D and NL boars. Some of these genes have already been shown to affect levels of androstenone, although most of them have not previously been described in this context. Both CYB5 and CYP17 were up-regulated in both D high (DH) and NL high (NLH) animals and this was confirmed by real-time PCR. Over-expression of CYB5 cDNA clones in high androstenone animals agrees with results previously seen in protein expression studies [13, 14]. Earlier studies on CYP17 have shown that polymorphisms in the coding region of the gene are not associated with androstenone levels [15], and that levels of the CYP17 protein are not correlated to levels of 16-androstene steroids in fat [14]. Notably, up-regulation of the gene expression of CYP17 in this study may be explained by its direct interaction with CYB5 in the andien-β synthase enzyme system. Differential expression of both genes was confirmed by real-time PCR in the D boars of our study.

StAR is involved in the transport of cholesterol from the outer to the inner mitochondrial membrane, where CYP11A1 converts cholesterol to pregnenolone. This has been described as the rate-limiting step in steroid hormone synthesis, after which pregnenolone may exit the mitochondria and can be metabolised to steroid hormones [21]. The cholesterol side chain cleavage enzyme system contains ferredoxin (FDX1) which functions as an electron transport system, transferring electrons from NADPH-dependent ferredoxin reductase to CYP11A1 [22, 23]. StAR was over-expressed in NLH, while FDX1 and CYP11A1 were over-expressed in both the DH and NLH groups. Our rcPCR results show that both StAR and CYP11A1 are highly up-regulated in the DH and NLH groups, even though the microarray study was not able to detect differences in the gene expression of StAR in DH. Based on results from the rcPCR study, StAR seems to particularly affect the androstenone level since this gene is up-regulated as much as five times in DH pigs and twelve times in NLH pigs. Emopamil-binding protein (EBP), sterol-C4-methyl oxidase-like isoform 1 (SC4MOL) and the cytochrome P450 enzyme lanosterol 14-alpha-demetylase (CYP51) are enzymes involved in the cholesterol biosynthesis pathway [24–26]. EBP was up-regulated in NLH, while SC4MOL and CYP51 were up-regulated in the DH. Isopentenyl-diphosphate delta isomerase (IDI1) catalyses the inter-conversion of isopentenyl diphosphate (IPP) to dimethylallyl diphosphate (DMAPP), which is the substrate for the reaction that ultimately results in cholesterol, and this gene was up-regulated in both DH and NLH. The over-expression of all these genes may indicate that this early stage of steroid hormone biosynthesis is important for the production levels of androstenone. However, since cholesterol ultimately metabolises into both androstenone and testosterone, it might be impossible to use these genes as markers for low androstenone without simultaneously affecting levels of testosterone. However, further studies are needed to address this hypothesis.

In addition to FDX1, two additional ferritin cDNA clones were significantly up-regulated in the NLH and DH groups: FTL and the ferritin heavy-chain (FTH1). The differential expression of FTL was also confirmed by real-time PCR. Ferritin is an iron storage protein and plays a central role in numerous essential cellular functions [27]. Ferritin may also affect levels of androstenone through the CYB5/CYP450 electron transfer as the haem-containing cytochrome P450s need to receive electrons, e.g. from interaction with CYB5 [28], to be activated. 5-aminolevulinate synthase (ALAS1) is the first and rate-limiting enzyme in the haem biosynthesis pathway, providing haem for e.g. CYP450s [29]. This gene is up-regulated in DH. The gene ontology results also suggest an important role for ferric iron/iron ion binding for the levels of androstenone in adipose tissue.

Cytochrome P450 c19 (CYP19) encodes the enzyme aromatase, which catalyses the synthesis of estrogens from androgens. The pig secretes unusually high levels of estrogens from the testes [30] and is the only mammal known to express functionally distinct isoforms of this gene [31]. Our results show an up-regulation of CYP19A2 isoform in both NLH and DH boars, whereas the CYP19A1 isoform was up-regulated in DH. Estrogens are shown to be positively correlated with androstenone levels in fat [32] and estrone sulfate levels in plasma have been found positively correlated with both plasma and fat levels of androstenone [33]. SMPD1 is involved in the conversion of sphingomyelin to ceramide, which has been shown to inhibit CYP19 activity through induction of transcription factors [34]. SMPD1 was down-regulated in DH animals in the microarray study, but we were not able to confirm this in the rcPCR analyses.

HSD17B4 was found up-regulated for both DH and NLH, in both microarray and rcPCR analyses. The 17β-HSDs catalyse the last step in the formation of androgens and estrogens, and the HSD17B4 inactivates hormones very efficiently [35]. This has been reported for the substrates 17β-estradiol and 5-androstene-3β,17β-diol [36] but the enzyme also accepts other substrates [37]. The enzymes have previously been assigned to porcine Leydig and Sertoli cells [38] and porcine tissues express HSD17B4 as a predominant dehydrogenase [39]. Progesterone, a metabolite of pregnenolone, has been shown to increase porcine HSD17B4 activity [39] and the progesterone receptor membrane component 1 (PGRMC1) was over-expressed in DH. Another 17β-HSD, isoform 11, also called dehydrogenase/reductase family member 8 (DHRS8) was down-regulated in DH boars.

AKR1C4 belongs to the cytosolic aldo-keto reductases that act as 3α-/3β-/17β-/20α-hydroxysteroid dehydrogenases (HSDs) in human [40, 41] and was highly up-regulated in DH and NLH pigs. This was also confirmed by rcPCR. Penning et al. [36] showed that all the isoforms AKR1C1-AKR1C4 could interconvert active androgens and oestrogens with their associated inactive metabolites, which prevents an excess of circulating steroid hormones and makes the steroids into substrates for conjugation reactions [40].

Nuclear receptor subfamily 5, group A, member 1 (NR5A1), also called steroidogenic factor 1 (SF-1), was over-expressed in DH boars. This transcription factor is involved in the regulation of numerous genes, including StAR [42], CYP11A [43], CYP17 [44] and CYP19 [45]. Also, nuclear receptor co-activator 4 (NCOA4), an androgen receptor (AR) activator often referred to as ARA70, was up-regulated in the DH group. The nuclear receptor co-activator 2 (NCOA2), another AR activator, was significantly over-expressed in NLH. Other genes regulating transcription and translation that were differentially expressed in this study include class III POU transcription factor (POU3F2), microtubule-associated protein light chain 3 isoform A (MAP1LC3A) and EEF1A1 which were up-regulated in the DH and NLH groups, eukaryotic translation elongation factor 1 beta (EEF1B2) which was up-regulated in DH boars, and nuclear receptor subfamily 3, group C, member 2 (NR3C2) which was down-regulated in DH pigs.

Conjugation reactions contribute to the levels and pattern of steroids present in the plasma circulation of the boar. The enzymes involved in conjugation reactions were first thought to strictly inactivate and eliminate the compounds by rendering them more water-soluble. However, the high levels of conjugated steroids that are present in the plasma of the boar may suggest that the biological role of these reactions is more complex [46, 47]. A common family of conjugation enzymes is the sulfotransferases. SULT2A1, also called hydroxysteroid sulfotransferase (HST), was up-regulated in both DH and NLH boars in the microarray experiment and subsequently confirmed by rcPCR. SULT2A1 is a key sulfotransferase enzyme in terms of the 16-androstene steroids, and its activity in the testis is negatively correlated to androstenone concentrations in fat [19, 48]. Sulfotransferase family 2 isoform B1 (SULT2B1) is also a HST enzyme, showing selectivity for the conjugation of 3β-hydroxysteroids [49]. The SULT2B1 was significantly up-regulated in NLH pigs. The over-expression of the sulfotransferase genes in this study does not correspond to the negative correlation to levels of androstenone previously reported for this enzyme family [48].

Another conjugation reaction is catalysed by the glutathione S-transferases (GSTs). They have been shown to bind hormones [50] and influence their transport, metabolism and action [51]. One isoform of these enzymes (GSTalpha) is shown to be active in the pig testis Leydig and Sertoli cells and positively regulated by both estradiol and testosterone [52]. Glutathione S-transferase omega (GSTO1) and glutathione S-transferase (MGST1) were over-expressed in DH. These enzymes have not previously been studied in connection with androstenone; however, no correlation was found between fat skatole levels and glutathione S-transferase activity [53].

D and NL show many similar molecular functions and biological processes in this study, although they also differ with respect to which genes are differentially expressed. [54]. This miscellany might be due to the generally higher androstenone levels observed in D compared to NL [8]. The variability in the potential for androstenone production or elimination may to some extent be explained by breed differences in age at sexual maturity [55]. Differences between D and NL have also been shown by differential expression of 3β-hydroxysteroid dehydrogenase and SULT2B1 proteins, both enzymes critical for steroid production [54]. The higher number of significant genes in D compared to NL may be explained by the more extreme androstenone values in this breed (See Additional file 9: Androstenone values), giving higher contrasts.

In addition to the genes listed, we also identified several highly differentially expressed cDNA clones with no homology to known human or mouse sequences. Further characterisation of these genes may uncover new and unexpected roles in association with androstenone. The over-expression of for example AKR1C4, SULT2A1 and SULT2B1 in high androstenone animals does not explain a role for these genes as inactivation enzymes and the role of these genes needs to be further studied in relation to androstenone levels. Additional studies of all the genes are necessary to see if their proteins show the same differential expression. Furthermore, identification of single nucleotide polymorphisms (SNPs) in the genes or in association with the genes can be used for breeding purposes. The exact function of genes from interesting gene ontology pathways, like ferric iron binding, iron ion binding, electron transport activities and oxidoreductase activities, also needs to be clarified.

Conclusion

Our study detected differentially expressed genes that are previously found to affect androstenone levels in boars, as well as genes from pathways not formerly described in this aspect. We confirm the involvement of CYP17 and CYB5 and detect a number of other genes involved in the steroid hormone pathway that seem to be essential for androstenone levels. Besides SULT2A1 we identified other conjugation enzyme genes that might be important, including SULT2B1, AKR1C4, GSTO1, MGST1 and HSD17B4. The genetic profiles identified should be further examined to clarify their potential as molecular markers for reduced boar taint.

Methods

Animals

D and NL boars from NORSVIN's three boar testing stations were included in this study. The D and NL boars were on average 143 and 156 days at 100 kg live weight, respectively, and were slaughtered on average 14 days later. Tissue samples from testicle were frozen in liquid N₂ and stored at -80°C until used for RNA isolation as described below. Fat samples were collected from the neck and stored at -20°C until used for androstenone measurements. Androstenone levels in fat were analysed by a modified time-resolved fluoroimmunoassay [56] using antiserum produced at the Norwegian School of Veterinary Science (NVH) [57]. A total of 1533 NL boars and 1027 D boars were analysed and statistical power calculations showed that selecting animals from each tail of the androstenone distribution was an effective way of obtaining high probability of finding differentially expressed genes with limited number of arrays. The power calculations suggested 30 arrays to be sufficient. Hence, for each breed, the 30 most extreme boars from each tail of the distribution were selected from all of the phenotyped animals (See Additional file 9: Androstenone values). In order to reduce family effects a maximum of two and three half sibs were chosen from NL and D, respectively. The same animals were used for verification of selected genes by rcPCR. Only D pigs were included in the real time PCR analysis.

Expression profiling using microarrays

The microarrays (UltraGAPS coated slides, Corning Incorporated, MA, USA) were produced at the Faculty of Agricultural Sciences, Aarhus University and contained 27,774 cDNA clones printed in duplicates [58]. Out of these, 26,877 were PCR products amplified from cDNA clones and 867 were control DNA fragments. Of the 26.877 cDNAs, 21.417 map to 15.831 human gene transcript IDs corresponding to roughly 1.35 cDNAs per gene transcript. The remaining 5.460 cDNAs were thus estimated to cover around 4.036 gene transcripts yielding 19.867 gene transcripts in total. The annotation was done using NCBI Reference Sequence (RefSeq) database Release 11. The cDNA clones originated from the Sino-Danish sequencing project, covering 0.66× of the pig genome [59]. A more detailed description of the microarray can be found at the NCBI Gene Expression Omnibus (GEO [60]) database using the platform accession number GPL3608.

Total RNA for microarray studies and real-time PCR were extracted from testicle tissue using the RNeasy midi kit according to the manufacturer's instruction (Qiagen, CA, USA). Quantities were measured using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, DE, USA) and qualities were examined by the 28S:18S rRNA ratio using the RNA 6000 Nano LabChip^® Kit on 2100 Bioanalyzer (Agilent Technologies, CA, USA). Aminoallyl-cDNA was synthesised from 20 μg of total RNA using the SuperScript indirect cDNA labelling system (Invitrogen Corporation, CA, USA) and labelled using ARES Alexa Fluor 488 or 594 labelling kit (Molecular Probes, OR, USA). Half of the samples were labelled with one dye and the other half with the other dye facilitating a direct balanced block hybridisation design, where the dye swap is balanced between the samples. As the animals are selected from a very large number of boars and the contrast between the two groups is extremely large, animals within each of the experimental group were treated as equal. We have prioritised a larger number of animals instead of hybridising all animals twice (with a dye swap), as we expect the biological variation to be greater than the technical variation. It may be noted that our experimental design eliminates the dye bias in the contrast 'high' versus 'low' androstenone, since SUM('high expressions')/30 minus SUM('low expressions')/30 is (30*H + 15*R + 15*G)/30 - (30*H + 15*R + 15*G)/30 = H - L, where H (L) denotes effect of high (low) androstenone, R (G) denotes the effect of red (green) dye and 30 is the number of arrays. Spike-in RNA from the Lucidea Universal ScoreCard (Amersham Biosciences) was added to the cDNA reactions. "Green" spike-in RNA was added to the samples labelled with Alexa-594 and "red" spike-in RNA was added to the samples labelled with Alexa-488. Purification of the amino-modified and fluorescently labelled cDNA was done using the NucleoSpin 96 Extract II PCR Clean-up kit (Macherey-Nagel, Düren, Germany), and a hybridization blocker (Invitrogen Corporation, CA, USA) containing polydA (Invitrogen Corporation, CA, USA) and Yeast tRNA (Invitrogen Corporation, CA, USA) was added. Each microarray was hybridised with one high and one low androstenone sample from the same breed, giving a total of 30 arrays for each breed. The high and low samples were paired randomly within each breed. Hybridisation was conducted in a Discovery XT hybridisation station (Ventana Discovery Systems, AZ, USA), followed by a manual wash and drying by centrifugation. More detailed descriptions of the microarray experiments are available at the GEO database through the series accession number GSE 7409.

The microarrays were scanned using a ScanArray Express scanner (Perkin Elmer Inc., MA, USA). Signal intensities were quantified using the ScanArray Express software and further analysis carried out in R version 2.2.1 [61] software package Linear Models for Microarray Analysis (Limma version 2.7.2) [62–64]. Mean foreground intensities were background corrected using the Edwards method [65] implemented in Limma using the median background intensities, and log-ratios were printtip loess normalised within arrays. The duplicate correlation function in Limma was used to consider the duplicate printing of each feature. As a quality check, MA-plots (M = log₂594/log₂488, A = (log₂594 + log₂488)/2), image plots and box plots were constructed using the Limma tools for visualisation of the data before and after normalisation. For assessing differential expression, Limma uses linear models in combination with an empirical Bayes method to moderate the standard errors of the estimated log-fold changes. The nominal p-values were corrected for multiple testing by false discovery rates (FDR) using Benjamini and Hochberg approach [66] and adjusted p-values < 0.01 were considered significant.

Since the Limma statistic provides a parametric test, it might be affected by outlier records. Therefore, we also conducted a more robust, but less powerful non-parametric test, namely Fisher's Sign Test (FST) [67].

The features of the arrays were mapped to a LocusLink identifier and an annotation package was built using the Bioconductor package AnnBuilder (version 1.9.14). Tests for significantly overrepresentation of GO terms (p < 0.01 and more than 10 significant genes) were conducted using the GOHyperG function of the Bioconductor package GOstats (version 1.6.0) [68, 69].

Quantitative rcPCR analysis and quantitative real-time RT-PCR analyses

A real competitive (rc) PCR gene expression analysis was used to verify some of the results from the microarray study. The method is build upon the MassARRAY methodology using the Quantitative Gene Expression (QGE) iPLEX system (Sequenom, CA, USA). Total RNA was isolated from testes by using the automatic DNA/RNA extractor (BioRobot M48 workstation; Qiagen; CA, USA) and first strand cDNA synthesis was conducted on 0.5 μg total RNA using SuperScript™-II Rnase H^- Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA). Purified total RNA was treated with TURBO DNA-free™ (Ambion, Huntingdon, UK) for removal of contaminating DNA. Assay designs for the genes included in this study (See Additional file 10: Gene transcripts included in the rcPCR analyses) were multiplexed into a single reaction using MassARRAY QGE Assay Design software (Sequenom, CA, USA). The competitors, a synthetic DNA molecule that matched the sequence of the targeted cDNA region at all positions except for one single base served as an internal standard for each transcript. A 10-fold dilution of competitor was initially used over a wide range of concentrations to determine an approximate equivalence point, followed by a 7-fold dilution of competitor from 4.04 × 10^-11 to 1.43 × 10^-19 to achieve more accurate results. The cDNA and competitor were co-amplified in the same PCR reaction with the conditions 95°C for 15 minutes, 45 cycles each of 95°C for 20 second, 56°C for 30 seconds and 72°C for 1 minute, and finally 72°C for 3 minutes. After a clean-up step to remove unincorporated nucleotides, the PCR products served as templates for the primer extension reaction. The iPLEX reaction cocktail mix and PCR conditions were done as described in the Sequenom application guide [70]. Parallel PCR-reactions were performed for all samples and each of the products was printed with 2 replicates on a SpectroCHIP. The primer extension reaction generates distinct mass signals for competitor and cDNA-derived products, and mass spectrometric analysis generated signals from which the peak areas were calculated.

Because of collaborative reasons between two labs, another verification method used for the D boars was quantitative real-time PCR. Gene specific primers for four selected genes were generated from the Primer 3 software (See Additional file 11: Real-time PCR primers). Porcine ACTB and 18S rRNA amplifier set probes were used as endogenous control for normalisation. The gene ACTB was detected by Taqman probe whereas other genes were detected by SYBR Green probes. The same Duroc RNA samples as the ones used for the microarray experiments were applied. The RNA was synthesised into first-strand cDNA using SuperScript II Reverse Transcriptase (Invitrogen Corporation, CA, USA). The real-time PCR reaction was composed of 5 μL of Taqman master mix, 2 μL cDNA, 0.3 μL of forward and reverse primers (10 μM) and 0.05 μL of probe (10 μM). The real-time cycler conditions were 50°C for 2 minutes, 95°C for 10 minutes and 40 cycles each of 95°C for 15 seconds and annealing/extension at 60°C for 1 minute. Each reaction was conducted in triplicate on each individual sample with ACTB and 18S rRNA amplified as internal control genes. The real-time PCR amplification was performed using ABI PRISM 7900 HT sequence detection system. A cycle threshold value (C_T) was recorded for each sample and a standard curve made from 4×, 2×, 1×, 0.5×, 0.125× and 0.0625× was used to calculate the relative mRNA levels.

In the rcPCR study, the gene expression levels of H and L androstenone boars were analysed using the software TITAN version 1.0–13 [71] that runs in the R statistical environment. The raw data from the Genotype Analyzer Software (Sequenom) was imported into TITAN and analysed using the default values of linear least squares polynomial regression with 4000 bootstrap replicates. The cDNA concentrations were corrected with respect to the housekeeping gene (HPRT), and p-values and confidence intervals for the fold changes were calculated. Details about the TITAN software are available at web [72].

Quantification of the real-time PCR amplification was performed using ABI PRISM 7900 HT sequence detection system. The standard curve method was used to calculate the relative mRNA levels.