C Y Chen

I Misztal

I Aguilar

† A Legarra

W M Muir

Effect of different genomic relationship matrices on accuracy and scale 1

Keywords: chicken, genetic evaluation, genomic prediction

Phenotypic data on BW and breast meat area were available on up to 287,614 broilers. A total of 4,113 birds were genotyped for 57,636 SNP. Data were analyzed by a single-step genomic BLUP (ssGBLUP), which accounts for all phenotypic, pedigree, and genomic information. The genomic relationship matrix (G) in ssGBLUP was constructed using either equal (0.5; GEq) or current (GC) allele frequencies, and with all SNP or with SNP with minor allele frequencies (MAF) below multiple thresholds (0.1, 0.2, 0.3, and 0.4) ignored. Additionally, a pedigree-based relationship matrix for genotyped birds (A 22 ) was available. The matrices and their inverses were compared with regard to average diagonal (AvgD) and off-diagonal (AvgOff) elements. In A 22 , AvgD was 1.004 and AvgOff was 0.014. In GEq, both averages decreased with the increasing thresholds for MAF, with AvgD decreasing from 1.373 to 1.020 and AvgOff decreasing from 0.722 to 0.025. In GC, AvgD was approximately 1.01 and AvgOff was 0 for all MAF. For inverses of the relationship matrices, all AvgOff were close to 0; AvgD was 2.375 in A 22 , varied from 11.563 to 12.943 for GEq, and increased from 8.675 to 12.859 for GC as the threshold for MAF increased. Predictive ability with all GEq and GC was similar except that at MAF = 0.4, they declined by 0.01 for BW and improved by 0.01 for breast meat area. Compared with BLUP, EBV in the ssGBLUP were, on average, increased by up to 1 additive SD greater with GEq and decreased by 2 additive SD less with GC. Genotyped animals were biased upward with GEq and downward with GC. The biases and differences in EBV could be controlled by adding a constant to GC; they were eliminated with a constant of 0.014, which corresponds to AvgOff in A 22 . Unbiased evaluation in the ssGBLUP may be obtained with GC scaled to be compatible with A 22 . The reduction of SNP with small MAF has a small effect on the real accuracy, but it may falsely increase the estimated accuracies by inversion.

INTRODUCTION

Genomic selection can be carried out either by estimating the effects of individual SNP or by using BLUP with a genomic relationship matrix, G (VanRaden et al., 2009b). For traits without large major genes, most methods provide similar accuracies [START_REF] Hayes | Invited review: Genomic selection in dairy cattle: Progress and challenges[END_REF]. Usually, not all animals have genotypes. A methodology exists to combine pedigree and genomic relationships [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF]. A single-step approach that uses the combined relationships allows for a simple evaluation that takes all available phenotypic, pedigree, and genomic information and has no limitations on models for analyses [START_REF] Vitezica | One-step vs. multistep methods for genomic prediction in presence of selection[END_REF]Chen et al., 2011).

One of the issues plaguing genomic predictions is bias. For instance, in dairy cattle the bulls evaluated genomically seem to have EBV biased upward ( Van-Raden et al., 2009a). Bias from genomic evaluation has been reported in other studies [START_REF] Vitezica | One-step vs. multistep methods for genomic prediction in presence of selection[END_REF][START_REF] Forni | Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information[END_REF]. With bias, an accurate comparison of young and older animals is difficult.

The bias can partially be due to a genomic relationship matrix that is suboptimal. Such a matrix can be constructed by assuming different gene frequencies and can possibly be shifted and scaled, and SNP with small minor allele frequencies (MAF) may be ignored ( Van-Raden, 2008). In Holsteins, [START_REF] Vitezica | One-step vs. multistep methods for genomic prediction in presence of selection[END_REF] found that the use of the equal allele frequencies resulted in the greatest realized accuracy and the smallest bias. In pigs, [START_REF] Forni | Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information[END_REF] found that the best matrix used current allele frequencies and was scaled similarly to the numerator relationship matrix for genotyped animals.

The first goal of this study was to evaluate the effect of different options to construct G on statistics [means of the diagonal (AvgD) and off-diagonal (AvgOff)] of G and its inverse, and the accuracy of EBV and biases. The second goal was to investigate the origin of biases with some G. Broiler chickens from a pure line with 1.5% animals genotyped for 57,636 SNP were used to evaluate the different corrections on G.

MATERIALS AND METHODS

Animal Care and Use Committee approval was not obtained for this study because the data were obtained from an existing database.

Data

Phenotypic data of broiler chickens from a single pure line across 4 generations were provided by Cobb-Vantress Inc. (Siloam Springs, AR). Traits analyzed included BW at 6 wk (BW, 100 g) and ultrasound area of breast meat (BM, cm 2 ). A total of 287,614 broilers were available for phenotypes. Body weight was recorded for all birds, whereas BM was recorded for only 69,057 broilers. There were 289,994 birds in the pedigree.

A total of 4,113 birds were genotyped for 57,636 SNP based on the SNP panel developed by [START_REF] Groenen | A highdensity SNP-based linkage map of the chicken genome reveals sequence features correlated with recombination rate[END_REF]. Descriptions of phenotypic records for genotyped birds and all birds are shown in Table 1.

The data set was split into training and validation data sets for the genetic prediction. The training data set consisted of 234,806 records in the first 3 generations, whereas the validation data set contained 806 genotyped birds in the last generation.

Model

In our previous study (Chen et al., 2011), the singletrait model was implemented for a similar data set. In the current study, the multiple-trait model was

y t = X t b t + Z t u t + W t mp t + e t ,
where t is 1 for BW and 2 for BM; y is the vector of observations for traits; b is the vector of fixed effects, including contemporary group (house-hatch) and sex; u is the vector of random additive genetic effects, combining polygenic and genomic breeding values; mp is the vector of random maternal permanent environmental effects; X, Z, and W are incidence matrices; and e is the vector of random residuals.

Effects of the multiple-trait model used in the present paper were the same as those used in the single-trait model by Chen et al. (2011), except that the maternal permanent environmental effect was also included for BM in the present paper. In Chen et al. ( 2011), the abnormality of accuracy obtained with a single-trait model may have been due to selection bias resulting from multiple-trait selection. For the next run of selection using the multiple-trait model, the anomaly was absent. Furthermore, the accuracy for a trait with missing records, BM, increased slightly with the multiple-trait model. Therefore, the multiple-trait model was used in the current study instead of the single-trait model.

For the regular BLUP analysis, the (co)variance matrix was assumed to be var ,

u mp e A D I Q I R             =             ⊗ ⊗ ⊗ 0 0 0 0 0 0
where I is an identity matrix, A is a numerator relationship matrix, and D, Q, and R are matrices of the additive, maternal permanent, and residual variances, respectively. In the single-step genomic BLUP (ssG-BLUP) proposed by [START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF], the A matrix was replaced by the H matrix with the following inverse [START_REF] Vitezica | One-step vs. multistep methods for genomic prediction in presence of selection[END_REF]:

H A G A - - - - = + -           1 1 1 22 1 0 0 0 ,
where H is a modified relationship matrix incorporating genomic information, A 22 1 -corresponds to the inverse of the numerator relationship matrix for genotyped birds, and G is a genomic relationship matrix. Genomic relationship matrices were constructed as in VanRaden ( 2008) by using either equal (GEq) or current (GC) allele frequencies:

G QQ = ∑ - ' ( ) , 2 1 p p j j
where Q is an incidence matrix for SNP effects with elements

q ij j j p if homozygous 11 p if heterozygous 12 or 21 p = - - - 0 2 1 2 2 2 j j if homozygous 22 ,             
for bird i and SNP j with allele frequency p j . These matrices were constructed with the MAF ranging from none to 0.4. The purpose of the thresholds was to determine the effect of eliminating SNP with low MAF on GEq and GC, and to examine the effect of the smaller SNP size.

Estimates of variance components were obtained using the multiple-trait model with the complete data set, which included records across 4 generations. Genetic evaluations were performed by modified BLUP90IOD [START_REF] Tsuruta | Use of the preconditioned conjugate gradient algorithm as a generic solver for mixed-model equations in animal breeding applications[END_REF][START_REF] Misztal | BLUPF90 and related programs (BGF90)[END_REF][START_REF] Vitezica | One-step vs. multistep methods for genomic prediction in presence of selection[END_REF] using BLUP (no genomic information) and ssGBLUP (with genomic information). Predictive ability, r u u mp e (ˆ, ) , + + was defined as the correlation between predicted breeding value and the observed phenotype, composed of the sum of the true breeding value, maternal permanent environmental effect, and residual, using the formula shown by [START_REF] Legarra | Performance of genomic selection in mice[END_REF]. Accuracy, or correlation between predicted and true breeding values, was calculated as

r u u r u u mp e h (ˆ, ) (ˆ, ) / , = + +
where h is the square root of heritability.

RESULTS AND DISCUSSION

Estimates of the additive variance (heritability) were 1.3 (0.23) for BW and 5.0 (0.25) for BM. Estimates of the genetic correlations between BW and BM were 0.4. Table 2 shows the AvgD and AvgOff of G constructed with equal or current allele frequencies and with different cutoffs for MAF. Averages for A 22 are presented for reference. With GEq, as MAF increased to 0.4, AvgD decreased from 1.37 to 1.02 and AvgOff decreased from 0.72 to 0.02. The differences from 1.0 and 0.0, respectively, decreased in GEq as the threshold was increased because the remaining SNP had average allele frequencies closer to 0.5. The differences between AvgD and AvgOff increased from 0.65 to 1.0. This difference, which was smaller than 1.0, was found by [START_REF] Forni | Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information[END_REF] to be closely associated with upward biases in estimates of the additive variance. With GC, AvgD was close to 1.0 and AvgOff was 0, regardless of the threshold. This was expected because the contribution of SNP to off-diagonals of G was proportional to the allelic frequencies, and they roughly canceled out. This was also the case with GEq for MAF >0.4, in which case the assumed frequency of 0.5 roughly corresponded to the true frequency. In addition, the differences between the averages were close to 1.0.

Table 3 shows the AvgD and AvgOff of the inverse of G, with values for the inverse of A 22 presented for reference. In G -1 , all AvgOff were close to 0 and AvgD were between 8.6 and 12.9. For comparison, in A 22 1 -, AvgOff were close to 0 and AvgD were 2.37. With GEq, the greatest of AvgD was at MAF extremes. With GC, AvgD increased with the MAF cutoff.

The diagonal of the inverse of the numerator relationship matrix A 22 has the following form:

a ii = 4/(4 -n p ) + n pr1 /2 + n pr2 /3,
where n p is the number of known parents of bird i, n pr1 is the number of progeny with the mate known, and n pr2 is the number of progeny with the mate unknown. Assuming that differences between AvgD in the inverses of A 22 and G are due to additional information in G relative to A 22 as well as to additional factors, we can write

g ii = a ii + q i /2 + e i ,
where q i is the extra information attributable to the genomic information in units of effective progeny, and e i is noise attributable to both incorrect scaling and the limited number of SNP.

When G is "adequately" scaled and derived from a large number of SNP, quantity q i can be estimated well because e i is likely close to 0. In separate experiments (results not shown), AvgD increased as the number of SNP was reduced. Thus, when the number of SNP used to create G decreased, the absolute value of e i increased. Assuming that GEq was scaled correctly, the increase in AvgD with the increasing threshold was likely due to extra noise in estimating G when the number of SNP was reduced. As the number of SNP was decreased, the numerical accuracy of G was decreased, creating "fake" relationships. Therefore, e i may have been reflecting information for false relationships when G was inaccurate (computed with too few SNP). [START_REF] Yang | Common SNPs explain a large proportion of the heritability for human height[END_REF] looked at underestimation of the additive variance by assuming that causal variants have small allele frequencies but that the SNP information is restricted to SNP with larger MAF. They used unrelated individuals and SNP information from 50 to 500 k and provided corrections for G based on MAF and the number of SNP. Corrections by [START_REF] Yang | Common SNPs explain a large proportion of the heritability for human height[END_REF] did not improve the accuracy in the study by [START_REF] Forni | Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information[END_REF].

Table 4 shows the predictive ability for the validation data with different G. The predictive abilities were very similar across all G. For equal allele frequencies, the predictive abilities reached the peak at the threshold of 0.1 to 0.3 for BW, and increased with the threshold for BM. Increasing the threshold decreased the amount of information in G but brought the average allele frequency close to 0.5. The somewhat greater accuracy for the cutoff of 0.4 in BM could be due to a major gene in the trait that has MAF close to 0.5 and then reduced noise in the estimation with a smaller number of SNP. Greater accuracies at larger MAF are in contrast with the assumptions of [START_REF] Yang | Common SNPs explain a large proportion of the heritability for human height[END_REF]. However, they looked for causative genes in a population of unrelated individuals, whereas G here reflects improved pedigrees of closely related individuals.

With the current allele frequencies, the accuracies were the same for the cutoff <0.4 for BW, with a peak for BM at the threshold of 0.4. With GEq, the accuracies could be greater at low thresholds because the off-diagonals of G were not affected by incorrect computations of contributions from SNP with small allele frequencies.

In the dairy cattle population (Aguilar et al., 2010), the accuracy was greater with GEq than with GC. The reverse was true here. In dairy cattle, the parent average for young bulls was biased even without the genomic information, with the bias attributable to the preferential treatment and specific population structure of the dairy cattle population. It is possible that this bias was partially compensated by bias introduced by GEq. The population structure in chickens is much simpler than that in dairy cattle. [START_REF] Moser | Accuracy of direct genomic values in Holstein bulls and cows using subsets of SNP markers[END_REF] looked at the accuracy of GEBV for several dairy traits by using subsets of SNP information selected by different methods. The subset of 5,000 SNP allowed >90% accuracy to be achieved compared with the complete SNP information. Although the greatest accuracy for each trait was obtained with trait-specific SNP, very few of the highest ranked SNP for 1 trait were highest ranked for other traits. The accuracy with 5,000 SNP selected by MAF alone was greater for cows but was reduced for bulls. In this study, we experienced a marginal decline in accuracy below 22,000 SNP for BW, but not for BM.

Table 5 shows the average EBV for the validation population and for the complete population. If we assume that BLUP provides unbiased EBV, then, on average, all EBV with GEq are biased upward and all Table 3. Means of the diagonal (AvgD) and off-diagonal (AvgOff) for the inverse of the pedigree-based relationship matrix A 22 and the genomic relationship matrix G for genotyped birds (n = 4,113) Predictive ability is defined as correlations between the predicted breeding value and the sum of the true breeding value, maternal effect, and residual.

EBV with GC are biased downward. The biases are approximately 2 additive SD for BW and 1 SD for BM. Additionally, on average, the difference in EBV between genotyped birds and all birds appeared biased upwards with GEq and downwards with GC. Thus, the biases seemed to be dependent on the structure of G. In particular, this could have been due to differences in AvgOff between G and A 22 1 -.

In the single-step method

var u H ( ) = ,
where . .

H A A A G A A A A A G GA A G = + - ( )           - - - - 11
Assume that the genomic relationship matrix is computed with an offset as G* = G + 11′α, where α is a constant and 1 is a vector of ones (that is, adding α to all elements of G). The new matrix H* will be defined as in [1]. This is equivalent to fitting a model where

u u i i i * , var( ) , = + = δ µ µ α
where δ i is 1 for genotyped birds (u g ) and 0 for the remaining birds; thus, In this case, u* contains the mean of all genotyped birds. That mean is automatically applied to the EBV of young genotyped birds. Therefore, if the genotyped birds are selected, their difference from the mean is added to all new genotyped birds, causing bias. Although the value α treated as variance causes regression of the estimate of μ toward zero, that regression is small for all but very small α because of the quantity of information that many genotyped birds contribute to a single variable μ.

To test the hypothesis that a mean of G would affect bias in ssGBLUP, different α were applied to G with the current allele frequencies and all SNP. Only results for positive α are presented because G with α < 0 was not positive definite (and does not make sense as the variance of μ). Average EBV for genotyped birds and all birds, and for their differences, are presented in Table 6. As α was increased to 0.05, the averages for all birds changed by about 1 SD for both BW and BM. Thus, the averages, and thus biases (compared with BLUP), were strongly influenced by AvgOff. All averages were similar to those with BLUP at α = 0.014. With this offset, AvgOff in A 22 and G + 11′α were identical. This suggests that the optimal G should have AvgD and AvgOff close to that of A 22 . Although similar AvgD -AvgOff in G and A 22 ensured unbiased estimates of the additive variances, identical AvgOff seemed to remove biases for the EBV of genotyped birds. The average differences between the EBV of genotyped birds and all birds as a function of a difference between the AvgOff elements of G and A 22 are shown in Figure 1 (using BW). The smallest bias was with the MAF cutoff of 0.4. In that case, both the AvgD and AvgOff of G were the closest to those in A 22 .

Adding α had a very small effect on the predictive abilities. The predictive ability for BW was 0.29 for any α. The predictive ability for BM was 0.36 for α < 0.015, and was decreased to 0.35 for larger α. Thus, the offset in G had a negligible effect on the accuracy of the validation population in this study.

The bias attributable to the offset in G mainly affects the genotyped birds, although it also affects the ungenotyped birds through the H matrix. The size of that bias depends on the strength of selection, the superiority of the genotyped as opposed to the ungenotyped birds, and the amount of information in older generations. In populations under weak or no selection, and when genotyped birds are a random sample from the population, the genotyped birds may have an average merit close to that of ungenotyped birds. Thus, any bias incurred by the offset of G is likely to be small or none. Under strong selection, the offset will bias the EBV of young genotyped birds. Additionally, it may reduce the accuracy of EBV of the older birds and, subsequently, their progenies if the genomic information is a large fraction of the total information. In this study, the offset affected only the bias and not the accuracy of the newest generation because the phenotypic information for older birds was large (about 100 times more records than genotypes). The stronger bias for BM suggests that selection for this trait was stronger than for BW.

Vitezica et al. ( 2010) looked at the biases and accuracy of BLUP, ssGBLUP, and a 2-step method with simulated data under weak and strong selection. The data consisted of 10 generations of animals with 10% of animals genotyped with G constructed using base allele frequencies. The EBV by BLUP was always unbiased. The accuracy by ssGBLUP was greater by the 2-step method under weak selection, but it was less under strong selection. After adjusting the offset in G, the accuracy of ssGBLUP increased, especially under strong selection, and exceeded that of the 2-step method [START_REF] Vitezica | Selection bias in genomic predictions of single-step and multi-step methods[END_REF]. They also provided a formal proof for the calculation of the offset and showed that it was the same as that developed in human studies by [START_REF] Powell | Reconciling the analysis of IBD and IBS in complex trait studies[END_REF]. Therefore, the correct offset of G is critical for small populations, especially under strong selection when the genomic information is a large fraction of the population and where progeny sizes are small. The use of adjusted G was also evaluated for Holsteins for type by using a large number of phenotypes and genotypes [START_REF] Tsuruta | Multipletrait genomic evaluation of linear type traits using genomic and phenotypic data in US Holsteins[END_REF]. The bias removed by adjusting the offset in G was small because of the large progeny sizes for bulls.

The population in this study consisted of a single line of a single breed. When 2 lines of chickens were evaluated simultaneously, even though using different G, the rankings within each line were almost identical but the means of the lines were different [START_REF] Simeone | Evaluation of a multi-line broiler chicken population using a single-step genomic evaluation procedure[END_REF]. This was due to different allele frequencies and different selection pressures in each line. In a joint evaluation of multiple lines or multiple breeds, sections of G corresponding to each line and each cross may require separate adjustments (e.g., as in [START_REF] Harris | Genomic predictions for New Zealand dairy bulls and integration with national genetic evaluation[END_REF].

Removing SNP with MAF ≤0.3 in a genomic evaluation of a commercial chicken population changed the realized accuracy of prediction only marginally. Those accuracies were slightly greater when assuming the current rather than the equal gene frequencies, especially when SNP with low MAF were retained. Removing alleles with low MAF increased the average elements of the inverse of G created with the current allele frequencies. This increase attributable to false relationships introduced by less accurate G caused inflated accuracies computed by inversion. In the genomic evaluation, differences in average EBV of birds with and without the genomic information varied with G and the threshold of MAF. These differences depended on the level of selection and the data structure, and they could be controlled by adding a constant to G constructed by assuming current allele frequencies. When AvgOff elements in A 22 and G were identical, the differences were the same as in BLUP. Assuming that BLUP is unbiased, an unbiased genomic evaluation could be obtained by adding a constant to G derived by assuming current allele frequencies.

Table 1 .

 1 Description of phenotypic records for genotyped birds and all birds

		No. of
	Item 1	records	Mean	SD
	Genotyped
	BW, 100 g	4,113	25.79	3.20
	BM, cm 2	3,923	43.98	5.75
	All
	BW, 100 g	287,614	26.25	4.91
	BM, cm 2	69,057	45.52	7.34

Table 2 .

 2 Means of the diagonal (AvgD) and off-diagonal (AvgOff) for the pedigree-based relationship matrix A 22 and the genomic relationship matrix G for genotyped birds(n = 4,113)

		No. of
	Item 1	SNP	AvgD	AvgOff	Diff 2
	A 22	-	1.004	0.014	0.990
	G
	GEq
	All SNP	57,636	1.373	0.722	0.651
	MAF0.1	41,749	1.197	0.368	0.829
	MAF0.2	32,885	1.121	0.217	0.904
	MAF0.3	22,619	1.062	0.101	0.961
	MAF0.4	11,253	1.020	0.025	0.995
	GC
	All SNP	57,636	1.017	0	1.017
	MAF0.1	41,749	1.015	0	1.015
	MAF0.2	32,885	1.013	0	1.013
	MAF0.3	22,619	1.012	0	1.012
	MAF0.4	11,253	1.008	0	1.008
	1 GEq = equal allele frequencies; MAF = minor allele frequencies;
	GC = current allele frequencies.

Table 4 .

 4 Predictive ability 1 with BLUP and single-step genomic BLUP (ssGBLUP) for BW and breast meat area (BM)

	Item 2	BW	BM
	BLUP	0.19	0.20
	ssGBLUP
	GEq
	All SNP	0.28	0.34
	MAF0.1	0.29	0.35
	MAF0.2	0.29	0.35
	MAF0.3	0.29	0.36
	MAF0.4	0.28	0.37
	GC
	All SNP	0.29	0.36
	MAF0.1	0.29	0.36
	MAF0.2	0.29	0.37
	MAF0.3	0.29	0.37
	MAF0.4	0.28	0.38
	1

Table 5 .

 5 Means of EBV with BLUP and single-step genomic BLUP (ssGBLUP) for genotyped birds (n = 4,113) and all birds (n = 289,994) for BW and breast meat area (BM)

		Genotyped	All
	Item 1	BW	BM	BW	BM
	BLUP	2.77	2.27	1.67	1.38
	ssGBLUP
	GEq
	All SNP	3.94	3.21	2.51	2.15
	MAF0.1	3.89	3.21	2.54	2.16
	MAF0.2	3.82	3.17	2.51	2.13
	MAF0.3	3.61	3.02	2.38	2.03
	MAF0.4	3.00	2.53	1.93	1.67
	GC
	All SNP	0.45	0.32	-0.05	-0.09
	MAF0.1	0.46	0.32	-0.05	-0.10
	MAF0.2	0.46	0.32	-0.05	-0.10
	MAF0.3	0.46	0.32	-0.05	-0.09
	MAF0.4	0.46	0.32	-0.04	-0.08
	1 GEq = equal allele frequencies; MAF = minor allele frequencies;
	GC = current allele frequencies.

Corresponding author: ching-yi_chen@newsham.com

Current address: Newsham Choice Genetics, 701 Crown Industrial Ct., Chesterfield, MO 63005. Received September 29, 2010. Accepted March 21, 2011.

at INRA Institut National de la Recherche Agronomique on July 24, 2013 www.journalofanimalscience.org Downloaded from

BW = BW at 6 wk; BM = ultrasound area of breast meat.

Diff = AvgD -AvgOff.

GEq = equal allele frequencies; MAF = minor allele frequencies; GC = current allele frequencies.

Different genomic matrices at INRA Institut National de la Recherche Agronomique on July 24, 2013 www.journalofanimalscience.org Downloaded from

comments by the 3 anonymous reviewers. We also thank Cobb-Vantress (Siloam Springs, AR) for access to the data for this study. This study was partially funded by the Holstein Association (Brattleboro, VT), Smithfield Premium Genetics (Rose Hill, NC), and Agriculture and Food Research Initiative (AFRI) grants 2009-65205-05665 and 2010-65205-20366 from the USDA National Institute of Food and Agriculture (NIFA) Animal Genome Program (Washington, DC). A. Legarra acknowledges funding from Agence Nationale de la Recherche (ANR, Paris, France) projects AMASGEN and SheepSNPQTL.

LITERATURE CITED


N Bottinelli
email: nicobottinelli86@hotmail.com

Y Capowiez

V Hallaire

J Ranger

P Jouquet

Inability of Near Infrared Reflectance Spectroscopy (NIRS) to identify belowground earthworm casts in no-tillage soil

Keywords:

Several studies have emphasised the ability of Near Infrared Reflectance Spectroscopy (NIRS) to identify surface earthworm casts in the field. However, less is known about casts deposited within the soil, which usually represent the majority found in the field. This study tested the ability of NIRS to identify belowground casts in agricultural systems. Casts and surrounding soils were sampled at depths of 20-30 cm in a loamy soil under no tillage for 12 years. To distinguish different types of cast, sizes and orientations relative to the horizontal plane were measured. NIRS analyses and analyses of carbon and nitrogen content were also performed to compare casts to surrounding soils. Casts were classified into 4 size classes, with no preferential orientation. Cast carbon and nitrogen content were not influenced by their size and did not differ from surrounding soils. PCAs performed on the NIRS data did not allow casts to be differentiated from surrounding soils, regardless of size class. However, soil aggregates were clearly differentiated probably due to their spatial distribution in the soil. Although this study did not identify specific NIRS signatures for casts, it shows the utility of this method to investigate the origin of the soil consumed by earthworms. In our case, NIRS analyses suggest that the high bulk density of the soil (1.42 g cm -3 ) forced ingestion by endogeic earthworms, simply to move around, without preferential selection for organic matter. Consequently, their casts were deposited a few mm from where they had ingested soil with similar organic matter quality.

Introduction

As an organised system, any modification of soil structure results, though not exclusively, from the activity of plant roots and soil fauna. Among the main regulators of soil structure are earthworms, often called 'soil engineers' due to their high abundance and activity in soils and their ability to create galleries and soil organo-mineral aggregates called 'casts' [START_REF] Lee | Soil fauna and soil structure[END_REF][START_REF] Jouquet | Soil invertebrates as ecosystem engineers. Intended and accidental effects on soil and feedback loops[END_REF]. Earthworms have been the focus of many studies since the pioneering work of [START_REF] Darwin | The Formation of Vegetable Mould Through the Action of Worms with Some Observations on Their Habits[END_REF]. Nevertheless, the influence of earthworms on soil functioning remains difficult to quantify, mainly because of the difficulty in assessing the dynamics of earthworm biostructures in the field. Recent studies suggest that the specific composition of organic matter in earthworm casts (e.g. resulting from the selective ingestion of organic matter followed by multiple transformations during transit) analysed by Near Infrared Reflectance Spectroscopy (NIRS) could allow significant progress on this topic [START_REF] Hedde | Specific functional signature in soil macro-invertebrate biostructures[END_REF][START_REF] Vélasquez | This ped is my ped: visual separation and near infrared spectra allow determination of the origins of soil macroaggregates[END_REF][START_REF] Jouquet | Relevance and limitations of biogenic and physicogenic classification: a comparison of approaches for differentiating the origin of soil aggregates[END_REF][START_REF] Zhang | Cast production and NIR spectral signatures of Aporrectodea caliginosa fed soil with different amounts of half-decomposed Populus nigra litter[END_REF][START_REF] Cécillon | Soil macroaggregate dynamics in a mountain spatial climate gradient[END_REF][START_REF] Jouquet | Utilization of near infrared reflectance spectroscopy (NIRS) to quantify the impact of earthworms on soil and carbon erosion in steep slope ecosystem A study case in Northern Vietnam[END_REF][START_REF] Jouquet | What are the limits of the drilosphere? An incubation experiment using Metaphire posthuma[END_REF][START_REF] Zangerlé | Do earthworms and roots cooperate to build soil macroaggregates? A microcosm experiment[END_REF].

NIRS is a rapid and reproducible analytical method that measures the intensity of the absorption of near-infrared light of a sample. The resulting spectra give unique chemical fingerprints, which have been used to differentiate earthworm-impacted soils from surrounding soil that has not been recently ingested, usually considered as reference. The use of NIRS to differentiate earthworm casts from surrounding soil has been used mainly for surface soils in the field [START_REF] Hedde | Specific functional signature in soil macro-invertebrate biostructures[END_REF][START_REF] Vélasquez | This ped is my ped: visual separation and near infrared spectra allow determination of the origins of soil macroaggregates[END_REF][START_REF] Jouquet | Relevance and limitations of biogenic and physicogenic classification: a comparison of approaches for differentiating the origin of soil aggregates[END_REF][START_REF] Jouquet | Utilization of near infrared reflectance spectroscopy (NIRS) to quantify the impact of earthworms on soil and carbon erosion in steep slope ecosystem A study case in Northern Vietnam[END_REF] or on earthworm casts produced in simplified systems in laboratory conditions [START_REF] Zhang | Cast production and NIR spectral signatures of Aporrectodea caliginosa fed soil with different amounts of half-decomposed Populus nigra litter[END_REF][START_REF] Jouquet | What are the limits of the drilosphere? An incubation experiment using Metaphire posthuma[END_REF][START_REF] Zangerlé | Do earthworms and roots cooperate to build soil macroaggregates? A microcosm experiment[END_REF]. Less is known, however, about the relevance of this method for differentiating subsurface earthworm casts in the field [START_REF] Cécillon | Soil macroaggregate dynamics in a mountain spatial climate gradient[END_REF], which can represent a large proportion of casts (especially for endogeic species) [START_REF] Whalen | Quantifying surface and subsurface cast production by earthworms under controlled laboratory conditions[END_REF]; Bottinelli et al., 2010b;[START_REF] Capowiez | Using X-ray tomography to quantify earthworm bioturbation non-destructively in repacked soil cores[END_REF], since their origins and dynamics are by nature more complex than those recreated in laboratory conditions or described on the soil surface.

The objective of this study was therefore to test the relevance of NIRS to differentiate belowground casts in the field. We hypothesised that NIRS signatures could differentiate earthworm casts from surrounding soil. Casts and surrounding soils were sampled at depths of 20-30 cm in soil under no tillage for 12 years. NIRS analyses and analyses of carbon and nitrogen content were performed to compare casts of different types (e.g., from different sizes and orientations) to surrounding soils lacking visible earthworm activity.

Material and methods

The study was performed in November 2011 at the Kerguéhennec Experimental Station located in western France (47 • 52 58 N, 2 • 43 59 W). The climate is temperate oceanic, with a mean annual precipitation of 890 mm. The soil is classified as a Humic Cambisol (FAO/WRB classification) with, at a depth of 20-30 cm, a loamy texture (17% clay, 42% silt, 42% sand), 2.0% of organic carbon content and a bulk density of 1.42 g cm -3 [START_REF] Viaud | Response of organic matter to reduced tillage and animal manure in a temperate loamy soil[END_REF], with fertilisation by poultry manure and no tillage. More details of management practices are described in [START_REF] Viaud | Response of organic matter to reduced tillage and animal manure in a temperate loamy soil[END_REF]. During the study, the field was planted with rapeseed as a cover crop in the crop sequence wheat/catch crop/maize. In May 2007, earthworms had been sampled using formaldehyde [START_REF] Bouché | Stratégies lombriciennes[END_REF] coupled with a hand sorting method [START_REF] Cluzeau | Indicateurs macrobiologiques des sols: les lombriciens-Méthode d'échantillonnage dans les agrosystèmes en zone tempérée[END_REF] and identified using [START_REF] Bouché | Lombriciens de France. Ecologie et systématique[END_REF] key in laboratory. Nine species were found: four endogeic species (Aporrectodea caliginosa, Allolobophora icterica, Allolobophora rosea and Allolobophora chlorotica typica), two anecic species (Lumbricus terrestris and Aporrectodea giardi) and three epigeic species (Dendrobaena mammalis, Lumbricus rubellus castaneus and L. rubellus). In total, earthworm abundance was 489 individuals.m -2 and dominant species were A. caliginosa (44%), L. terrestris (20%) and A. icterica (18%).

Soil sampling and preparation

Because roots of rapeseed were abundant close to the surface which impeded the excavation of intact blocks, sampling was performed in the 20-30 cm soil layer. Twenty-three soil blocks (15 cm long × 10 cm wide) were collected manually from a single soil profile 5 m long. Blocks were oriented, labelled and carefully transported to the laboratory. Thereafter, blocks were air-dried and photographed with a Nikon D700 camera with a resolution of 30 m per pixel. On an image of each block, clusters of casts were meticulously outlined and manually labelled (Fig. 1). In total, 109 clusters were identified for analysis.

Types of casts

To identify different types of casts, the diameter of individual casts and the orientation of cast clusters (i.e., the angle between a cluster's primary axis and a line parallel to the x-axis of the image) were measured with imageJ software. It is likely that different cast sizes and orientation of cast clusters could be related to different species or development stages of earthworms that produced them. Five individual casts were randomly selected per cluster and their size measured. The size distribution was bimodal, with two peaks from 1.4-1.8 mm to 1.8-2.4 mm (Fig. 2), which represented 85% size (mm) Individual casts (%) 1.0 1.4 1.8 2.2 2.6 3.0 0 5 10 15 20 25 of all casts. The orientation of clusters showed no particular distribution (Fig. 3), demonstrating that the orientation of clusters was random. Thus, cast types were classified only by four cast sizes (1-1.4, 1.4-1.8, 1.8-2.4 and 2.4-3 mm) in the subsequent chemical analyses.

Chemical analyses

Chemical analyses were measured on cast clusters and surrounding soils (e.g. soil collected few mm from casts and lacking visible earthworm activity with the naked eye). The carbon and nitrogen contents of casts and surrounding soils were measured in three replicates for the four size classes with a CHN analyser (Elementar Vario EL III). NIRS analyses were performed on casts and adjacent surrounding soil on soil blocks without any pre-treatment ('unground' samples) or extracted from blocks and manually ground to a fine powder (∼0.1-0.5 g < 200 m) ('ground' samples). Because some clusters were too small to provide enough material or because surrounding soil without earthworm features was not clearly identified, only 62 and 55 clusters were analysed in unground and ground samples, respectively. Samples were scanned with a spectrometer (Matrix-F, Bruker) in the 1100-2500 nm spectral range. The reflectance measurements were made at 0.5 nm intervals. Reflectance (R) was converted to absorbance (A) using the equation: A = log(1/R) and further transformed to second derivative according to general procedures recommended to remove baseline shifts and separate overlapping absorption [START_REF] Cécillon | Soil macroaggregate dynamics in a mountain spatial climate gradient[END_REF]. Finally, average values were calculated for 20 nm intervals in order to reduce the number of variables processed.

Statistical analysis

All statistical calculations were carried out using R software. Differences in carbon and nitrogen contents between casts and surrounding soils were tested using multiple comparison Kruskal-Wallis tests with R's 'pgirmess' library. Principal component analyses (PCA) were performed from the NIRS data with R's 'ade4' library. Data corresponded to a matrix of 68 wavelengths × 124 objects for the unground samples and 68 wavelengths × 110 objects for the ground samples. Differences in PCA coordinates between casts and surrounding soils were tested using permutation tests. Differences were evaluated at the 0.05 probability level of significance.

Table 1

Mean (standard deviation) carbon and nitrogen contents of casts in the 1-1.4, 1.4-1.8, 1.8-2.4 and 2.4-3 mm size classes and surrounding soils.

Size (mm)

Type C (%) N (%)

1.0-1.4 Cast 1.9 (0.5) 0.2 (0.1) Soil 1.5 (0.4) 0.2 (0.0)

Results

The carbon and nitrogen contents of casts and surrounding soils were not statistically different regardless of cast size (Table 1), and values were similar among the 4 sizes considered (P > 0.05).

From the PCA of NIRS spectra, only 2 and 4 groups were observed in unground and ground samples, respectively (Fig. 4). For unground samples, separation mainly occurred along the first axis of the PCA, which explained 13% of the variability. For ground samples, aggregates were separated on the first and second axes, which explained 16 and 12%, respectively. The groups separated, however, referred neither to aggregate origin or size, but seemed to result from natural soil heterogeneity. In the PCA, NIRS spectra of casts were always similar to those of their surrounding soils (P > 0.05 for all sizes), whether unground or ground. (1) 1-1.4 mm, (2) 1.4-1.8 mm, (3) 1.8-2.4 mm and (4) 2.4-3 mm. Lines connect each cast cluster with its surrounding aggregates. Tags with group names represent the barycentre of each group.

C1 S1 C2 S2 C3 S3 C4 S4 C1 C2 C3 C4 S1 S2 S3 S4 Eigenvalue Eigenvalues (a) Unground samples Axis 1 = 13% Axis 2 = 8% -3 0 1 2 3 -2 -1 C1 S1 C2 S2 C3 S3 C4 S4 C1 C2 C3 C4 S1 S2 S3 S4 Eigenvalues Eigenvalues Axis 1 = 16% Axis 2 = 12% (b) Ground samples 0 1 2 3 4 -4 -3 -2 -1 -3 0 1 2 3 -2 -1

Eigenvalues

Discussion

Origin of belowground earthworm casts

Casts found within the soil are thought to be associated mainly with the casting activity of endogeic species. However, laboratory studies have shown that anecic species can also deposit casts within the soil without lining their burrows [START_REF] Curry | Cast production and soil turnover by earthworms in soil cores from South Australian pastures[END_REF][START_REF] Felten | Earthworm burrowing behaviour in 2D terraria with single-and multi-species assemblages[END_REF][START_REF] Capowiez | Using X-ray tomography to quantify earthworm bioturbation non-destructively in repacked soil cores[END_REF]. To differentiate anecic from endogeic origins of casts, several criteria could be considered. First, burrows created by anecic species are mainly vertically to sub-vertically oriented, but those of endogeic species have no preferential orientation [START_REF] Felten | Earthworm burrowing behaviour in 2D terraria with single-and multi-species assemblages[END_REF][START_REF] Capowiez | Using X-ray tomography to quantify earthworm bioturbation non-destructively in repacked soil cores[END_REF]; therefore, clusters of casts may have a preferential orientation according to the type of earthworm. Second, adult anecic species are larger than adult endogeic species [START_REF] Bouché | Lombriciens de France. Ecologie et systématique[END_REF][START_REF] Pérès | Earthworm activity and soil structure changes due to organic enrichments in vineyard systems[END_REF], meaning that their casts may be larger, too. Finally, anecic species select their food, which is richer in organic matter than that ingested by endogeic species [START_REF] Curry | The feeding ecology of earthworms -a review[END_REF]. As a result, casts produced by anecic species generally have more organic matter than surrounding soils [START_REF] Hedde | Specific functional signature in soil macro-invertebrate biostructures[END_REF][START_REF] Jouquet | Chemical and physical properties of earthworm casts as compared to bulk soil under a range of different land-use systems in Vietnam[END_REF], which is not the case for endogeic species, depending on soil properties [START_REF] Barois | Ecology of earthworm species with large environmental tolerance and/or extended distributions[END_REF]Bottinelli et al., 2010a).

In the present study, the orientation of cast clusters seemed to be random. Cast size ranged from 1 to 3 mm, with most ranging from 1.8 to 2.4 mm. Finally, the carbon and nitrogen contents were similar between casts and their surrounding soils, suggesting that earthworms did not select their food. All these results strongly suggest that casts found in the 20-30 cm soil layer were all produced by endogeic species, presumably A. caliginosa and/or A. icterica, which had previously been shown to represent 64 and 27% of the total abundance of endogeic earthworms, respectively. Because A. caliginosa and A. icterica are nearly a similar size at the study site, their casts probably have the same size. Hence, the 4 size classes of casts more likely reflect the development stage of these earthworms instead of the species.

NIRS as a tool to differentiate endogeic earthworm casts in crop fields

Grinding soil aggregates increased the number of groups differentiated from 2 to 4, confirming the utility of extracting and grinding soil aggregates before measuring their NIRS signatures, as suggested by [START_REF] Brunet | Determination of carbon and nitrogen contents in Alfisols, Oxisols and Ultisols from Africa and Brazil using NIRS analysis: effects of sample grinding and set heterogeneity[END_REF]. It is likely that reflectance is influenced by the size of aggregates and that this effect is reduced by grinding. Differentiation of groups seemed to result from soil heterogeneity. NIRS fingerprints are known to be plot-dependent at ecosystem scales [START_REF] Vélasquez | This ped is my ped: visual separation and near infrared spectra allow determination of the origins of soil macroaggregates[END_REF][START_REF] Cécillon | Soil macroaggregate dynamics in a mountain spatial climate gradient[END_REF], but at smaller scales, our study confirms the high spatial heterogeneity within soil [START_REF] Tisdall | Organic matter and water-stable aggregates in soils[END_REF][START_REF] Young | Can there be a contemporary ecological dimension to soil biology without a habitat? Discussion[END_REF][START_REF] Young | Microbial distribution in soils: physics and scaling[END_REF] and highlights the usefulness of NIRS to describe it [START_REF] Patzold | Soil heterogeneity at the field scale: a challenge for precision crop protection[END_REF][START_REF] Wielopolski | Soil analysis based on samples withdrawn from different volumes: correlation versus calibration[END_REF].

Regardless of treatment, NIRS did not distinguish casts and their surrounding soils in this study, suggesting that NIRS is not always useful for differentiating recently emitted casts from their surrounding soils in the field. The similar NIRS signatures and carbon and nitrogen contents of casts and surrounding soils highlight that their quantity and quality of soil organic matter was similar. These results may be explained by (i) the high soil bulk density (1.42 g cm -3 ), which may have forced earthworms to ingest soil simply to move around [START_REF] Whalley | Biological effects of soil compaction[END_REF] and did not allow them to select their food. As a result, casts were deposited a few mm from where they were ingested and had similar soil organic matter content and quality as surrounding soils. In addition, (ii) the surrounding soil may not have represented a true control, instead consisting mainly of weathered casts not discernible with the naked eye, since this soil is highly bioturbated (Bottinelli et al., in press). Finally, (iii) earthworms may have influenced the NIRS signature of the surrounding soil by leaching nutrients, which may have extended the limit of the drilosphere around the casts [START_REF] Jouquet | What are the limits of the drilosphere? An incubation experiment using Metaphire posthuma[END_REF]. Hence, the absence of an influence of casting activity on NIRS fingerprints was not completely surprising in this context. These results, however, contradict previous studies performed in laboratory conditions [START_REF] Zhang | Cast production and NIR spectral signatures of Aporrectodea caliginosa fed soil with different amounts of half-decomposed Populus nigra litter[END_REF][START_REF] Zangerlé | Do earthworms and roots cooperate to build soil macroaggregates? A microcosm experiment[END_REF] that found specific NIRS signatures for casts. Differences between these experiments and ours may be explained by the fact that the laboratory casts (i) were enriched in organic matter since endogeic earthworms were fed with leaf litter and (ii) were compared to a non-ingested soil, sampled from microcosms without earthworms [START_REF] Jouquet | What are the limits of the drilosphere? An incubation experiment using Metaphire posthuma[END_REF]. It is evident that these laboratory conditions did not reflect the complex conditions of our field study. Surprisingly, our results did not confirm those of [START_REF] Cécillon | Soil macroaggregate dynamics in a mountain spatial climate gradient[END_REF], who identified specific NIRS signatures of belowground casts when compared to physicogenic aggregates sampled in the Ah horizons of three mountain forest ecosystems. These different results may be due to (i) higher soil organic matter content and macroporosity of the soils sampled and (ii) different earthworm species producing belowground casts.

NIRS was unable to identify casts produced by juveniles or adults of A. caliginosa or A. icterica. In addition, the carbon and nitrogen contents of casts were not statistically different. Hence, our study does not confirm the ability of NIRS to differentiate casts produced by different endogeic species, as shown by [START_REF] Zangerlé | Do earthworms and roots cooperate to build soil macroaggregates? A microcosm experiment[END_REF] in laboratory conditions.

Conclusion

In conclusion, it appears that the usefulness of using NIRS to differentiate cast aggregation may be context-dependent. Nonetheless, though NIRS could not differentiate earthworm casts from their surrounding soil environment, it pinpointed strategies that earthworms, probably the endogeic A. caliginosa and/or A. icterica may use to move, feed and produce casts in the field.

Fig. 1 .

 1 Fig. 1. Example of a vertically oriented soil block sampled for analyses. Clusters of casts are outlined in white.

Fig. 2 .Fig. 3 .

 23 Fig.2. Size distribution of casts sampled in no-tillage with poultry manure fertilisation between 20 and 30 cm depth.

EigenvaluesFig. 4 .

 4 Fig. 4. PCA of NIRS spectra of casts (C) and surrounding soils (S) of (a) unground and (b) ground samples. Numbers indicate the four size classes of casts considered:(1) 1-1.4 mm, (2) 1.4-1.8 mm, (3) 1.8-2.4 mm and (4) 2.4-3 mm. Lines connect each cast cluster with its surrounding aggregates. Tags with group names represent the barycentre of each group.

Acknowledgements

We thank the personnel of the Brittany Chamber of Agriculture, in particular Djilali Heddadj and Patrice Cotinet, who helped perform this work on the Kerguéhennec research site. We are grateful to Patrick Lavelle for constructive discussion, Lauric Cécillon for his help in data analyses and Michael and Michelle Corson for useful English corrections on the manuscript.