{'fraction_non_alphanumeric': 0.030934866273688313, 'fraction_numerical': 0.0326399376431042, 'mean_word_length': 4.816803627089827, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ovine HSP90AA1 Expression Rate Is Affected by Several SNPs at the Promoter under Both Basal and Heat Stress Conditions

The aim of this work was to investigate the association between polymorphisms located at the HSP90AA1 ovine gene promoter and gene expression rate under different environmental conditions, using a mixed model approach. Blood samples from 120 unrelated rams of the Manchega sheep breed were collected at three time points differing in environmental conditions. Rams were selected on the basis of their genotype for the transversion G/C located 660 base pairs upstream the gene transcription initiation site. Animals were also genotyped for another set of 6 SNPs located at the gene promoter. Two SNPs, G/C−660 and A/G−444, were associated with gene overexpression resulting from heat stress. The composed genotype CC−660-AG−444 was the genotype having the highest expression rates with fold changes ranging from 2.2 to 3.0. The genotype AG−522 showed the highest expression levels under control conditions with a fold change of 1.4. Under these conditions, the composed genotype CC−601-TT−524-AG−522-TT−468 is expected to be correlated with higher basal expression of the gene according to genotype frequencies and linkage disequilibrium values. Some putative transcription factors were predicted for binding sites where the SNPs considered are located. Since the expression rate of the gene under alternative environmental conditions seems to depend on the composed genotype of several SNPs located at its promoter, a cooperative regulation of the transcription of the HSP90AA1 gene could be hypothesized. Nevertheless epigenetic regulation mechanisms cannot be discarded.


Introduction
Current concern about global warming and its effects over agricultural and livestock production systems have opened novel scientific opportunities to study the adaptation of organisms to new and harsher environmental conditions. In this context, the study of the genetic basis of traits linked with adaptation and fitness has great importance. Heat is one of the main sources of stress which has an important impact in livestock production. The genetic variability underlying animal's thermo tolerance could be exploited in livestock breeding programs to achieve animals that could cope with the effects of heat stress over productive and functional traits. Among the livestock animals, sheep (Ovis aries), is one of the oldest domesticated species [1]. It is widely distributed throughout the world due to its high plasticity and adaptability to withstand poor nutrient diets and tolerance to extreme climatic conditions [2], [3]. Sheep is thus an interesting biological material to study the genetic basis of thermo-tolerance. There is some literature about heat stress effects over physiological and productive traits in cattle [4], [5], [6] and sheep [7], [8], [9], [10]. Also, at the molecular level, genes involved in the heat stress response have been described [11], [12], [13], [14]. Among them, those encoding heat shock proteins have been the most studied. However, in sheep there are few works regarding this topic.
Heat shock proteins (HSPs) [15], play a fundamental role in the maintenance of cellular homeostasis, under both physiological and stress conditions [16]. HSPs are organized into several families according to their molecular size (kDa). HSPs are highly conserved across species [17], particularly the 90 kDa heat shock protein (HSP90) [18]. In eukaryotes there are two major isoforms of HSP90 constituted by gene duplication [19]: the inducible form, HSP90a and the constitutive form, HSP90b. The HSP90AA1 ovine gene (DQ983231) which encodes the HSP90a protein has been sequenced, mapped and characterized in sheep by Marcos-Carcavilla et al. [20]. In their work 34 polymorphisms (12 in the coding region, 14 in the promoter - Figure S1-and 8 in the intron 10) were detected. Further on, also a new INDEL (insertion/ deletion) was observed at the promoter [21] ( Figure S1). The transversion G/C located at position 2660 in the gene promoter was associated with resistance/susceptibility to scrapie [22], with sperm DNA fragmentation in rams [23] and with the adaptation pattern of different sheep breeds to the thermal conditions in where they are reared [24]. In this last study, this polymorphism was associated with differences in the transcription rate of the HSP90AA1 gene, being either the causal mutation which putatively modifies a transcription factor binding site or in linkage disequilibrium with the causal mutation. However, the study was based on a limited number of animals and on standard basic statistical methods used to analyze quantitative real-time PCR (qPCR) data and to asses differences in expression levels of alternative genotypes under control and heat stress conditions.
In general, qPCR has provided a powerful tool for quantifying gene expression. Nevertheless, it makes necessary to carefully consider some technical and analytical factors to ensure reproducible and accurate measurements and not lead to misinterpretations [25]. However, commonly, several of these essential procedures have been widely ignored. Those technical factors include the initial sample amount, RNA recovery and RNA purity and integrity [26], [27] among others. Some factors considered as analytical are the selection of the suitable housekeeping gene(s) (HK), the experimental design [28], the statistical method used, etc.
Traditional statistical analyses have been restricted to pair-wise comparisons of treatments in which Cq (quantification cycle) values of GOIs (genes of interest) were previously normalized using standard HKs. This kind of approach does not allow to include technical and biological effects having influence over gene expression data. The joint analysis of GOIs and HKs data can lead to a better partition of such sources of variation [29] and allows checking HK stability and subsequent normalization of GOIs. Mixed model methodology makes possible this kind of approach, giving the possibility of including systematic and random effects, and interactions among them. They constitute a powerful tool in qPCR analyses including more than two treatments and multiple experimental factors [30], [31].
The objectives of this study were to 1) confirm gene expression differences observed for alternative genotypes of the G/C 2660 transversion of the HSP90AA1 gene promoter using more animals and a wider range of climatic conditions that those in [12]; 2) study the effect over the expression levels of another six polymorphisms (A/C 2601 , G/A 2528 , G/T 2524 , A/G 2522 , G/T 2468 and A/ G 2444 ) selected on the basis of their population genotype frequencies (Table S1) and relative position regarding SNP G/ C 2660 ; and 3) to simultaneously select the best HK and analyze expression data using a mixed model statistical approach that includes technical and biological sources of variation.

Ethics Statement
The current study was carried out under a Project License from the INIA Scientific Ethic Committee. Animal manipulations were performed according to the Spanish Policy for Animal Protection RD 53/2013, which meets the European Union Directive 86/609 about the protection of animals used in experimentation. We hereby confirm that the INIA Scientific Ethic Committee, which is the named IACUC for the INIA, specifically approved this study.
Animals belong to an artificial insemination centre, were raised in small groups in different barns and fed according to their necessities.

Linkage Disequilibrium Analysis and Haplotype Determination
Animal material, nucleic acid isolation, DNA amplification and SNPs genotyping. Peripheral whole blood samples were collected from 103 animals of the Manchega Spanish sheep breed in order to analyse linkage disequilibrium among the 7 SNPs of interest located at the HSP90AA1 promoter.
Animals were grouped in 48 parent-offspring trios. Trios consist of 10 sires, 48 dams and 1 offspring per pair (all females). Three of the dams were also daughters from another trio. Genomic DNA was extracted from lymphocytes according to the salting out procedure [32]. The polymerase chain reaction was performed from 100 ng of genomic DNA using CERTAMP complex amplifications kit chemistry (Biotools, Madrid, Spain) with specific primers (Forward: 59CGAGGCTCTGGCAGGCACTTGTTG39 and Reverse: 59 GCCGCCGTTCCCA GCCCTACCT 39). A 499 bp fragment of the promoter containing SNP 2660 and 6 more SNPs (2601, 2528, 2524, 2522, 2468, 2444) was obtained. The resulting PCR fragment was purified with ExoSAP-IT (USB Corporation, OH, USA) and sequenced with specific primers (shown above).
Linkage disequilibrium estimation. PLINK software [33] was used to estimate linkage disequilibrium among all pairs of the 7 SNPs measured as r 2 , the squared correlation based on genotypic allele counts [34]. Hardy-Weinberg equilibrium exact test and observed and expected heterozygosities for each SNP were also calculated using PLINK.

Expression Analysis
Animal material. In order to confirm the association of the HSP90AA1 polymorphism (G/C 2660 ) previously associated with the adaptation to different thermal conditions in sheep previously described [24], 428 unrelated rams of Manchega Spanish sheep breed were genotyped (same protocol and primers as described in [24]). All animals belonged to an artificial insemination centre, and therefore they were reared under the same environmental and management conditions. A total of 120 out of 428 rams were selected based on their genotype: 40 CC 2660 , 40 GC 2660 and 40 GG 2660 . Genomic DNA from these 120 animals was used to genotype the previosly defined 499 bp amplicon of the HSP90AA1 promoter. Genotype frequencies are shown in Table 1.
Peripheral whole blood samples from the 120 rams were collected in 3 time points, corresponding to different climatic conditions in a dry region of central Spain (Ciudad Real). The 3 time points were in March, when environmental temperature conditions are mild, and in July and August when heat stress temperatures occur. Hereafter, we will refer to the March collection as the control. The temperature humidity index (THI) equation proposed by Marai et al. [8] was used as another indicator of thermal stress. This index combines both temperature and relative humidity. The enviromental parameters for the 3 points in time are shown in Table 2.
Total RNA isolation and cDNA synthesis. Total RNA was isolated from 9 ml of whole blood using the LeukoLock kit (Ambion, Inc., TX, USA), following manufacturers instructions. RNA concentration was determined using a NanoDrop ND-1000 UV/Vis spectrophotometer (Nanodrop Technologies, Inc., DE, USA). Degradation of RNA samples was assessed with the Agilent 2100 bionalyzer (Agilent Technologies Hewlett-Packard-Str.8 76337 Waldbronn, Germany) in RNA Nano Chips, following manufacturers instructions. RIN (RNA Integrity Number) values were obtained. cDNA was synthesized using the ImProm-II TM Reverse Transcription System (Promega Corp., WI, USA).
Quantitative reverse transcription polymerase chain reaction (qRT-PCR). qRT-PCR was performed on all samples collected. Three HKs were tested, MDH1, SDHA and HSP90AB1. MDH1 and SDHA became the most stable HK pair for the heat stress response in sheep under similar conditions [38]. Also the HSP90AB1 gene was included as HK candidate since its expression is ubiquitous, less inducible and more constitutive than that of the HSP90AA1 gene [39], [40]. Primers were designed with NetPrimer software (Biosoft International, CA, USA), and are listed in Table 3 together with amplicon sizes and CG content. Primers were designed avoiding possible genomic DNA amplifications. In silico specificity of the amplicons was screened by BLAST searches. qPCR amplification reactions were performed from 100 ng of cDNA using LightCyclerH 480 SYBR Green I Master kit (Roche, Switzerland). Reactions were run in triplicate on a LightCyclerH 480 (Roche, Switzerland) following manufacturer's cycling parameters. Dissociation curves were performed for each gene to check primer specificity and to confirm the presence of a unique PCR product. The corresponding mRNA levels were measured and analyzed by their Cq.
To estimate PCR efficiencies, standard curves based on 6 serial dilutions (1/20 from a departure concentration of 50 ng/ml) of a cDNA stock (a cDNA mixture of more than 121 samples accounting for the 3 genotypes and the 3 time points) were performed. Efficiencies (E) were calculated from the slope of curves as in Rasmussen and coworkers' [41]. Estimated E for each gene are shown in Table 3.

Statistical Procedures
Statistical analysis of RIN values. A mixed model was fitted by using the MIXED procedure of the SAS statistical package [42] for determining factors affecting RIN values. RIN value of all samples were included as a dependent variable. Fixed effects included were genotype G/C 2660 (G) -3 levels: CC, GC and GG -; date of collection (D) -3 levels: Control, July and August -; group of sample processing (GP) -4 levels corresponding to the barn where a group of animals were located and sampledand the interaction date of collection x group of sample processing (DxGP) were included as fixed effects. The barn needs to be included because it is related to the period of time between samples collection and processing. The animal (A) was included as random effect. Goodness of fit statistics AIC (Akaike's Information Criterion) and BIC (Schwarz's Bayesian Criterion) were used as criteria for model selection. A type III fixed effects test was used to determine significance of the effects included in the model. P,0.05 was established as threshold for statistical significance.
HK selection. HK selection among HSP90AB1, MDH1 and SDHA genes followed the strategy from Serrano et al. [38]. As amplification efficiencies of some genes were ,2 (,100%), Cq data were transformed using the equation proposed by Steibel et al. [29] to rescale Cq values. The equation of the mixed model used was the following: where y oijkmr is the transformed Cq data of the j th gene, from the r th well, in the k th plate, collected from de m th animal under the i th treatment; M o is the fixed effect of the o th genotype; T i is the fixed effect of i th treatment; G j is the fixed effect of the j th gene; P k is the effect of the k th plate; b(RG) imnj is the interaction between the RIN value of the mi th sample and the j th gene, where b is the regression coefficient of RIN x gene variable on Cq; S im is the random effect of the biological sample (S im *N(0,s 2 s )); A m is the random effect of the animal from where samples were collected (A m *N(0,s 2 A )); MTG oij is the random interaction effect among the o th genotype, the i th treatment and the j th gene (MTG oij *N(0,s 2 MTG )); e oijkmr is the random residual. Gene specific residual variance (heterogeneous residual) was fitted to the gene by treatment effect (e oijkmr *N(0,s 2 eij )). Expression stability values were obtained by calculating the Mean Square Error (MSE), which was defined as in [38].
Analysis of expression results. Statistical analysis of gene expression was carried out following the method proposed by Steibel et al. [29]. As amplification efficiencies for HSP90AA1, HSP90AB1 and MDH1 genes were ,2, Cq data were transformed as aforementioned. The mixed model fitted was: where effects were as in model 1, except that in this case the MTG factor was included in the model as fixed effect and the residual variance was heterogeneous for the gene effect) .
To test differences, diff GOI , in the expression rate of alternative genotypes and to obtain fold change (FC) values from the estimated MTG differences, the approach suggested in [29] was used. Significance of diff GOI estimates was determined with the t statistic.
Also asymmetric 95% confidence intervals (up and low) were calculated for each FC value by using the standard error (SE) of diff GOI : Only contrasts between genotypes expression data that remained significant after the Holm-Bonferroni correction and with a FC .1 are going to be discussed. Supplementary Tables  (Tables S2, S3 and S4) show estimates, standard errors, FC and confidence intervals (FC up -FC low ) of significant contrasts between genotypes in each treatment. Also FCs are graphically represented in Figures 1, 2 and 3 where segments indicate 95% confidence interval. AvT-5daysBC = Average mean temperature of the five days previous to collection (uC).
MaT-5daysBC = Average maximum temperature of the five days previous to collection (uC). Rh-5daysBC/100 = Average relative humidity of the five days previous to collection (%)/100. Rhmax-5daysBC/100 = Average maximum relative humidity of the five days previous to collection (%)/100.
THIavr-5daysBC = THI calculated with the average temperature and relative humidity of the 5 days before collection. THImax-5daysBC = THI calculated with the average maximum temperature and maximum relative humidity of the 5 days before collection.

Linkage Disequilibrium Analysis
Results from Hardy-Weinberg equilibrium exact test and expected and observed heterozygosities are shown in Table S5. No deviations from the Hardy-Weinberg equilibrium were observed for any of the SNPs genotyped on population composed by trios. The average expected and observed heterozygosities were 0.256 and 0.306, respectively. The SNPs with the lowest allele frequencies were A/G 2522 and A/G 2444 (0.018 and 0.053, respectively).
Results from linkage disequilibrium analyses are shown in Table  S6. The analysis revealed the existence of two blocks of SNPs segregating jointly. One block was constituted by the SNPs at positions 2660 and 2528 (r 2 = 0.95), and the other was integrated by the SNPs mapped at 2601, 2524 and 2468 (r 2 between them ranged from 0.92 to 1.00). Linkage disequilibrium between other SNP pairs showed values lower than 0.06, ranging from 0.001 to 0.059.

Statistical Analysis of RIN Values
The model including the interaction DxGP as a fixed effect and the animal as a random effect, showed the lowest values for the goodness of fit criteria (AIC and BIC). Estimated animal and residual variance were 0.11 and 2.25, respectively. Type III fixed effects test showed a highly significant (p,0.0001) effect of DxGP on RIN values but no significant effect were observed for G, D and GP on the trait.  Thus, as RIN values only depend on the order in which samples were processed after their collection, it can be included as a systematic effect in the statistical model used to analyse expression data.

Best HKs
Given the linkage disequilibrium results for the 7 SNPs, the ''genotype'' included as an effect in the mixed model to test the stability of genes was G/C 2660 -A/C 2601 -A/G 2522 -A/G 2444 . Table 4 shows MSE values obtained for each gene within treatments and across genes. HSP90AB1 was in all cases the most stable gene, followed by HSP90AA1. Therefore, HSP90AB1 was selected as the only HK to normalize the expression results of HSP90AA1. The highest stability values for all genes corresponded to samples collected in August and lowest stability values corresponded to the control samples.

Environmental Conditions and Statistical Analysis of Gene Expression
As it is shown in Table 2 the maximum temperatures for days of samples collection in July, 35.0uC, and August, 34.4uC, exceeded the sheep thermoneutral zone, but this is not the case for the average temperatures, 26.8uC and 24.7uC, respectively. For both time points, average and maximum THI values occurred in the zones of severe and extreme heat stress [43]. The THImax was one unit higher in August than in July. Also in August the THImax, calculated over the five days before collection, was two  For samples collected in August, CC 2660 -CC 2601 -GG 2522 -AG 2444 showed the highest expression rate (differences in FC from 2 to 3) when comparing with other 8 genotypes. Seven of these 8 genotypes were GG 2444 , highlighting the importance of this position for the expression rate of the gene under heat stress conditions independently of the other three SNPs. Two significant contrasts pointed out the effect of SNP 2660 in terms of expression efficiency. CC 2660 showed differences in FC of 1.27 and 1.31 when comparing with GC 2660 and GG 2660 , respectively. It is important to note that contrasts showing differences in FC from 2.3 to 3.0 were those with wider confidence intervals which indicate higher estimated standard errors. This is due to the fact that in many of these contrasts, composed genotypes are at very low frequencies. In some cases only one animal exhibit the genotype compared. For comparisons where FC differences ranged between 1.2 and 1.3, confidence intervals were narrower due to the higher frequencies of the genotypes compared.
For samples collected in July, no comparisons between alternative genotypes were significant after the Holm-Bonferroni correction.
For the 5 significant contrasts of control samples the loose of the effect over the expression rate of the SNP 2444 in the gene   30). Finally, the SNP 2601 did not seem to be involved in the regulation of the basal expression of the gene. In this case, differences in the confidence intervals were smaller since genotypes compared had higher frequencies.
2) Genotype A/C 2601 -A/G 2522 -A/G 2444 . To determine the impact of SNP 2660 over the expression rate of the gene under the different environmental conditions previously described [24], we have tested a genotype of three SNPs excluding SNP 2660 . Table S3 and Figure 2, show results for contrasts among A/ C 2601 -A/G 2522 -A/G 2444 genotypes under the different climatic conditions considered.
Regarding the high decrease in the magnitude of the FC differences for the 4 significant comparisons, we can confirm the critical influence of SNP 2660 over the gene expression rate. Surprisingly, higher differences in FC were detected for contrasts in control than in August, revealing a more important effect of these SNPs on the basal expression of the HSP90AA1 gene under control conditions.
In samples collected in August, the two significant contrasts with the highest expression rate implied the AG 2444 genotype in all cases. Again, AG 2444 showed higher expression rate (FC = 1.3-1.4) over GG 2444 under heat stress conditions. Genotypes at 2601 and 2522 did not show any clear effect. No significant contrast was found among genotypes from samples collected in July.
Once more, under mild temperatures the SNP 2444 lost its effect over the HSP90AA1 expression which appeared under heat stress conditions. As in the case above mentioned, AG 2522 showed a positive effect over the expression of the gene comparing with GG 2522 (FC = 1.37-1.39).
3) Genotype G/C 2660 . Table S4 and Figure 3, show results for contrasts among G/C 2660 genotypes under the different climatic conditions considered.
Contrasts of this genotype showed lower differences in FC than that observed for the previous 2 sets of genotypes. Two comparisons were significant in samples collected in August. Differences in FC were 1.22 and 1.20 for the contrasts CC 2660 vs. Thus, when considering only SNP 2660 , results indicated no differences in HSP90AA1 expression rate across treatments, but the existence of such differences across genotypes.
In order to understand better the results obtained, contrasts involving SNP 2601 and SNP 2522 , were carried out ( Figure S2 and Figure S3). In the first analysis, alternative genotypes of A/C 2601 -A/G 2522 were compared. Only one significant contrast was found under control conditions. A FC value of 1.4 was observed for the contrast CC 2601 -AG 2522 vs. CC 2601 -GG 2522 , confirming the effect of A/G 2522 over the basal expression of the gene under mild temperatures. SNP 2601 , did not seem to play any clear role in the HSP90AA1 expression rate.
In the second additional analysis, genotypes of G/C 2660 -A/ G 2522 were compared to elucidate the relevance of SNP 2522 under heat stress conditions. In this case, significant contrasts were found for control samples and those collected in August. The importance of SNP 2660 in August was again revealed. CC 2660 had higher expression rate than GC 2660 and GG 2660 (FC = 1.24 and 1.27, respectively). Under control conditions the effect of SNP 2522 was clear (FC = 1.28 at least). Changes in the behavior of G/C 2660 under control conditions were observed as well. In this case, GC 2660 was superior to CC 2660 and GG 2660 as it was showed in previous analyses.

Discussion
In this study, we aimed to confirm the results obtained previously [24], that showed an association of the SNP 2660 of the ovine HSP90AA1 gene promoter with the expression levels of this gene under different environmental temperatures, using a higher amount of samples. We also aimed to increase the information available of the biological process underlying this type of stress conditions. With these purposes we have studied new polymorphisms found at the promoter region, which would affect the expression rate of the gene not only in heat stress events as it was firstly thought, but also modulating its basal expression.

RIN Effect
Best conservation and minimum degradation processes are critical points when sampling commercial livestock animals for expression studies. The degree of RNA degradation in the samples affects gene expression measurements. We have established that RIN values depended neither on the source of biological sample (the animal) nor on the environmental conditions surrounding samples collection. The only factor having a significant effect on RIN values was the period of time occurring between blood extraction and blood processing in LeukoLock platforms for each time point (DxGP). The higher this period of time was, the higher the RNA was degraded (lower values of RIN). Therefore, we proposed to include RIN values as a fixed effect or as a covariate in the statistical model used to analyze expression differences. In fact, the same results were obtained using both approaches (data not shown).
RIN values affect Cq of samples depending on amplicon size [44]. The higher the amplified DNA fragment is, the higher the probability is to be broken down. The length of the amplified product was more correlated with RIN values than expected. Amplicon sizes were in the range of 70 to 250 bp (Table 3) for which Fleige and Pfaffl [26] indicates a more or less independence of qPCR products and RNA quality. Furthermore, DNA CG content did not seem to affect RNA stability as it has been previously described, where lower CG degree content was correlated with higher RIN values [45]. SDHA amplicon was the one with the highest CG content (65%) but was the most affected by RNA degradation. Only for MDH1, which has the lowest CG content and a high RIN effect, this relationship seems to be true. These results reveal that the effect of RNA integrity over both the GOI and the HK should be taken into account in expression analyses.

HK Selection
A crucial aspect revealed in this work, is the need to test the stability of the candidate HKs and the GOI simultaneously. MDH1 and SDHA were previously selected [38] among 16 candidates tested, as the most stable pair in similar conditions to those evaluated here. HSP90AA1 was not included in that experiment. In the present work, we have verified that the GOI, HSP90AA1, is much more stable than the two previously selected HKs, MDH1 and SDHA. Therefore, none of them can be used to normalize the GOI expression data. The constitutive counterpart of the HSP90AA1 gene, HSP90AB1, showed the best stability values within and between treatments (Table 4), and it was chosen as the HK to normalize the expression data of the HSP90AA1. Differences in the stability of both chaperone genes might be due to the effect of the SNPs existing at the promoter of the HSP90AA1 [13] [40] and to the more inducible behavior of this gene.

qPCR Experimental Design and Statistical Methods for Expression Data Analyses
When a great number of samples and treatments are included in a qPCR study the experimental design is important since qPCR plates have a limited capacity (96 or 384 wells). In our design, plates contained a randomized set of animals, treatments, genotypes, RIN values and genes to avoid estimation biases. The repetition of one or more samples in all plates connects the plate's system allowing to remove technical nuisance from this source of variability and to compare results from all plates. We have confirmed that the plate effect is an important source of variability since differences in Cq among plates can reach values up to 1.4.
Traditional statistical methods to analyze qPCR data was restricted to pair-wise comparisons of treatments in which expression data from GOIs are previously normalized with one or more HKs. This kind of approach does not include systematic nor random effects and their interactions that could affect expression results. In the linear mixed model used in this study, GOI and HKs data are simultaneously analyzed.This model includes different sources of biological and technical variation (i.e. plate, RIN, genotypes, genes, and interactions among them) as fixed or random effects. Fitting this model let us checking HK stability, normalization of GOI data with the most stable HK(s) and test the linear hypothesis of the existence of different expression levels of the HSP90AA1 gene depending on the genotype of the mutations located at its promoter and on diverse environmental conditions.

Environmental Conditions and Gene Expression
Sheep are believed to be one of the most resistant species to climatic extremes, especially to high environmental temperatures. Environmental conditions in Ciudad Real often exceed sheep thermo neutral zone which is comprised between 5uC and 25uC [46]. As expected, expression results differed between heat stress and control conditions. However, unexpectedly, differences in expression rate among genotypes were observed in samples collected in August but not in July. The scarce differences in climatic parameters existing between August and July collects did not explain the observed differences between these time points in terms of FC. The higher THImax values at collection time and during 5 days before collection in August than in July could be the clue to such differences. Other environmental factors here unknown such as wind speed, number of hours over the comfort temperature, insulation, etc. included in Fanger's comfort equation [47] would have also contributed in such differences.
Significant differences among A/C 2601 -A/G 2522 -A/G 2444 and G/C 2660 among genotypes, found for August and control samples but not for July ones would be explained by the existence of a transition in the expression state of the gene between the basal transcription and the heat stress response. Changes in the expression ranking of G/C 2660 observed between control and August samples, and also for other SNPs in a less evident way, would support this hypothesis. Also, since the heat stress response is not a permanent state, in terms of gene expression, even when heat shock conditions are still present, acclimatization processes cannot be discarded as possible source of differences found in samples collected in July and August [48].

Expression Analysis and Genotype Comparisons
Our initial hypothesis was that differences in the expression rate of the gene with different G/C 2660 genotypes would be observed only under heat stress conditions [24]. Surprisingly, after including a higher amount of samples and a set of 6 additional polymorphisms located also at the HSP90AA1 promoter, the existence of expression differences under non heat stress situations was confirmed as well. Thus, polymorphisms located in the promoter of HSP90AA1 affect not only its expression rate as response to heat shock but also its basal transcription levels.
Differences in the expression rate found for the contrasts among alternative genotypes for the SNPs studied here suggest that the transcription of this gene may be multiply regulated by cross-talk of various transcription factors, as it was pointed out for this gene in human [39]. Although much of the heat-induced gene expression can be explained by HSF1, a perfect correlation between its binding and induction has not been found [12]. Signal transduction cascades activated by p53, Jak and Ras pathways via HSF1 binding to the heat-shock response element (HSE) and integrating to modulate HSP transcription have been reported [49]. Additional positive or negative factors may modulate the transcriptional induction of HSF1-bound genes. Moreover, eukaryotic gene expression is tightly regulated at many levels, and can vary its regulation complexity [50]. The core promoter (TATA box, initiator -INR-and downstream promoter element -DPE-), is the essential part. Next to the core, proximal enhancers as cis-control elements (i.e. CCAAT box, GC box, B recognition element (BRE) and STRE elements) might be acting. Upstream, distal enhancers (hormone responsive elements -HRE-and nuclear factor element -NFE-) and a huge diversity of regulators that recruit a cascade of more transcription factors contribute to gene transcription regulation [51], [52]. The SNPs studied in this work are located enough upstream to the beginning of the transcription initiation to consider them as binding sites of these co-regulators, or distal enhancers that do not directly activate the transcription of the gene but modulate its expression.
The role of the SNP G/C 2660 in the transcription of the gene under heat stress has been confirmed through analysis in which only this mutation is tested. However, results from the analyses of composed genotypes G/C 2660 -A/C 2601 -A/G 2522 -A/G 2444 and A/C 2601 -A/G 2522 -A/G 2444 revealed a cooperative relationship among several SNPs in terms of transcription efficiency. Thus, alternative genotypes of SNP 2660 -SNP 2444 seem to affect the expression of the gene in response to heat stress and those of SNP 2660 -SNP 2522 the basal transcription of HSP90AA1, which may occur under climatic conditions comprising comfort temperatures. Under heat stress conditions, the superiority of CC 2660 over GC 2660 and GG 2660 and of GC 2660 over GG 2660 indicated an additive effect for this mutation. However, for the control samples, the GC 2660 genotype was superior to CC 2660 and GG 2660 . The effect of the SNP 2444 was less clear due to the low frequencies of the AG 2444 and AA 2444 genotypes; however, for the SNP located at 2522, more clear conclusions can be extracted.
Several putative TFs have been predicted (Table S7 and Table  S8) for the presence of C and A at SNPs 2660 and 2444, respectively. Some TFs that could co-activate gene expression as distal enhancers only with CC 2660 were NFI/CTF (Nuclear factor I or CCAAT box-binding transcription factor) and VDR (Vitamin D receptor) together with RxR-alpha (Retinoid X receptor Alpha). The last two TFs form a heterodimer which attracts a complex of co-activators proteins. This complex links the heterodimer to the initiation complex formed at the TATA box, promoting the transcription machinery [53]. Both TFs bind putatively at the sequence around the SNP 2660 , and VDR only when this position is C. For AG 2444 a heat shock element that could bind a heat shock factor was predicted for the presence of the A nucleotide [39].
The presence of an INDEL of two adenines (AA) located at position -704 in the promoter [21] completely linked, at least in this breed with SNP 2660 , must also be considered in the expression regulation of the gene under heat stress conditions. Thus, CC 2660 animals are also homozygous for the AA insertion (AA/AA), animals with GC 2660 are heterozygous ins/del (AA/2) and animals with GG 2660 are homozygous for the AA deletion (2/2). This INDEL (AA) was located within a putative glucocorticoid receptor (GR) transacting factor binding site. The AA deletion (2/ 2) created a GR transcription site. It has been pointed out that glucocorticoids can suppress the heat shock response in stressed cells by inhibiting the action of the heat shock factor 1 (HSF1) [54]. Therefore, this mutation would be the responsible of the expression differences observed for SNP 2660 . Because of the high linkage disequilibrium between the SNPs at 2660 and 2528 (R 2 = 0.95) the possible effect of the SNP 2528 over the transcription rate of the gene under heat stress conditions is masked by the first, and therefore no conclusions could be extracted from this position.
Under control conditions, A/G 2522 seems to have a predominant effect over the transcription rate of the gene, being AG 2522 (AA 2522 was not found in these samples) superior than GG 2522 . Due to the proximity of this mutation to the SNPs located at 2528 (5pb) and 2524 (1pb) TF binding sites were predicted for a sequence containing the three SNPs (Table S7 and Table S8). Several putative TF binding sites linked to the presence of adenine at 2522 and thymine at 2524 were found. Among them, the stress response element (STRE) [39] and the JunD (functional component of the AP1 -activator protein 1-transcription factor complex) related with transcription coactivator activity, oxidative stress response [55] and spermatogenesis [56] seemed to be closer to the HSP90AA1 functions. The TF c-Fos stimulates transcription of genes containing AP-1 regulatory elements and was predicted for the sequence AtagTcA for the SNPs at 2528, 2524 and -522. In our samples animals with AG 2522 were always CC-601 -TT 2524 -TT 2468 and in most cases (70%) CC-601 -AG 2528 -TT 2524 -TT 2468 . Two putative TF binding sites implying the presence of cytosine were found for A/C 2601 , HES-1 (hairy and enhancer of split-1) [57], which can act as a repressor or activator, and USF1 (Upstream stimulatory factor 1) [58], that has been found to be involved in the stress-activated signaling cascade [59] and in the cessation of Sertoli cell proliferation and differentiation to spermatozoids [60]. For the SNP 2468 one interesting homolog of the human ZNF395 binding Sp1 was found for maize [61] linked to the response to oxidative stress (Table S7).
Most stress-response genes are regulated in a concordant manner with respect to transcript levels and translational efficiency. A strong overall correlation has been observed between transcriptional/translational induction of genes and induction of the corresponding proteins [62]. During environmental stress in fission yeast, most mRNAs are regulated both at transcription and translation level but only up-regulated mRNAs showed a strong correlation with protein expression, while down-regulated mRNAs showed no such correlation [62]. Therefore changes in the expression rate of the HSP90AA1 here observed as a function of environmental temperatures and genotypes at the SNPs located in its promoter will be accompanied by changes in the amount of protein produced. Thus, those genotypes which showed higher expression levels under heat shock also will display higher protein amounts. Despite the low magnitude of the changes in HSP90AA1 expression rate observed, even a small proportion may be significant as HSP90 is one of the most abundant proteins in most cells [39]. Higher amounts of HSP90a protein would increase its capacity to exert its protective role over the effects caused by heat stress at the cellular level. Effects of changes in the protein amount as consequence of several stress sources must be studied in tissues in which HSP90a predominates, such as brain and testis [39]. In this context, our results for two of the known functions of the HSP90a, refolding proteins with aberrant conformation [39] and spermatogenesis and meiotic progression in testis [63], confirm this hypothesis. GG 2660 was related with higher sperm DNA fragmentation values in rams under heat stress conditions [23] and also to lower scrapie incubation period in sheep (missfolding stress) [22].
Future studies will focus on testing TF interaction at binding sites where polymorphisms of the HSP90AA1 promoter are located, by employing in vitro techniques of EMSA (Electrophoretic Mobility Shift Assay). In addition, an epigenetic regulation of the HSP90AA1 expression cannot be discarded. Methylation of CpG islands located at gene promoters is a well known mechanism of expression regulation and gene silencing [64]. A CpG island has been predicted for the promoter of this gene and some of the polymorphisms here analyzed are susceptible to be methylated (SNPs located at 2660, 2601, 2528 and 2522). Therefore, future research will be also focused in the study of the methylation pattern of alternative genotypes of the SNPs containing cytosine by bisulfite sequencing techniques. Table S1 SNP frequencies of the polymorphisms found at the HSP90AA1 promoter as described in Marcos-Carcavilla's work [24]. (XLSX)