An F2 design resource population was developed between 2000 and 2006 (Guo et al., 2009) as follows: two White Duroc sires and 17 Erhualian dams were mated to produce F1 animals, from which 9 F1 boars and 59 F1 sows were intercrossed (avoiding full-sib mating) to produce 967 F2 males and 945 F2 females (in total n = 1,912) in six batches. All the F2 animals were kept under standard indoor conditions at the experimental farm of Jiangxi Agricultural University (China). Then the F2 piglets were weaned at 46 days, and males were castrated at 90 days. At 240 ± 6 days of age, 1,030 F2 animals including 549 gilts and 481 barrows were slaughtered at 70–120 kg live weight.

Genomic DNA was isolated from ear tissue with a standard phenol/chloroform extraction method. All DNA samples were diluted to a final concentration of 50 ng/μl in 96-well plates. In total, 933 F2 were genotyped with the Illumina PorcineSNP60 BeadChip on an iScan System (Illumina, United States) following the manufacturer’s protocol (Ramos et al., 2009). Quality control procedures were implemented by PLINK (version 1.07) (Chang et al., 2015). Briefly, SNPs with unspecific positions on the genome build 10.2, a call rate lower than 90%, and a minor allele frequency (MAF) lower than 1% were eliminated, and animals with a missing typing rate higher than 10% were also removed. In total, 374 phenotypes were measured on the individuals of the F2 population, including carcass traits, reproductive traits, immune traits, meat traits, growth traits, and epigenetic traits (see Additional file 1: Supplementary Table S1). These 374 traits were then divided into three groups according to their heritability, i.e., 68 traits with high heritability ( h 2 > 0.4 ), 148 traits with a moderate heritability ( 0.2 < h 2 < 0.4 ), and 158 with low heritability ( h 2 < 0.2 ).

We used the GEMMA software to calculate the substitution effects of 14,320,159 SNPs on the 374 traits included in the standard linear model (Zhou and Stephens, 2012). Sex was included as a fixed effect, and heritability was estimated by using the − l m m procedure implemented in GEMMA. Population stratification was adjusted by including a genomic relationship matrix. Briefly, the model was as follows: y = W a + X β + u + e   ; u ∼ M V N n ( 0 , σ u 2 K ) , e ∼ M V N n ( 0 , σ e 2 I n ) where y is an n element vector of phenotypic values, all the traits were normalized before calculation so that the substitution effects were comparable among all the phenotypes, W is a design matrix of covariates, a is a vector of fixed effects, X is a vector of genotypes at each locus, β is the effect size of SNPs, and u is the vector of random effects following a multivariate normal distribution M V N n ( 0 , σ u 2 K ) , e is the vector of errors following M V N n ( 0 , σ e 2 I n ) , σ u 2 and σ e 2 are polygenic variance and residual variance, respectively, which are estimated based on the REML average information (AI) algorithm. K is a known kinship matrix estimated from genome sequence variants, and I n being an n × n identity matrix.

Three genotypes “AA,” “Aa,” and “aa” were assumed each locus and were represented by 0, 1, and 2, respectively, with p and q the frequencies of alleles “A” and “a,” respectively. Assuming that the effect value of this locus is estimated as β (with no dominance), the additive genetic variance can be expressed as 2 p q β 2 (Park et al., 2011). In the group of traits with high heritability, all the loci for each phenotype were put together and ranked by additive genetic variance from large to small. And all the ordered loci were equally divided into four groups. For each group, the proportion of the sum of the additive genetic variances of all loci to the total additive genetic variance (or variance-ratio thereafter) was calculated, equal to a 1 , b 1 , c 1 , and d 1 , respectively (Subsequently called variance ratio). Similarly, the same method was used for the groups of traits with a moderate heritability and a low heritability, resulting in a 2 , b 2 , c 2 , d 2 , and a 3 , b 3 , c 3 , d 3 , respectively. Therefore, the four expected variances in each of the three groups can be expressed as:

Group of traits with high heritability:   δ 1 2 = a 1 V g γ M ; δ 2 2 = b 1 V g γ M ; δ 3 2 = c 1 V g γ M ; δ 4 2 = d 1 V g γ M

Group of traits with a moderate heritability:   δ 1 2 = a 2 V g γ M ; δ 2 2 = b 2 V g γ M ; δ 3 2 = c 2 V g γ M ; δ 4 2 = d 2 V g γ M

Group of traits with a low heritability:   δ 1 2 = a 3 V g γ M ; δ 2 2 = b 3 V g γ M ; δ 3 2 = c 3 V g γ M ; δ 4 2 = d 3 V g γ M Where V g and M is the additive variance and the number of markers, respectively. γ is set to 0.25.

The linear model for genomic prediction was as follows: y = Xb + Bg + e Where n individuals and m SNPs were assumed. Thus, y is a n × 1 vector of phenotypes recorded; b is the vector of fixed effects; g is a m × 1 vector of additive SNP effects; e is a vector of residual errors; X is the design matrix for fixed effects; and B is standardized design matrix for additive SNP effects (coded as 0 for genotype “AA,” 1 for “Aa” and 2 for “aa,” respectively).

In this study, the prior distribution with four zero-mean normal distributions was written as a function of prior distributions of SNP variance as determined above: π ( g j ) =   γ ϕ ( g j | 0 , δ 1 2 ) + γ ϕ ( g j | 0 , δ 2 2 ) + γ ϕ ( g j | 0 , δ 3 2 ) + γ ϕ ( g j | 0 , δ 4 2 ) ; γ = 0.25 (1) π ( g j ) is the univariate normal distribution, the effects of SNPs are obtained by using the Iterative Conditional Expectation algorithm (Meuwissen et al., 2009). In brief, assume that E ( g j | y − j ) is estimated, the current effects of all the other SNPs are used to calculate the y − j as follows: y − j = y − X b − ∑ k ≠ j B k g k where B k is a vector from the k t h column of B, the expectation of SNP effect, E ( g j | y − j ) , is then estimated by a Bayesian model in the next round: E ( g j | y − j ) = ∫ − ∞ + ∞ g j f ( g j | y − j ) d g j = ∫ − ∞ + ∞ g j f ( y − j | B j g j , I δ e 2 ) π ( g j ) d g j f ( y − j ) = ∫ − ∞ + ∞ g j f ( y − j | B j g j , I δ e 2 ) π ( g j ) d g j ∫ − ∞ + ∞ f ( y − j | B j g j , I δ e 2 ) π ( g j ) d g j (2) f ( y − j ) is a marginal distribution function of y − j and can be calculated using the law of total cumulance: ∫ − ∞ + ∞ f ( y − j | B j g j , I δ e 2 ) π ( g j ) d g j . Calculating f ( y − j | B j g j , I δ e 2 ) is computationally demanding because it is a multivariate normal density, which involves calculating the determinant and the inverse of the variance-covariance matrix for the data y − j . Therefore, we simplified the derivation by using a univariate normal densities   f ( Y | g j , δ 2 ) to replace f ( y − j | B j g j , I δ e 2 ) , where Y = ( B j ' B j ) − 1 B j ' y − j and δ 2 = ( B j ' B j ) − 1 δ e 2 ; details of the derivation process are as follows: f ( y − j | B j g j , I δ e 2 ) ∝ e x p [ − ( y − j − B j g j ) , ( y − j − B j g j ) 2 δ e 2 ] = exp ( − y − j , y − j − 2 B j , y − j g j + B j , B j g j 2 2 δ e 2 ) = exp [ − g j 2 − 2 ( B j , B j ) − 1 B j , y − j g j + ( B j , B j ) − 1 y j , y − j 2 ( B j , B j ) − 1 δ e 2 ] = e x p { − [ g j − ( B j , B j ) − 1 B j , y − j ] 2 − [ ( B j , B j ) − 1 B j , y − j ] 2 + ( B j , B j ) − 1 y j , y − j 2 ( B j , B j ) − 1 δ e 2 } ∝ e x p { − [ g j − ( B j , B j ) − 1 B j , y − j ] 2 2 ( B j , B j ) − 1 δ e 2 } ∝ e x p [ − ( g j − Y ) 2 2 δ 2 ] ∝ f ( Y | g j , δ 2 )

Therefore, the equation E ( g j | y − j ) can be written as: E ( g j | y − j ) = ∫ − ∞ + ∞ g j f ( Y | g j , δ 2 ) π ( g j ) d g j ∫ − ∞ + ∞ f ( Y | g j , δ 2 ) π ( g j ) d g j (3) the numerator of Eq. 3 can be broken down into four terms combined with Eq. 1 as follows: γ ∫ − ∞ + ∞ g j f ( Y | g j , δ 2 ) ϕ ( g j | 0 , δ 1 2 ) d g j + γ ∫ − ∞ + ∞ g j f ( Y | g j , δ 2 ) ϕ ( g j | 0 , δ 2 2 ) d g j + γ ∫ − ∞ + ∞ g j f ( Y | g j , δ 2 ) ϕ ( g j | 0 , δ 3 2 ) d g j + γ ∫ − ∞ + ∞ g j f ( Y | g j , δ 2 ) ϕ ( g j | 0 , δ 4 2 ) d g j (4)

The first term of Eq. 4 can be derived as follows: γ ∫ − ∞ + ∞ g j f ( Y | g j , δ 2 ) ϕ ( g j | 0 , δ 1 2 ) d g j =   γ ∫ − ∞ + ∞ g j 1 δ 2 π e x p [ − ( Y − g j ) 2 2 δ 2 ] 1 δ 1 2 π e x p [ − g j 2 2 δ 1 2 ] d g j = γ 2 π ∫ − ∞ + ∞ g j δ δ 1 2 π e x p [ − ( Y − g j ) 2 2 δ 2 − g j 2 2 δ 1 2 ] d g j = γ 2 π ∫ − ∞ + ∞ g j δ δ 1 2 π e x p [ − ( g j 2 − 2 Y δ 1 2 δ 2 + δ 1 2 g j + Y 2 δ 1 2 δ 2 + δ 1 2 ) ( δ 2 + δ 1 2 ) 2 δ 2 δ 1 2 ] d g j =   γ 2 π ∫ − ∞ + ∞ g j δ δ 1 2 π e x p [ − ( g j − Y δ 1 2 δ 2 + δ 1 2 ) 2 ( δ 2 + δ 1 2 ) 2 δ 2 δ 1 2 − Y 2 2 ( δ 2 + δ 1 2 ) ] d g j = γ 2 π e x p [ − Y 2 2 ( δ 2 + δ 1 2 ) ] 1 δ 2 + δ 1 2 ∫ − ∞ + ∞ g j δ δ 1 δ 2 + δ 1 2 2 π e x p [ − ( g j − Y δ 1 2 δ 2 + δ 1 2 ) 2 2 ( δ δ 1 δ 2 + δ 1 2 ) 2 ] d g j

Here, the last term of this formula equals Y δ 1 2 δ 2 + δ 1 2 as it can be taken as calculating the expected value of g j in the normal distribution with mean Y δ 1 2 δ 2 + δ 1 2 , and variance δ 2 δ 1 2 δ 2 + δ 1 2 . Therefore, the first term of Eq. 4 can be written as follows: γ 2 π exp [ − Y 2 2 ( δ 2 + δ 1 2 ) ] 1 δ 2 + δ 1 2 Y δ 1 2 δ 2 + δ 1 2

Here, the derivation process for the remaining terms of Eq. 4 was the same as for this term, and therefore, the final form of the numerator of Eq. 3 is: γ 2 π exp [ − Y 2 2 ( δ 2 + δ 1 2 ) ] 1 δ 2 + δ 1 2 Y δ 1 2 δ 2 + δ 1 2 + γ 2 π exp [ − Y 2 2 ( δ 2 + δ 2 2 ) ] 1 δ 2 + δ 2 2 Y δ 2 2 δ 2 + δ 2 2 + γ 2 π exp [ − Y 2 2 ( δ 2 + δ 3 2 ) ] 1 δ 2 + δ 3 2 Y δ 3 2 δ 2 + δ 3 2 + γ 2 π exp [ − Y 2 2 ( δ 2 + δ 4 2 ) ] 1 δ 2 + δ 4 2 Y δ 4 2 δ 2 + δ 4 2 (5) there is no g j in the integrand of the denominator in Eq. 3 compared to that of the numerator. Therefore, it should calculate the cumulative probability from − ∞ to + ∞ , but not calculate the expected value, and this value is 1. Thus, the denominator in Eq. 3 can be written as: γ 2 π exp [ − Y 2 2 ( δ 2 + δ 1 2 ) ] 1 δ 2 + δ 1 2 + γ 2 π exp [ − Y 2 2 ( δ 2 + δ 2 2 ) ] 1 δ 2 + δ 2 2 + γ 2 π exp [ − Y 2 2 ( δ 2 + δ 3 2 ) ] 1 δ 2 + δ 3 2 + γ 2 π e x p [ − Y 2 2 ( δ 2 + δ 4 2 ) ] 1 δ 2 + δ 4 2 (6) thus, the final form for Eq. 3 is derived: E ( g j y − j ) = Y 2 δ 2 + δ 1 2 + exp [ Y 2 2 ( δ 2 + δ 1 2 ) − Y 2 2 ( δ 2 + δ 2 2 ) ] δ 2 + δ 1 2 δ 2 + δ 2 2 Y δ 2 2 δ 2 + δ 2 2 + exp [ Y 2 ( δ 2 + δ 1 2 ) − Y 2 ( δ 2 + δ 3 2 ) ] δ 2 + δ 1 2 δ 2 + δ 3 2 Y δ 3 2 δ 2 + δ 3 2 + exp [ Y 2 2 ( δ 2 + δ 1 2 ) − Y 2 2 ( δ 2 + δ 4 2 ) ] δ 2 + δ 1 2 δ 2 + δ 4 2 Y δ 4 2 δ 2 + δ 4 2 1 + exp [ Y 2 2 ( δ 2 + δ 1 2 ) − Y 2 2 ( δ 2 + δ 2 2 ) ] δ 2 + δ 1 2 δ 2 + δ 2 2 + exp [ Y 2 2 ( δ 2 + δ 1 2 ) − Y 2 2 ( δ 2 + δ 3 2 ) ] δ 2 + δ 1 2 δ 2 + δ 3 2 + exp [ Y 2 2 ( δ 2 + δ 1 2 ) − Y 2 2 ( δ 2 + δ 4 2 ) ] δ 2 + δ 1 2 δ 2 + δ 4 2

Here, the fixed effects are estimated at each iteration by the formula: b ^ = ( X ′ X ) − 1 X ′ ( y − B g ^ ) . Convergence of solutions at the t th iteration was judged based on the formula ( G t − G t − 1 ) ′ ( G t − G t − 1 ) ( G t ) ′ G t < 10 − 8 , where G = ( b ^ ' g ^ ′ ) ′ . It ends at the iteration when all the SNPs have been calculated once.

In the following analysis, we used GBLUP (Meuwissen et al., 2001; VanRaden, 2008), SSgblup (Legarra et al., 2009; Christensen and Lund, 2010), FMixFN, MIX (Xavier et al., 2019), BayesR (Kemper et al., 2015), BayesA, and BayesB with respective model fittings to compare their performance, the variance components were pre-estimated using the mixed model. The details of these analyses were as follows:

GBLUP: GBLUP was used to estimate the effects of the markers by BLUP, assuming that each marker explains an equal proportion of the total genetic variance. The software GEMMA was used to implement the GBLUP calculation process (Zhou and Stephens, 2012).

SSgblup: Single-step genomic BLUP (SSgblup), which was developed by Aguilar et al. (2010) and Christensen and Lund (2010), opened the way to perform genomic prediction using phenotype, pedigree, and genomic information simultaneously on both genotyped and non-genotyped individuals via a combined relationship matrix (H). Implementation of SSgblup is completed by the R package “Hiblup” (https://hiblup.github.io/).

MIX: the MIXTURE model assumed that the marker effects came from a mixture of two distributions: one distribution with large variance (accommodating large marker effects) and one with small variance (accommodating small marker effects). The distribution to which the marker belongs is sampled from the Bernoulli distribution. The variances of the two distributions underlying the mixture are estimated using a noninformative chi-square distribution. Implementation of MIX is completed by the R package “VIGoR” (https://cran.r-project.org/web/packages/VIGoR/index.html).

BayesR: BayesR starts the hierarchical model and poses a mixture of four zero-mean normal distributions as a conditional prior for a specific SNP effect. We use BayesR software to implement the calculation process (https://cnsgenomics.com/software.html).

BayesA: BayesA assumes that the distribution of SNP effects follows a Student’s t-distribution. Mathematically, it is assumed that each SNP effect comes from a normal distribution but σ 2 can be varied among the SNPs because the t-distribution is not easy to incorporate into a prediction of the marker. A scaled inverted chi distribution, X 2 ( ν , S ) is usually used as prior for the variance components.

BayesB: The prior distribution of BayesB is a mixture distribution with some SNPs with zero effects and the rest with a t-distribution, and the prior hypothesis of the SNP with non-zero effect is the same as BayesA. The implementation of BayesA and BayesB is completed by the R package “BGLR” (Perez and de los Campos, 2014). All MCMC sampling was run for 50,000 cycles, and the first 20,000 cycles were discarded as burn-in for BayesR, BayesA, and BayesB.

To test the performance of FMixFN in terms of predictive ability and calculation time, we did the following. Firstly, the variance ratio of the prior distribution of FMixFN was assumed to be random, then two traits were selected to verify the unbiasedness, and the compatibility of the variance ratio estimated in the F2 population. The specific assumptions of the variance ratio are shown in Additional file 1: Supplementary Table S2, the first one is to average all variances, i.e. to set the classification with the largest variance ratio at 50%, and the second was to centralize all variances, i.e., the classification with the largest variance ratio is assumed to be 90%. Secondly, we compared the predictive ability and calculation time of FMixFN and other mainstream genomic selection methods and selected two phenotypes with different genetic structures and different heritabilities from each group for cross-validation, one trait is controlled by numerous polygenic genes and the other one is controlled by several loci with large variances. Supplementary Figure S1 (see Additional file 2: Supplementary Figure S1) shows that traits 1, 3, and 5 were controlled by SNPs with large variances, and traits 2, 4, and 6 were controlled by many SNPs, each with a very small effect. After quality control, the remaining number of individuals with the six traits were 839, 832, 834, 840, 784, and 838, respectively, with 33,901, 33,893, 33,891, 33,891, and 33,894 SNPs, respectively, the heritability of each trait was 0.600, 0.560, 0.380, 0.369, 0.145, and 0.107, respectively. Details on the number of individuals, number of SNPs, and heritability estimates are shown in Additional file 1: Supplementary Table S3. In addition, we also evaluated the stability of FMixFN by using data from the duroc experiment and Asian rice experiments, respectively, more specific information from this population is also shown in Additional file 1: Supplementary Table S3. Finally, to demonstrate that FMixFN can perform genomic prediction analysis based on large-scale sample data with no data overflow error, our study simulated 20 sets of sample data using QMsim software (Sargolzaei and Schenkel, 2009), which contains 10,000, 20,000, ……, 190,000, and 200,000 individuals, respectively. Each set of data was obtained through eight generations of mating, combining genotype and phenotype data from generations 3–8, and determining the number of individuals per generation by parameter setting. Genomic information of each individual was set with 10 chromosomes, each chromosome is set 100 cM long and including 101 markers and 100 QTLs, respectively, with a marker mutation rate of 2.5 and QTL mutation rate of 3. Genomic prediction by a replicated training-testing method was used to evaluate the predictive results. Cross-validation of nine replicates was performed. All individuals were randomly and evenly divided into nine groups. In each replicate, one of the groups was selected as the testing data set while the remaining eight groups were used as the training data set, and the results of each cross-validation are shown in Additional file 1: Supplementary Table S4. Predictive ability is defined as the correlation between GEBV and the phenotypes adjusted for the covariates ( y − X u ^ ) (Meuwissen et al., 2001).

In this study, all traits of the F2 population were divided into three groups based on the heritability of the traits: high, moderate, and low. For the group of traits with high heritability, the calculated a 1 , b 1 , c 1 , and d 1 were equal to 0.8752, 0.0958, 0.0256, and 0.0032, respectively. For the group of traits with a moderate heritability, the calculated a 2 , b 2 , c 2 , and d 2 were equal to 0.8367, 0.1246, 0.0342, and 0.0043, respectively. And for the group of traits with low heritability, the calculated a 3 , b 3 , c 3 , and d 3 were equal to 0.8225, 0.1413, 0.0324, and 0.0036, respectively. Those parameters were composed in the procedure of FMixFN, as FMixFN starts running, the program determines which group of variances is calculated based on the heritability of the experimental trait.

In this study, we selected phenotype 3 and phenotype 4 to verify the unbiasedness of the variance ratio of the prior distribution. When the variance ratio was assumed to be 0.5, 0.25, 0.125, and 0.125, the predictive ability of phenotype 3 and phenotype 4 are 0.4773 and 0.4911, respectively. When the variance ratio is assumed to be 0.9, 0.005, 0.045, 0.005, the predictive ability of each of these phenotypes are 0.4789 and 0.4911, respectively. In contrast, the predictive ability of these two phenotypes estimated by using the original parameters is 0.4787 and 0.4913, respectively.

The predictive ability of each of the six F2 traits under the six predictive methods is shown in Figure 1. For phenotypes 1, 3, and 5, the predictive ability with BayesR, BayesA, and BayesB was, respectively, 0.0377, 0.0308, and 0.0374 higher than that of FMixFN, and the predictive ability of FMixFN was slightly better than of GBLUP by 0.0045, 0.0013, and 0.0103, respectively. For those three phenotypes, there is almost no difference between SSgblup and FMixFN in predictive ability, and FMixFN performs better than SSgblup for phenotype5. For phenotypes 2, 4, and 6, FMixFN performed best for phenotype 4, with a predictive ability 0.0129 higher that of BayesR, BayesA, and BayesB. For phenotype 2, the predictive ability of the five software was similar and was highest with BayesA but only 0.0053 higher than that of FMixFN. For phenotype 6, FMixFN ranked second in the predictive ability, just 0.0027 lower than that of SSgblup. It was worth mentioning that the predictive ability of FMixFN was slightly better than that of other ICE-based Bayesian mixture regression (MIX) by 0.0221, 0.0116, 0.0801, and 0.0263 for phenotype 1, 4, 5, and 6, respectively. The specific information of the predictive ability was also shown in Table 1. Table 2 reports the predictive ability performances of FMixFN and other methods using the Duroc and rice datasets. In the Duroc population, the prediction accuracies were 0.3655, 0.3300, 0.3998, 0.3476, and 0.3589 for FMixFN, GBLUP, BayesR, BayesA, and BayesB, respectively. From the mean value, we found that FMixFN performed slightly worse than BayesR, but outperformed GBLUP, BayesA, and BayesB. In general, the MCMC-based Bayes genome selection algorithm showed some advantages in the traits controlled by several major QTLs, which explained a large proportion of phenotypic variance in some SNPs, while FMixFN performs better than GBLUP. FMixFN is slightly better than some other mainstream methods when traits follow a polygenic model. In this study, we measured the calculation speed of five methods as the average time necessary for the first cross-validation of the six traits. As shown in Figure 2A, the average calculation time was 0.54, 0.37, 0.29, 30.2, 42, and 50.5 min for FMixFN, GBLUP, MIX, BayesR, BayesA, and BayesB, respectively.

The computational time per iteration of FMixFN increases almost linearly as the number of individuals increases in the reference group, as shown in Figure 2B. Data reading time also increased linearly as the amount of sample data increased. Through simulation studies, we also found that FMixFN can calculate GEBV for 200,000 large samples without data overflow. Data overflow usually occurred in exponential functions, where data overflows or underflows could occur when the exponential part of the exponential function was very large or very small. The exponential part of Eq. 5 and Eq. 6 was exp [ − Y 2 2 ( δ 2 + δ i 2 ) ] , in which Y = ( B j , B j ) − 1 B j , y − j , B j was the j column of the accompanying matrix of SNPs. When B j growth unlimited, the limit of exp [ − Y 2 2 ( δ 2 + δ i 2 ) ] is negative infinity, and data underflows will occur. In this study, we found FMixFN could calculate GEBV for 200,000 samples without data overflow. So FMixFN has the potential for use on large-scale sample data.

In conclusion, when there are fewer individuals in the reference group, the computational speeds for ICE-based FMixFN are on the same order of magnitude with GBLUP, and they were much faster than the MCMC-based Bayesian methods. When the number of individuals in the reference dimension of the kinship matrix exceeds hundreds of thousands or even millions, the process to inverse the matrix becomes very difficult for the direct method of genomic selection. FMixFN has excellent computational efficiency and can handle large-scale sample data.