Study system
We used Arctic charr of the so called ‘brown’ morph from lake Vatnshlíðarvatn, Iceland (65°51´693 N, -19°63’ 536 W). Vatnshlíðarvatn is a physically simple, small (70 ha) and shallow (average depth 2 - 3 m) lake (48). Arctic charr is the only fish species in the lake and has diverged into two weakly divergent morphs (silver and brown; (48,75)). The brown morph shows extensive variation in egg size and egg diameter, and standard length of embryos at H and FF are strongly correlated (N = 12 families, H: R2 = 0.96, N offspring = 201, FF: R2 = 0.97, N offspring = 168; both P <0.001, Leblanc et al. unpublished data).
The clutches and rearing protocols used in this study are described in (29). In short, mature females (N=7) and males (N= 5) were caught using gill nets in September 2014. Fertilised eggs were allowed to water-harden before transport. Each female was mated to a single male (i.e. full sib families), with a subset of (N = 3) females sharing the same male (i.e. half-sib families). This design causes variation in the relatedness of offspring and does not allow us to fully disentangle direct genetic from maternal effects but was used to minimize unsuccessful crosses. We henceforth use the term ‘family effects’. All adults were euthanized by lethal cranial concussion in the field after stripping and before transport back to Verið aquaculture facilities. Age of females was estimated by reading otoliths, as described by Tsinganis (76).
Fish were grown under the ethics permit of Verið aquaculture station and in accordance with Icelandic law. Experiments in this study were conducted within the ethical framework of the 3Rs (Replacement, Reduction and Refinement) as follows: 1) alternatives to the experimental procedures carried out by this research study were not available and thus the use of fish in this study was considered essential; 2) only the minimal number of fish required for the experiments were used, and all experiments were conducted on early-life stages and ended prior to independent feeding; and 3) animal welfare levels were high as experiments were conducted in an experienced aquaculture station with excellent expertise in fish husbandry practices.
Embryos were reared in the laboratory in darkness in mesh cages in a shelf incubator system with a constant flow of 95% recycled water, at a temperature of 3.5-4.5ºC. Before eye stage, a subset of at least 15 eggs were measured per family to determine mean family egg size. Developmental timing was tracked with an accumulative temperature estimate (degree days, DD; (77)). Embryos were randomly collected at: 1) H, when individuals have emerged from the egg but still rely on nutrition from the yolk sac (at 461 ± 3.11 DD); and 2) FF, when individuals begin feeding (at 658 ± 12.58 DD, see Table 2 for sample summary). Samples for the gene expression (N = 82) and craniofacial shape (N = 94) analyses were collected at the same time for each family (Table 2). All embryos were euthanized using 600ppm of 2-phenoxyethanol (78), a widely approved method (79). All individuals (eggs, or left-side of H and FF embryos) were digitally photographed (Canon EOS 650D) with a mm scale, and eggs measured four times each to obtain average diameter as a measure of egg size, and once for measurement of standard length, SL (for H and FF stages; (80)), to the nearest 0.01mm, using the program Fiji (81).
Gene expression data collection
Gene expression data was taken from Beck et al. (29), which uses offspring from the same dataset to characterise inter-individual variation in the expression of genes related to growth and skeletogenesis. A candidate gene approach using genes found to exhibit expression differences due to either family, developmental stage or offspring size in a previous study, was selected over an analysis of the whole transcriptome in order to maximise the number of individuals per family (Beck et al., in prep.; (82)). Briefly, a total of 14 candidate genes were chosen: six target genes involved in promoting trophic skeletogenesis (Timp2, Ets2, Sparc, Ctsk, Mmp2 and Mmp9) (83) and display differential expression in divergent Arctic charr morphs from lake Þingvallavatn (33); whilst the remaining eight candidate genes (Star, Igf1, Igf2, Gr, Mtor, Sgk1, Rictor and Ghr1) were chosen using a literature search based on evidence of their involvement in promoting growth, especially during early embryonic development (84–92).
Total RNA was extracted from the whole embryo at H using TRI reagent (Sigma–Aldrich, St Louis, MO) after homogenizing tissues in the Bead Beater (Biospec). The same RNA extraction methods were applied to offspring at FF. However, due to their larger size, individuals were decapitated behind the pectoral fin and only the RNA extracted from the head was included in this study as most of the phenotypic variation in Arctic charr is related to trophic morphology. As in Beck et al. (29), RNA was precipitated using isopropanol, washed with ethanol and air‐dried. The RNA pellet was resuspended in RNase‐free water and treated with DNase I (New England, Biolabs, Ipswich, MA) to remove any contaminating DNA. A subset of extracted RNAs were electrophoresed on agarose gels to test RNA quality. Single stranded cDNA was synthesized from 1 µg of total RNA, using the High Capacity cDNA Reverse Transcription kit (Applied Biosystems, Foster City, CA) according to the manufacturer's protocol. cDNA was consequently diluted in nuclease‐free water in preparation for qPCR.
Primers (29) were designed using an assembled Arctic charr transcriptome (93) and exon boundaries mapped to Salmo salar orthologs from the salmonid species database (94). Primers spanned at least one exon boundary and were selected based upon their short amplicon size (<250bp). RT-qPCR was performed on 96-well PCR plates on a 7500 Real-Time PCR System (Applied Biosystems) using 2x Fermentas Maxima SYBR Green qPCR Master Mix (Fisher Scientific, Pittsburgh, PA, USA), with a final reaction volume of 10 μl. For family-level variation we had 4-8 individuals (replicates) per family, and a total of 7 families that represent biological replicates. Each biological sample was run in technical duplicates with a non-template control in each run for each gene. Primer efficiencies were calculated (29) and both RT-qPCR and differential mRNA gene expression calculations for each target gene performed according to (33), using two validated early developmental Arctic charr reference genes, Actb and Ef1a (95). Reference genes were normalised by randomly selecting one individual within each developmental stage as a calibrator sample. Relative expression quantities (RQ) were calculated according to (33).
Staining of craniofacial elements and digitization
For analyses of morphology, individuals were fixed in 4% paraformaldehyde (PFA) and stained with Alizarin red for bone, and Alcian blue for cartilage, using a modified acid-free double staining protocol (46,96), see Appendix S1. Based on a previous study (56), as well as preliminary trials, we selected craniofacial elements that were clearly visible at both H and FF stages (Fig. 1; for an overview on all craniofacial structures in developing embryos and time of appearance (56)).
All samples were stained simultaneously within each developmental stage and the same staining solutions were used across both developmental stages to ensure as much consistency as possible. Once stained, individuals were placed in a petri-dish containing 2% and 3% transparent methylcellulose (Appendix S1) for H and FF, respectively, to allow easy maneuverability during photography and to ensure that the embryos lay flat. Individuals were photographed ventrally using a HD digital microscope camera (LEICA MC170 HD) mounted on a stereomicroscope (LEICA M165 C), with a scale set for each photograph. Each individual was photographed twice, each time removing and re-positioning the fish to account for measurement error due to placement. The two sets of photos per individual were duplicated so each individual was represented by a total of four photos: two photos to account for placement error, and two photos to account for digitizing error. Photos of disfigured individuals (N = 8) were removed and remaining photos randomized before digitization. To capture the ventral craniofacial shape of developing Arctic charr and to enable comparisons between H and FF stages, we chose 17 fixed landmarks and 54 sliding semi-landmarks (to increase effect sizes (97)) that were homologous between developmental stages, and placed on ventral surfaces of individuals in tpsDig2 v.2.31 ((98); Fig. 1). Landmarks were placed on trophic structures, such as the hyoid arch, Meckel’s cartilage and ethmoid plate. Semi-landmarks were placed at set equal distances along a curve between two fixed landmarks, using tpsDig2. The relative measurement error due to the differences in placement of individuals for photos, as well as digitising error, was evaluated for each developmental stage using Procrustes ANOVA (99) in MorphoJ (100), before taking the average of all photos per individual to remove such measurement errors.
Statistical analyses
Gene expression
Analyses were performed on log2 transformed RQ values and differences in gene expression across stages and families can be found in (29).
Geometric morphometrics for craniofacial shape
All morphometric analyses were conducted in R v.3.3.2 (101) using the R package geomorph v.3.0.3 (102). Procrustes shape residuals were obtained using a generalized Procrustes analysis (103), which optimally superimposes landmarks according to location, size and orientation. The resulting Procrustes residuals were then used for analysis of object symmetry (99), using the function bilat.symmetry to average landmarks across the line of symmetry to remove variation due to side, as well as averaging our replicates.
All analyses were conducted both with and without the effect of offspring size. For size-free shapes, residuals from the relationship between offspring size and morphology were obtained using procD.lm. The influence of developmental stage (H and FF) and family effects, as well as their interaction, on offspring craniofacial shape was determined using ANOVA from procD.lm (Model 1). As procD.lm uses Type1 sum-of-squares, the position of the last variable was alternated in sequential regressions to obtain an unbiased estimate of the proportion of shape variance attributable to each variable (104). A PCA was conducted on Procrustes shape residuals to visualize shape changes across families and developmental stage by plotting the first two axes from the PCA, whilst also showing morphological variation at the extremes of those axes in relation to the average morphology using the function plotTangentSpace. The first two PC axes explaining the most variation in craniofacial shape were then correlated with individual offspring size (standard length, mm), family and developmental stage (Model 1b), and significant relationships plotted. Data on egg size was only available as family means (Table 2). Therefore, family and mean egg size were not included within the same model. For Model 2 we tested whether female mean egg size had an effect on offspring craniofacial shape across developmental stages.
Further analyses were then conducted within each developmental stage: Model 3 examined the effect of offspring size on morphology, with family included as a random effect when the inclusion of both family and offspring size were not confounded. In cases where offspring size and family were confounded, they were tested separately; andfinally, Model 4 determined the extent that mean female egg size had on craniofacial shape of offspring.
Gene expression and craniofacial shape
To enable comparison of gene expression and craniofacial shape, we used family means within each data set (gene expression and shape at H and FF stages; N= 7) to determine the extent of their covariance. Using family means was necessary as data could not be collected on the same individuals for both gene expression and shape. Two-block partial least squares (PLS; (105)) analysis was used to assess whether craniofacial shape (Procrustes residuals) covaried with the relative expression of genes related to skeletogenesis and/or growth within each developmental stage (due to the significant difference in gene expression between H and FF (29)).
A singular value decomposition is obtained from the two covariance matrices (craniofacial shape and relative gene expression), whereby the resulting PLS axes are uncorrelated, and where the first axis explains maximal covariation between both blocks (gene expression and shape; (104,105)). The amount of covariation between the two blocks was measured using a multivariate correlation coefficient (rPLS: (106)), with associated P-values based on 10,000 permutations under the null hypothesis of independence between both blocks of variables.
Individual PC plots for both gene expression and craniofacial shape at each developmental stage were used to facilitate our understanding of how differences in gene expression levels (for both sets of genes) may correspond with differences in craniofacial shape. These PC plots will be displayed with gene loadings from the PLS analyses to determine the strength of genes in any covariance between gene expression and craniofacial shape.