Genome-wide SNP analysis reveals an increase in adaptive genetic variation through selective breeding of coral

,


Introduction
Coral reefs globally are undergoing significant degradation due to mass coral mortality driven by high sea surface temperatures (Heron et al., 2017;Hughes et al., 2018).Increasingly warm seawater temperatures and losses in coral cover are also leading to changes in fundamental reproductive processes (Baums et al., 2019;Hughes et al., 2019;Shlesinger & Loya, 2019), with major implications for ecosystem resilience and recovery in the aftermath of severe, acute events such as mass coral bleaching.To safeguard the persistence and resilience of coral reefs, assisted evolution interventions to increase the heat tolerance of corals at the early life-history stages as well as in adults are being tested in the laboratory and the field (National Academies of Sciences and Medicine, 2019;Reef Restoration and Adaptation Program: Intervention Technical 2019).These interventions may also protect populations of high importance.Assisted evolution interventions include selective breeding within species, hybridisation between species, and microbial manipulations (Chakravarti, Beltran, & Oppen, 2017;Chan, Peplow, Menéndez, Hoffmann, & van Oppen, 2019;Damjanovic, Blackall, Webster, & Oppen, 2017;Morgans, Hung, Bourne, & Quigley, 2019;Quigley, Bay, & van Oppen, 2019;van Oppen, Oliver, Putnam, & Gates, 2015).The knowledge gained from experimentation of genetic interventions will facilitate optimisation of breeding designs, reduce the incidences and impacts of trait trade-offs and the help to determine the probabilities of success and timeframes associated with genetic gains.
Improved forecasting of future coral reef health based on the scope and rates of recovery is urgently needed.This will rely on estimating rates of acclimatisation and adaptation to predict how populations will respond to a suite of selective pressures.The incorporation of the effects of selection and adaptation to understand population evolutionary trajectories has been applied to the conservation of species like Eltham's butterfly and bighorn sheep (Creech et al., 2017;Roitman et al., 2017).Corals exist close to their thermal limits (reviewed in Drury, 2019) but exhibit a wide range of phenotypes in response to heat stress (Ainsworth et al., 2015;Kenkel et al., 2013), suggestive of a high adaptive capacity predicated on genetic diversity across populations and species (Matz, Treml, Aglyamova, & Bay, 2018).There is evidence that corals have increased their heat tolerance by ˜0.5°C in the last decade, which suggests substantial capacity to adapt to increasing temperatures (Sully, Burkepile, Donovan, Hodgson, & van Woesik, 2019).Part of this heat tolerance can be attributed to the diversity and relative abundance of dinoflagellate symbionts (family Symbiodiniaceae) that inhabit coral tissues (Baird, Bhagooli, Ralph, & Takahashi, 2009;reviewed in Quigley, Baker, Coffroth, Willis, & van Oppen, 2018).In particular, changes in symbiont community composition can attribute 1.0-1.5°C of increased tolerance to adults (Berkelmans & van Oppen, 2006) and 26x increased tolerance in juvenile corals relative to warm-warm and cool-warm crosses (Quigley, Randall, van Oppen, & Bay, 2020) and has been found to explain up to 24% of bleaching variability (Mizerek, Baird, & Madin, 2018).The genomic architecture of heat tolerance is often polygenic (Bay & Palumbi, 2015;Jin et al., 2016).Multiple candidate genes of varying effect sizes have been examined, including heat shock proteins and genes involved in immunity (Louis, Bhagooli, Kenkel, Baker, & Dyall, 2017), although some target genes may more broadly be associated with stress tolerance generally and not heat tolerance specifically.Therefore, to understand corals adaptive capacity, comprehensive knowledge is needed concerning the variation in alleles at particular loci and the topology of polymorphic genomic regions.This knowledge can then be used to understand their influence in driving the emergence of different phenotypes that ultimately enable adaptation in populations over ecological and evolutionary time frames.
Evolutionary genetics describes gene frequency changes across populations and over time.The analysis of genetic diversity in wild and artificially-bred populations is central to understanding the rates and potential of adaptation.Shallow and deep whole genome sequencing is now widely available and relatively cost effective, including sequencing whole genome single-nucleotide polymorphisms (SNPs) (Davey et al., 2011;Helyar et al., 2011) and are increasingly illuminating the ecological and evolutionary mechanisms coral use to respond to their environment (Dixon et al., 2015;Fuller et al., 2019;Palumbi, Barshis, Traylor-Knowles, & Bay, 2014).SNPs are particularly valuable for identifying common variants underlying trait variation and can be applied to genetic conservation practices through breeding programs that utilize quantitative genetic principles, breeding values, and narrow and broad-sense heritability (h 2 , H 2 ) (Visscher, Hill, & Wray, 2008).SNPs are intrinsically linked to the h 2 given the proportionality to the product of SNP quantity and effect size (Holland et al., 2019).This marker-assisted selection approach identifies markers ultimately underlying the heritability estimates, although single alleles may explain only a small proportion of the measured heritable variation (Manolio et al., 2009).Specifically, SNPs provide the marker data used as inputs for relatedness models, that when combined with phenotypic data, are used to calculate additive genetic effects and variances and finally, heritability and the potential for selection (O'hara, Cano, Ovaskainen, Teplitsky, & Alho, 2008).
Here we apply high-density genome-wide marker sequencing to samples collected from selectively bred reproductive crosses to elucidate how artificially produced corals may evolve under heat stress.Selective breeding of corals from historically warmer reefs with others sourced from cooler reefs and then exposing juveniles to temperature tolerant symbionts like Durusdinium trenchii significantly improved coral fitness when exposed to 31°C (e.g.increased survival, growth and bleaching tolerance (Quigley, Randall, van Oppen, & Bay, 2020)), even upon inspection of a limited number of crosses produced from only five parents.Here we link those traits measured in the five reproductive crosses to their underlying allele frequency changes and genomic contribution (h 2 ) and identify the putative causative variants separating the five families to explain and potentially link these heritability estimates at both ambient and elevated temperatures to explain the genetic variation in these three quantitative traits.

Spawning, larval rearing, settlement and symbiont infection
Gravid adult corals of the species Acropora spathulata were collected from one reef in the far north of the Great Barrier Reef (Tijou, far northern GBR; 13°10'44.0"S,143°56'54.6"E,permit G16/38488.1)and one reef in the central GBR (Backnumbers,central GBR;18°30'49.8"S,147°09'10.7"E,permit G12/35236.1)and brought back to the National Sea Simulator (SeaSim) for coral spawning in November 2017 as outlined in (Quigley et al., 2020).Briefly, on the night of spawning the separated eggs and sperm from three far northern colonies were mixed with those of three central colonies, resulting in larvae from 30 distinct familial crosses that were raised in 27.5°C ultrafiltered, flow-through seawater (see Table 1 for pedigree and parental colony designations).Only five families were available at the time of larval settlement and used to produce recruits (hereafter referred to as juveniles given measurements were taken after 70 days of development).Here we focus on those five familial crosses sequenced at both the larval stage (27.5°C) and the juvenile stage (27.5°C and 31°C).These crosses include three produced from two parents sourced from the warm far northern reef (WW1, WW2, WW3), one cross with a warm dam and cool sire (WC) and one with a cool dam and warm sire (CW).Aposymbiotic juveniles were infected with Symbiodinium tridacnidorum ,Cladocopium goreaui , and Durusdinium trenchii with cultured material obtained from the AIMS Symbiont Culture Facility (further details in Quigley et al., 2020).Full details of the collection, environmental parameters of each location, permits, spawning, larval and juvenile rearing, and symbiosis establishment can be found in Quigley et al., 2020.Briefly, symbiosis establishment in coral juveniles was performed by adding 1x10 5 ml -1 of one of the three cultured symbiont strains separately to each juvenile tank (n = 3 replicate tanks per culture), suspending the flow to allow for contact between corals and symbionts, and repeating one further time until infection was confirmed with microscopy.

Larval sampling and trait measurements in juveniles
Acropora spathulata larvae were sampled into 100% ethanol after rearing at 27.5°C ˜2.5 days after fertilization.DNA from tissue samples were extracted using a modified SDS protocol previously used for sequencing of single coral eggs (Quigley, Willis, & Bay, 2017).Juvenile growth (percent change in area), bleaching (percent change in coloration), and survival (alive or dead) were quantified after 70 days of experimental exposure to 27.5°C and 31°C, and a detailed assessment of this physiological data can be found in (Quigley et al., 2020).
Sampling, DNA extraction, library preparation, sequencing Individual larvae were genotyped (total n = 68) across five families representing three population crosses (n = WW: 34, WC: 15, CW: 19).All samples were sequenced on the Illumina HiSeq2500 and genotyped using proprietary DArT-seq TM technology at Diversity Arrays Technology (Canberra, Australia).Digestion and ligation reactions were carried out with the PstI and HpaII restriction enzymes following (Jaccoud, Peng, Feinstein, & Kilian, 2001;Kilian et al., 2012).DNA purification and library preparation followed DArT proprietary methods (Kilian et al., 2012).As part of the DArT pipeline, raw sequence data was filtered based on Q-scores and filtered against viral and bacterial databases to remove contamination.The restriction site fragments generated were used to call single nucleotide polymorphisms (SNPs) with the KDCompute pipeline and DArTsoft14 algorithm (Diversity Arrays Technology, http://www.kddart.org/kdcompute.html).SNPs were referenced against the Acropora digitifera genome (Shinzato et al., 2011) and further quality controlled and filtered to produce a reduced dataset of only high quality loci using both proprietary DART software and the following criteria, which included filtering loci for average repeatability of alleles at the locus ('filter.repavg',>99%), collapsing duplicated sequences ('keep monomorphs'), individuals with missing data ('filter.callrate',>25%), filtering by Minor Allele Frequencies (< 0.02), coverage (< 8 read depth), and filtering loci out of Hardy-Weinberg Equilibrium.

Statistical Analyses
Multivariate statistics SNP metrics were recalculated in R for post-filtered data using the 'DArTR package' (Gruber, Unmack, Berry, & Georges, 2018).PCoAs were constructed from individual larval SNP genotypes with Gower PCoA ordination and Euclidean distances using the packages 'Adegenet' and 'DArTR' (Gruber et al., 2018;Jombart, 2008).Clusters of genetically related individuals were inferred using Discriminant Analysis of Principal Components (Jombart, Devillard, & Balloux, 2010) constructed with and without the between-and withingroup variance components from PCoA outputs.Multiple loading thresholds were explored, including the default value of the third quartile (75%) of x-values (Grünwald, Kamvar, & Everhart, Jombart & Collins, 2015).Diversity metrics were calculated using 'DArTR' following (Sherwin, Chao, Jost, & Smouse, 2017).These analyses are similar to classic population genetic analyses that measure differentiation via F st , in which high allele frequency differentiation suggests population differentiation although they can be fraught with false positives (Rajora, 2019).
Given the influence of close relationships in STRUCTURE analyses (Anderson & Dunham, 2008), DAPC was used as an alternative method to infer admixture (Grünwald et al., 2010).Hence, "true k " was calculated based on Bayesian Information Criteria (BIC, Jombart et al., 2010), in which the largest drop in BIC was at k = 5.

Single Nucleotide Polymorphism annotation
SNPs significantly contributing to the PCoA separation between each family cross were identified, translated to nucleotide sequences by searching and aligning sequences to the Acropora digitiferagenome (v.OIST v1.1) using custom scripts.All nucleotide matches were downloaded and re-formatted into a searchable BLAST database using the "blastdbcmd" command.Nucleotide matches were blasted to annotated proteins (OIST v.0.9, see reference for full instructions on annotations).Predicted proteins were filtered for only protein matches with E-values <1e-4 were retained (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM= blastx&PAGE TYPE=BlastSearch&LINK LOC=blasthome).

Modelling SNP frequency distributions
The distributions of expected Hardy-Weinberg (HW) equilibrium genotype frequencies were simulated using a custom R script (Chivers, 2011) for three population sizes of the SNP profiles from single larval samples (number of larvae in each population, n = 34, 19, 15).A null distribution of allele frequencies in a population that is not evolving follows HW proportions.SNP profiles were constructed from single larval samples for all families as well as bioinformatically pooled by geographic source of parental corals (WW, WC, CW).Parameters included were number of parents, total population sizes, number of genotypes, and number of generations to simulate in randomly mated populations.Expected HW genotype frequencies were statistically compared to the experimentally observed genotype frequencies from the families and a multidimensional, distribution-free Kolmogorov-Smirnov Two-Sample test was used to assess statistical significance of differences between paired continuous two-dimensional distributions using the package 'Peacock.test'(Fasano & Franceschini, 1987).

Narrow-sense heritability
Physiological data on survival, bleaching, and growth (Quigley et al., 2020) were combined with pedigree information (Supplementary Table 1) to calculate narrow-sense heritability estimates for each trait at 27.5°C and 31°C.Relatedness between individual juveniles (sibships and family groups) are calculated from the pedigree input as a random effect.Narrow-sense heritability estimates (h 2 ) were derived from additive genetic variance calculated from the 'animal model' (Kruuk, 2004) using the package 'MCMCglmm' (Hadfield, 2010), in which the coefficient of relatedness between individual juveniles is parameterized as a random effect.The 'animal model' is a type of specialized mixed-effect model that calculates narrow-sense heritability by parameterizing relatedness between measured individuals using random effects (Kruuk, 2004).Hence, this analysis calculates the relative contribution of genetic and environmental influence on phenotypic variance of each trait (Wilson et al., 2010).This analysis builds upon that performed in Quigley et al., 2020 by: 1) using the underlying physiological data and 2) extending the mixed effects models used in that study that parameterized the influence of population of origin (WW, WC, CW) and parental cross (WW1, WW2, WW3, WC, CW) on survival, growth, and bleaching.Three Symbiodiniaceae strains were used for inoculation of juveniles.For clarity, heritability estimates for only Cladocopium , the punitively principal symbiont of Acropora spathulata on the GBR (sensu van Oppen, Palstra, Piquet, & Miller, 2001), is presented in the main results and the additional strains in the Supplementary Information.Models of h 2 of juvenile survival at both temperatures were run with time as a fixed factor, pedigree and individual as random factors, using binary logistic regression (categorical distribution), 150 × 10 6 iterations, 0.1% burn-in of total iterations and at multiple thinning levels of 7500 (27.5°C:Symbiodinium tridacnidorum , Cladocopium goreaui , andDurusdinium trenchii , 31°C: C. goreaui and D. trenchii ) or 37,500 (31°C: S. tridacnidorum ).Models of h 2 of juvenile bleaching and growth at both temperatures were run using 1.5 × 10 6 iterations, a thinning level of 1,500 and a burn-in of 10% of the total iterations.Fixed and random factors were parameterized in the survival models with a Gaussian distribution.
A non-informative flat prior specification was used, following an inverse gamma distribution (Wilson et al., 2010).Assumptions of chain mixing, normality of posterior distributions and autocorrelation were met.The posterior heritability was calculated by dividing the model variance attributed to relatedness by the sum of additive and residual variance.Given differences in iterations and thinning levels across models, a scaling factor (x' = x -x min /x max -x min ) was applied to the posterior distributions of each trait using the "rescale" function in the 'plottrix' package (Lemon, 2006) to visualize each heritability estimate equally to ease interpretation.Specifically, each posterior was re-scaled to the minimum (y = 0.00598) and maximum (y = 57.42421)y-density values for the S. tridacnidorum treatment at 31°C.

Interpopulation differentiation
After filtering, 9,031 high quality SNPs remained from the larval samples collected at ambient temperature across 68 individually genotyped larvae from the five families (n = 7 -20 for each).The five families separated in multidimensional space, in which genetic variation explained 43% of the variation between multilocus genotypes from the five families (PC1-PC3: 25.8; 13.4; 4.1%) (Figure 1A).Families that shared the same sire (Tijou 2A) roughly clustered together in Principle Components Analysis (PCoA) space (WW1 and CW with WW2), whereas the other two families were more divergent.
Analysis of population structure using Discriminant Analysis of Principle Components (DAPC) further confirmed this pattern of clustering between WW1 and CW along the first two discriminant functions.The 95% confidence ellipses demonstrate that the dispersion of multilocus genotypes also varied across the families, where WW2 shows the greatest and WW1 shows the least variance.Cluster membership probabilities calculated from DAPC analysis also showed flat distributions in between-group and within-group allele variances, demonstrating that admixture has occurred, especially between WW3, WC, WW1, and to a lesser extent with WW2, and particularly CW.

Allele frequency profiles
The majority of loci were fixed at either 0 or 1.The average distribution of loci in population purebred crosses was significantly different compared to interpopulation crosses, with a greater abundance of loci at intermediate frequencies in the interpopulational families (compared to null distributions; Kolmogorov-Smirnov Two-Sample test,p = 2.1e-26 and 1.3e-21, Figure 2).This suggests that these interpopulation hybrids have significantly greater genetic diversity than within population purebreds.When individual purebred families were analysed, the majority of the on average reduced diversity (e.g.lower abundance of loci at intermediate frequencies) was attributed to crosses WW1 and WW2, whereas WW3 demonstrated a similarly greater abundance of loci at intermediate frequencies as the interpopulational crosses (Supplementary Figure 1).WW1 and WW2 exhibited significantly different distributions compared to the distribution of WC (p = 5.5e-15 and 2.4e-4, respectively) and CW (p = 1.2e-11 and 6.7e-4, respectively) (Figure 2).In contrast, the distribution of loci in the WW3 cross was not significantly different compared to WC (p =8.7e-2) or CW (p = 0.77).All families differed significantly from expected, modelled HW frequencies when modelled again populations of the same size (p = WW and WC: 1.03e-45; CW: 2.2e-46).Observed heterozygosity (H o ), and alpha Diversity ( 0 D α and 2 D α at q = 0 and 2) was highest in the interpopulation crosses (WC, CW) compared to the purebreds, with the exception of WW2, which was also high (Table 1).Alpha Entropy at q = 1 was the highest in WC, WW2 and lowest in WW3.

QTL analysis
SNPs were correlated with the five families reared under ambient temperature.The number of SNPs associated with these five families varied depending on threshold cut-offs ranging from 2,224 SNPs (75% quartile) to 87 SNPs at the stricter threshold value (99% quartile).Twenty-six of the 87 SNPs (30%) could not be assigned annotations within the A. digitifera genome, a common issue in coral and symbiont 'omics (Figure 3 inset, sensu (Barshis et al., 2013;Peng et al., 2010;Shinzato et al., 2011;Sogin, Putnam, Anderson, & Gates, 2016)).
In the 87 SNPs that contributed to the separation of crosses at ambient temperatures, we identified proteins involved in a range of biological processes, mainly immunity and stress, cell functioning and metabolism, and calcification (Figure 3 inset).Stress related proteins included lysosomal-trafficking regulator-like proteins, CEPU-1-like protein, E3 ubiquitin-protein ligase RNF213-like, spondin, NFX1-type zinc finger-containing protein, NACHT, MAP/microtubule affinity and other processes, NLRC3-like proteins.Proteins involved with collagen production varied across families, as well as 2 dTDP-glucose 4,6-dehydratase-like and two sodium bicarbonate transporter-like proteins.A host of kinase related proteins, including those involved in cell functioning and transcription were detected.
PCoA loadings defined cut-offs for significant SNPs and hence assigned proteins.Homeodomain-interacting proteins, NFX1-type zinc finger-containing proteins and src proteins (Figure 3A), which are broadly associated with transcription, immunity and stress, and cell functioning were broadly positively related to the separation of families WW2/WW1/CW and WW3 and WC along PCoA axis 1 (Figure 1A).Alternatively, differences in cell functioning proteins were more negatively associated with their separation.Proteins contributing to the separation of WW2/WW3, WW1 and CW, and WC along PCoA2 were associated with responses to UV radiation, and MAP proteins (associated with cell functioning and immunity/stress, respectively) whereas multiple proteins associated with immunity and transcription were negatively associated with their separation.

Narrow-sense heritability
Population of origin (WW, WC, CW), family cross identity (WW1, WW2, WW3, WC, CW) and symbiont type significantly influenced juvenile responses to growth, bleaching and survivorship at 27.5°C and 31°C (full summary in Quigley et al., 2020).Specifically, juveniles with at least one parental coral sourced from the northern GBR survived significantly better, grew more, and bleached less at 31°C compared to other crosses, especially if infected with the Symbiodiniaceae symbiont, D. trenchii .Physiological data for juvenile survival, growth, and bleaching were previously presented in Quigley et al., 2020.These raw data were used to calculate Bayesian narrow-sense heritability (h 2 ) estimates for survival, bleaching, and growth demonstrate that only survival of corals infected withC.goreaui was strongly affected by temperature (Figure 4, Supplementary Table 1).Bleaching responses were generally associated with low heritability estimates at both temperatures, although juveniles associated with C. goreaui exhibited more moderate heritability estimates.Heritable genotypic variation across families contributed little to differences in growth rates at 27.5 and 31°C (Table 2).With the other three Symbiodiniaceae taxa, h 2 estimates across the three traits varied in juveniles from the five families infected at 27.5 and 31°C (Supplementary Figure 1, Supplementary Table 1).

Selective breeding produces increased occurrence of alleles with intermediate frequencies on which natural selection can act
Selective breeding is one genetic intervention strategy that may quickly increase adaptive genetic variation in corals, as such facilitating adaptation to increasing sea surface temperatures (Chan, Hoffmann, & van Oppen, 2019;van Oppen et al., 2015).This approach has been used across a range of aquaculture and mariculture species to improve a number of commercially important traits.Commercially important species like oysters and mussels share similar life-history attributes with corals, including the production of numerous, potentially low-quality gametes, and high levels of genetic diversity (e.g.R-strategist) (Ellegren & Galtier, 2016).In those systems, only a few generations of selective breeding has resulted in 30% higher growth under elevated pCO 2 conditions compared to wild populations (Parker et al., 2012) and 50% higher growth in redclaw crayfish (Stevenson, Jerry, & Owens, 2013) without significantly eroding genetic diversity compared to wild populations after pooling breeding combinations (O'Connor, Dove, & Knibb, 2016).Selective breeding over only one generation in corals has shown that significant increases in heat tolerance of larvae (Dixon et al., 2015) and juveniles (Quigley et al., 2020) are possible using these methods.Lessons learnt here may therefore provide insights into the mechanisms underpinning the success of these techniques in a wild marine context.
We found that alleles at intermediate frequencies were at greater abundances than expected under HWE in WC and CW crosses relative to the on average distribution of loci in WW crosses, suggesting an overall increase in genetic diversity resulting from selective breeding of corals from different regions of the GBR.Distributions of allele frequencies typically follow that most loci nearly reach fixation at either 0 or 1, where alleles at intermediate frequencies are less common but are important given an increased abundance of intermediate alleles are the raw material for selection (Jombart & Collins, 2015).Hence, genetic material on which selective adaptive processes can operate likely exists at greater abundances in the WC and CW relative to the WW crosses.Hybridisation may lead to increased performance/fitness relative to the parental generation, with this increased performance having been linked to dominant, regulatory, single quantitative trait nucleotides (Jakobson & Jarosz, 2019).Increases in genetic diversity, hybrid vigour and genetic rescue are well-known but inconsistent features of intra-and inter-specific hybridization (Chan, Hoffmann & van Oppen, 2019;Flowers et al., 2019;Hazzouri et al., 2019;Weeks et al., 2011) but has not yet been demonstrated in the selective breeding of corals.For example, it was unknown whether the breeding of divergent populations would cause a decrease in genetic diversity due to processes like genome incompatibility (Hogenboom, 1975).These results support that even in a small number of crosses and contributing parents, genetic diversity improves.Observed heterozygosity and alpha Diversity ( 0 D α and 2 D α at q = 0 and 2) was highest in the interpopulation crosses (WC, CW).The 0 D α metric is sensitive to the presence of rare alleles, which is important given breeding across populations may initially introduce novel mutants (Sherwin et al., 2017), suggesting that WC, CW, and WW2 had an increased occurrence of rare variants. 1H α , also known as the Shannon Information criteria, is informative as a natural measure of evolvability, and was highest in WC and WW2.Finally, preliminary analysis of adult colonies shows Tijou corals had on average higher observed heterozygosity compared to Backnumbers corals (H o = 0.065 vs. 0.0217, unpublished data), suggesting that there could have been a risk of reduced diversity in offspring.Importantly, we did not observe a loss of genetic diversity by crossing these two divergent populations.
Finally, we found that interpopulation selective breeding significantly changed the resulting allele frequencies of offspring compared to modelled HWE distributions.This is perhaps unsurprising given assumptions of random mating, no gene flow, and infinite population sizes were not met.However, this suggests the influence of breeding on the genetic architecture of the resulting F1 generation.Therefore, the interbreeding of even a small number of corals from different reefs across the GBR may result in extensive introgression and therefore accelerate the potential for adaptation to warming, although the production of F2 generations would be needed to confirm this.
We also demonstrate that variants associated with immune responses, growth and cellular operating may be re-arranged during breeding but are maintained within the next generation.This suggests that the important functional diversity (i.e., at stress tolerance genes) associated with focal populations can be maintained in the breeding process.In response to heat stress, corals alter their gene and protein expression patterns, as reflected in changes in their structural lipids, metabolism, and immune responses (Barshis et al., 2013;Sogin et al., 2016).Differences in proteins associated with collagen production and sodium bicarbonate transport were important in differentiating the five families produced in our study.Collagen is important for the production of the extracellular matrix, required for multicellularity and the spatial organization of functional units of cells (Helman et al., 2008) whereas bicarbonate transporters are pivotal for coral calcification and hence growth (Zoccola et al., 2015).Basic cellular functioning also potentially varied through the differences in dTDP-glucose 4,6-dehydratase-like proteins detected and their involvement in the non-oxidative pentose phosphate pathway (Buerger, Wood-Charlson, Weynberg, Willis, & van Oppen, 2016;Yuyama, Watanabe, & Takei, 2011), critical for glucose utilization.Differences between crosses in these foundational processes like cellular organization and biomineralization therefore suggests that even breeding across relatively few individuals has the potential to substantially create distinct genetic combinations.Protein NLRC3-like, which has been previously implicated in acroporid immune suppression in response to heat stress by acting in Toll-like receptor modification (Zhou et al., 2017), also varied significantly across families.Other proteins involved in immunity and stress were also detected (lysosomal-trafficking regulatorlike proteins), and these have also been linked to the mounting of innate immune responses through Toll-like receptor activity in mice (Westphal et al., 2017).Additional immunity related proteins were detected, including CEPU-1-like protein (Spaltmann & Brummendorf, 1996), E3 ubiquitin-protein ligase RNF213like (Iguchi et al., 2019), NACHT (Hamada et al., 2013), and spondin (Palmer & Traylor-Knowles, 2012).Interestingly, NFX1-type zinc finger-containing protein was found here and downregulated in resistant corals exposed to disease (Polato et al., 2010).Expression of nascent polypeptide-associated complex subunit alpha protein (Bellantuono, Hoegh-Guldberg, & Rodriguez-Lanetty, 2012) was similarly downregulated in resistant corals, suggesting that these proteins provide an important role in protective immune responses.
Genotyping individual aposymbiotic larva from the five families reared under ambient conditions (27.5°C) provides foundational knowledge as to how selective breeding influences underlying genetic architecture.It also sheds light on the underlying molecular origins and mechanisms of heritability, a long-standing goal of quantitative genetics (Jakobson & Jarosz, 2019).Broad and narrow-sense heritability has been quantified for a range of traits, but the underlying mechanisms have rarely been described (Dixon et al., 2015;Dziedzic, Elder, & Meyer, 2017).The majority of alleles fixed at the extremes of allele frequency distributions ("Ushaped"sensu Hill, Goddard, & Visscher, 2008) are likely driven by the small number of parents from which the larvae in each cross were derived (5 unique parental colonies used across the families, Supplementary Table 1), but may also suggest selection against heterozygotes during the aquarium rearing period.Interpopulation hybrids displayed greater genetic diversity relative to within population purebreds, which may be the result of the varying effects of selection on purebreds versus hybrids in the aquarium environment.U-shaped distributions may arise under strong cases of artificial selection (e.g.aquariums) combined with rare mutations (Hill et al., 2008).

Functional variation associated with selectively bred families
How does genomic variation lead to phenotypic differences between corals?The location and effect size of the SNP difference are important to determining its eventual phenotypic effect, and simplistically, differences in coding vs. non-coding regions are predicted to cause phenotypic effects through a variety of mechanisms (Cavallo & Martin, 2005).Assigning function to SNP differences is challenging given the majority of SNPs detected by association studies are non-coding (Nishizaki & Boyle, 2017), whereas the majority of key changes are coding (Cavallo & Martin, 2005), setting up a situation of difficult detection and classification.Analysis of population structure using DAPC links the genomic patterns seen in the multilocus genotypes with the underlying biological processes quantified in the heritability models.Using this approach, we identified alleles contributing to the separation of these selectively bred families, revealing that breeding of the selected populations targets changes to the immunity and stress responses and growth, likely important processes in survival generally.Assigning differences in SNPs to phenotypic differences between individuals will be key to understanding and increasing thermal tolerance for intervention methods.
Selective breeding influences the level of admixture and population discontinuity Admixture events affect members of species, populations and individuals differently (Lawson, Van Dorp, & Falush, 2018).We saw this in the extent of admixture and its associated variance across the five families, in which some families exhibited very little admixture relative to others (especially CW).Although this can be somewhat dependent on the k structure of the model used, both patterns were explored independently using two techniques (DAPC and PCoA) and in conjunction with AIC, confirmed statistically the likelihood of population differentiation.Therefore, this may suggest that the shared ancestry of the colony sourced from Backnumbers2 may be limited between the other Backnumbers and Tijou corals or whose origins were from few or divergent founders (Lawson et al., 2018).This would suggest that adult colony Backnumbers2 is not highly related to other Backnumbers or Tijou colonies.Furthermore, PCoA and DAPC both demonstrated that the five families separate out in multidimensional space given the magnitude of allelic covariance between individuals.Irrespective of population labels, DAPC analysis also recapitulated the number of selectively bred families produced, although interestingly, the purebred families were not assigned to single population clusters but instead retained hybrid clustering structure in which the proportion of ancestry was shared between multiple two to three clusters simultaneously.The discontinuity between populations was also surprising, as demonstrated by the reduced spread of individuals between clusters, especially in WW1, suggestive that selective breeding produces discrete differences in the underlying genomic architecture of F1 individuals, even in populations likely undergoing some underlying level of gene flow (Lukoschek, Riginos, & van Oppen, 2016).

Heritability
Heritability estimates provide information as to the extent to which phenotypic trait variation has a genetic basis and the potential for adaptive change in these traits (Visscher et al., 2008).SNPs are intrinsically linked to narrow-sense heritability given h 2 is proportional to the product of the number of SNPs and their effect size (Holland et al., 2019).We found a high level of variation in survival, bleaching, and growth attributed to narrow-sense heritability (h 2 ) when pooled at the level of population cross, which varied in juveniles infected with C. goreaui across the two temperatures tested here (see Supplementary Information for further details on S. tridacnidorum and D. trenchii ).
Interestingly, additive genetic variance was only important in influencing survival at 31°C when infected with C. goreaui but not at 27.5°C.Hence, it is likely that there is significant selection on survival at high temperatures forC.goreaui , consistent with the generally overall poorer performance of these symbionts at both temperatures (S. tridacnidorum ) or at higher temperatures in coral juveniles (C.goreaui ) (see Table 6.1 in Quigley, Baker, Coffroth, Willis, & van Oppen, 2018).The shapes of the posterior heritability distributions are informative for generating hypotheses concerning the underlying genomic architecture associated with selective breeding.Bimodal heritability estimates can be driven by single large effect eQTLs (as seen generally for traits with high heritability), whereas multiple loci of small effect size generally result in unimodal distributions (Rudra et al., 2017).The influence of cis and trans co-regulation of the traits may also influence the underlying distribution (Yang et al., 2014).Bimodal heritability (often with wide credibility intervals) estimates can also be driven by high underlying within-population variability even with large sample sizes (mirrored in the high phenotypic variability between genotypes, Quigley et al., 2020), and have been recorded across a range of disease-related expression of genes (Yang et al., 2014).High heritability at ambient temperatures may also be common (Kronenberg et al., 2019), and potentially indicative of core traits under continual selective pressures.Finally, it is likely that not all the SNPs associated with these three phenotypes have been captured although narrow-sense heritability should incorporate additively all the common SNPs associated (Holland et al., 2019) and it will be important to use these methods under a variety of experimental settings to fully develop models of adaptive selection.

Conclusion
The identification of potential allele frequency differences among familial crosses can link specific genotypes to important fitness traits.Although only a small number of crosses are presented in this pilot study produced from the mixing of five adult colonies, our findings suggest that selective breeding has measurable effects on the genomic diversity of the resulting offspring.This was apparent in the significant differences in genetic variants and phenotypes associated with stress tolerance when exposed to varying temperatures.These results are in line with studies demonstrating the influence of recent immigrants having significant influence on allele frequency shifts.h 2 (mean ± 95% Bayesian Credibility Intervals) h 2 (mean ± 95% Bayesian Credibility Intervals) Bleaching 0.13 (0.02 -0.34) 0.15 (0.03 -0.34)Growth 0.01 (0.00-0.76) 0.01 (0.00 -0.67)

Figure 1 .
Figure 1.Genetic differential among population crosses.Single Nucleotide Polymorphism (SNP) cluster analysis for 68 individually genotyped Acropora spathulata larvae produced from five reproductive crosses (red: WW1, WW2, WW3), warm-cool (yellow: WC), cool-warm (blue: CW).(A) Principal Coordinates Analysis and(B) Discriminant Analysis of Principal Components (DAPC) of genetic structure across the SNP dataset.Point colors correspond to five genetic crosses.(C) Admixture inference using DAPC atk = 5.Each vertical bar represents an individually genotyped larva, with colors indicative of the membership probability for that population into each of the assigned genetic clusters.

Figure 3 .
Figure 3. Variation in functional proteins as a result of selective breeding of corals.(A) Relative contribution of each protein to PCoA loadings (PC1 variation).Significantly different annotated proteins across the five reproductive crosses in PCoA-1 (A) and PCoA-2 (B).Proportion of proteins per functional category of proteins (inset).Colors correspond to the functional categories of proteins.

Figure 4 .
Figure 4. Narrow-sense heritability (h 2 ) of survival, bleaching and growth under ambient (blue: 27.5°C) and heat stress (red: 31°C) temperature conditions in juvenile corals infected with Cladocopium goreaui .Dashed lines correspond to mean h 2 of the distribution.