BACKGROUND
it is now feasible to carry out large scale characterization of transcription activity using microarrays.
this technology has opened new avenues to characterize and to dissect the transcriptome's genetic basis.
it is a complementary approach to the classical ascertainment of the genetic architecture of complex traits, such as quantitative trait loci studies.
there is now overwhelming evidence that the levels of mrna are affected by a number of environmental, physiological and genetic factors, much the same as for any other quantitative, complex trait  <cit> .
the extent and influence of each factor is, however, controversial and unknown to a large extent.
some authors reported that sex was far more important than the genetic line or age in drosophila  <cit> , whereas individual variability have been found to be very large in fundulus fish  <cit> .
the study of the genetic basis of variability at the transcriptome level is particularly relevant because, it has been claimed, evolution and thus phenotypic variability is due primarily to regulatory rather than structural mutations  <cit> .
in many microarray studies so far, the goal has been to compare two physiological statuses, e.g., disease vs. healthy, by analyzing a single tissue in several individuals.
some studies have also compared different breeds or strains, again focussing on one or very few tissues in order to gain insight as to how much phenotypic variability correlates with differences at the transcriptome level; see for instance reiner-benaim et al.
<cit> .
nevertheless, the number of studies comparing different breeds is not very large yet, particularly in livestock.
for instance, only four published studies have compared the global transcriptome of at least two porcine breeds  <cit>  and one compared two divergent lines  <cit> , whereas over a hundred dealing with porcine microarrays are indexed in pubmed <cit>  as of november  <dig>  thus, although a great deal is known about the transcriptome changes brought about by diseases, say cancer, much less is known about the relevance of genetic or environmental factors.
moreover, an organism is made up of many different tissues and organs, each with its own genetic programme.
how different are these genetic programmes?
how relevant are they in terms of total transcriptome variability?
what are the relative contributions of breed or sex when assessed across tissues?
divergent livestock breeds offer extremely valuable genetic material to study these issues.
in order to contribute to answering these questions, here we study the variability in a large number of tissues pertaining to a reduced number of individuals.
by allocating more experimental resources to a much larger number of tissues than usual, we aim at better characterizing the transcriptome variability across the whole organism.
although some studies have compared the transcriptome across different tissues  <cit> , the contribution of each effect  has not been quantified.
the variability between tissues is not as well studied as that within tissues and some tissues remain poorly characterized, all the more in livestock species.
the transcriptome of some tissues like liver and muscle have been analyzed extensively, well over a thousand citations in pubmed <cit>  as of november  <dig>  some other tissues have been only rarely studied despite their known physiological importance: only four microarray studies are reported for pineal gland,  <dig> for the hypophysis.
certainly, phenotypic differences in a given tissue can actually be due to changes that occur in a different organ.
for instance, the fatness differences between individuals, a trait of utmost relevance in animal breeding, are more likely to be caused by genetic signals that originate in the hypothalamus or in the pineal gland rather than in the fat tissue itself.
the aim of this work was to study the global transcripome across tissues in two highly divergent pig breeds.
we report a mixed model analysis of  <dig> tissues pertaining to four pigs, two large white  and two iberian , one male and one female per breed.
these two breeds are phenotypically extreme for most traits of economic relevance, e.g., growth, fatness, reproductive performance.
the tissues studied were olfactory bulb, hypothalamus, pineal gland, adenohypophysis, neurohypophysis, adrenal cortex, adrenal medulla, thyroid gland, diaphragm, biceps femoris muscle, back fat tissue, abdominal fat tissue, stomach, liver, ileum and whole blood.
RESULTS
transcriptome diversity across tissues and samples
clustering is a useful starting exploratory tool to visualize highly dimensional data, and has been widely used to microarray data since the seminal paper of eisen and cols.
<cit> .
here we applied two clustering methdos, the classical one based on the upgma criterion  <cit> , and the neighbor-joing  clustering.
in both cases, we used the distance one minus the squared correlation  between the samples, after normalizing the raw data with the rma procedure  <cit> , as detailed in material and methods.
results are drawn in figure  <dig>  where it can be seen that samples were clearly grouped by tissue, next by breed in both trees.
this was neatly observed for ileum, liver, thyroid gland, adeno and neurohypohysis and olfactory bulb.
muscle samples were clustered by tissue  but less clearly within each muscle.
as for fat, the similarity was larger between tissues than between breeds, and samples of both back and abdominal fat origins were clustered together.
the same was observed between cortex and medulla from adrenal gland.
in this case, contamination between both tissues cannot be ruled out because of the irregular limits of the medulla that make not easy to separate that region neatly from the cortex collecting rapidly enough amount of tissue for analysis.
other authors have described previously contamination of medulla in the cortex sample when mechanical separation is performed <cit> .
thus, this resemblance was not completely unexpected.
the only outlier sample seemed to be the pineal gland of the large white male , which clustered with the rest of hypothalamus microarrays.
here contamination can be discarded in all likelihood because the two regions, hypothalamus and pineal gland, are in distinct areas of the brain.
however, the pineal gland works in harmony with the hypothalamus.
the former produces melatonin, which directly influences the function of various brain centers, including the hypothalamus.
in stomach, less clearly in blood, samples were grouped by sex rather than by breed.
as expected, both upgma and nj methods provided identical clustering at the first coalescence level.
however, there were some interesting differences at higher levels.
nj identified four groups of tissues.
the first group comprised brain tissues, including hypophysis; the second group, thyroid and adrenal endocrine glands; the third, muscle and fat tissues; whereas stomach, ileum, blood and liver were in the last group.
in contrast, upgma clustered blood, liver and muscles in separate groups, while brain tissues were grouped together with fat, thyroid and adrenal glands; blood and liver were in distinct groups.
it is well known that nj have better properties than upgma when reconstructing evolutionary trees because it does not assume a constant evolutionary rate  <cit> .
but is not evident how to translate this advantage in relation to the performance for classifying transcriptomes.
the upgma method is equivalent to clustering algorithms implemented in popular microarray packages  <cit> , whereas nj have not been used so far, to our knowledge, in microarray studies.
here we found that, overall, nj provided more biologically sound results than upgma: the first two groups in the nj-tree comprised brain tissues and endocrine glands; muscle and fat tissues have a predominantly structural function in the body, while the rest of tissues  have an important metabolic role.
in contrast, upgma results were not always so easily interpretable, e.g., it grouped hypophysis closer to fat tissues than to other brain organs like hypothalamus or the olfactory bulb.
the trees depicted in figure  <dig> were obtained with the original data .
a better alternative to study relationships between tissues is to consider the probe × tissue solutions obtained from model , as these solutions are corrected for 'noisy' factors such as sex, breed and global tissue effects.
thus we constructed a nj-tree of tissue transcriptomes using the 1-rt <dig> distance .
this tree is shown in figure  <dig>  although the pattern was similar to the nj-tree in figure  <dig>  this tree provided a clearer picture of relationships between tissues.
for instance, it emphasized that blood is the most distant tissue, corroborated also by the fact that blood harbored the largest number of extreme probes .
we defined a tissue's extreme probe as a probe for which all four mrna levels of the tissue were either smaller or larger than the remaining mrna levels .
in contrast, the two adrenal tissues and both fats were the closest pair of tissues.
the brain is not a uniform organ, and this well known fact  <cit>  was corroborated by a relevant variability within the different brain tissues sampled.
olfactory bulb was neighbor to hypothalamus.
the two pituitaries, adenohypohysis and neurohypohysis, exhibited neighbor but distinct transcriptomes.
this result agrees with the fact that both organs carry out very different physiological functions and have separate embryonic origins from the ectoderm layer .
in fact, it was more interesting to note that neurohypophysis and hypothalamus were relatively distant.
a close relationship might have been expected because neurohypophysis consists primarily of axonic projections of the hypothalamus.
the pineal gland exhibited an intermediate transcriptome between hypophyses and hypothalamus – olfactory bulb.
a the complete list is in supplementary table s <dig>  see material and methods for definition of extreme probe.
interestingly, there was a relation between embryonic origin and clustering.
most tissues with the same embryonic origin clustered together in the nj-tree .
all brain tissues are of ectodermic origin, as is hypophysis, and these tissues were clustered together in the first group of tissues.
the second group comprised the adrenal and thyroid glands, tissues with double embryonic origin.
the adrenal medulla derives from ectoderm , while the cortex develops from mesoderm.
the thyroid gland develops from two distinct embryonic lineages: follicular cells  and parafollicular c-cells  and are of endodermal and ectodermal neural crest origins, respectively  <cit> .
mesoderm gives origin to fat and muscle.
liver, and stomach and ileum epithelia are derived from endoderm.
finally, blood, the most separate and outlier tissue, is of mesoderm origin.
thus, it seems that embryonic development leaves a standing footprint on the transcriptome in each of the tissues.
functional analysis of extreme genes across embryo origins and tissues
we sought to investigate more in detail the genetic basis of the tissue arrangement by embryonic layers, and to pick up genes that can be differentially expressed in concerted action within each embryonic layer.
to that end, we obtained the extreme gene probes differentiating the ectoderm tissues , mesoderm , and endoderm .
an extreme probe for each embryo layer was defined as for individual tissues, i.e., a probe whose all expression levels for that group of tissues were either smaller or larger than the levels for the rest of tissues .
we excluded thyroid and adrenal glands for being mixed tissues and blood, for outlier.
the complete list in additional file  <dig>  a pie-plot with differential gene ontology  frequencies obtained with onto-express  <cit>  is in figure  <dig>  a wide variety of go were over represented in each layer; this fact, together with a large percentage of unknown or unannotated genes makes it somewhat difficult to interpret the results in detail.
neverheless, some interesting results appeared.
for instance, genes with ontologies like nervous development and ion transport – clearly related to central nervous system – were more frequent than expected within ectoderm extreme genes, as were genes involved in development .
in particular, some of the most overexpressed genes like ina , astn  or napil <dig>  are involved in neuron development, whereas others like gabrg <dig>  or syt <dig>  are involved in normal neuron functioning.
the case of nap1l <dig>  is particularly interesting.
although, according to pubmed, 'the function of this family member is unknown, mouse studies suggest that it represents a class of tissue-specific factors interacting with chromatin to regulate neuronal cell proliferation'.
our results strongly suggest that nap1l <dig> is involved in neuronal system and that the function is maintained across species.
in parallel to results for ecctoderm tissues, the most significant enriched ontology was muscle development in the mesoderm .
here, some of the most extreme genes were hox genes involved in limb development and myogenesis .
in particular, it has been suggested that prrx <dig> and hoxd <dig> are also involved in the urogenital tract development, also of mesoderm origin.
the two top extreme genes for mesoderm tissues gamma-sarcoglycans  are also involved in muscle development.
as for the endoderm , an excess of genes involved in general metabolism was observed, although we also found a significant enrichment in nervous development genes .
two of the most overexpressed extreme genes, hnf4a  and foxa <dig>  are known to be involved in liver development, and it is suspected that hnf4a can play a role also in intestine development, which is supported by our data.
intriguingly, some of the most extreme overexpressed genes were involved in the complement pathway , may be as a consequence of the defense function of the ileum.
finally, when all extreme genes across the three embryo layers were considered jointly , the most significant ontologies were nervous system development, cell differentiation, cell adhesion and multicellular organism development, which again suggests a sound biological interpretation.
we also obtained the extreme probes for each of the individual tissues.
the main over represented biological process per tissue is in table  <dig>  and a complete list of the genes is in additional file  <dig>  results confirmed, overall, previous biological knowledge.
brain tissues  exhibited a large number of genes in signalling pathways .
the highly oxidative muscle diaphragma was enriched in genes involved in atp synthesis through citrate cycle and oxidative phosphorylation, whereas biceps femori  extreme genes were often related to the production of atp by glycolysis and sugar pathway.
the liver extreme genes were in the urea cycle pathway and bile acid and amino acid metabolisms; blood genes pertained to hematopoietic processes and signaling pathway systems.
mixed model analysis
mixed model analyses were highly successful in modelling transcriptome variability .
for instance, for model  the fraction of unexplained variance was very small,  <dig> - hp <dig> - hpt <dig> =  <dig> , and ~15% in the rest of models.
this means that the set of expression levels can be very accurately modelled as a random normal variate in the log scale.
note that models  were highly parsimonious: a few parameters, i.e., the variance components and the fixed effects tissue, sex and breed were needed to model the data.
by far, the largest variance component was that of the probe, which explained overall at least 85% of the total variability .
the interaction probe × tissue  accounted for 11% of phenotypic variance , whereas the total variances explained by probe × breed or probe × sex were marginal, only about 1% .
these were the global variance estimates, when all probes were considered jointly.
it does not follow that sex or breed were not associated to changes in gene expression.
the relative importance of sex or breed did increase when we restricted the analysis to a data subset.
for instance, the percentage of total variance explained by breed or sex increased to 25% or 48% when we used the  <dig> most differentially expressed genes between breeds or between sexes, respectively .
symmetrically, the relevance of probe effect decreased for these data subsets.
differential gene expression
at a false discovery rate fdr <  <dig> , we identified  <dig> and  <dig> genes differentially expressed between sexes and breeds, respectively.
the complete lists together with tentative annotations  <cit>  are in tables  <dig> and  <dig>  two of the genes, orthologous to human genes in the hla-a and hla-b complex, were differentially expressed simultaneously between breeds and between sexes.
however, there was no correlation overall between tissue, breed or sex effects.
this was evident from the pattern in figure  <dig>  which displays the relation between zb and zs scores  and between zb scores and the standard deviation of the probe × tissue solutions .
if there were a relation, we would see the dots around a diagonal rather than close to the axes.
thus, each of the factors studied here, sex and breed, influence the transcriptome through different genetic programmes.
z – scores correspond to male minus female contrasts.
z – scores correspond to large white minus iberian contrasts.
we performed a clustering of probes and samples using only the differentially expressed genes from tables  <dig> and  <dig>  results were drawn in figure  <dig>  where the discrimination between sexes  and breeds  is neat and clearly visible.
this clustering contrasts with that in figure  <dig> where all probes were employed.
the discrepancy occurs because the expression pattern of most genes is primarily affected by the tissue and not so much by sex or breed, as discussed.
however, as figure  <dig> clearly shows, this does not mean that a careful selection of probes does not allow us to discriminate samples between according to sex or breed.
the list of differentially expressed genes between sexes comprised eight x or y linked genes, based on the localization of orthologous human genes .
note that three x-linked genes were up-regulated in males , a phenomenon also observed by yang et al.
<cit>  although for different genes.
male up-regulated genes were predominantly involved in spermatogenesis , cell proliferation, migration, and differentiation , as well as genes genes from the histocompatiblity complex .
female up-regulated genes had functions related to positive regulation of the transforming growth factor beta receptor signaling pathway , proteins up-regulated in gonadotropinomas  or affecting the circadian rhytm .
up-regulation of clock in females is consistent with evidence showing that this gene is related to the cyclic nature of the reproductive process in females, both in mammal and in non mammal species  <cit> .
we found a good concordance with the list of differentially expressed genes obtained with the same microarray in lung and mesenteric lymph node by our group in a separate experiment .
thirteen genes  and ten genes  were also in table  <dig>  despite the fact that corresponded to different tissues, treatments and breeds.
yang et al.
<cit>  compared microarrays of  <dig> female and  <dig> male mice in four tissues , the largest study to date of sex differentially expressed genes.
they reported  <dig> genes that were differentially expressed in all four tissues, four of them  were also found here; ddx3y and eif2s3y were the most differentially expressed genes found here and in  <cit> .
it is interesting to note that these authors also reported that immune genes, e.g., histocompatiblity genes, were overexpressed in males with respect to females, in agreement again with our results  <cit> .
differential expression between sexes for sscy genes can be considered as validated.
we also validated probe , annotated as clock gene  <cit> , by real time quantitative rt-pcr  in backfat samples of  <dig> pigs.
qrt-pcr results clearly confirmed the microarray data .
the average relative expression was  <dig>  ±  <dig>  in females, whereas we were not able to reliably detect expression in males at 1/ <dig> cdna dilution because we were below the limit of detection.
to verify whether this gene product is actually expressed in males, we repeated the assay at 1/ <dig> cdna dilution.
we could detect low levels of expression in four out of the  <dig> males studied.
iberian male  <dig>  that with largest relative quantification , exhibited nonetheless a expression level ~ <dig> fold lower than in females.
sex: f, female; m, male
rq : relative quantification rq, average of three replicates each, at 1/ <dig> cdna dilution.
the results are given in terms of sample number  <dig>  used as calibrator, i.e., rq = 1
rq : relative quantification rq, average of three replicates each, at 1/ <dig> cdna dilution, only analyzed for males
bld, below limit of detection
na, not available
the list of probe sets corresponding to differentially expressed genes at fdr <  <dig>  between lw and ib breeds is in table  <dig>  the thirty three probes corresponded to  <dig> genes and  <dig> hypothetical proteins;  <dig> were over expressed in large white  and  <dig> over expressed in iberian .
the lw-biased list comprised genes involved in atp biosynthetic process, nervous system development, lipid metabolism, rna processing, and control of growing and cellular division.
the ib-biased genes affected the regulation of cell growth, rna splicing factor activity or muscle development.
these results hint that transcriptome differences between breeds seem to be manifold.
no particular single biological process was predominantly affected.
it is to be noticed, nonetheless, that there were several genes involved in nervous development  or growth and cell cycle .
as there is currently very little information on the transcriptome differences between pig breeds, it is not possible to know how general these results will be.
some recent results have compared muscle in duroc vs. taoyuan breeds  <cit> , in landrace vs. tongcheng  <cit>  and liver and muscle in duroc vs. pietrain  <cit> .
while these studies have focussed on a particular tissue and with different techniques and microarrays, all suggest that breeds differ in genes related to a variety of functional categories.
in some cases, the microarray included predominantly genes involved in muscle structure and development or energy metabolism  <cit>  so their results may not be representative of the whole genome.
it should be noted that the differentially expressed genes reported here  were those that exhibited largest differential expression across all tissues.
this follows from the specification of models , where all probes were analyzed simultaneously and where the tissue effect was included in the models as just an additional effect.
thus, it is important to bear in mind that other genes could show larger differential expressions in a specific tissue than those listed in tables  <dig> and  <dig>  the experimental design of this study did not allow us to characterize with reasonable fdr the specific expression of a gene in a single tissue.
nevertheless, we can get a glimpse of the effect of tissue in differential expression by subdividing the tissues in groups.
we reasoned that the probe × tissue interaction should be larger between tissues that clustered far apart, as their genetic programmes are more different that tissues that are very similar, e.g., back and abdominal fats.
we defined the following tissue groups base on the nj-tree shown in figure 2: brain , endocrine , structural  and metabolic tissues .
to a large extent, these groups reflect also a common embryonic origin .
we reanalyzed each of the groups with models  <dig> and  <dig>  obtaining new zb and zs scores for each of the tissue groups.
table  <dig> presents the correlation between the z-scores obtained with all tissues analyzed simultaneously and each of the groups.
in order to focus on the most relevant genes, the correlations shown were those obtained with the  <dig> genes with largest absolute z-scores.
the pattern shown was highly illuminating.
for sex, there seemed to be little interaction, as zs-scores were highly correlated across groups of tissues.
correlation coefficients were always >  <dig> .
in contrast, correlations between zb scores were much more variable and, importantly, much lower overall.
thus, whereas the correlation between brain and metabolic tissues was  <dig>  for sex z-scores, it was only  <dig>  for breed z-scores.
this means that most differentially expressed genes between sexes were shared across tissues whereas this was much more unlikely when comparing two distant breeds like iberian and large white.
it is interesting to note that the clusters in figure  <dig> also hinted that heterogeneity was larger within breed differentially expressed genes than within sex differentially expressed genes.
note that the color patterns of sex differentially expressed genes were more uniform across samples than for breed differentially expressed genes.
besides, clustering was much 'flatter' – i.e., more uniform – within sexes than within breeds.
this result may have important implications.
it suggests that physiological changes responsible for breed differences have targeted differentially the transcriptome across tissues.
not all tissues have been equally affected.
it remains to be studied on which of the tissues the effect of breed differentiation has been the largest.
our current data set does not allow us to respond to this question accurately: the probe × tissue heritabilities  were very similar in all tissue groups and no clear pattern emerged from figure  <dig>
the tissues included in each group are listed in table  <dig>  the z-scores correspond to the  <dig> largest differentially expressed genes when all tissues are jointly analyzed.
upper diagonal, sex differentially expressed genes.
lower diagonal, breed differentially expressed genes.
functional annotation of sex and breed differentially expressed genes
we carried a go automatic annotation  <cit>  with the  <dig> most differentially expressed genes between sexes and  <dig> between breeds.
the corresponding plots are in figure  <dig>  in neither case, sex or breed, was a given biological process go clearly over represented .
as a result, the discussion should be considered as tentative or provisional.
nevertheless, we found defense genes to be over represented for both sex and breed specific genes.
for sex, there was a significant excess of transcription and translation related genes.
some of these genes are sex linked , as mentioned above, so an over representation of this ontology is not unexpected.
the go biological processes were more scattered for breed than for sex .
interestingly, there were more genes involved in nervous system development and cell differentiation than expected among breed specific genes.
this might provide some clues as to what are the primary changes exerted by selection and breeding in the organism' transcriptome but further work is needed to get a definitive answer.
global discussion and
CONCLUSIONS
this study shows that a much larger fraction of transcriptome variability is due to tissue differences rather than to sex or breed.
our results agree with recent literature  <cit>  that showed that transcriptional co-expression is indicative of gene function and that tissues with similar function tend to cluster together.
an additional relevant observation from the nj-clustering  was that embryonic development seems to leave an enduring footprint on the transcriptome: tissues with similar developmental origin tend to cluster closer than those with different embryonic origins.
as clustering was carried out with all probes simultaneously, this means that the footprint extends to the majority of genes.
we further characterized the go biological process of the extreme genes between tissues or groups of tissues sharing a common embryonic origin .
we found that, often, a sound biological interpretation was possible.
for instance, when all extreme genes across the three embryo layers were considered jointly , the most frequent and significant ontologies were related to development.
we also indentified several genes involved in tissue development that have a distinct expression pattern in tissues with a common embryonic origin.
this distinct expression pattern is maintained beyond embryo development.
a key element of our analysis was mixed model methods.
although well known to statisticians and animal breeders  <cit> , these techniques are not so widely employed to dissect microarray data.
nevertheless, its use have been advocated and successfully employed in a number of studies, e.g.,  <cit> .
we report that mixed modeling explained a large part of total variability, with a very parsimonious parameterization.
a relevant advantage over standard methods is that all data are analyzed simultaneously as opposed to more utilized methods like anova – based techniques or bioconductor's package limma  <cit> ; as a consequence, reduced standard errors are expected.
in a second step, one can use the solutions from the models to carry out clustering or differential expression studies, as we did here.
although the variance explained by the effect of tissue was much larger  than that explained by sex or breed, which were comparable , we were able to identify genes with differential expression between breeds and between sexes.
noticeably, different genes were affected by each factor, i.e., the genes primarily affected by sex were different than those affected by breed or tissue.
are genes whose expression differ between breeds differentially expressed across all tissues or is there any interaction tissue × gene?
our experimental design does not allow to resolve this issue with enough power.
however, the results  hint that interaction can be important for differentially expressed genes between breeds but not so much between sexes.
resolving this question is of the utmost interest, both to understand the profound influences of artificial selection on the organisms and to propose novel more efficient strategies for animal breeding.
although more work is needed to calibrate the actual relevance of breed or genetic differences to the pig's transcriptome , it seems difficult that the importance of breed or sex differences exceeds that of tissue.
thus, the argument that regulatory rather than structural mutations are a more important source of phenotypic variability  <cit>  needs then to be considered in perspective.
it follows from our study that differences within the same individual largely exceeds those between individuals.

 it is now feasible to carry out large scale characterization of transcription activity using microarrays. this technology has opened new avenues to characterize and to dissect the transcriptome's genetic basis. it is a complementary approach to the classical ascertainment of the genetic architecture of complex traits, such as quantitative trait loci studies. there is now overwhelming evidence that the levels of mrna are affected by a number of environmental, physiological and genetic factors, much the same as for any other quantitative, complex trait. the extent and influence of each factor is, however, controversial and unknown to a large extent. some authors reported that sex was far more important than the genetic line or age in drosophila, whereas individual variability have been found to be very large in fundulus fish. the study of the genetic basis of variability at the transcriptome level is particularly relevant because, it has been claimed, evolution and thus phenotypic variability is due primarily to regulatory rather than structural mutations. in many microarray studies so far, the goal has been to compare two physiological statuses, e.g., disease vs. healthy, by analyzing a single tissue in several individuals. some studies have also compared different breeds or strains, again focussing on one or very few tissues in order to gain insight as to how much phenotypic variability correlates with differences at the transcriptome level; see for instance reiner-benaim et al.. nevertheless, the number of studies comparing different breeds is not very large yet, particularly in livestock. for instance, only four published studies have compared the global transcriptome of at least two porcine breeds and one compared two divergent lines, whereas over a hundred dealing with porcine microarrays are indexed in pubmed as of november thus, although a great deal is known about the transcriptome changes brought about by diseases, say cancer, much less is known about the relevance of genetic or environmental factors. moreover, an organism is made up of many different tissues and organs, each with its own genetic programme. how different are these genetic programmes? how relevant are they in terms of total transcriptome variability? what are the relative contributions of breed or sex when assessed across tissues? divergent livestock breeds offer extremely valuable genetic material to study these issues. in order to contribute to answering these questions, here we study the variability in a large number of tissues pertaining to a reduced number of individuals. by allocating more experimental resources to a much larger number of tissues than usual, we aim at better characterizing the transcriptome variability across the whole organism. although some studies have compared the transcriptome across different tissues, the contribution of each effect has not been quantified. the variability between tissues is not as well studied as that within tissues and some tissues remain poorly characterized, all the more in livestock species. the transcriptome of some tissues like liver and muscle have been analyzed extensively, well over a thousand citations in pubmed as of november some other tissues have been only rarely studied despite their known physiological importance: only four microarray studies are reported for pineal gland, for the hypophysis. certainly, phenotypic differences in a given tissue can actually be due to changes that occur in a different organ. for instance, the fatness differences between individuals, a trait of utmost relevance in animal breeding, are more likely to be caused by genetic signals that originate in the hypothalamus or in the pineal gland rather than in the fat tissue itself. the aim of this work was to study the global transcripome across tissues in two highly divergent pig breeds. we report a mixed model analysis of tissues pertaining to four pigs, two large white and two iberian, one male and one female per breed. these two breeds are phenotypically extreme for most traits of economic relevance, e.g., growth, fatness, reproductive performance. the tissues studied were olfactory bulb, hypothalamus, pineal gland, adenohypophysis, neurohypophysis, adrenal cortex, adrenal medulla, thyroid gland, diaphragm, biceps femoris muscle, back fat tissue, abdominal fat tissue, stomach, liver, ileum and whole blood.  transcriptome diversity across tissues and samples clustering is a useful starting exploratory tool to visualize highly dimensional data, and has been widely used to microarray data since the seminal paper of eisen and cols.. here we applied two clustering methdos, the classical one based on the upgma criterion, and the neighbor-joing clustering. in both cases, we used the distance one minus the squared correlation between the samples, after normalizing the raw data with the rma procedure, as detailed in material and methods. results are drawn in figure where it can be seen that samples were clearly grouped by tissue, next by breed in both trees. this was neatly observed for ileum, liver, thyroid gland, adeno and neurohypohysis and olfactory bulb. muscle samples were clustered by tissue but less clearly within each muscle. as for fat, the similarity was larger between tissues than between breeds, and samples of both back and abdominal fat origins were clustered together. the same was observed between cortex and medulla from adrenal gland. in this case, contamination between both tissues cannot be ruled out because of the irregular limits of the medulla that make not easy to separate that region neatly from the cortex collecting rapidly enough amount of tissue for analysis. other authors have described previously contamination of medulla in the cortex sample when mechanical separation is performed. thus, this resemblance was not completely unexpected. the only outlier sample seemed to be the pineal gland of the large white male, which clustered with the rest of hypothalamus microarrays. here contamination can be discarded in all likelihood because the two regions, hypothalamus and pineal gland, are in distinct areas of the brain. however, the pineal gland works in harmony with the hypothalamus. the former produces melatonin, which directly influences the function of various brain centers, including the hypothalamus. in stomach, less clearly in blood, samples were grouped by sex rather than by breed. as expected, both upgma and nj methods provided identical clustering at the first coalescence level. however, there were some interesting differences at higher levels. nj identified four groups of tissues. the first group comprised brain tissues, including hypophysis; the second group, thyroid and adrenal endocrine glands; the third, muscle and fat tissues; whereas stomach, ileum, blood and liver were in the last group. in contrast, upgma clustered blood, liver and muscles in separate groups, while brain tissues were grouped together with fat, thyroid and adrenal glands; blood and liver were in distinct groups. it is well known that nj have better properties than upgma when reconstructing evolutionary trees because it does not assume a constant evolutionary rate. but is not evident how to translate this advantage in relation to the performance for classifying transcriptomes. the upgma method is equivalent to clustering algorithms implemented in popular microarray packages, whereas nj have not been used so far, to our knowledge, in microarray studies. here we found that, overall, nj provided more biologically sound results than upgma: the first two groups in the nj-tree comprised brain tissues and endocrine glands; muscle and fat tissues have a predominantly structural function in the body, while the rest of tissues have an important metabolic role. in contrast, upgma results were not always so easily interpretable, e.g., it grouped hypophysis closer to fat tissues than to other brain organs like hypothalamus or the olfactory bulb. the trees depicted in figure were obtained with the original data. a better alternative to study relationships between tissues is to consider the probe × tissue solutions obtained from model, as these solutions are corrected for'noisy' factors such as sex, breed and global tissue effects. thus we constructed a nj-tree of tissue transcriptomes using the 1-rt distance. this tree is shown in figure although the pattern was similar to the nj-tree in figure this tree provided a clearer picture of relationships between tissues. for instance, it emphasized that blood is the most distant tissue, corroborated also by the fact that blood harbored the largest number of extreme probes. we defined a tissue's extreme probe as a probe for which all four mrna levels of the tissue were either smaller or larger than the remaining mrna levels. in contrast, the two adrenal tissues and both fats were the closest pair of tissues. the brain is not a uniform organ, and this well known fact was corroborated by a relevant variability within the different brain tissues sampled. olfactory bulb was neighbor to hypothalamus. the two pituitaries, adenohypohysis and neurohypohysis, exhibited neighbor but distinct transcriptomes. this result agrees with the fact that both organs carry out very different physiological functions and have separate embryonic origins from the ectoderm layer. in fact, it was more interesting to note that neurohypophysis and hypothalamus were relatively distant. a close relationship might have been expected because neurohypophysis consists primarily of axonic projections of the hypothalamus. the pineal gland exhibited an intermediate transcriptome between hypophyses and hypothalamus – olfactory bulb. a the complete list is in supplementary table s see material and methods for definition of extreme probe. interestingly, there was a relation between embryonic origin and clustering. most tissues with the same embryonic origin clustered together in the nj-tree. all brain tissues are of ectodermic origin, as is hypophysis, and these tissues were clustered together in the first group of tissues. the second group comprised the adrenal and thyroid glands, tissues with double embryonic origin. the adrenal medulla derives from ectoderm, while the cortex develops from mesoderm. the thyroid gland develops from two distinct embryonic lineages: follicular cells and parafollicular c-cells and are of endodermal and ectodermal neural crest origins, respectively. mesoderm gives origin to fat and muscle. liver, and stomach and ileum epithelia are derived from endoderm. finally, blood, the most separate and outlier tissue, is of mesoderm origin. thus, it seems that embryonic development leaves a standing footprint on the transcriptome in each of the tissues. functional analysis of extreme genes across embryo origins and tissues we sought to investigate more in detail the genetic basis of the tissue arrangement by embryonic layers, and to pick up genes that can be differentially expressed in concerted action within each embryonic layer. to that end, we obtained the extreme gene probes differentiating the ectoderm tissues, mesoderm, and endoderm. an extreme probe for each embryo layer was defined as for individual tissues, i.e., a probe whose all expression levels for that group of tissues were either smaller or larger than the levels for the rest of tissues. we excluded thyroid and adrenal glands for being mixed tissues and blood, for outlier. the complete list in additional file a pie-plot with differential gene ontology frequencies obtained with onto-express is in figure a wide variety of go were over represented in each layer; this fact, together with a large percentage of unknown or unannotated genes makes it somewhat difficult to interpret the results in detail. neverheless, some interesting results appeared. for instance, genes with ontologies like nervous development and ion transport – clearly related to central nervous system – were more frequent than expected within ectoderm extreme genes, as were genes involved in development. in particular, some of the most overexpressed genes like ina, astn or napil are involved in neuron development, whereas others like gabrg or syt are involved in normal neuron functioning. the case of nap1l is particularly interesting. although, according to pubmed,'the function of this family member is unknown, mouse studies suggest that it represents a class of tissue-specific factors interacting with chromatin to regulate neuronal cell proliferation'. our results strongly suggest that nap1l is involved in neuronal system and that the function is maintained across species. in parallel to results for ecctoderm tissues, the most significant enriched ontology was muscle development in the mesoderm. here, some of the most extreme genes were hox genes involved in limb development and myogenesis. in particular, it has been suggested that prrx and hoxd are also involved in the urogenital tract development, also of mesoderm origin. the two top extreme genes for mesoderm tissues gamma-sarcoglycans are also involved in muscle development. as for the endoderm, an excess of genes involved in general metabolism was observed, although we also found a significant enrichment in nervous development genes. two of the most overexpressed extreme genes, hnf4a and foxa are known to be involved in liver development, and it is suspected that hnf4a can play a role also in intestine development, which is supported by our data. intriguingly, some of the most extreme overexpressed genes were involved in the complement pathway, may be as a consequence of the defense function of the ileum. finally, when all extreme genes across the three embryo layers were considered jointly, the most significant ontologies were nervous system development, cell differentiation, cell adhesion and multicellular organism development, which again suggests a sound biological interpretation. we also obtained the extreme probes for each of the individual tissues. the main over represented biological process per tissue is in table and a complete list of the genes is in additional file results confirmed, overall, previous biological knowledge. brain tissues exhibited a large number of genes in signalling pathways. the highly oxidative muscle diaphragma was enriched in genes involved in atp synthesis through citrate cycle and oxidative phosphorylation, whereas biceps femori extreme genes were often related to the production of atp by glycolysis and sugar pathway. the liver extreme genes were in the urea cycle pathway and bile acid and amino acid metabolisms; blood genes pertained to hematopoietic processes and signaling pathway systems. mixed model analysis mixed model analyses were highly successful in modelling transcriptome variability. for instance, for model the fraction of unexplained variance was very small, - hp - hpt =, and ~15% in the rest of models. this means that the set of expression levels can be very accurately modelled as a random normal variate in the log scale. note that models were highly parsimonious: a few parameters, i.e., the variance components and the fixed effects tissue, sex and breed were needed to model the data. by far, the largest variance component was that of the probe, which explained overall at least 85% of the total variability. the interaction probe × tissue accounted for 11% of phenotypic variance, whereas the total variances explained by probe × breed or probe × sex were marginal, only about 1%. these were the global variance estimates, when all probes were considered jointly. it does not follow that sex or breed were not associated to changes in gene expression. the relative importance of sex or breed did increase when we restricted the analysis to a data subset. for instance, the percentage of total variance explained by breed or sex increased to 25% or 48% when we used the most differentially expressed genes between breeds or between sexes, respectively. symmetrically, the relevance of probe effect decreased for these data subsets. differential gene expression at a false discovery rate fdr <, we identified and genes differentially expressed between sexes and breeds, respectively. the complete lists together with tentative annotations are in tables and two of the genes, orthologous to human genes in the hla-a and hla-b complex, were differentially expressed simultaneously between breeds and between sexes. however, there was no correlation overall between tissue, breed or sex effects. this was evident from the pattern in figure which displays the relation between zb and zs scores and between zb scores and the standard deviation of the probe × tissue solutions. if there were a relation, we would see the dots around a diagonal rather than close to the axes. thus, each of the factors studied here, sex and breed, influence the transcriptome through different genetic programmes. z – scores correspond to male minus female contrasts. z – scores correspond to large white minus iberian contrasts. we performed a clustering of probes and samples using only the differentially expressed genes from tables and results were drawn in figure where the discrimination between sexes and breeds is neat and clearly visible. this clustering contrasts with that in figure where all probes were employed. the discrepancy occurs because the expression pattern of most genes is primarily affected by the tissue and not so much by sex or breed, as discussed. however, as figure clearly shows, this does not mean that a careful selection of probes does not allow us to discriminate samples between according to sex or breed. the list of differentially expressed genes between sexes comprised eight x or y linked genes, based on the localization of orthologous human genes. note that three x-linked genes were up-regulated in males, a phenomenon also observed by yang et al. although for different genes. male up-regulated genes were predominantly involved in spermatogenesis, cell proliferation, migration, and differentiation, as well as genes genes from the histocompatiblity complex. female up-regulated genes had functions related to positive regulation of the transforming growth factor beta receptor signaling pathway, proteins up-regulated in gonadotropinomas or affecting the circadian rhytm. up-regulation of clock in females is consistent with evidence showing that this gene is related to the cyclic nature of the reproductive process in females, both in mammal and in non mammal species. we found a good concordance with the list of differentially expressed genes obtained with the same microarray in lung and mesenteric lymph node by our group in a separate experiment. thirteen genes and ten genes were also in table despite the fact that corresponded to different tissues, treatments and breeds. yang et al. compared microarrays of female and male mice in four tissues, the largest study to date of sex differentially expressed genes. they reported genes that were differentially expressed in all four tissues, four of them were also found here; ddx3y and eif2s3y were the most differentially expressed genes found here and in. it is interesting to note that these authors also reported that immune genes, e.g., histocompatiblity genes, were overexpressed in males with respect to females, in agreement again with our results. differential expression between sexes for sscy genes can be considered as validated. we also validated probe, annotated as clock gene, by real time quantitative rt-pcr in backfat samples of pigs. qrt-pcr results clearly confirmed the microarray data. the average relative expression was ± in females, whereas we were not able to reliably detect expression in males at 1/ cdna dilution because we were below the limit of detection. to verify whether this gene product is actually expressed in males, we repeated the assay at 1/ cdna dilution. we could detect low levels of expression in four out of the males studied. iberian male that with largest relative quantification, exhibited nonetheless a expression level ~ fold lower than in females. sex: f, female; m, male rq: relative quantification rq, average of three replicates each, at 1/ cdna dilution. the results are given in terms of sample number used as calibrator, i.e., rq = 1 rq: relative quantification rq, average of three replicates each, at 1/ cdna dilution, only analyzed for males bld, below limit of detection na, not available the list of probe sets corresponding to differentially expressed genes at fdr < between lw and ib breeds is in table the thirty three probes corresponded to genes and hypothetical proteins; were over expressed in large white and over expressed in iberian. the lw-biased list comprised genes involved in atp biosynthetic process, nervous system development, lipid metabolism, rna processing, and control of growing and cellular division. the ib-biased genes affected the regulation of cell growth, rna splicing factor activity or muscle development. these results hint that transcriptome differences between breeds seem to be manifold. no particular single biological process was predominantly affected. it is to be noticed, nonetheless, that there were several genes involved in nervous development or growth and cell cycle. as there is currently very little information on the transcriptome differences between pig breeds, it is not possible to know how general these results will be. some recent results have compared muscle in duroc vs. taoyuan breeds, in landrace vs. tongcheng and liver and muscle in duroc vs. pietrain. while these studies have focussed on a particular tissue and with different techniques and microarrays, all suggest that breeds differ in genes related to a variety of functional categories. in some cases, the microarray included predominantly genes involved in muscle structure and development or energy metabolism so their results may not be representative of the whole genome. it should be noted that the differentially expressed genes reported here were those that exhibited largest differential expression across all tissues. this follows from the specification of models, where all probes were analyzed simultaneously and where the tissue effect was included in the models as just an additional effect. thus, it is important to bear in mind that other genes could show larger differential expressions in a specific tissue than those listed in tables and the experimental design of this study did not allow us to characterize with reasonable fdr the specific expression of a gene in a single tissue. nevertheless, we can get a glimpse of the effect of tissue in differential expression by subdividing the tissues in groups. we reasoned that the probe × tissue interaction should be larger between tissues that clustered far apart, as their genetic programmes are more different that tissues that are very similar, e.g., back and abdominal fats. we defined the following tissue groups base on the nj-tree shown in figure 2: brain, endocrine, structural and metabolic tissues. to a large extent, these groups reflect also a common embryonic origin. we reanalyzed each of the groups with models and obtaining new zb and zs scores for each of the tissue groups. table presents the correlation between the z-scores obtained with all tissues analyzed simultaneously and each of the groups. in order to focus on the most relevant genes, the correlations shown were those obtained with the genes with largest absolute z-scores. the pattern shown was highly illuminating. for sex, there seemed to be little interaction, as zs-scores were highly correlated across groups of tissues. correlation coefficients were always >. in contrast, correlations between zb scores were much more variable and, importantly, much lower overall. thus, whereas the correlation between brain and metabolic tissues was for sex z-scores, it was only for breed z-scores. this means that most differentially expressed genes between sexes were shared across tissues whereas this was much more unlikely when comparing two distant breeds like iberian and large white. it is interesting to note that the clusters in figure also hinted that heterogeneity was larger within breed differentially expressed genes than within sex differentially expressed genes. note that the color patterns of sex differentially expressed genes were more uniform across samples than for breed differentially expressed genes. besides, clustering was much'flatter' – i.e., more uniform – within sexes than within breeds. this result may have important implications. it suggests that physiological changes responsible for breed differences have targeted differentially the transcriptome across tissues. not all tissues have been equally affected. it remains to be studied on which of the tissues the effect of breed differentiation has been the largest. our current data set does not allow us to respond to this question accurately: the probe × tissue heritabilities were very similar in all tissue groups and no clear pattern emerged from figure the tissues included in each group are listed in table the z-scores correspond to the largest differentially expressed genes when all tissues are jointly analyzed. upper diagonal, sex differentially expressed genes. lower diagonal, breed differentially expressed genes. functional annotation of sex and breed differentially expressed genes we carried a go automatic annotation with the most differentially expressed genes between sexes and between breeds. the corresponding plots are in figure in neither case, sex or breed, was a given biological process go clearly over represented. as a result, the discussion should be considered as tentative or provisional. nevertheless, we found defense genes to be over represented for both sex and breed specific genes. for sex, there was a significant excess of transcription and translation related genes. some of these genes are sex linked, as mentioned above, so an over representation of this ontology is not unexpected. the go biological processes were more scattered for breed than for sex. interestingly, there were more genes involved in nervous system development and cell differentiation than expected among breed specific genes. this might provide some clues as to what are the primary changes exerted by selection and breeding in the organism' transcriptome but further work is needed to get a definitive answer. global discussion and  this study shows that a much larger fraction of transcriptome variability is due to tissue differences rather than to sex or breed. our results agree with recent literature that showed that transcriptional co-expression is indicative of gene function and that tissues with similar function tend to cluster together. an additional relevant observation from the nj-clustering was that embryonic development seems to leave an enduring footprint on the transcriptome: tissues with similar developmental origin tend to cluster closer than those with different embryonic origins. as clustering was carried out with all probes simultaneously, this means that the footprint extends to the majority of genes. we further characterized the go biological process of the extreme genes between tissues or groups of tissues sharing a common embryonic origin. we found that, often, a sound biological interpretation was possible. for instance, when all extreme genes across the three embryo layers were considered jointly, the most frequent and significant ontologies were related to development. we also indentified several genes involved in tissue development that have a distinct expression pattern in tissues with a common embryonic origin. this distinct expression pattern is maintained beyond embryo development. a key element of our analysis was mixed model methods. although well known to statisticians and animal breeders, these techniques are not so widely employed to dissect microarray data. nevertheless, its use have been advocated and successfully employed in a number of studies, e.g.,. we report that mixed modeling explained a large part of total variability, with a very parsimonious parameterization. a relevant advantage over standard methods is that all data are analyzed simultaneously as opposed to more utilized methods like anova – based techniques or bioconductor's package limma; as a consequence, reduced standard errors are expected. in a second step, one can use the solutions from the models to carry out clustering or differential expression studies, as we did here. although the variance explained by the effect of tissue was much larger than that explained by sex or breed, which were comparable, we were able to identify genes with differential expression between breeds and between sexes. noticeably, different genes were affected by each factor, i.e., the genes primarily affected by sex were different than those affected by breed or tissue. are genes whose expression differ between breeds differentially expressed across all tissues or is there any interaction tissue × gene? our experimental design does not allow to resolve this issue with enough power. however, the results hint that interaction can be important for differentially expressed genes between breeds but not so much between sexes. resolving this question is of the utmost interest, both to understand the profound influences of artificial selection on the organisms and to propose novel more efficient strategies for animal breeding. although more work is needed to calibrate the actual relevance of breed or genetic differences to the pig's transcriptome, it seems difficult that the importance of breed or sex differences exceeds that of tissue. thus, the argument that regulatory rather than structural mutations are a more important source of phenotypic variability needs then to be considered in perspective. it follows from our study that differences within the same individual largely exceeds those between individuals. 
as an organism is made up of many different tissues and organs, each with its own genetic programme, it is pertinent to ask: how relevant is tissue in terms of total transcriptome variability? in order to gain insight on these issues, we conducted microarray expression profiling of different tissues from four animals of two extreme pig breeds, large white and iberian, two males and two females. nevertheless, we were able to confidently identify genes with differential expression across tissues between breeds and between sexes. mixed model analysis and neighbor – joining trees showed that tissues with similar developmental origin clustered closer than those with different embryonic origins. artificial selection has resulted in animal breeds with extreme phenotypes. embryonic development leaves an enduring footprint on the transcriptome. often a sound biological interpretation was possible for overrepresented gene ontology categories within differentially expressed genes between groups of tissues. interaction with tissue can be important for differentially expressed genes between breeds but not so much for genes whose expression differ between sexes.
 artificial selection has resulted in animal breeds with extreme phenotypes. as an organism is made up of many different tissues and organs, each with its own genetic programme, it is pertinent to ask: how relevant is tissue in terms of total transcriptome variability? which are the genes most distinctly expressed between tissues? does breed or sex equally affect the transcriptome across tissues?  in order to gain insight on these issues, we conducted microarray expression profiling of different tissues from four animals of two extreme pig breeds, large white and iberian, two males and two females. mixed model analysis and neighbor – joining trees showed that tissues with similar developmental origin clustered closer than those with different embryonic origins. often a sound biological interpretation was possible for overrepresented gene ontology categories within differentially expressed genes between groups of tissues. for instance, an excess of nervous system or muscle development genes were found among tissues of ectoderm or mesoderm origins, respectively. tissue accounted for ~ times more variability than sex or breed. nevertheless, we were able to confidently identify genes with differential expression across tissues between breeds and between sexes. the genes primarily affected by sex were overall different than those affected by breed or tissue. interaction with tissue can be important for differentially expressed genes between breeds but not so much for genes whose expression differ between sexes.  embryonic development leaves an enduring footprint on the transcriptome. the interaction in gene × tissue for differentially expressed genes between breeds suggests that animal breeding has targeted differentially each tissue's transcriptome. 