Multi-trait GWAS using imputed high-density genotypes from whole-genome sequencing identifies genes associated with body traits in Nile tilapia

Abstract

Background

Body traits are generally controlled by several genes in vertebrates (i.e. polygenes), which in turn make them difficult to identify through association mapping. Increasing the power of association studies by combining approaches such as genotype imputation and multi-trait analysis improves the ability to detect quantitative trait loci associated with polygenic traits, such as body traits.

Results

A multi-trait genome-wide association study (mtGWAS) was performed to identify quantitative trait loci (QTL) and genes associated with body traits in Nile tilapia (Oreochromis niloticus) using genotypes imputed to whole-genome sequences (WGS). To increase the statistical power of mtGWAS for the detection of genetic associations, summary statistics from single-trait genome-wide association studies (stGWAS) for eight different body traits recorded in 1309 animals were used. The mtGWAS increased the statistical power from the original sample size from 13 to 44%, depending on the trait analyzed. The better resolution of the WGS data, combined with the increased power of the mtGWAS approach, allowed the detection of significant markers which were not previously found in the stGWAS. Some of the lead single nucleotide polymorphisms (SNPs) were found within important functional candidate genes previously associated with growth-related traits in other terrestrial species. For instance, we identified SNP within the α1,6-fucosyltransferase (FUT8), solute carrier family 4 member 2 (SLC4A2), A disintegrin and metalloproteinase with thrombospondin motifs 9 (ADAMTS9) and heart development protein with EGF like domains 1 (HEG1) genes, which have been associated with average daily gain in sheep, osteopetrosis in cattle, chest size in goats, and growth and meat quality in sheep, respectively.

Conclusions

The high-resolution mtGWAS presented here allowed the identification of significant SNPs, linked to strong functional candidate genes, associated with body traits in Nile tilapia. These results provide further insights about the genetic variants and genes underlying body trait variation in cichlid fish with high accuracy and strong statistical support.

Background

Tilapia is one of the most important fish species cultivated in the world, and is currently farmed in more than 125 countries. Total farmed finfish production reached 54.3 million tons globally in 2018, and Nile tilapia (Oreochromis niloticus) represented 8.3% of this volume [1]. Tilapia is generally sold as whole fish or fillets, making body traits, such as body and fillet weight, among the most economically important traits for this species. In fact, body size traits represent the primary breeding objective in genetic improvement programs for tilapia and other aquaculture species [2]. The most important body traits in Nile tilapia are body weight measured at a specific age (e.g. body weight at harvest), fillet weight or fillet yield (fillet weight/body weight). These traits show heritability values ranging from 0.06 to 0.48, when using pedigree-based estimates [3,4,5,6,7,8,9]. Previous studies have estimated high values of genetic correlations between harvest weight and fillet weight (> 0.96) and moderate to high values between harvest weight and fillet yield (0.21 to 0.74) [7, 9, 10], suggesting that is not possible to improve fillet traits independently of body weight [11]. Although, previous reports have also identified negative or null genetic correlation between harvest weight and fillet yield [12], which suggests the importance of assessing these relationships on each particular population. Other body traits which have been proposed as selection criteria to generate more profitable commercial fish populations, are reduced waste (sum of bones, viscera, head, and fins) and carcass weight, due to their higher heritability values, less correlation to body weight, compared to fillet weight, and null or even favourable impact on fillet yield [13, 14].

The availability of a chromosome-level reference genome assembly [15] and high-throughput whole-genome sequencing (WGS) methods [16, 17], have allowed for the assessment of genetic variation of different Nile tilapia populations at a genome-wide level and the recent development of single nucleotide polymorphism (SNP) panels [18, 19]. The availability of Nile tilapia SNP panels made it possible to use modern molecular breeding approaches; including mapping of quantitative trait loci (QTL) through genome-wide association studies (GWAS), marker-assisted selection (MAS) and genomic selection [20, 21]. The GWAS approach evaluates the association between genotypes and phenotypes, with both sources of information available for a large number of individuals. This method captures the linkage disequilibrium (LD) between markers and causative mutations that tend to be inherited together across generations [22]. GWAS has been applied to provide insights into both the genetic architecture and loci underpinning the genetic variation of growth-related traits in different finfish species, including Atlantic salmon and catfish [23,24,25,26], using high-density SNPs arrays (ranging from 108 K to 218 K SNPs) and, more recently, Nile tilapia by using a medium-density (50 K) SNP array [20]. These studies revealed the polygenic nature of growth-related traits and identified some genes harboring significant SNPs, which are well-known to be involved in growth and bone development, including meprin A subunit beta-like (MEP1A), fibroblast growth factors (FGF), disintegrin and metalloproteinase domain 12 (ADAM12), myosin light chain kinase (MYLK) and transforming growth-factor beta receptor type 3 (TGFBR3).

The use of ultra-high-density SNPs or WGS can improve the accuracy and power of GWAS to detect QTLs associated with complex traits [27,28,29,30]. Although the cost of WGS is rapidly decreasing, it is still expensive to sequence all available phenotyped individuals in a GWAS design. To solve this, genotype imputation to WGS data can be successfully implemented to detect putative causal loci in a cost-efficient manner. Previous studies using imputed genotypes from WGS for GWAS have been reported in cattle [27, 28], pigs [29, 30] and sheep [31]. In addition, new strategies such as multi-trait GWAS (mtGWAS) analysis are required to increase the power to detect QTL through GWAS [32]. mtGWAS improves the power of GWAS through the incorporation of summary information contained in the output of single-trait GWAS (stGWAS). Thus, mtGWAS jointly exploits information from genetically correlated traits to increase statistical power, due to fact that the true SNP effects and their estimated error may be correlated across traits. For instance, multi-trait approaches have been implemented in pertinent software, e.g. MTAG v0.9.0 [33], and successfully applied to boost the discovery of genetic variants associated with important traits in humans [34,35,36].

To the best of our knowledge, no previous studies have shown the use of imputation to high-density SNP genotypes, in a combination with mtGWAS, to uncover putative causative genetic variants associated with body traits in aquaculture species. The objective of this study was to use mtGWAS and high-density SNP genotypes to increase the accuracy and power to identify both QTLs and genes associated with eight body traits in Nile tilapia.

Results

Descriptive statistics, quality control and genetic parameters

A total of 1309 animals averaging 370 days-old were phenotyped and genotyped. Average, standard deviation, minimum and maximum phenotypic values for average daily gain (ADG), body weight at harvest (BWH), waste weight (WW), head weight (HW), gutted head-on weight (HON), body length at harvest (BLH), fillet weight (FW) and fillet yield (FY) are reported in Table 1. The coefficient of variation ranged between 6.86 and 27.47%, with the lowest and the highest values calculated for trait FY and FW, respectively.

Table 1 Descriptive statistics for phenotypic values of body traits recorded in a breeding Nile tilapia population

For WGS, the call-rate parameter excluded the highest number of SNPs (~ 12 million), whereas MAF discarded ~ 7.8 million and ~ 253 K SNPs, for WGS and imputed WGS data, respectively. The HWE filter discarded a low number of markers, ~ 1.8 million for WGS and 79 K for imputed WGS data, respectively. After quality control applied to the 50 K SNP chip, 5905, 4114 and 3665 SNPs were removed by HWE, MAF and genotyping call-rate filters, respectively, 29,587 SNPs remained for subsequent analyses. After applying sample call-rate, all samples in both WGS and 50 K SNP chip were retained (Supplementary Table 1).

Heritability estimates calculated using the SNP-based genomic-relationship matrix (GRM) constructed with about 1 million markers ranged from 0.21 to 0.45 for the body traits analyzed here, with the lowest and the highest value determined for FY and HW, respectively (Table 2). The correlation of SNP effects among all body traits analyzed here ranged from 0.20 to 1.00, with small values only reported for correlations between FY and the rest of the traits (Fig. 1).

Table 2 Genetic parameters and comparison of association results between single- and multi-trait GWAS for Nile tilapia

Fig. 1

Correlation of SNP effects (standard error) among eight body traits in Nile tilapia. ADG: average daily gain (g); BWH: body weight at harvest (g); WW: waste weight (g); HW: head weight (g); HON: gutted head-on weight (g); BLH: body length at harvest (cm); FW: fillet weight (g); FY: fillet yield (%)

Comparison between single-trait and multi-trait GWAS

The average gain in statistical power for mtGWAS compared to stGWAS was assessed by the increase in the mean χ² statistic. Thus, we  calculated how much larger the stGWAS sample size would be expected, to be equivalent to the increase observed in χ² statistic. We found that the mtGWAS analysis corresponded to gains equivalent to increase in the original sample size from 13 to 44%. These values corresponded to an increase in sample size from 1309 for stGWAS to a value ranging from 1474 to 1890 for mtGWAS (Table 2). For instance, the number of SNP surpassing the Bonferroni corrected significance threshold for stGWAS and mtGWAS, respectively, was: 1 and 1359 for ADG, 1 and 1209 for BWH, 1 and 1347 for WW, 0 and 1595 for HW, 1 and 1138 for HON, 0 and 827 for BLH, 1 and 833 for FW, and 1 and 1920 for FY. In addition, the maximum -log(p-value) increased from 7.52 to 14.58 for ADG, from 7.63 to 14.39 for BWH, from 7.45 to 14.60 for WW, from 5.71 to 14.39 for HW, from 7.45 to 13.00 for HON, from 5.63 to 17.15 for BLH, from 7.59 to 17.75 for FW, and from 8.50 to 11.62 for FY, when comparing stGWAS against mtGWAS (Table 2).

The stGWAS identified a single significant genomic region on LG16, in position 4,178,535 base pairs (bp), associated with ADG, BWH, WW, HON and FW, and a significant SNP on LG07, in position 16,847,179 bp, for FY (Supplementary Fig. 1). When combining the summary statistics of all body traits, using mtGWAS, we identified several novel genomic regions associated with different traits. The number of SNPs surpassing the genome-wide significance threshold ranged from 827 to 1920 depending on the trait analyzed, with the lowest and the highest number of significant variants associated with BLH and FW (Table 3). The greatest number of significant variants were located on LG03 and LG12 for all traits, except FW where most of the variants were located on LG13 (Fig. 2). The location of significant variants on different chromosomes, and representation of several loci, suggest that these body traits are under polygenic control.

Table 3 Genomic regions and the closest candidate genes for the top five lead SNPs associated with body traits based on multi-trait GWAS in Nile tilapia

Fig. 2

Manhattan plot for multi-trait GWAS (mtGWAS) for eight body traits in Nile tilapia. Manhattan plots of SNPs associated with: a Average daily gain. b Body weight at harvest. c Waste weight. d Head weight. e Gutted head-on weight. f Body length at harvest. g Fillet weight. h Fillet yield. The x-axis presents genomic coordinates along chromosomes 1–23 in Nile tilapia. On the y-axis the negative logarithm of the SNPs associated p-value is displayed. The dashed black line represents the genome-wide significance threshold after Bonferroni correction (−log₁₀ (p-value > 7.21e-8)

Most of the lead SNPs were on LG01, LG03 and LG12 for ADG, BWH, WW, HW, HON and BLH. Some variants were common between body traits, such as two SNPs at positions 24,557,870 and 24,557,984 on LG12, that were the most significant SNPs (p-value < 9.893E-14) common in ADG, BWH, WW, HW, and HON. The lead SNPs for FW and FY were found on LG04 and LG13, and none of those were identified in other body traits (Table 3).

Candidate genes

The full list of genes located within 100 kb upstream and downstream of the lead SNP is available in additional file (Supplementary Table 2). Some lead SNPs for ADG, BWH, WW, HON, BLH are close to candidate genes, including collagen type IV alpha 1 chain (COL4A1) and growth differentiation factor 6 (GDF6) on LG16 and LG22, respectively, and ankyrin repeat and SOCS box containing 2 (ASB2) associated with BWH and HON, located on LG19. The genes intercepted by lead SNPs, located in exonic or intronic regions are shown in Table 4. Some of these genes have been associated with body traits in previous studies. For FW, the gene A disintegrin and metalloproteinase with thrombospondin motifs 9 (ADAMTS9), located in LG05, was intercepted by a SNP in an exon region at 29,062,243 bp. Two lead SNPs for WW, located on LG09, at positions 14,670,077 and 14,674,835 bp, intercepted introns of the gene solute carrier family 4 member 2 (SLC4A2). Intronic regions of α1,6-fucosyltransferase (FUT8) and the heart development protein with EGF like domains 1 (HEG1), located on LG15 and LG16, were intercepted by lead SNPs associated with ADG and FY, respectively. Two SNPs within nucleoporin 107 (NUP107), located on LG17, were associated with both BWH and HON, on positions 19,609,147 and 19,612,729 bp, respectively, with the first SNP hitting an intronic region and the second one located in an exon region. Others genes such as Coiled-Coil Domain Containing 102A (CCDC102A), SLIT-ROBO Rho GTPase Activating Protein 1 (SRGAP1), MutS Homolog 6 (MSH6) Myosin VI (MYO6), Myosin XVI (MYO16), and Kinectin 1 (KTN1) were intercepted by one or more lead SNPs, but no clear evidence of a close association with body size and growth-related traits has been reported.

Table 4 Genes intercepted by a lead SNP associated with body traits based on multi-trait GWAS in Nile tilapia

Discussion

We found moderate to high heritability values for ADG, BWH, WW, HW, HON, BLH, FW and FY, which is consistent with previous estimates for Nile tilapia calculated using pedigree and genomic methods [8, 9, 20, 21]. The additive genetic variance and heritability estimated for BWH using genotypes imputed to high-density genotypes increased about 15% in comparison to the value previously estimated for the same population using a 50 K SNP panel [20].

The use of genomic information can help in the identification of QTLs controlling complex traits which are economically important for aquaculture purposes, such as growth-related traits. Previous studies have identified loci and candidate genes associated with growth-related traits in aquaculture species [20, 23, 24, 26, 37, 38]. However, similar to what we found when using stGWAS (Supplementary Fig. 1), few or no markers surpassed the genome-wide significance threshold, or represented a small proportion of genetic variance for all body traits studied here. No studies have found evidence of major QTLs for growth-related traits, and GWAS signals were moderate even when a relatively large sample size (> 4600 animals) and more than 100 K markers were used, as in the case of GWAS for body weight in Atlantic salmon [23].

To increase the statistical power, in order to detect genetic association between SNPs and traits of interest, recent studies have used mtGWAS, which can leverage multiple summary statistics from GWAS performed on the same trait with different measures or different traits with a high genetic correlation among them [33, 39, 40]. We combined the use of genotypes imputed to high-density and the mtGWAS approach implemented in MTAG software to increase the statistical power and accuracy of QTL detection [33]. The imputation proceeded from a medium-density (50 K) SNP panel to high-density, where the markers from the reference dataset were previously selected based on quality control, and an expected accuracy of imputation higher than 0.80. The mtGWAS increases statistical power by using information from different traits that are genetically correlated with each other [33]. Here, the correlation of the overall SNP effects ranged from 0.86 to 1.00, except for the correlation between FY and all of the other traits, which ranged from 0.20 to 0.47 (Fig. 1), and the samples were overlapped for all traits. The better resolution of the genotypes imputed to high-density, combined with the power of the mtGWAS approach, lead to the detection of several novel significant markers not previously found when using stGWAS.

A difference in the number of significant SNPs between stGWAS and mtGWAS is expected given the substantial increase in statistical power which has been documented for the mtGWAS approach. However, it has also been shown that original associations detected by single-trait GWAS can disappear when running multi-trait GWAS. For instance, in the paper describing the application of mtGWAS [33], the increase of significant lead SNPs was from two up to four times higher when comparing mtGWAS against stGWAS. Nevertheless, there were also SNPs associated in the stGWAS analyses which were not found to be associated when running a multi-trait GWAS. If the SNP association is not confirmed by the mtGWAS, we may assume that the previous association identified by the stGWAS is spurious and interpretations on these unconfirmed associations have to be taken with caution.

We found numerous significant markers associated with body traits, dispersed in almost all linkage groups (LG; Fig. 2), probably due to the polygenic architecture of these traits in Nile tilapia. However, a major common association peak on LG12 was found for all traits analyzed, except for FW where the major peak was found on LG13; suggesting that part of the genetic variation that affects body traits might be explained by loci on these linkage groups. No gene was intercepted by the two most significant lead SNPs in this region, but a nearby gene on LG12, hydhroxysteroid 17-Beta Dehydrogenase 4 (HSD17B4), a possible regulator of muscle development in Berkshire pigs, has been reported to play an important role during the early stages of myogenesis when the expression of its mRNA is significantly high [41].

Some lead SNPs identified in this study were located close or intercepted several strong functional candidate genes associated with body and growth-related traits in previous studies. For instance, strong functional candidate genes were found in windows within 100 kb downstream and upstream from the lead SNP, such as COL4AI, located in LG16, associated with different body traits, including ADG, BWH, WW, HON and BLH. In catfish COL4A1 was identified within QTLs associated with body length and body length of the fish without the head. Collagen is an important component of the extracellular matrix of cartilage and bone, playing a key role in skeletal development [42]. We also found GDF6, located in LG22, which was associated with different traits including ADG, BWH, WW, HON, BLH. Mutations in GDF6 in zebrafish is related with reduced eye size and different skeletal defects [43]. In a study aimed to compare the orthologous sequences from 14 species (including human, mice, livestock, fugu, and zebrafish), the GDF6 gene was found to control developmental patterning of skeletal joints [44]. Inactivation of the GDF6 gene can cause defects in the joints, ligaments, and cartilage formation in mouse [45]. In addition, the ASB2 gene, located in chromosome 19, was associated with BWH and HON. In Atlantic salmon, the ASB2 gene is not involved in muscle differentiation but may play an important role in growth inhibition. The high expression of ASB2 observed in skeletal muscle of fasting fish is strongly downregulated in response to feeding [46].

We also found strong candidate genes intercepted by lead SNPs that may contribute to a better understanding of the biological mechanisms controlling body traits in Nile tilapia. Growth is considered a continuous function during the life of an animal and ADG is an important trait which can be targeted to select for rapid growth. ADG was previously reported as a selection criteria in a breeding program for Nile tilapia, which has applied selection for at least five generations [47]. We found a lead SNP associated with ADG on LG15, located in an intronic region of the FUT8 gene, which has been associated with ADG (from birth to six months-age) in a sheep population from Iran [48]. In mice, the disruption of the FUT8 gene induces severe growth retardation and early mortality during postnatal development [49,50,51,52]. The insulin-like growth factor binding protein-3 (IGFBP-3) has growth inhibitory effects, and the alteration in the function of low-density lipoprotein receptor-related protein-1 (LRP-1) is a result of the loss of core fucosylation that might cause an elevated serum concentration of IGFBP-3 in FUT8-null mice [52]. The loss of function of FUT8 has also been reported to be related to down-regulation of transforming growth factor-beta 1 (TGF-β1) receptor and epidermal growth factor (EGF) receptor, proteinase-activated receptor and integrin activity, which contributes to emphysema-like changes in the lung, and growth retardation in FUT8-null mice [51].

Two lead SNPs associated with BWH and HON were found on LG17, in an intronic and exonic region of the NUP107 gene which plays an important role in the development of vertebrate embryos. The zygotic deficiency of NUP107 in zebrafish embryos can result in loss of pharyngeal skeletons, degeneration of intestinal and retinal epithelia, and implications in cartilage and bone formation [53]. In senescent fibroblasts of humans and organs of aged mice, a decreased level of NUP107 is suggested to be involved in the hypo-responsiveness to growth [54].

The waste weight is the sum of the weight of the head, viscera, bones and fins, and has been suggested as an alternative phenotypic record to improve fillet yield through the application of various index (e.g. fillet to waste ratio). However, based on simulated data of ten generations of selection using real genetic parameters of five farmed fish populations, direct selection on fillet yield was generally the best approach to improve the trait [13]. The potential limitation for selecting against waste weight is the probability of decreasing the volume of essential organs in the visceral cavity. A negative genetic correlation (< − 0.52) between fillet yield and head, and bone development has been reported in rainbow trout [55]. We found a lead SNP that intercepts the SLC4A2 gene, a strong biological candidate for waste weight in Nile tilapia. The loss of function of this gene causes emaciation and achlorhydric condition [56], generating severe growth retardation, reduced osteoclast numbers and/or a reduction in osteoclast activity, resulting in osteopetrosis in mice [57]. Osteopetrosis is a skeletal disorder that can affect humans and animals, characterized by the formation of overly dense bones [57, 58]. In Red Angus cattle, a deletion mutation in SLC4A2 is associated with an osteopetrosis phenotype [58].

Fillet traits are key economic characteristics for aquaculture species and new insights regarding the underlying genetic variants controlling them can help to enhance yield. We found two lead SNPs associated with FW and FY, intercepting an exonic and intronic region of genes ADAMTS9 and HEG1, respectively. The ADAMTS9 gene is highly expressed during embryo development and continues to be expressed in adult tissues of mice [59, 60]. A significant expression of ADAMTS9 during skeletal development of mouse was suggested by Jungers et al. (2005), including mandible, ossification centers, initial condensation of mesenchyme to form the cartilage centers, perichondrium around formed cartilage, proliferative zones of cartilage and long bones. Skeletal development may be correlated with organic growth [61]. ADAMTS9 is responsible for the regulation of the epidermal growth factor receptor (EGFR) and TGF-β1 [62, 63]. Tang et al. (2019) [64] identified a 22-bp indel in ADAMTS9 associated with chest width, chest width index and chest circumference index, and a14-bp indel associated with height across the hip in cashmere goats, which suggests that ADAMTS9 is a functional molecular marker that can be used to improve growth traits in goats [62].

A lead SNP associated with FY intercepted the HEG1 gene, located in LG16. The HEG1 gene was initially reported to be responsible for regulating zebrafish heart growth and the development of heart and blood vessels development [65]. Recently the HEG1 gene was identified as one of several novel genes associated with human skeletal muscle growth, exhibiting a significant correlation with the percentage of change in lean mass [66, 67]. In a comparative transcriptomic analysis aimed to identify differentially expressed genes related to product performance and meat quality from the longissimus dorsi in sheep, Cheng et al. (2020) [68] identified six differentially expressed genes, including HEG1 [66].

Some genes such as MSH6, SRGAP1, MYO6, MYO16, KTN1 and other molecules presented in Table 4, are intercepted by lead SNPs and thought to be involved in different biological functions and conditions, such as some types of cancer [69, 70], hearing loss [71] and schizophrenia [72]. The mutation of MSH6, for example, may increase the risk of developing colorectal carcinomas [73, 74], and MYO16 appears to have an important role in neural development and the function of the nervous system [75]. The functional relationship between these genes and the variation in growth-related traits in Nile tilapia is unclear. Thus,  the function of the identified genes and their potential relationship with body traits in Nile tilapia must be better characterized.

Conclusions

We used dense genotypic information to refine association mapping analysis for body traits in Nile tilapia and found that mtGWAS provided substantial improvements in the number of significant SNPs identified when compared to stGWAS. These results confirm the increase of statistical power to identify trait-specific genetic associations in multi-trait analysis. Interestingly, we found several lead SNPs within or nearby genes related to cartilage, bone, skeletal growth and development in humans, mice, livestock and aquaculture species. These results can provide further knowledge and a better understanding of genetic variants and genes underlying complex body traits in Nile tilapia.

Material and methods

Animals and phenotypes

We used a total of 1309 phenotyped animals from 72 families (mean = 18, minimum = 7, and maximum = 25 animals per family) belonging to a breeding nucleus owned by Aquacorporación International group (ACI), Costa Rica. More details about the breeding program, the origin of the Nile tilapia population and production conditions are described in detail in previous studies [18, 20, 76]. Briefly, a mating design of two dams per sire was used to produce the 72 full-sib families. The eggs of each full-sib family were incubated and reared in separate hapas until individual tagging by using PIT (passive Integrated Transponder)-tags at an average weight and age of 13 g (SD = 8 g) and 104 days (SD = 18 days), respectively. After tagging, the fish were grown in excavated ponds for about 370 days until harvest. All animals were slaughtered by hypothermia in ice slurry at commercial processing plant, and different body traits were measured at harvest time: body weight at harvest (BWH in g), fillet weight (FW in g), waste weight (WW in g = BWH – FW), head weight (HW in g), gutted head-on weight (HON in g = BWH – Viscera), body length at harvest (BLH in cm), average daily gain (ADG in g = (BWH - body weight at tagging)/(age at harvest - age at tagging)), and fillet yield (FY in % = FW/BHW*100).

Genotypes and imputation to whole-genome sequences

Genomic DNA was extracted and purified from 1309 fin clip samples using the DNeasy Blood & Tissue Kit (QIAGEN) according to the manufacturer’s protocol. The samples were genotyped using a 50 K SNP Illumina BeadChip [18] and filtered using departure from Hardy-Weinberg Equilibrium (HWE, p-value 10^− 6), minor allele frequency (MAF < 0.01), and a genotyping call-rate for SNPs and samples of < 0.90. After quality control 29,587 SNPs and 1309 samples were retained.

Initially, 26.6 million non-redundant SNP variants were identified through Illumina HiSeq 2500 re-sequencing performed in 143 animals from the breeding nucleus owned by Aquacorporación International group (ACI), Costa Rica [18]. Quality control of WGS genotypes was performed using the following thresholds: HWE (p-value 10^− 8), MAF < 0.01 and call-rate for SNPs < 0.80. A total of 5,011,051 SNPs were retained after applying the filters described above. In order to estimate the overall accuracy of imputation and remove the variants with low imputation accuracy we used a five-fold cross validation scheme. Briefly, the 143 animals with data from the WGS-derived genotypes were randomly divided into five exclusive reference sets (80% of animals genotyped with ~ 5 million SNPs) and the remaining animals were used as the validation set (20% of animals with medium-density genotypes). The accuracy of imputation was estimated as the correlation between true and imputed genotypes (R² value). A total of 1,324,420 SNPs with R² value higher than 0.80 were used as the final ultra-dense SNP panel for imputation. The 143 re-sequenced animals and 1,324,420 SNPs were used as a reference dataset to impute the 1309 animals with medium-density SNP genotypes using the software FImpute v. 3.0 [77]. A post-imputation quality control excluded SNPs with MAF < 0.05 and HWE p-value < 10^− 8, resulting in a total of 992,494 SNPs available for downstream analyses.

Single-trait genome-wide association

The single-trait genome wide association analyses (stGWAS) were performed using the mlma option of the software GCTA v. 1.24 [78], which was used to apply the following linear mixed model:

$$ {\mathrm{y}}_{\mathrm{i}\mathrm{j}}=\upmu +{\mathrm{b}}_1\ast {\mathrm{a}\mathrm{ge}}_{\mathrm{j}}+{\mathrm{b}}_2\ast {\mathrm{SNP}}_{\mathrm{i}}+{\mathrm{a}}_{\mathrm{i}\mathrm{j}}+{\mathrm{e}}_{\mathrm{i}\mathrm{j}} $$

(1)

where y_ij is the phenotypic value of the j-th animal, μ is the fixed effect of the overall mean, b₁ and b₂ are the regression coefficients for age and the allele substitution effect for SNP, respectively, age_j is the age covariate of the j-th animal and SNP_i is the i-th SNP genotype of animal j, coded as 0, 1 and 2 for genotype A₁A₁, A₁A₂ and A₂A₂, respectively, a_ij is the random polygenic effect of the j-th animal $\sim \mathrm{N}\left(0,\mathbf{G}{\upsigma}_{\mathrm{a}}^2\right)$, with G representing the genomic relationship matrix (GRM) calculated using the imputed genotypes and ${\upsigma}_{\mathrm{a}}^2$ the genetic variance [78, 79], and e_ij is the random residual effect $\sim \mathrm{N}\left(0,\mathbf{I}{\upsigma}_{\mathrm{e}}^2\right)$, with I representing an identity matrix and ${\upsigma}_{\mathrm{e}}^2$ the residual variance. The GRM is calculated based on the relationship from a genome-wide sample of SNPs obtained by using a common-sense weighting scheme [78]. The GRM restricted maximum likelihood (GREML) [78] implemented in GCTA was used to estimate the genetic and residual variances. Heritability (h²) was calculated as h² = ${\upsigma}_{\mathrm{a}}^2$ /$\left({\upsigma}_{\mathrm{a}}^2+{\upsigma}_{\mathrm{e}}^2\right).$ For each SNP, the allele substitution effect and its p-value were also estimated using GCTA.

Multi-trait genome-wide association

The summary statistics from stGWAS were used as input for the multi-trait analysis of GWAS (mtGWAS) performed using the software MTAG v0.9.0 [33]. In MTAG, the SNP effect estimated for each trait can be improved when different traits that are correlated are included in the analysis. This multi-trait approach can increase the power to detect loci in any of the traits assessed. The first step of MTAG is to filter variants based on discarding non common SNPs, duplicated SNPs, or SNPs with strand ambiguity. In our study, out of the 992,494 SNPs available after imputation and initial quality control, a total of 183,401 SNPs with strand ambiguity were filtered out. The remaining 809,093 SNPs were used for mtGWAS analyses. A bivariate linkage disequilibrium (LD) score regression was used, thus summary statistics do not need to come from independent samples [33]. The MTAG output consists of a file per trait with updated results of SNP effects and p-values from a mtGWAS, which can be interpreted in the same way as stGWAS. Significance thresholds were determined for both single-trait and mtGWAS using Bonferroni correction (0.05/ number of SNPs).

To calculate how much larger the stGWAS sample size would have to be to give the same mean χ² statistics than mtGWAS, the following equation was used [33]:

$$ {\mathrm{N}}_{\mathrm{GWAS}\ \mathrm{equivalent}}={\mathrm{N}}_{\mathrm{GWAS}}\frac{\overline{\upchi_{\mathrm{mtGWAS}}^2}}{\overline{\upchi_{\mathrm{stGWAS}}^2}} $$

(2)

where, ${\upchi}_{\mathrm{mtGWAS}}^2$ and ${\upchi}_{\mathrm{stGWAS}}^2$ are the mean χ² statistic for mtGWAS and stGWAS results, respectively, and N_GWAS is the number of actual sample size in stGWAS (1309 animals).

Identification of QTL and candidate genes

The most significant SNP per chromosome per each trait, detected using mtGWAS, was selected as the lead SNP and furtherly used to search for candidate genes based on proximity to the variant. Genes located within 100 kb upstream and downstream of the lead SNP were considered putative candidate genes associated with the trait. The gene search was performed using BLAST (Basic Local Alignment Search Tool) against the latest version of the Oreochromis niloticus reference genome (O_niloticus_UMD_NMBU [15]), which is publicly available at NCBI (GenBank assembly accession GCA_001858045.3).