A novel role for pigment genes in the stress response in rainbow trout (Oncorhynchus mykiss)

In many vertebrate species visible melanin-based pigmentation patterns correlate with high stress- and disease-resistance, but proximate mechanisms for this trait association remain enigmatic. Here we show that a missense mutation in a classical pigmentation gene, melanocyte stimulating hormone receptor (MC1R), is strongly associated with distinct differences in steroidogenic melanocortin 2 receptor (MC2R) mRNA expression between high- (HR) and low-responsive (LR) rainbow trout (Oncorhynchus mykiss). We also show experimentally that cortisol implants increase the expression of agouti signaling protein (ASIP) mRNA in skin, likely explaining the association between HR-traits and reduced skin melanin patterning. Molecular dynamics simulations predict that melanocortin 2 receptor accessory protein (MRAP), needed for MC2R function, binds differently to the two MC1R variants. Considering that mRNA for MC2R and the MC1R variants are present in head kidney cells, we hypothesized that MC2R activity is modulated in part by different binding affinities of the MC1R variants for MRAP. Experiments in mammalian cells confirmed that trout MRAP interacts with the two trout MC1R variants and MC2R, but failed to detect regulation of MC2R signaling, possibly due to high constitutive MC1R activity.

ASIP : To sequence the full CDS for ASIP, the primers (Suppl. Table 1), 7663 and 7664 were designed using an EST sequence available at NCBI (BX309917.3). The PCR protocol consisted of the following: 1) an initial denaturation of 2 min at 92 o C; 2) 45 cycles of 30 sec at 92 o C, 30 sec at 56 o C, and 1 min at 68 C. The PCR product was purified and sequenced by BigDye Terminator v3.1 kit (Applied Biosystems), using 7663 and 7664 primers. In this case Herculase (Agilent Technologies) was used for amplification.

MRAP and MRAP2:
For amplifying MRAP, the primers 7806 and 7813 were designed (Suppl. Table 1) based on the salmon contig 107111 provided by the International Collaboration to Sequence the Atlantic Salmon Genome (ICSASG). The PCR protocol using Herculase (Agilent Technologies) consisted of the following steps: 1) an initial denaturation of 10 min at 92 o C; 2) 45 cycles of 30 sec at 92 o C, 30 sec at 58 o C, and 1 min at 68 o C. The PCR product was purified and sequenced by BigDye Terminator v3.1 kit (Applied Biosystems). To amplify MRAP2, primer 7771 was designed based on the salmon contig 30619, while primer 7860 was based on contig 25970 (Suppl. Table 1). Both contigs were provided by ICSASG. The PCR protocol using Herculase (Agilent Technologies) consisted of the following steps: 1) an initial denaturation of 10 min at 92 o C, 2) 45 cycles of 30 sec at 92 o C, 30 sec at 58 o C, and 1 min at 68 o C. The PCR product was purified and sequenced by BigDye Terminator v3.1 kit (Applied Biosystems).

MC2R:
Full length MC2R nucleotide is available (EU119870) in GenBank. Primers 7650 and 7651 (Suppl. Table 1) were designed to amplify the complete coding region. A previously described PCR protocol was followed 2 .

Tissue specific gene expression
Quantitative gene expression was assessed by real competitive PCR (rcPCR) 3 as described by the manufacturer. The primers and competitors were designed for the standard homogenous MassEXTEND (hME) stop reaction, although iPLEX enzyme and iPLEX protocol were used for the extension PCR, as they increase the efficiency in the analysis of high multiplexes. The relative gene expression was determined using specific PCR primers to co-amplify cDNA, and competitors of known concentrations. The amplification was followed by a second extension PCR and scoring on the MassARRAY system. Furthermore, 2 technical PCR replicates and 2 printing replicates were made to ensure the reproducibility. Samples were run on four different plates with 2 fish from each of the four groups, LR vs. HR and spotted vs. non-spotted on each plate. Skin and interrenal tissue samples were collected from 8 LR and 8 HR fish, and 8 spotted and 8 non-spotted fish from a commercial fish farm in Valdres, Norway (Valdres Ørretoppdrett, Valdres, Norway, stock is AquaGen ASA, www.aquagen.no). All these fish were non-stressed, i.e. sampled directly from their rearing tank or from social isolation. Expression of MC1R (paralogs 1 and 2) was measured in one multiplex using allelotyping 4 with ACTB and GPDH as housekeeping genes (for primer sequence see Suppl. Table 2). ASIP and MC2R were amplified in another multiplex (for primer sequence see Suppl. Table 3), together with the housekeeping genes ACTB, GPDH, S18 and RS11. Purified total RNA was treated with TURBO DNA-free™ (Ambion) for removal of contaminating DNA. First strand cDNA synthesis was conducted using SuperScript™-II RNase H, Reverse Transcriptase (Invitrogen), oligo-dT primer T270, and 5 μg of total RNA from each tissue was used as template.

Effect of cortisol on ASIP expression
In total, 16 LR trout were isolated in 125 l compartments in glass observation aquaria. The fish were allowed to acclimate to the experimental set-up for 10 days. During acclimation they were fed pelleted trout feed (1% of body weight) once daily, between 12.00 and 16.00, by dropping pellets one by one into the aquarium. Surplus feed was removed. On day 11, 8 fish were injected with high dose of cortisol (84 mg/kg body weight) prepared in coconut oil:palm oil, (60:40) ratio, in the peritoneal cavity. Another 8 fish were injected with the vehicle solution to serve as control. After 14 days post injection, fish were anesthetized and blood samples were collected from the caudal vein. Skin and interrenal tissue samples were collected for expression analyses. The multiplex shown in Suppl. Table 3 was used as described in the paragraph above.

Correlation between MC2R, MRAP and MRAP2 expression.
The correlation of MC2R mRNA level with MRAP and MRAP2 mRNA level in trout was analyzed by real time PCR. Head kidney tissue from 8 non-stressed HR and 8 non-stressed LR trout was collected in RNAlater. Total RNA was isolated, and qtRT-PCRs were carried out using a Roche LC480 light cycler (Roche Diagnostics). Reaction volumes were 10 µl, and included Light cycler ® 480 SYBR Green I Master (Roche Diagnostics GmbH), primers (5 µM each, for primer sequence refer to Suppl. Table 4) and cDNA. Cycling conditions were as follows: 10 min at 95°C, 42 cycles of 10 sec at 95°C, 10 sec at 60°C and 10 sec at 72°C followed by melt curve analysis. All reactions were run in duplicate and controls without DNA template were included to verify the absence of contamination. Relative gene expression data was calculated from qtRT-PCR raw data using the formula:ICE Cp /GOIE Cp = Expression of GOI in ratio to IC, where IC is internal control (ACTB), GOI is gene of interest, E is priming efficiency, and Cp is the crossing point. E values were calculated for each qtRT-PCR reaction using LinRegPCR software, version 11.30.0 4 .

Statistical analysis.
Since two housekeeping genes were common in both multiplexes for the real competitive PCR analysis, the data were controlled for possible multiplex effects. The correlation between AF and ACTB was 0.78 (p<0.001) and between ACTB and GDPH was 0.90 (p<0.001), indicating few multiplex effects in the estimation of the expression level. As the other housekeeping genes showed high variability between samples, expression levels were normalised using the average of GPDH and ACTB. Hence normalised expression levels were obtained as zig = yignfi where yig is the observed log10 concentration of sample i and gene g, and nfi is the average of the log10 concentration of AF and ACTB for sample i (i=1,…,64, g=1,…,4). For each gene and tissue, a linear model including group (HR/LR and non-spotted/spotted), plate and the interaction between group and plate. More precisely the model used for both data sets was zijk = mean + groupj + platek + groupj x platek + errorijk (i=1,…,64, j=1,2, k=1,…,4). All these analyses were done using Matlab. Gene expression levels were then normalized with LR (or spotted) average = 1, and fold changes were compared statistically using t-test. The latter approach was also used for data from real time PCR, which was utilized when paralog specific analysis was not required. Behavior of spotted and non spotted fish was compared using t-test for the parameters total feed consumed (expressed as % of individual body weight) during the observation period, and time spent moving in the acute stress test. Plasma cortisol samples of spotted and non-spotted trout, stressed and non-stressed, were analyzed with two-way ANOVA (n=8 pr group) indicating an effect of stress at the level of p= 0.01 and an effect of group (spotted vs non-spotted) at p=0.05. Correlations were analyzed statistically by linear regression.

MC1R:
Sequencing of MC1R revealed two distinct variants of this gene that showed 97 % identity (100 % coverage) at the nucleotide level. We assume that these two variants represent two paralogs of the MC1R gene in rainbow trout, and the paralogs were termed MC1R_paralog1 (Accession#: FN821693) and MC1R_paralog 2 (Accession#: FN821694), respectively. Both paralogs contain an open reading frame that encodes a protein consisting of 338 amino acids. The MC1R paralogs differed only at one amino acid (aa) position. While paralog 1 has Methionine (Met) in aa-position 23, paralog 2 has Valine (Val) in the corresponding position. Alignment of the predicted full-length amino acid sequence of rainbow trout MC1R together with other known vertebrates suggests that trout MC1R has highest sequence identity with Dicentrarchus labrax (79 % identity, 94 % coverage). The high number of synonymous (24) compared to non-synonymous (1) nucleotide differences between the two paralogs indicates a strong selection against non-synonymous mutations. This is a bit suprising given the existence of two paralogs of this gene, and the relatively large amount of MC1R variation found in other species e.g. humans 5 . The apparent selection pressure, and the fact that the two paralogs are almost equally expressed in the experimental groups, may indicate that both paralogs are functional genes.

ASIP:
The putative ASIP sequence contained an open reading frame of 122 amino acids (Accession#: FN821692). The overall amino acid sequence of trout agouti is 55 % identical to Takifugu rubripes (100 % coverage). No indication of any paralogous ASIP sequence was observed. There is one potential SNP (c.181A>C) observed, leading to replacement of Thr with Pro. The SNP did not follow any of the phenotypic groups.

MRAP and MRAP2:
The putative MRAP sequence (FR837908) contained an open reading frame of 78 amino acids, which is shorter than reported in other organisms. The overall amino acid sequence showed 72% identity (69 % coverage) to Meleagris gallopara and Gallus gallus. MRAP2 (FR837909) contained an open reading frame of 238 amino acids. The overall amino acid sequence was 48 % identical (91 % coverage) to Monodelphis domestica.

MC2R:
No indication of any paralogous MC2R sequence was observed. Several SNPs were observed, but neither followed any particular phenotypic group.

Modeling of MC1R_paralog 2
A structural model of MC1R_ paralog 2 was made using I-TASSER. The seven conserved transmembrane helices (TM1-7) are clearly visible ( Figure S1). The residue displayed in ball stick corresponds to the L176M residue, and is located in the 4th transmembrane domain (TM4). Residue 176 in rainbow trout MC1R corresponds to human MC1R residue 170. This residue is Leu in most bony fishes and birds, while Val or Ala is most common in mammals. However, in mouse Met is found in the corresponding position. Both the location and the lack of conservation of residue 176 might suggest that it is not crucial in ligand binding. However, this transmembrane domain may be involved in MRAP binding in the endoplasmic reticulum and subsequent membrane trafficking 6 .

Regression between the relative expression of Melanocortin 2 Receptor and (A) Melanocortin 2
Receptor Accessory Protein (MRAP) and (B) MRAP2 in head kidney tissue from high-(HR) and low-responsive (LR) trout lines. Gene expression data are presented as fold change with LR average = 1. MRAP produced a significant common regression line with MC2R (difference between slopes p=0.62, deviation from linearity p=0.62). Separate regression lines were also significant for both HR and LR, with linear regression p=0.02 and p=0.03, respectively. MRAP2 showed no such relationship (two separate non-significant regression lines, difference between slopes p=0.04).