Baptiste Monsion

Re ´my Froissart

Yannis Michalakis

Ste ´phane Blanc

Large Bottleneck Size in Cauliflower Mosaic Virus Populations during Host Plant Colonization

The effective size of populations (Ne) determines whether selection or genetic drift is the predominant force shaping their genetic structure and evolution. Despite their high mutation rate and rapid evolution, this parameter is poorly documented experimentally in viruses, particularly plant viruses. All available studies, however, have demonstrated the existence of huge within-host demographic fluctuations, drastically reducing Ne upon systemic invasion of different organs and tissues. Notably, extreme bottlenecks have been detected at the stage of systemic leaf colonization in all plant viral species investigated so far, sustaining the general idea that some unknown obstacle(s) imposes a barrier on the development of all plant viruses. This idea has important implications, as it appoints genetic drift as a constant major force in plant virus evolution. By co-inoculating several genetic variants of Cauliflower mosaic virus into a large number of replicate host plants, and by monitoring their relative frequency within the viral population over the course of the host systemic infection, only minute stochastic variations were detected. This allowed the estimation of the CaMV Ne during colonization of successive leaves at several hundreds of viral genomes, a value about 100-fold higher than that reported for any other plant virus investigated so far, and indicated the very limited role played by genetic drift during plant systemic infection by this virus. These results suggest that the barriers that generate bottlenecks in some plant virus species might well not exist, or can be surmounted by other viruses, implying that severe bottlenecks during host colonization do not necessarily apply to all plant-infecting viruses.

Introduction

The main forces driving changes of the frequencies of alleles within populations are selection (deterministic) and genetic drift (stochastic). While both forces obviously act concomitantly in most cases, the balance of their respective action is governed by a number of factors, one of the most important being the effective population size (Ne) [START_REF] Nei | The bottleneck effect and genetic variability in populations[END_REF][START_REF] Gherman | Population bottlenecks as a potential major shaping force of human genome architecture[END_REF][START_REF] Kimura | The mutation load in small populations[END_REF]. This parameter, commonly defined as the size of an ideal population which would drift at the same rate as the observed population [START_REF] Wright | Evolution in Mendelian populations[END_REF], is widely investigated in the literature both theoretically [START_REF] Criscione | Effective sizes of macroparasite populations: a conceptual model[END_REF][START_REF] Nei | Genetic drift and estimation of effective population size[END_REF][START_REF] Tajima | Note on genetic drift and estimation of effective population size[END_REF][START_REF] Crow | The effective number of a population with overlapping generations: a correction and further discussion[END_REF] and experimentally [START_REF] Araki | A potential bias in the temporal method for estimating Ne in admixed populations under natural selection[END_REF][START_REF] Wang | Estimation of effective population sizes from data on genetic markers[END_REF][START_REF] Johnson | Temporal changes in allele frequencies and low effective population size in greater prairie-chickens[END_REF] in a broad variety of organisms. Ne significantly informs whether the evolution of a given population might be better described by a deterministic or a stochastic evolution model [START_REF] Rouzine | Transition between stochastic evolution and deterministic evolution in the presence of selection: general theory and application to virology[END_REF][START_REF] Jain | Deterministic and stochastic regimes of asexual evolution on rugged fitness landscapes[END_REF][START_REF] Wahl | Models of experimental evolution: the role of genetic chance and selective necessity[END_REF]. In brief, when Ne is large, competition between genetic variants is fully acting, with no or little interference of random processes, and selection shapes the genetic composition of corresponding populations. Conversely, in populations with a small Ne, various processes resulting in stochastic ''sampling'' of genetic variants that will engender the next generation are prevalent and counter the effect of selection.

Experimental evolution using viruses as biological models has developed enormously during the past two decades. A series of converging reports on viruses of bacteria, animals, and plants confirmed that protocols maintaining constant large population sizes lead to selection of the fittest variants, thus augmenting the mean fitness of viral populations [START_REF] Weinreich | Rapid evolutionary escape by large populations from local fitness peaks is likely in nature[END_REF][START_REF] Willi | Genetic rescue in interconnected populations of small and large size of the self-incompatible Ranunculus reptans[END_REF][START_REF] Clarke | Genetic bottlenecks and population passages cause profound fitness differences in RNA viruses[END_REF]. On the other hand, the repeated experimental imposition of severe bottlenecks (temporarily reducing Ne) is associated with relaxed selection and strong drift, dramatically reducing the mean fitness in the population [START_REF] Chao | Fitness of RNA virus decreased by Muller's ratchet[END_REF][START_REF] De La Iglesia | Fitness declines in tobacco etch virus upon serial bottleneck transfers[END_REF][START_REF] Novella | Exponential fitness gains of RNA virus populations are limited by bottleneck effects[END_REF][START_REF] Yuste | Drastic fitness loss in human immunodeficiency virus type 1 upon serial bottleneck events[END_REF]. Viruses are intuitively perceived as having extremely large population sizes during host invasion. While this is true for the census population size, i.e. the number of viral genomes accumulating in single hosts, it does not indicate what fraction thereof is actually actively replicating and propagating, or is efficiently progressing into new organs and tissues, yielding subsequent generations upon host colonization. Estimating viral Ne during systemic infection of a host has proved experimentally difficult, as illustrated by the impressive number of studies published, and still debated, on HIV [START_REF] Kouyos | Stochastic or deterministic: what is the effective population size of HIV-1?[END_REF]. Consequently, virus species for which Ne has been formally evaluated during invasion of various organs and tissues of the host remain extremely few in both plants and animals. The question of whether deterministic or stochastic models better explain the evolution of viral populations is thus pending further investigation of Ne.

In plants, several specific barriers could impose severe bottlenecks on populations of invading viruses. The level at which the existence of such bottlenecks has been most documented is the colonization of new leaves, not only because of the obvious practical ease of comparing the genetic content of different leaves from the same host, but also because this level results from longdistance movement of the virus population, which is loaded into the vascular system from source leaves and unloaded into sink leaves (for review, see [START_REF] Waigmann | Viral transport in plants[END_REF]). As discussed later, other phenomena could also generate population bottlenecks, but virus movement within the host plant through plasmodesmata is most often considered as a major putative obstacle to exponential expansion.

While not directly estimating Ne, two experimental evolution studies demonstrated the existence of viral population bottlenecks during leaf colonization. With a mixture of 12 engineered genetic variants of Cucumber mosaic virus (CMV) co-inoculated into tobacco plants, Li and Roossinck [START_REF] Li | Genetic bottlenecks reduce population variation in an experimental RNA virus population[END_REF] showed that diversity decreased randomly but continuously in the viral population, as increasing numbers of variants were lost when systemic infection progressed into upper, newly-formed leaves. Even more illustrative was the spatial distribution of the genetic diversity of a Plum pox virus (PPV) population, maintained in a perennial host tree for over 13 years [START_REF] Jridi | Distinct viral populations differentiate and evolve independently in a single perennial host plant[END_REF]. While a large number of related PPV variants could be distinguished and detected in various organs and tissues, distinct subpopulations were shown to be isolated in different branches. Extreme population bottlenecks were further evidenced when the virus progressed into newly formed leaves, which were all colonized by one single viral genotype.

The effective size of virus populations during systemic infection of host plants has been evaluated more formally in two instances. Populations of Wheat streak mosaic virus (WSMV) invading a wheat tiller [START_REF] French | Evolution of Wheat streak mosaic virus: dynamics of population growth within plants may explain limited variation[END_REF], and of Tobacco mosaic virus (TMV) invading a tobacco leaf [START_REF] Sacristan | Estimation of population bottlenecks during movement of tobacco mosaic virus in tobacco plants[END_REF], were shown to be founded by as few as 4 and 2-20 virus particles, respectively, hence again resulting in severe bottlenecks.

All the studies cited above converged to the conclusion that the population size of plant viruses fluctuates dramatically and can be temporarily remarkably small, i.e. of the order of one to a few genome units founding the population that subsequently develops to billions of genomes in each systemically infected leaf. Such a demographic regime suggests that genetic drift is a major force in plant virus evolution, as proposed and discussed by several authors [START_REF] Li | Genetic bottlenecks reduce population variation in an experimental RNA virus population[END_REF][START_REF] Jridi | Distinct viral populations differentiate and evolve independently in a single perennial host plant[END_REF][START_REF] French | Evolution of Wheat streak mosaic virus: dynamics of population growth within plants may explain limited variation[END_REF][START_REF] Sacristan | Estimation of population bottlenecks during movement of tobacco mosaic virus in tobacco plants[END_REF][START_REF] Garcia-Arenal | Variation and evolution of plant virus populations[END_REF]. The fact that extreme population bottlenecks are consistently described for unrelated virus species infecting monoor di-cotyledonous hosts appealingly suggests that all plant viruses might be subject to the same phenomenon [START_REF] Li | Genetic bottlenecks reduce population variation in an experimental RNA virus population[END_REF][START_REF] Sacristan | Estimation of population bottlenecks during movement of tobacco mosaic virus in tobacco plants[END_REF][START_REF] Garcia-Arenal | Variation and evolution of plant virus populations[END_REF][START_REF] French | Population structure within lineages of Wheat streak mosaic virus derived from a common founding event exhibits stochastic variation inconsistent with the deterministic quasi-species model[END_REF], perhaps related to the unavoidable physical barriers that hamper the systemic movement of all viruses in plants. However, this tempting generalisation definitely requires closer inspection as, if proven true, it would illustrate the impossibility of viral adjustments and trade-offs on important traits requiring large Ne. In contrast, a single counter example would demonstrate that, in some virusplant associations, such barriers might not exist or might be surmounted by the virus. Further research on more diverse virus species, with different replication strategies and different life cycles, is still needed in order to assume extremely small Ne as a general rule during host plant colonization.

In this report, we assessed the importance of bottlenecks during systemic host colonization in populations of Cauliflower mosaic virus (CaMV), a DNA virus whose biological properties differ largely from those of the RNA viruses investigated in the studies cited above. Monitoring the frequency of several engineered allelic variants within leaves of single host plants revealed remarkably small stochastic fluctuations in the genetic structure of numerous CaMV populations, even over considerable periods of time. Beyond demonstrating that CaMV populations are not subject to intense genetic drift, and hence do not undergo severe demographic bottlenecks, the slight stochastic fluctuations detected were exploited to infer the effective size of CaMV populations upon systemic leaf colonization. Our estimates consistently indicate that several hundreds of genome units founded the viral populations in all leaves analyzed (originating from different plants). This value, which is ,100-fold higher than values previously reported for other plant viruses, demonstrates that extremely small viral Ne during host plant colonization is not a general rule and opens up the possibility of trade-offs on viral traits that directly or indirectly depend on Ne.

Materials and Methods

Engineered CaMV variants

The six plasmids (pCa-VIT1 to pCa-VIT6) used in this study to generate the six CaMV allelic variants have been described in detail and characterized previously [START_REF] Monsion | Quantitative Single-letter Sequencing: a method for simultaneously monitoring numerous known allelic variants in single DNA samples[END_REF].

All are infectious full-length clones of the CaMV Cabb-S isolate [START_REF] Franck | Nucleotide sequence of cauliflower mosaic virus DNA[END_REF], where a specific genetic marker (a dsDNA non-coding sequence of 40 bp) has been inserted between CaMV ORFs II and III. These markers should not affect any viral function, as noncoding sequences between CaMV ORF do not affect translation. When inoculated individually into turnip plants, each CaMV clone (CaMV-VIT1 to CaMV-VIT6) induces symptoms similar to CaMV wild type, and all six genetic markers have been shown to be stably maintained in the viral genome even after three successive passages in plants [START_REF] Monsion | Quantitative Single-letter Sequencing: a method for simultaneously monitoring numerous known allelic variants in single DNA samples[END_REF].

Plant growth conditions and inoculation with mixed CaMV-VIT populations

All turnip plants (Brassica rapa var. ''Just Right'') were maintained in an insect-proof greenhouse under controlled conditions: 25/19uC day/night with a photoperiod of 16/8 h day/night.

Plants infected for 21 days with one of the CaMV-VIT1 to -VIT6 clones (one plant per clone) were used to prepare virus particle-enriched fractions as previously described [START_REF] Monsion | Quantitative Single-letter Sequencing: a method for simultaneously monitoring numerous known allelic variants in single DNA samples[END_REF]. Equal volumes of each of the six virus particle-enriched fractions were pooled to produce a mixed inoculum (designated Mix6VIT) containing all six allelic variants of CaMV. Fifty healthy young plantlets, at the three-leaf stage, were mechanically infected by rubbing 20 ml of Mix6VIT solution onto the entire surface of the two larger leaves, previously powdered with abrasive carborundum. First symptoms indicative of CaMV infection appeared on non-inoculated leaves within 7 to 9 days, and all plants proved systemically infected at 11 days post-infection (dpi).

Author Summary

Whether selection or stochastic genetic drift is the major force driving the evolution of a virus depends largely on the size of the viral population, with the former being predominant in large populations and the latter taking over when population sizes are transiently or durably reduced. This question has been intensively debated in both plant and animal viruses, as demographic fluctuations throughout viral life cycles are poorly understood. In plant viruses, an extremely small population size-down to a few founder genome units colonizing each leaf-has been formally estimated in two instances, and all other virus species investigated so far have consistently been shown to undergo extreme demographic bottlenecks during systemic invasion of their host. This situation conveys the general idea that all viruses are confronted with ''universal barriers'' in plants, imposing repeated transient decreases in their population size, thus making genetic drift a major constant driver of their evolution. Here, using the example of Cauliflower mosaic virus, we mitigate this general idea by showing that at least one virus species can overcome such putative limiting barriers and massively invade leaves with hundreds to thousands of founding genome units.

Harvesting of leaf samples

At 13 dpi, i.e. very soon after systemic infection had developed, the inoculated leaves were discarded and all large expanded leaves were collected. In all cases, two to three young expanding or emerging leaves in the centre of the rosette were left unsevered in order to allow continuous growth of the plants for further sampling as described below. The viral population extracted from each plant at this early stage was considered as the initial population.

Thirty-two days later (at 45 dpi) the same plants had grown continuously and produced 10 to 15 newly expanded leaves. At this stage, one single leaf from each plant, randomly chosen between the 5 th and the 12 th newly formed leaves, was collected. The viral population extracted from each single leaf at this late stage was considered as the final population. All leaf samples were stored at 220uC until further viral DNA purification.

The rationale of this sampling protocol is further discussed below.

Purification of viral DNA and genetic composition analysis

Purification of viral DNA from the harvested leaf samples, as well as the analysis of the genetic structure of the corresponding viral genome populations were performed exactly as described previously, using the QSS (Quantitative Single-letter Sequencing) method [START_REF] Monsion | Quantitative Single-letter Sequencing: a method for simultaneously monitoring numerous known allelic variants in single DNA samples[END_REF]. Briefly, QSS allows the simultaneous quantification of numerous allelic variants in a single DNA sample. After a PCR amplification step using a pair of primers flanking the markercontaining region, the PCR product is submitted to single-letter sequencing primed with a fluorescently-labelled oligonucleotide, located upstream of the markers' position. The resulting monochromatic electropherogram exhibits numerous specific diagnostic peaks, attributable to specific variants, signifying their presence/absence in the DNA sample. Finally, peak fluorescence can be quantified and used to estimate the frequency of the corresponding variant in the DNA population. The accuracy and reproducibility of the QSS method have been fully evaluated and shown to be equivalent or higher than that of competing technologies (including those based on real-time PCR) for quantifying variants with a relative frequency above 5% in the DNA population [START_REF] Monsion | Quantitative Single-letter Sequencing: a method for simultaneously monitoring numerous known allelic variants in single DNA samples[END_REF].

Statistical analysis

Changes in the relative frequencies of the six alleles (-VIT1 to -VIT6) between the initial and final populations, sampled as described above, were precisely monitored in 50 replicate test plants.

In order to calculate the number of founder genomes initiating each final population in a single leaf (N), all parameters accounting for the transition from p (initial frequency of a given marker) to p9 (final frequency of the same marker) must be evaluated. For any given plant the difference between p and p9 (Dp) may be potentially attributed to genetic drift and to selection:

Dp~Dp DriftzDp Selection ð1Þ
There is no a priori reason to assume that any selection affecting the frequency of the markers is heterogeneous across plants (e.g. favouring a marker in a plant and selecting against it in another plant). Moreover, given the way our markers were constructed, there is no a priori reason to expect directional selection in favour of any of them. For the moment, we will thus assume that selection is negligible, and the markers effectively neutral, and will provide further arguments supporting this assumption at the end of this section.

We used two methods to estimate N, both based in the change of genetic variance between the two sampling events.

The first method directly tracks changes in variance. The variance of the estimates of p and p9 between plants can be written as:

Var p 0 ð Þ~Var p ð ÞzVar Drift ð Þ ð2Þ
The variance due to drift is equal to:

Var Drift ð Þ~p 1{p ð Þ=N ð3Þ
Where p is the frequency of a given marker in initial populations.

From Equations 2 and 3 we obtain:

N~p 1{p ð Þ Var p 0 ð Þ{Var p ð Þ ð4Þ
The second technique is based on Fst statistics. Fst statistics were introduced by Wright [START_REF] Wright | Evolution in Mendelian populations[END_REF] and represent a way to partition genetic variance within and between populations. One way to express Fst is:

Fst~H T {H S ð Þ =H T ð5Þ
where H T is gene diversity assuming all populations form a single large population and H S is the average gene diversity within each population. In our case, each plant represents a population. Gene diversities express the probability to randomly draw two different alleles and are thus equal to 1{ P i p 2 i , where p i is the frequency of each allele at the subdivision level under consideration. In our case H T represents the gene diversity obtained after calculating the average frequency of each marker across all plants, while H S represents the average across all plants of within plants gene diversities.

Using standard population genetics theory (e.g. [START_REF] Whitlock | Temporal fluctuations in demographic parameters and the genetic variance among populations[END_REF]) it can be shown that, for a haploid such as a virus,

F 0 ST ~1 N z 1{ 1 N F ST ð6Þ
where the prime denotes sampling at two different points in time.

Rearranging Equation 6 yields an expression that can be used to estimate N:

N~1 {F ST F 0 ST {F ST ð7Þ
We will refer to the estimates obtained by Equations 4 and 7 as N V and N F respectively.

To provide a confidence interval for these estimates we used a resampling technique. We bootstrapped over plants: for each bootstrap we randomly drew a sample of 50 plants with replacement, calculated Var(p9), Var(p) and p, and from that estimated N V . We also calculated H T and H S , and from that estimated N F . We repeated this procedure 1,000 times and constructed a distribution of N V and N F . The 95% confidence intervals correspond to the 2.5 and 97.5 percentiles of these distributions. In several bootstraps, the final variance was smaller than the initial variance, or the final F ST was smaller than the initial F ST , yielding negative values of N V or N F . In such cases, our estimation method does not work. Such cases would correspond to conditions leading to convergence of marker frequencies across plants, i.e. conditions where drift is negligible, and thus N can be regarded as infinite. We conservatively discarded such cases from the calculation of the upper limit of the confidence interval.

The average changes in frequency of the markers proved very small (see Table 1). Neutrality tests applied to VIT1, 3 and 4, demonstrated that only VIT1 and 3 might be slightly selected for and against, respectively (data not shown). Nevertheless, both these variants were used for estimating N because, while directional selection would lead to a situation where marker frequencies would be more similar across plants at the second sampling event than at the first, yielding smaller final variance and F ST , our data show the opposite, strongly contradicting the directional selection hypothesis. Another a posteriori justification is that while the neutrality tests indicate that one marker is slightly selected for, another is slightly selected against and another is effectively neutral, the N V estimates obtained form all three are of the same order of magnitude. These arguments suggest that selection, if present, can be neglected relative to drift in our experiment.

Results

Genetic structure of initial CaMV populations in whole plants

Following systemic spread of infection in plants inoculated with a mixture of 6 CaMV genetic variants (Mix6VIT), populations of viral genomes in each plant were purified from all fully expanded and systemically infected leaves, excluding inoculated leaves. DNA samples from each of 50 replicate plants were considered as the initial populations and were submitted to QSS analysis to detect markers in CaMV genomes and quantify their relative frequency as described in Materials and Methods. The mean frequency of each marker and the variance among the 50 CaMV populations collected at this stage is given in Table 1 (full dataset with details for each plant is provided in Table S1).

The six CaMV-VIT genotypes were found in widely unequal ratios in all 50 cases, with mean relative frequencies being close to 0 for CaMV-VIT5 and -VIT6, and up to around 40% for CaMV-VIT1 and -VIT3. However, a striking observation was that each marker was found with a very similar frequency in all 50 replicate CaMV populations (Table S1), as demonstrated by the very low standard deviation among repeats (Table 1). This result clearly indicates that little stochastic variation is generated by the inoculation process, or between the inoculation process and the sampling of initial populations. The frequency differences between markers being very consistent in all 50 replicates, they are certainly due to a non-stochastic phenomenon (discussed below).

Genetic structure of final CaMV populations in single leaves

At late infection stages, 32 days after collecting the initial populations, a single systemically infected leaf per individual plant, randomly chosen between the 5 th and the 12 th new expanded leaves, was harvested. The CaMV populations extracted from these sampled leaves are considered as ''final populations,'' and their genetic structure (Table 1 andS1) was evaluated in exactly the same manner as that of the initial populations.

Remarkably, in all 50 repeats, the final population resembled the corresponding ancestor population despite the different sampling process (individual leaves for final populations versus pools of leaves for initial populations), demonstrating that very little change had occurred over this considerable time period, whatever the position or age of the analyzed leaf (see E(Dp) values [Table 1], and compare p and p9 values [Table S1] for each of the 50 plants). Consequently, as between initial populations, there were only small variations between the final populations, although the standard deviation was slightly higher in the latter case [Table 1].

Taken together, these results argue for the absence of large stochastic variations in the genetic structure of a CaMV population upon progression of the infection into newly formed leaves. This suggests that the effective size (Ne) of CaMV populations during host systemic colonization is likely to be surprisingly large (compared to data previously published on other plant viruses) as specifically evaluated in the next section.

Estimation of bottleneck size in CaMV populations during systemic leaf colonization

Because the average relative frequency of markers CaMV-VIT2, -VIT5 and -VIT6 was close to, or even below, the limit of accuracy of the QSS method [START_REF] Monsion | Quantitative Single-letter Sequencing: a method for simultaneously monitoring numerous known allelic variants in single DNA samples[END_REF], their quantification was considered poorly reliable, and hence these markers were not used for further analysis. We thus used only the three markers -VIT1, -VIT3 and -VIT4 to estimate the average number of genomes founding the population in each leaf (N), as described in the Materials and Methods.

The estimates of N from these three markers yielded remarkably high values, corresponding to several hundreds of viral genomes. Two different statistical methods were applied to the data set and provided very consistent results, both for the observed bottleneck size and for the limit of the 95% confidence interval (compare N V in Table 2, and N F in Table 3). The lower limit of the confidence interval was .100 in all cases, whereas the higher limit reached thousands.

Discussion

This report evaluates the effective size of CaMV populations during systemic invasion of plant leaves. Several previous studies have used distinct experimental protocols to tackle similar The genetic structure of the 50 sampled CaMV populations is determined by the different markers present in each, as well as by their respective relative frequency. For reasons of clarity, only mean frequency values and standard deviations among the 50 repeats are shown here. The full dataset with details of the 50 CaMV populations is available (Table S1).

a Sampling protocol for initial and final populations is described in Materials and Methods.

b Mean relative frequency of each marker over the 50 replicates. questions with various viral species; the specifics and rationale of the protocol used here are discussed below. To eliminate variations in the genetic structure of viral populations that could be related to the inoculation process, the virus population present in the whole plant (excluding the inoculated leaves) soon after systemic spread of infection was considered as the starting point of the experiment, putative subsequent changes thus occurring only under the influence of within-plant processes. At this initial stage, the virus population uploaded into, and circulating within, the vascular system is most likely best represented by the overall content of the systemically infected leaves, as they have either received viruses from the vasculature, exported viruses into it, or both. Expanded infected leaves were therefore harvested for initial analysis, carefully preserving two to three young newly expanding leaves on the still-growing plant. The infected plants were then left to grow for a period of 32 days, during which the virus population successively colonized 10 to 15 emerging and expanding new leaves. Regardless of where the virus population originates from during this process (vascular system and roots at the beginning, then increasing numbers of leaves later on), analysis of single leaf (between positions 5 and 12 above the initial harvest) contents collected 32 days later should reveal the existence of any putative bottlenecks at any stage of the systemic infection. Indeed, final CaMV populations result largely from the successive leaf-to-leaf passages that occur sequentially when young sink leaves are infected, become sources and subsequently export virus into new sinks.

All genetic variants (CaMV-VIT1 to CaMV-VIT6) were coinoculated at similar locations and at the same time point onto the two first true leaves of turnip plantlets. The reason why the relative proportions of the six variants are highly unequal later in systemically infected plants was mostly that unequal proportions were already present in the initial inoculum, Mix6VIT (further discussed in Table S2). Previous studies, based on co-infection of turnip plants by two distinct CaMV variants with seemingly equal growth rate, have compared the variant ratios in the initial inoculum and in resulting systemically infected plants. No differences were observed when concentrated virus particles were used for inoculation [START_REF] Zhang | Competition between isolates and variants of cauliflower mosaic virus in infected turnip plants[END_REF], whereas stochastic fluctuations were detectable when inoculum consisted in viral DNA prepared from infectious clones [START_REF] Riederer | The mode of cauliflower mosaic virus propagation in the plant allows rapid amplification of viable mutant strains[END_REF], presumably due to the lower infectivity of DNA preparations, engendering a stochastic founding effect in the latter case. These two studies indicate that indeed the inoculation process could induce unwanted fluctuations in repeated inoculations and prompted us to use virus particles enriched preparation for the purpose of our experiment (albeit development of systemic infection in between the two observation time points chosen in our protocol should not be affected by inoculation variations). Consistently, all 50 test plants in our study contained similar ratios of the six variants, suggesting very little stochastic variation during the inoculation process.

The values of p9 were estimated from individual leaf samples, each collected on a different test plant at a random position (between leaf position 5 and 12 above the initial harvest point). The remarkable observation that only relatively small variations in p9 are recorded in leaves from different plants (see standard deviation in Table 1) strongly indicates that very little variation occurs between leaves of the same plant. We therefore conclude that a similarly large founding population of several hundred viral genomes colonizes each newly formed leaf in plants systemically infected by CaMV.

The number of founder CaMV genomes in systemically infected leaves appears to be at least 100-fold higher than that determined in previously published studies for various RNA viruses [START_REF] Jridi | Distinct viral populations differentiate and evolve independently in a single perennial host plant[END_REF][START_REF] French | Evolution of Wheat streak mosaic virus: dynamics of population growth within plants may explain limited variation[END_REF][START_REF] Sacristan | Estimation of population bottlenecks during movement of tobacco mosaic virus in tobacco plants[END_REF]. Although there are clear differences in the inoculation and sampling protocols described in different studies, the presence of severe bottlenecks in WSMV [START_REF] French | Evolution of Wheat streak mosaic virus: dynamics of population growth within plants may explain limited variation[END_REF], TMV [START_REF] Sacristan | Estimation of population bottlenecks during movement of tobacco mosaic virus in tobacco plants[END_REF], CMV [START_REF] Li | Genetic bottlenecks reduce population variation in an experimental RNA virus population[END_REF] and PPV [START_REF] Jridi | Distinct viral populations differentiate and evolve independently in a single perennial host plant[END_REF], and their obvious absence in CaMV, can hardly be explained by artifacts of experimental designs. Indeed, at least in the present study and in three of the cited examples [START_REF] Li | Genetic bottlenecks reduce population variation in an experimental RNA virus population[END_REF][START_REF] Jridi | Distinct viral populations differentiate and evolve independently in a single perennial host plant[END_REF][START_REF] French | Evolution of Wheat streak mosaic virus: dynamics of population growth within plants may explain limited variation[END_REF], viruses moved not only from inoculated to systemic tissues, but also from systemically infected tissues to newly formed leaves or tiller. As mentioned in the Introduction, the previous repeated demonstration that different viruses undergo extreme bottlenecks during systemic infection of their host plant suggested that unavoidable barriers (e.g. connections between cells, etc.) may exist similarly for all plant viruses. The results presented here clearly demonstrate that such putative limiting barriers can be surmounted by some viruses and that the size of the viral population circulating in planta might thus be directly or indirectly controlled and needs to be evaluated for each virus individually. Whether the absence of severe bottlenecks demonstrated here for CaMV is common, and whether this phenomenon is related to the biology of the host plant, that of the virus, or more intricately on specific virus-host associations remains an open question.

The situation described here for CaMV re-opens the question as to what phenomena actually generate bottlenecks during virus infection. A reasonable hypothesis explaining the presence or absence of severe demographic bottlenecks would be the regulation of multiple infection of cells by several genomes of the viral population, in other words, regulation of the multiplicity of infection (MOI) of host cells. Indeed, it was proposed recently that bottlenecks seen in a PPV population infecting a Prunus tree could stem from the fact that viral genomes cannot secondarily invade tissues that have already been infected by a closely related genome [START_REF] Jridi | Distinct viral populations differentiate and evolve independently in a single perennial host plant[END_REF], thus preventing extensive mixing of genetic variants within host. Consistently, two isogenic PPV genomes, but expressing fluorescent proteins of different colours, revealed their mutual exclusion within infected tissues, with apparently rare coinfection of single cells [START_REF] Dietrich | Fluorescent labelling reveals spatial separation of potyvirus populations in mixed infected Nicotiana benthamiana plants[END_REF], a phenomenon also reported for other virus species, such as WSMW [START_REF] French | Evolution of Wheat streak mosaic virus: dynamics of population growth within plants may explain limited variation[END_REF], TMV [START_REF] Giritch | Rapid high-yield expression of full-size IgG antibodies in plants coinfected with noncompeting viral vectors[END_REF], and a few other RNA viruses [START_REF] Takahashi | Analysis of the spatial distribution of identical and two distinct virus populations differently labeled with cyan and yellow fluorescent proteins in coinfected plants[END_REF]. This phenomenon could logically induce a very low MOI, in turn engendering bottlenecks (and small Ne), due to competition for host ''territories'' between variants of the same population. Accordingly, under this hypothesis, the absence of severe bottlenecks found for CaMV would suggest a higher MOI that is totally consistent with the remarkably high within-host recombination rate described for this virus [START_REF] Froissart | Recombination every day: abundant recombination in a virus during a single multi-cellular host infection[END_REF], recombination being possible solely in multiply infected cells. We believe that explaining, at least in part, variations of within-plant bottlenecks in different virus species by their different capacity to multiply infect host cells is a very appealing hypothesis that will stimulate further novel research, the ''natural'' MOI during infection of a multi-cellular host being virtually unknown in viruses of both animals and plants.

c

  Standard deviation of marker frequencies across the 50 replicates. d Mean difference between initial and final frequencies in each replicate, and for each marker. All values are expressed as percent of the viral genome population. doi:10.1371/journal.ppat.1000174.t001

Table 1 .

 1 Genetic structure of initial and final CaMV populations.

		Initial populations a	Final populations a	E(Dp) d
		p b 610 2 (SD c . 10 2 )	p 9 b 610 2 (SD c . 10 2 )
	VIT1	41.6 (3.52)	45.8 (4.27)	4.2
	VIT2	4.2 (1.90)	3.3 (1.55)	20.9
	VIT3	42.5 (4.19)	40.9 (5.08)	21.6
	VIT4	10.2 (2.88)	9.2 (3.19)	21.0
	VIT5	1.7 (1.89)	1.9 (1.61)	0.2
	VIT6	1.7 (0.98)	3.6 (1.70)	1.9

Table 2 .

 2 Number of CaMV genomes founding the population in individual leaves, estimated through the analysis of the variance of the frequency of markers.Only markers VIT1, VIT3, and VIT4 were used for estimating Nv, for reasons explained in the text. b N v (and 95% confidence interval) is calculated according to Equation4(see Materials and Methods), where p is estimated from values given in Table1. doi:10.1371/journal.ppat.1000174.t002

					95% confidence
	Marker a	Var(p9)610 4	Var(p)610 4	N V b	interval
	VIT1	18.13	12.39	423.00	180-5034
	VIT3	25.79	17.60	298.45	122-3219
	VIT4	10.19	8.30	484.61	166-5193

a

Table 3 .

 3 Number of CaMV genomes founding the population in individual leaves, estimated through Fst statistics.Only markers VIT1, VIT3, and VIT4 were used for estimating N F , for reasons explained in the text. b N F was estimated from Fst values using Equation7(see Materials and Methods). The confidence interval, given in parentheses, was estimated by bootstrapping over plants. doi:10.1371/journal.ppat.1000174.t003

		Initial populations a	Final populations a
	H T	0.63530089	0.61474821
	H S	0.6314719	0.6093363
	Fst	0.006027	0.008803
	N F = 358 (184-1909) b
	a

PLoS Pathogens | www.plospathogens.org

October 2008 | Volume 4 | Issue 10 | e1000174

Acknowledgments

We are grateful to one anonymous referee for suggesting the Fst estimation method. We thank three referees for helpful comments and constructive criticism. Turnip seeds were kindly provided by the Takii Seed Company (Japan). We are thankful to the CaGeTE group for critical reading and discussion of the manuscript and to Jean-Luc Macia for all greenhouse work, plant production, and maintenance. YM acknowledges support from CNRS and IRD.

This work was funded by an INRA grant JE77, and a PhD grant for BM from Region Languedoc-Roussillon and INRA SPE. YM acknowledges support from CNRS and IRD.

Supporting Information

 S2Analysis of the initial inoculum Mix6VIT. a The mixture Mix6VIT was prepared and analysed by QSS as described in Materials and Methods. b p 0 = mean relative frequency determined from five independent repeats of QSS measurement. c Standard deviation calculated from the five repeated QSS analysis. All values are expressed as percent of the viral genome population. The distribution of the relative frequencies of all VIT1-6 variants in the initial Mix6VIT inoculum was globally similar to that in the 50 inoculated plants (VIT1<VIT3. VIT4.VIT2.VIT5<VIT6). Nevertheless, apart from VIT1 and VIT3 which increased by about 10% in frequency from inoculum to infected plants, all other variants equally decreased by approximately 5% (compare p 0 values in Table S2 and p values in Table 1). This phenomenon remains unclear and could be due to several different explanations, as for example: i) selection acting specifically at inoculation, could favour or disfavour some of the variants, ii) an undetermined threshold effect at the inoculation step could have positively and negatively affected the most and less frequent variants, respectively; iii) each variant in Mix6VIT originating from a different plant extract, they might have been differentially infectious due to unwanted and unequal damages of virus particles during extraction. Nevertheless, it is important to note that these considerations concern only the inoculation step, as only very minute changes in the mean frequency of all markers were detected later, over the 32 days separating the initial and final populations (see E(Dp) values in Table 1). Found at: doi:10.1371/journal.ppat.1000174.s002 (32 KB DOC)

Author Contributions


A Legarra

I Aguilar

Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information

Keywords: best linear unbiased predictor, genomic selection, single nucleotide polymorphism, genetic evaluation

Currently, genomic evaluations use multiple-step procedures, which are prone to biases and errors. A single-step procedure may be applicable when genomic predictions can be obtained by modifying the numerator relationship matrix A to H = A + A Δ , where A Δ includes deviations from expected relationships. However, the traditional mixed model equations require H -1 , which is usually difficult to obtain for large pedigrees. The computations with H are feasible when the mixed model equations are expressed in an alternate form that also applies for singular H and when those equations are solved by the conjugate gradient techniques. Then the only computations involving H are in the form of Aq or A Δ q, where q is a vector. The alternative equations have a nonsymmetric left-hand side. Computing A Δ q is inexpensive when the number of nonzeros in A Δ is small, and the product Aq can be calculated efficiently in linear time using an indirect algorithm. Generalizations to more complicated models are proposed. The data included 10.2 million final scores on 6.2 million Holsteins and were analyzed by a repeatability model. Comparisons involved the regular and the alternative equations. The model for the second case included simulated A Δ . Solutions were obtained by the preconditioned conjugate gradient algorithm, which works only with symmetric matrices, and by the bi-conjugate gradient stabilized algorithm, which also works with nonsymmetric matrices. The convergence rate associated with the nonsymmetric solvers was slightly better than that with the symmetric solver for the original equations, although the time per round was twice as much for the nonsymmetric solvers. The convergence rate associated with the alternative equations ranged from 2 times lower without A Δ to 3 times lower for the largest simulated A Δ . When the information attributable to genomics can be expressed as modifications to the numerator relationship matrix, the proposed methodology may allow the upgrading of an existing evaluation to incorporate the genomic information.

IntrODuCtIOn

Availability of dense molecular markers of type SNP led to the recent introduction of the genome-wide or genomic selection evaluation models. Those models are most often based on the simultaneous estimation of SNP marker effects a. Differences among methods are mostly on the a priori distribution of marker effects a [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Gianola | Genomic-assisted prediction of genetic value with semiparametric procedures[END_REF]. Efficient procedures exist for the computation of a, even for large data sets [START_REF] Legarra | Technical note: Computing strategies in genome-wide selection[END_REF].

The genomic evaluation is currently implemented as a multistep procedure. For example, an implementation for US dairy cattle [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Vanraden | Invited Review: Reliability of genomic predictions for North American Holstein bulls[END_REF] requires 3 steps: a) regular evaluation by the animal model, b) estimation of genomic effects for a relatively small number of genotyped animals, and c) estimation of genomic breeding values by a selection index. The elements in the index include a parent average or PTA from step a), genomic solutions from step b), and a parent average or PTA computed based on genotyped ancestors. Weights in the index are functions of heritability and accuracy. The marker-assisted selection program in France simultaneously fits QTL and polygenic effects, with weights depending on associated variance components [START_REF] Guillaume | Short communication: correlations of marker-assisted breeding values with progeny-test breeding values for eight hundred ninety-nine French Holstein bulls[END_REF].

Advantages of the multistage procedure include no change to the regular evaluations and simple steps for predicting genomic values for young genotyped animals. Disadvantages are requirements for parameters in steps b) and c) such as prior variances and weights, and loss of accuracy and biases attributable to selection. Whereas the model in a) uses the information on all animals and can be multitrait, the model in b) is equivalent to a single-trait sire model for a highly selected set of sires. Incorrect parameters in b) and c) can result in unexpected changes for high-reliability bulls. [START_REF] Neuner | Strategies for estimating genetic parameters in marker-assisted best linear unbiased predictor models in dairy cattle[END_REF] claimed that problems associated with the multistep procedure reduce its benefits, especially for cows.

VanRaden (2008) investigated 2 options for step b): "nonlinear," based on estimating effects attributable to SNP markers with a prior mixture distribution for those effects, and "linear," based on prior normal distribution for SNP markers. The latter is equivalent to using mixed model equations with a genomic relationship matrix. For most dairy traits, predictions based on the estimation of marker effects with nonlinear predictions were practically equivalent to linear predictions and thus to predictions with BLUP using a genomic relationship matrix [START_REF] Cole | Distribution and location of genetic effects for dairy traits[END_REF][START_REF] Vanraden | Invited Review: Reliability of genomic predictions for North American Holstein bulls[END_REF]. Therefore, using the genomic relationship matrix results in little or no loss of accuracy.

One way to simplify the multistep procedure is by incorporating the genomic information into step a), resulting in a single-step procedure. This could be accomplished by modifying the numerator relationship matrix A in that evaluation to include the genomic information. Such modifications are presented and discussed by [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] in a companion paper.

Assume that such a modification is known and that it involves relatively few elements of A. The mixed model equations require A -1 , which is very easy to create for large populations because of its sparsity and its special structure [START_REF] Henderson | A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values[END_REF]. However, obtaining the inverse of the modified matrix is likely to be impossible in general for large populations. This is not only because the cost of inversion is high, but also because A is dense and thus too large to store for large pedigrees. Thus, an approach using a modified A is of little value unless a feasible computing approach is available. The purpose of this study is to develop an efficient computing strategy to obtain solutions to mixed model equations in which the numerator relationship matrix is modified by a known matrix accounting for the genomic information.

materIaLS anD metHODS

Assume regular mixed model equations as used in a traditional genetic evaluation, for simplicity with only a single random effect:

y = Xb + Zu + e,
where y is a vector of records, b is a vector of fixed effects, and u is a vector of animal effects. Under a polygenic infinitesimal model of inheritance, and X and Z are appropriate incidence matrices.

Assume that the numerator relationship can be modified to account for genomic information:

H = A + A Δ ,
where A Δ is a matrix that can be stored explicitly, and H is the new modified matrix. In the simplest case, a genomic relationship matrix G replaces the numerator relationship matrix for the genotyped animals. Let indices 1 and 2 refer to ungenotyped and genotyped animals, respectively. Then . [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] proposed several H based on the partition of animals into several groups, including ungenotyped and genotyped animals. Although their different H are more complex than in the simple case, most quantities can be computed efficiently without any steps involving large matrix multiplications. Therefore, for simplicity of presentations, the following computing formulas assume the simple case above.

Solving Algorithm

Assume that G and A 22 are available. Temporarily assume that H is positive definite. The regular mixed model equations are ¢ Assume that the system of equations is solved using an algorithm that does not require the elements of LHS explicitly but only its product by a vector, say LHS q, as in the preconditioned conjugate gradient (PCG) iteration on data [START_REF] Tsuruta | Use of the preconditioned conjugate gradient algorithm as a generic solver for mixed-model equations in animal breeding applications[END_REF]. Then

¢ ¢ ¢ ¢ é ë ê ê ê ù û ú ú ú é ë ê ê ê ù û ú ú ú = ¢ ¢ é ë ê ê ê ù û ú ú ú - X X X Z Z X Z Z +

LHS q

X Xq + X Zq

Z Xq + Z Zq + H q c c + c 1 2 1 2 -1 2 1 2 3 = ¢ ¢ ¢ ¢ é ë ê ê ê ù û ú ú ú = é ë ê ê ê a ù ù û ú ú ú , with q q q c Z Xq Z Zq c H q RHS r r 1 2 2 1 2 3 -1 2 1 2 = é ë ê ê ê ù û ú ú ú = ¢ + ¢ = = é ë ê ê ê ; ; ; a ù ù û ú ú ú .
However, H -1 can be computed only for small populations; furthermore, H might be singular or close to singularity. [START_REF] Henderson | Applications of Linear Models in Animal Breeding[END_REF][START_REF] Henderson | Best linear unbiased prediction using relationship matrices derived from selected base populations[END_REF] and [START_REF] Harville | Extension of the Gauss-Markov theorem to include the estimation of random effects[END_REF] described an unsymmetric set of mixed model equations in which only H, not necessarily of full rank, is required: The new formulas do not include H -1 but include A Δ c 2 , Ac 2 , A Δ r 2 , and Ar 2. For the simplistic H, the first term can be computed directly at a low cost. The second term can also be computed inexpensively following the algorithm by [START_REF] Colleau | An indirect approach to the extensive calculation of relationship coefficients[END_REF] see Appendix A), which uses only the pedigree information and is completed in the amount of time proportional to the number of animals. The same algorithm also can be used to compute Ar 2 . Selected elements of A can be computed recursively, for example, by using the algorithm by Aguilar and Misztal (2008).

¢ ¢ ¢ ¢ é ë ê ê ê ù û ú ú ú é ë ê ê ê ù û ú ú ú = ¢ ¢ é ë ê ê ê ù û ú ú X X

More Complicated Models

Assume a multiple trait model, possibly with effects such as random regression or maternal. The regular mixed model equations for such models can be presented as

      + Ä é ë ê ê ê ù û ú ú ú é ë ê ê ê ù û ú ú ú ¢ é ë ê ê ê ù û ú ú ú G H u Z R y 0a -1 -1 -1 ˆ,
where parts not listed (…) are due to effects other than ˆ.

u By expanding the unsymmetric model by [START_REF] Henderson | Applications of Linear Models in Animal Breeding[END_REF] to multiple traits, the quantities needed for the iterations become

LHS q c G Ac G A c Iq M 1 0 0 2 = Ä + Ä + é ë ê ê ê ù û ú ú ú 2 2 D with RHS r G A G A r M 1 0 0 2 = Ä + Ä ( ) é ë ê ê ê ù û ú ú ú ∆ ,
where quantities c 1 and r 1 are now associated with all effects other than the additive.

Nonsymmetric Solvers

The presented system of equations is nonsymmetric and the matrix H may be semipositive definite. The PCG algorithm [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF] is applicable only to symmetric systems of equations. Therefore, it is important to find a suitable conjugate-gradient type algorithm and ensure that it would converge even with a poorly conditioned H. [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF] and [START_REF] Van Der Vorst | Iterative Krylov Methods for Large Linear Systems[END_REF] reviewed and presented several algorithms for solving the linear systems of equations. Based on their studies, the standard algorithm for solving sparse systems with nonsymmetric LHS is bi-conjugate gradient stabilized (Bi-CGSTAB; Van der Vorst, 1992; see Appendix B). This algorithm requires 2 LHS times a vector products per round as opposed to just one with PCG. When that product uses the majority of the computing time, Bi-CGSTAB is about twice as expensive as PCG per round of iteration.

Choice of Preconditioners

In initial tests (results not reported), Bi-CGSTAB converged very quickly with the unsymmetric equations for small models, but not for large ones. This was traced to large off-diagonal elements of the unsymmetric equations. The standard way in conjugate-gradient types of algorithms to improve convergence is by choice of a preconditioner M, which approximates LHS but is easily invertible (Van der Vorst, 2003). Then the system of equations solved is equivalent to M -1 LHS w = M -1 RHS, which has better numerical properties than the original system. The preconditioner is never used explicitly, but only in multiplications with a vector.

Assuming a diagonal preconditioner,

M L HS D 0 0 D 1 2 - - = ( ) = é ë ê ê ê ù û ú ú ú 1 1 diag ,
the LHS for regular equations after preconditioning is

M LHS D X X DX Z D Z X D Z Z A 1 1 2 2 -1 - = ¢ ¢ ¢ ¢ + ( ) é ë ê ù û ú é ë ê ê ê ê ù û ú ú ú ú 1 α .
The symmetry can be partially restored with a modified preconditioner:

M D 0 0 D A M 1 2 -1 -= é ë ê ê ê ù û ú ú ú 1 . Then M M LHS M D 1 X X D 1 X Z D 2 I A -1 A Z X D 2 I A -1 A - = ¢ ¢ + ae è ç ç ç ö ø ÷ ÷ ÷ ¢ é ë ê ê ù û ú ú + 1 ∆ ∆ ∆ ae è ç ç ç ö ø ÷ ÷ ÷ ¢ + é ë ê ê ù û ú ú é ë ê ê ê ê ê ù û ú ú ú ú ú = - + Z Z A -1 M LHS 0 0 D 2 A -1 a 1 A A Z X D 2 A -1 A Z Z ∆ ∆ ¢ ¢ é ë ê ê ê ê ù û ú ú ú ú .
When the genomic information is missing (A Δ = 0), the preconditioned left-hand side of the unsymmetric system of equations is the same as with the preconditioned regular equations. With the genomic information, the off-diagonal elements are likely to be small for small A Δ . The cost of the extra preconditioning is low because the product D 2 A -1 q, where q is a vector, can be done sequentially as D 2 (A -1 q).

Data

The data set included 10.5 million final scores on 6.2 million Holsteins as used for the recent genetic evaluation by the Holstein Association. Analyses were by a repeatability animal model. Two sets of mixed model equations were considered: regular and unsymmetric. For the second set, the genomic information was simulated for 5,000 randomly chosen animals as random numbers from the uniform distribution from 0 to b, where b was set to 0.0, 0.01, 0.03, and 0.05. For b = 0 there was no adjustment (A Δ = 0). Only positive adjustments were included to avoid some elements of H being negative. Solving algorithms were PCG (for the regular equations only) and Bi-CGSTAB. The first algorithm used a diagonal preconditioner. The second algorithm used the modified preconditioner because no convergence was achieved with the diagonal preconditioner. In all cases, the stopping criterion was set at 10 -12 . Computing was by the regular and modified program 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] and was carried out on an Opteron system running at 3 GHz.

reSuLtS anD DISCuSSIOn

The purpose of testing with the simulated genomic changes was to evaluate the computing feasibility of the method, and especially the robustness of the computing methodology. The results presented for the unsymmetric equations are only with the modified preconditioner. The Bi-CGSTAB diverged with the regular preconditioner and large data sets although it converged with small data sets. This is because products of A were very large for rows corresponding to popular bulls as all elements of A are positive; those products with A -1 are small because of cancellations; a contribution to a parent by a progeny in A -1 is proportional to [… 1.0 … -0.5 … -0.5 …], which sums to 0.

Table 1 shows the number of rounds and computing time with PCG and Bi-CGSTAB for the regular and unsymmetric equations and with varying magnitudes of simulated changes. For the regular equations, Bi-CGSTAB was slightly faster but took twice the computing time (26 vs. 13 s). Figure 1 shows the convergence pattern for the regular equations. Whereas the pattern for PCG shows small fluctuations, the pattern for Bi-CGSTAB has more abrupt changes. Some differences in the number of rounds to convergence may be due to differences in the convergence criteria. However, the differences in solutions were very small (correlations >0.99999).

For the unsymmetric equations with no simulated changes, the number of rounds approximately doubled and the computing time increased by 30% (from 26 to 34 s). Adding small simulated changes (b = 0.01) increased the computing time per round by 10% (from 34 to 37 s) and slightly deteriorated convergence. The number of rounds increased by about 30% when changes were increased to b = 0.03 and again by 10% when changes were increased by b = 0.05. Figure 2 shows the convergence pattern for the modified equations and b = 0.0 and 0.03. Much larger fluctuations than with the regular equations were observed, which may have been due to a more complex preconditioner. For a multipletrait random regression model, Aguilar et al. (2008) observed much larger fluctuations in the convergence pattern with a block-diagonal preconditioner as compared with a diagonal one. Additional computations will be necessary in practical applications of the method with the real genomic relationship matrix. For simple H, additional steps include the multiplications of G-A 22 and A by a vector. The last one can be done efficiently using the algorithm of [START_REF] Colleau | An indirect approach to the extensive calculation of relationship coefficients[END_REF] in linear time (see Appendix A). The cost of this algorithm is equal to scanning the pedigree file twice and is small, especially with pedigrees in memory. [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] presented formulas for more realistic H and also computing details for a product of that H by a vector. With such a product, the only components that cannot be computed in linear but rather in quadratic time (matrix-vector multiplication) are those corresponding to G and possibly those due to A 22 . If A 22 needs to be available explicitly, it can be computed by the method of Aguilar and Misztal (2008). When applied to 17 million Holsteins, that method calculated about 80,000 inbreeding coefficients/s. Assuming that computing one relationship costs no more than computing one inbreeding coefficient, on average, the computation of A 22 for 20,000 animals would take 40 min. Alternatively, A 22 can be computed by the repeated applications of the algorithm of [START_REF] Colleau | An indirect approach to the extensive calculation of relationship coefficients[END_REF], in which the vector to multiply by would contain one 1.0 and zeros elsewhere.

When the number of genotyped animals is very high, say >50,000, storage and computations with matrix G and possibly A 22 can be quite involving. A few choices may be applicable. First, some computations may easily be done in parallel. Current computers routinely include 4 processors (cores) per processor module, and computers with 4 modules are readily available. Second, some elements in A Δ may be very small or unimportant and could be neglected. Neglecting small elements in the computation of sparse inverse for the purpose of calculating accuracies reduced the computations by 50 times while retaining high precision [START_REF] Thompson | Calculation of prediction error variances using sparse matrix methods[END_REF]. Finally, genotypes of some animals may be unimportant and do not have to be included.

In summary, we have demonstrated that mixed model equations with small modifications to the numerator relationship matrix can be solved efficiently by conjugate-gradient type algorithms. Only a few modifications may be required for existing programs using the PCG algorithm.

APPENDIX A

Below we show how to create a product Aq, where A is the numerator relationship matrix and q is a vector. The recurrence equation for the additive effect is

a = Pa + ϕ,
where a is a vector of animals ordered from oldest to youngest, ϕ is a diagonal matrix of Mendelian samplings, and P is a matrix relating animals to their parents; this matrix has at most 2 elements per row, both equal to 0.5. Following [START_REF] Quaas | Additive genetic model with groups and relationships[END_REF], Var(a) = A = (I -P) -1 D(I -P) -1′ , where D = var(ϕ). [START_REF] Colleau | An indirect approach to the extensive calculation of relationship coefficients[END_REF] showed that the product of A by a vector, for example, v = Aq = (I -P) -1 D(I -P) -1 ′q = (I -P) -1 D [(I -P) -1 ′q], can be solved in linear time. In particular, quantities r = (I -P) -1 ′q and v = (I -P) -1 Dr can be obtained by solving (I -P)′r = q and (I -P)v = Dr, each one in a single sweep because (I -P) is triangular. The scalar formulas are r i = r i + q i ; r si = r si + r i /2; r di = r di + r i /2; i = n, ..., 1 v i = d i r i + (v si + v di )/2, i = 1, ..., n, where s i and d i are positions of the sire and dam of animal i, respectively. The [START_REF] Colleau | An indirect approach to the extensive calculation of relationship coefficients[END_REF] algorithm can be used to compute products of sections of matrices. For instance, the products below show how to compute A 12 q, A 22 q, A 21 q, or A 22 q:

  s α , where A is the numerator relationship matrix based on pedigree. Furthermore,

  RHS using the usual notation, where LHS and RHS are the left-and right-hand side, and w = b u ˆˆ.

  preconditioned conjugate gradient; Bi-CGSTAB = bi-conjugate gradient stabilized.2 Changes in relationships simulated from uniform (0, b) distribution for 5,000 randomly selected animals.

Figure 1 .

 1 Figure 1. Convergence rate for the preconditioned conjugate gradient (PCG) and bi-conjugate gradient stabilized (Bi-CGSTAB) algorithms with the symmetric system of equations.

aCKnOWLeDGmentS

  Discussions with Rohan Fernando (Department of Animal Science, Iowa State University, Ames), Jeff O'Connell (University of Maryland School of Medicine, 4653 JoiNt geNoMiC aND PHeNotyPiC eValUatioN Journal of Dairy Science Vol. 92 No. 9, 2009

Figure 2 .

 2 Figure 2. Convergence rate for the bi-conjugate gradient stabilized (Bi-CGSTAB) algorithm with the unsymmetric system of equations and either no (b = 0) or a middle level (b = 0.03) of simulated changes attributable to the genomic relationship matrix.

Table 1 .

 1 The number of rounds (computing time per round in seconds) for different computing algorithms and different magnitudes of modification to the numerator relationship matrix

	Solving algorithm 1

JoiNt geNoMiC aND PHeNotyPiC eValUatioNJournal of Dairy Science Vol. 92 No. 9, 2009

Journal of Dairy Science Vol. 92 No. 9, 2009 MiSztal et al.

APPENDIX B

The pseudo-program below implements the Bi-CG-STAB (Van der Vorst, 1992) for a system of equations Ax = b with M being a preconditioner. The major expenses in the algorithm are products of A by a vector, possibly followed by products of M -1 , but only if M is of complex structure.

Compute r (0) = b -Ax (0) for some initial guess