The elegance of the Breeder’s Equation

# The elegance of the Breeder’s Equation

03 May 2022

The Breeder’s Equation is a well-known expression in quantitative genetics that is widely used to understand how genetic gain can be achieved. Here, we will focus on a few statistical aspects of this equation and the role they play in helping us to maximize genetic gain. An expression of this equation is:

x

where,

is the genetic gain per period of time (often in units/year),

is the narrow- or broad-sense heritability,

is the selection differential (in units),

is the period or cycle (often in years).

Let’s start with . This represents the genetic gain on an annual basis, and its units will depend on the target trait. Often, we use yield/year but can also correspond to an index from a trait composite with different weights. In our discussion, we will assume that we have estimated breeding values (EBVs) for a group of individuals that will be parents to constitute the next generation (breeding population) or for commercial use (production population). Hence, the focus is on additive effects. Alternatively, we could use total genetic values (TGVs). These are common in plant breeding, and correspond to effects that quantify the total genetic worth of an individual (often a clone), and therefore it contains both additive and non-additive effects.

The most important aspect of the above equation is the heritability (). If it is exactly one, then we have the perfect situation of ‘what we see is what we get’. That is, phenotypic selection of the top individuals will yield the best genotypes. However, this is never the case, and we have a range of heritabilities often from as low as 0.1 up to 0.8. For simplicity, we will consider the following simple heritability expression:

Narrow-sense:    / +

where, , and are the additive, total genetic and residual variances, respectively.

Here, anything aiming to reduce will increase heritability and, in turn, produce greater genetic gains. Statistically, there are multiple options to consider to reduce . For example, the use of a sound and optimized experimental design will help to eliminate bias and/or have better control of background noise. Also, on the linear mixed model analysis side, we can use, for example, spatial analyses on field trials, or inclusion of relevant nuisance random or fixed effects, to better model the data.

In practical terms there are many other options to reduce this background noise. One is to have better and more careful phenotypic measurements. This implies clear, consistent and well-defined protocols, but also, as in some cases, the use of more expensive measurement procedures and well-trained people. For example, use of sophisticated laboratory techniques to measure protein or sugar content. In field trials, good selection and preparation of experimental sites is also critical to ensure more homogeneous conditions for plants.

All of the above elements affect equivalently EBVs or TGVs, but for the former the inclusion of pedigree is important to improve the accuracy of these random effects. This pedigree should be of top quality, as any mistakes or inconsistencies (e.g., unrecorded pollen contamination) will result in increased background noise.

Additionally, the heritability from the Breeder’s Equation can take many definitions. An important one is associated with the use of the genotype’s mean based on several phenotyped replications (as done with clones), or multiple measurements over time (as done on trees and dairy cows). To illustrate this, consider the following expressions:

Narrow-sense:     / +

where and are the mean broad- and narrow-sense heritabilities, respectively; and is the replication (or number of repeated measurements). Note that the replication reduces the magnitude of the background noise in a direct way; therefore, any efforts to increase replication will increase heritability (unless  is zero). This implies, that in plant breeding for example, having an adequate number of replications will increase genetic gain. Defining the optimal number of replications or measurements over time of an organism is another statistical aspect to plan carefully, and is associated with good experimental design practices and helped with simulations.

The differential of selection, , has obvious implications. Its expression is: , where is the mean of the selected individuals and is the overall mean of the population under selection. Hence, the further away the selected individuals are from the overall mean, the greater the potential for large genetic gains. This seems a simple strategy that can occur in at least the two ways described below.

One way is to select very few (say < 5) genotypes. But this is difficult for most breeding programs, as there is a restriction on a minimum number of future progenitors to generate the next generation; and they also need to present a wide range of genetic diversity to ensure the long-term success of the program.

Alternatively, the other strategy is to play a numbers game, where we have a very large pool (thousands or millions) of genotypes to select from. Large numbers are required in order to increase the chance of finding an outstanding combination of alleles. Note that this, in order to be successful, still requires a non-zero heritability. This strategy has several side effects, such as high operational costs, large experimental study size, and in general complex logistics. Nevertheless, in practical terms we should always aim at having as many individuals as possible to select from.

The last component is , which is the period, cycle or generation interval. Depending on the organism it might or might not be possible to modify this component, as some breeding programs are limited on the reproduction cycle, and/or the time required to collect the phenotypic data. For example, milk production of sires requires to phenotype female offspring. Similarly with trees, a wait of several years is required to perform crosses and collect seeds on mature trees. But, in some cases, it might be possible to induce early maturation with hormones, thus reducing this cycle for all or some of the individuals. Alternatively, selecting genotypes early (with or without phenotypic data) speeds up this cycle, and it has been a successful strategy in perennial plants. However, this may imply a loss in the accuracy of the EBVs (or TGVs) and/or a reduction in heritability.

The emergence of genomic selection (GS) has also added interesting alternatives to the components of the Breeder’s Equation. For example, under the numbers game mentioned above, it is possible to genotype thousands of individuals for which we obtain genomic predictions based on a previously trained model. This will clearly increase the pool of individuals for selection, but at an important economic cost associated with genotyping. Additionally, genomic predictions can be implemented on a pool of individuals that are still immature, and this set can be used as soon as possible as parents, even before having phenotypic data associated with these genotypes.

Genomic prediction models are not perfect, and they add an extra complication to the genetic gain calculations. In order to understand this, we will use the expression of indirect genetic gain below:

x x    x

where,

is the genetic gain per period of time for the target (unmeasured) trait,

is the heritability of the target trait,

is the heritability of the indirect trait,

is the additive correlation between the target and indirect trait,

is the selection differential based on the indirect trait,

is the period or cycle.

This expression is modified by considering that GS generates an indirect predicted genetic value for each individual based on molecular data that is not perfectly accurate (with accuracy defined as the correlation between the true and the predicted EBVs). As we are dealing with the same trait,  = , and we write:

= x x

where, refers to the genetic gains based on the genomic prediction model available, and is the accuracy of this model (as defined before).

As is less than one, then we always have a loss of genetic gain for our trait, and this is going to depend on how far from one we are. Most common accuracy values range between 0.2 to 0.4; hence, there could be a considerable potential loss of genetic gain. Nevertheless, as indicated earlier, this can be overcome by genotyping a greater number of individuals, or by reducing the cycle given the availability of predictions at early stages, among many options not described here.

Statistically, there are many possible ways to increase the accuracy of genomic models. Some include larger training populations, more SNPs markers on the genotyping panel, close relatedness between training and evaluation populations, better or more appropriate genomic techniques for fitting models (such as GBLUP or Bayes B), to name a few. It is not our objective to focus on the above or other aspects, but it is important to mention that there are many alternatives to increase the accuracy of our models, and they will have a relevant effect on increasing genetic gains.

The last aspect that it is important to mention about our Breeder’s Equation is that exploiting its benefits depends on costs and careful economical evaluations. Most options to increase genetic gain are associated with greater costs. For example, in some countries it is cheaper to improve the quality of the phenotyping rather than increasing the number of genotypes to evaluate. All of these, and other, elements need to be considered specifically for the reality of each commercial breeding program and they are often complemented with detailed simulation studies. But, as shown here, the Breeder’s Equation presents an easy way to enumerate the elements that are worth assessing so we can make the most of our resources and statistical tools available.

Dr. Salvador Gezan is a statistician/quantitative geneticist with more than 20 years’ experience in breeding, statistical analysis and genetic improvement consulting. He currently works as a Statistical Consultant at VSN International, UK. Dr. Gezan started his career at Rothamsted Research as a biometrician, where he worked with Genstat and ASReml statistical software. Over the last 15 years he has taught ASReml workshops for companies and university researchers around the world.

Dr. Gezan has worked on agronomy, aquaculture, forestry, entomology, medical, biological modelling, and with many commercial breeding programs, applying traditional and molecular statistical tools. His research has led to more than 100 peer reviewed publications, and he is one of the co-authors of the textbook Statistical Methods in Biology: Design and Analysis of Experiments and Regression.

Kanchana Punyawaew

01 March 2021

Linear mixed models: a balanced lattice square

This blog illustrates how to analyze data from a field experiment with a balanced lattice square design using linear mixed models. We’ll consider two models: the balanced lattice square model and a spatial model.

The example data are from a field experiment conducted at Slate Hall Farm, UK, in 1976 (Gilmour et al., 1995). The experiment was set up to compare the performance of 25 varieties of barley and was designed as a balanced lattice square with six replicates laid out in a 10 x 15 rectangular grid. Each replicate contained exactly one plot for every variety. The variety grown in each plot, and the coding of the replicates and lattice blocks, is shown in the field layout below:

There are seven columns in the data frame: five blocking factors (Rep, RowRep, ColRep, Row, Column), one treatment factor, Variety, and the response variate, yield.

The six replicates are numbered from 1 to 6 (Rep). The lattice block numbering is coded within replicates. That is, within each replicates the lattice rows (RowRep) and lattice columns (ColRep) are both numbered from 1 to 5. The Row and Column factors define the row and column positions within the field (rather than within each replicate).

### Analysis of a balanced lattice square design

To analyze the response variable, yield, we need to identify the two basic components of the experiment: the treatment structure and the blocking (or design) structure. The treatment structure consists of the set of treatments, or treatment combinations, selected to study or to compare. In our example, there is one treatment factor with 25 levels, Variety (i.e. the 25 different varieties of barley). The blocking structure of replicates (Rep), lattice rows within replicates (Rep:RowRep), and lattice columns within replicates (Rep:ColRep) reflects the balanced lattice square design. In a mixed model analysis, the treatment factors are (usually) fitted as fixed effects and the blocking factors as random.

The balanced lattice square model is fitted in ASReml-R4 using the following code:

> lattice.asr <- asreml(fixed = yield ~ Variety,
random = ~ Rep + Rep:RowRep + Rep:ColRep,
data=data1)

The REML log-likelihood is -707.786.

The model’s BIC is:

The estimated variance components are:

The table above contains the estimated variance components for all terms in the random model. The variance component measures the inherent variability of the term, over and above the variability of the sub-units of which it is composed. The variance components for Rep, Rep:RowRep and Rep:ColRep are estimated as 4263, 15596, and 14813, respectively. As is typical, the largest unit (replicate) is more variable than its sub-units (lattice rows and columns within replicates). The "units!R" component is the residual variance.

By default, fixed effects in ASReml-R4 are tested using sequential Wald tests:

In this example, there are two terms in the summary table: the overall mean, (Intercept), and Variety. As the tests are sequential, the effect of the Variety is assessed by calculating the change in sums of squares between the two models (Intercept)+Variety and (Intercept). The p-value (Pr(Chisq)) of  < 2.2 x 10-16 indicates that Variety is a highly significant.

The predicted means for the Variety can be obtained using the predict() function. The standard error of the difference between any pair of variety means is 62. Note: all variety means have the same standard error as the design is balanced.

Note: the same analysis is obtained when the random model is redefined as replicates (Rep), rows within replicates (Rep:Row) and columns within replicates (Rep:Column).

### Spatial analysis of a field experiment

As the plots are laid out in a grid, the data can also be analyzed using a spatial model. We’ll illustrate spatial analysis by fitting a model with a separable first order autoregressive process in the field row (Row) and field column (Column) directions. This is often a useful model to start the spatial modeling process.

The separable first order autoregressive spatial model is fitted in ASReml-R4 using the following code:

> spatial.asr <- asreml(fixed = yield ~ Variety,
residual = ~ar1(Row):ar1(Column),
data = data1)

The BIC for this spatial model is:

The estimated variance components and sequential Wald tests are:

The residual variance is 38713, the estimated row correlation is 0.458, and the estimated column correlation is 0.684. As for the balanced lattice square model, there is strong evidence of a Variety effect (p-value < 2.2 x 10-16).

A log-likelihood ratio test cannot be used to compare the balanced lattice square model with the spatial models, as the variance models are not nested. However, the two models can be compared using BIC. As the spatial model has a smaller BIC (1415) than the balanced lattice square model (1435), of the two models explored in this blog, it is chosen as the preferred model. However, selecting the optimal spatial model can be difficult. The current spatial model can be extended by including measurement error (or nugget effect) or revised by selecting a different variance model for the spatial effects.

#### References

Butler, D.G., Cullis, B.R., Gilmour, A. R., Gogel, B.G. and Thompson, R. (2017). ASReml-R Reference Manual Version 4. VSN International Ltd, Hemel Hempstead, HP2 4TP UK.

Gilmour, A.R., Thompson, R. & Cullis, B.R. (1995). Average Information REML, an efficient algorithm for variance parameter estimation in linear mixed models. Biometrics, 51, 1440-1450.

09 March 2022

Meta analysis using linear mixed models

Meta-analysis is a statistical tool that allows us to combine information from related, but independent studies, that all have as an objective to estimate or compare the same effects from contrasting treatments. Meta-analysis is widely used in many research areas where an extensive literature review is performed to identify studies that had a similar research question. These are later combined using meta-analysis to estimate a single combined effect. Meta-analyses are commonly used to answer healthcare and medical questions, where they are widely accepted, but they also are used in many other scientific fields.

By combining several sources of information, meta-analyses have the advantage of greater statistical power, therefore increasing our chance of detecting a significant difference. They also allow us to assess the variability between studies, and help us to understand potential differences between the outcomes of the original studies.

The underlying premise in meta-analysis is that we are collecting information from a group of, say n, studies that individually estimated a parameter of interest, say . It is reasonable to consider that this parameter has some statistical properties. Mainly we assume that it belongs to a Normal distribution with unknown mean and variance. Hence, mathematically we say:

In meta-analysis, the target population parameter θ can correspond to any of several statistics, such as a treatment mean, a difference between treatments; or more commonly in clinical trials, the log-odds ratio or relative risk.

There are two models that are commonly used to perform meta-analyses: the fixed-effect model and the random-effects model. For the fixed-effect model, it is assumed that there is only a single unique true effect our single θ above, which is estimated from a random sample of studies. That is, the fixed-effect model assumes that there is a single population effect, and the deviations obtained from the different studies are only due to sampling error or random noise. The linear model (LM) used to describe this process can be written as:

where is each individual observed response, is the population parameter (also often known as  μ, the overall mean), and is a random error or residual with assumptions of . The variance component is a measurement of our uncertainty in the information (i.e., response) of each study. The above model can be easily fitted under any typical LM routine, such as R, SAS, GenStat and ASReml.

For the random-effects model we still assume that there is a common true effect between studies, but in addition, we allow this effect to vary between studies. Variation between these effects is a reasonable assumption as no two studies are identical, differing in many aspects; for example, different demographics in the data, slightly differing measurement protocols, etc. Because, we have a random sample of studies, then we have a random sample of effects, and therefore, we define a linear mixed model (LMM) using the following expression:

where, as before, is each individual observed response, is the study-specific population parameter, with the assumption of and is a random error or residual with the same normality assumptions as before. Alternatively, the above LMM can be written as:

where and is a random deviation from the overall effect mean θ with assumptions .

This is a LMM because we have, besides the residual, an additional random component that has a variance component associated to it, that is or . This variance is a measurement of the variability ‘between’ studies, and it will reflect the level of uncertainty of observing a specific  . These LMMs can be fitted, and variance components estimated, under many linear mixed model routines, such as nlme in R, proc mixed in SAS, Genstat or ASReml.

Both fixed-effect and random-effects models are often estimated using summary information, instead of the raw data collected from the original study. This summary information corresponds to estimated mean effects together with their variances (or standard deviations) and the number of samples or experimental units considered per treatment. Since the different studies provide different amounts of information, weights should be used when fitting LM or LMM to summary information in a meta-analysis, similar to weighted linear regression. In meta-analysis, each study has a different level of importance, due to, for example, differing number of experimental units, slightly different methodologies, or different underlying variability due to inherent differences between the studies. The use of weights allows us to control the influence of each observation in the meta-analysis resulting in more accurate final estimates.

Different statistical software will manage these weights slightly differently, but most packages will consider the following general expression of weights:

where is the weight and is the variance of the observation. For example, if the response corresponds to an estimated treatment mean, then its variance is , with MSE being the mean square error reported for the given study, and the number of experimental units (or replication).

Therefore, after we collect the summary data, we fit our linear or linear mixed model with weights and request from its output an estimation of its parameters and their standard errors. This will allow us to make inference, and construct, for example, a 95% confidence interval around an estimate to evaluate if this parameter/effect is significantly different from zero. This will be demonstrated in the example below.

### Motivating example

The dataset we will use to illustrate meta-analyses was presented and previously analysed by Normand (1999). The dataset contains infromation from nine independent studies where the length of hospitalisation (measured in days) was recorded for stroke patients under two different treatment regimes. The main objective was to evaluate if specialist inpatient stroke care (sc) resulted in shorter stays when compared to the conventional non-specialist (or routine management) care (rm).

The complete dataset is presented below, and it can also be found in the file STROKE.txt. In this table, the columns present for each study are the sample size (n.sc and n.rm), their estimated mean value (mean.sc and mean.rm) together with their standard deviation (sd.sc and sd.rm) for both the specialist care and routine management care, respectively.

### Statistical analyses

We will use the statistical package R to read and manipulate the data, and then the library ASReml-R (Butler et al. 2017) to fit the models.
First, we read the data in R and make some additional calculations, as shown in the code below:

 STROKE <- read.table("STROKE.TXT", header=TRUE)
 STROKE$diff <- STROKE$mean.sc - STROKE$mean.rm STROKE$Vdiff <- (STROKE$sd.sc^2/STROKE$n.sc) + (STROKE$sd.rm^2/STROKE$n.rm) STROKE$WT <- 1/(STROKE$Vdiff)

The new column diff contains the difference between treatment means (as reported from each study). We have estimated the variance of this mean difference, Vdiff, by taking from each treatment its individual MSE (mean square error) and dividing it by the sample size, and then summing the terms of both treatments. This estimate assumes, that for a given study, the samples from both treatments are independent, and for this reason we did not include a covariance. Finally, we have calculated a weight (WT) for each study as the inverse of the variance of the mean difference (i.e., 1/Vdiff).

We can take another look at this data with these extra columns:

The above table shows a wide range of values between the studies in the mean difference of length of stay between the two treatments, ranging from as low as −71.0 to 11.0, with a raw average of −15.9. Also, the variances of these differences vary considerably, which is also reflected in their weights.

The code to fit the fixed-effect linear model using ASReml-R is shown below:

 library(asreml) meta_f<-asreml(fixed=diff~1,                weights=WT,                family=asr_gaussian(dispersion=1),                data=STROKE)

In the above model, our response variable is diff, and the weights are indicated by the variate WT. As the precisions are contained within the weights the command family is required to fix the residual error (MSE) to exactly 1.0, hence, it will not be estimated.

The model generates output that can be used for inference. We will start by exploring our target parameter, i.e. θ, by looking at the estimated fixed effect mean and its standard error. This is done with the code:

 meta_effect <- summary(meta_f, coef=TRUE)$coef.fixed Resulting in the output: The estimate of θ is equal to −3.464 days, with a standard error of 0.765. An approximate 95% confidence interval can be obtained by using a z-value of 1.96. The resulting approximate 95% confidence interval [−4.963;−1.965] does not contain zero. The significance of this value can be obtained by looking at the approximated ANOVA table using the command:  wald.asreml(meta_f) Note that this is approximated as, given that weights are considered to be known, the degrees of freedom are assumed to be infinite; hence, this will be a liberal estimate. The results from this ANOVA table indicate a high significance of this parameter (θ) with an approximated p-value of < 0.0001. Therefore, in summary, this fixed effect model analysis indicates a strong effect of the specialised care resulting in a reduction of approximately 3.5 days in hospitalisation. However, as indicated earlier, a random-effects model might seem more reasonable given the inherent differences in the studies under consideration. Here, we extend the model to include the random effect of study. In order to do this, first we need to ensure that this is treated as a factor in the model by running the code:  STROKE$study <- as.factor(STROKE$study)_f) The LMM to be fitted using ASReml-R is:  meta_r<-asreml(fixed=diff~1,  random=~study,  weights=WT,  family=asr_gaussian(dispersion=1),  data=STROKE) Note in this example the only difference from the previous code is the inclusion of the line random=~study. This includes the factor study as a random effect. An important result from running are the variance component estimates. These are obtained with the command:  summary(meta_r)$varcomp

In this example, the variance associated with the differences in the target parameter (θ) between the studies is 684.62. When expressed as a standard deviation, this corresponds to 26.16 days. Note that this variation is large in relation to the scale of the data, reflecting large differences between the random sample of studies considered in the meta-analysis.

We can output the fixed and random effects using the following commands:

 meta_effect <- summary(meta_r, coef=TRUE)$coef.fixed BLUP <- summary(meta_r, coef=TRUE)$coef.random

Note that now that our estimated mean difference corresponds to −15.106 days with an standard error of 8.943, and that the approximate 95% confidence interval [−32.634;2.423] now contains zero. An approximated ANOVA based on the following code:

 wald.asreml(meta_r)

results in the output:

We have a p-value of 0.0912, indicating that there is no significant difference in length of stay between the treatments evaluated. Note that the estimates of the random effects of study, also known as BLUPs (best linear unbiased predictions) are large, ranging from −45.8 to 22.9, and widely variable. The lack of significance in the random-effects model, when there is a difference of −15.11 days, is mostly due to the large variability of 684.62 found between the different studies, resulting in a substantial standard error for the estimated mean difference.

In the following graph we can observe the 95% confidence intervals for each of the nine studies together with the final parameter estimated under the Random-effects Model. Some of these confidence intervals contain the value zero, including the one for the random-effects model. However, it can be observed that the confidence interval from the random-effects model is an adequate summarization of the nine studies, representing a compromising confidence interval.

An important aspect to consider is the difference in results between the fixed-effect and the random-effects model that are associated, as indicated earlier, with different inferential approaches. One way to understand this is by considering what will happen if a new random study is included. Because we have a large variability in the study effects (as denoted by ), we expect this new study to have a difference between treatments that is randomly within this wide range. This, in turn, is expressed by the large standard error of the fixed effect θ, and by its large 95% confidence interval that will ensure that for ‘any’ observation we cover the parameter estimate 95% of the time. Therefore, as shown by the data, it seems more reasonable to consider the random-effects model than the fixed-effect model as it is an inferential approach that deals with several sources of variation.

### Summary

In summary, we have used the random-effects model to perform meta-analysis on a medical research question of treatment differences by combining nine independent studies. Under this approach we assumed that all studies describe the same effect but we allowed for the model to express different effect sizes through the inclusion of a random effect that will vary from study to study. The main aim of this analysis was not to explain why these differences occur, here, our aim was to incorporate a measure of this uncertainty on the estimation of the final effect of treatment differences.

There are several extensions to meta-analysis with different types of responses and effects. Some of the relevant literature recommended to the interested reader are van Houwelingen et al. (2002) and Vesterinen et al. (2014). Also, a clear presentation with further details of the differences between fixed-effect and random-effects models is presented by Borenstein et al. (2010).

Dataset: STROKE.txt
R code: STROKE_METAA.R

### References

Borenstein, M; Hedges, LV; Higgins, JPT; Rothstein, HR. 2010. A basic introduction to fixed-effect and random-effects models for meta-analysis. Research Synthesis Methods 1: 97-111.

Butler, DG; Cullis, BR; Gilmour, AR; Gogel, BG; Thompson, R. 2017. ASReml-R Reference Manual Version 4. VSNi International Ltd, Hemel Hempstead, HP2 14TP, UK.

Normand, ST. 1999. Meta-analysis: Formulating, evaluating, combining, and reporting. Statistics in Medicine 18: 321-359.

van Houwelingen, HC; Arends, LR; Stignen, T. 2002. Advanced methods in meta-analysis: multivariate approach and meta-regression. Statistics in Medicine 21: 589-624.

Vesterinen, HM; Sena, ES; Egan, KJ; Hirst, TC; Churolov, L; Currie, GL; Antonic, A; Howells, DW; Macleod, MR. 2014. Meta-analysis of data from animal studies: a practical guide. Journal of Neuroscience Methods 221: 92-102.

Salvador Gezan is a statistician/quantitative geneticist with more than 20 years’ experience in breeding, statistical analysis and genetic improvement consulting. He currently works as a Statistical Consultant at VSN International, UK. Dr. Gezan started his career at Rothamsted Research as a biometrician, where he worked with Genstat and ASReml statistical software. Over the last 15 years he has taught ASReml workshops for companies and university researchers around the world.

Dr. Gezan has worked on agronomy, aquaculture, forestry, entomology, medical, biological modelling, and with many commercial breeding programs, applying traditional and molecular statistical tools. His research has led to more than 100 peer reviewed publications, and he is one of the co-authors of the textbook “Statistical Methods in Biology: Design and Analysis of Experiments and Regression”.

Dr. Andrew Illius and Dr. Nick Savill

20 July 2022

Quantification holds the key to controlling disease

Background

Andrew Illius with Nick Savill have, since 2018, studied the epidemiology and control of maedi-visna virus (MV) in sheep and have been looking at understanding and finding ways of controlling this incurable disease. Accessing published data and with the use of Genstat, they aimed to find ways of controlling MV.

When one of your sheep gets diagnosed with an incurable disease, you have to worry. How quickly will symptoms develop, and welfare and productivity suffer? And how soon will it spread throughout the rest of the flock? The disease in question is maedi-visna (MV, see Box 1), notorious for its impact in Iceland, where the disease was first described; extreme measures over 20 years were required before it was finally eliminated. Culling seropositive animals is the main means of control. For the farmer, the crucial question is whether living with the disease would be more expensive than trying to eradicate it. We are addressing such questions by analysing data from long-term experiments.

#### 1            MV – the tip of an iceberg?

Putting aside for a moment MV’s fearsome reputation, the way the pathogen works is fascinating. The small ruminant lentiviruses (SRLV, family retroviridae) are recognised as a heterogeneous group of viruses that infect sheep, goats and wild ruminants. Lentiviruses target the immune system, but SRLV does not target T-cells in the manner of immune deficiency lentiviruses such as HIV. Instead, SRLV infects monocytes (a type of white blood cell) which infiltrate the interstitial spaces of target organs (such as the lungs, mammary glands, or the synovial tissue of joints) carrying proviral DNA integrated into the host cell genome and hence invisible to the immune system. Virus replication commences following maturation of monocytes into macrophages, and the ensuing immune response eventually shows up as circulating antibodies (termed seroconversion). But it also causes inflammation that attracts further macrophages, slowly and progressively building into chronic inflammatory lesions and gross pathology. These take years to present clinical symptoms, hence the name lentivirus (from the Latin lentus, slow). By the time clinical signs become evident in a flock, the disease will have become well-established, with perhaps 30-70% of the flock infected. That is why MV is called one of the iceberg diseases of sheep – for every obviously affected individual, there are many others infected, but without apparent symptoms.

A large body of research into the pathology, immunology and molecular biology of small ruminant lentiviruses (SRLV) exists, as might be expected given its economic significance, welfare implications and its interest as a model for HIV. The main route of transmission of the virus is thought to be horizontal, via exhaled droplets of the highly infectious fluid from deep in the lungs of infected animals, suggesting a risk from prolonged close contact, for example in a sheep shed. But despite all the research into disease mechanisms, we were surprised to find that there has been almost no quantitative analysis of SRLV epidemiology, nor even an estimation of the rate of SRLV transmission under any management regime. So, our first foray into the data aimed to rectify this

#### Quantification to the rescue

We found an experiment published in 1987 with excellent detail on a five-year timecourse of seroconversions in a small infected sheep flock, and a further trawl of the Internet netted a PhD thesis that built on this with a focussed experiment. Karianne Lievaart-Peterson, its author, runs the Dutch sheep health scheme and became a collaborator. We also worked with Tom McNeilly, an immunologist at the Moredun Research Institute.

Nick Savill, a mathematical epidemiologist at Edinburgh University, did the hard work of developing and parameterising a mathematical model based on infectious disease epidemiology and a priori known and unknown aspects of SRLV biology. The model determines the probability of a susceptible ewe seroconverting when it did, and of a susceptible ewe not seroconverting before it was removed from the flock or the experiment ended. The product of these probabilities gives the likelihood of the data given the model. The model was prototyped in Python and then written in C for speed.

The striking result of this research is that MV is a disease of housing. Even brief periods of housing allow the virus to spread rapidly, but transmission is negligible between sheep kept on pasture So, although individual sheep never recover from the disease, it could be eliminated from flocks over time by exploiting the fact that transmission of the virus is too slow between grazing sheep to sustain the disease.

Our second striking result suggests the disease is unlikely to be spread by newly-infected animals, contrary to general expectation. We estimated that the time between an animal being infected and becoming infectious is about a year. This delay, termed epidemiological latency, is actually longer than the estimated time delay between infection and seroconversion.

We can now begin to see more clearly how disease processes occurring in the infected individual shape what happens at the flock, or epidemiological, level. It seems that, after a sheep becomes infected, the disease slowly progresses to the point when there is sufficient free virus to be recognised by the immune system, but then further development of inflammatory lesions in the lungs has to take place before there are sufficient quantities of infective alveolar macrophages and free virus for transmission by the respiratory route. There follows a further delay, perhaps of some years, before the disease has advanced to the stage at which symptoms such as chronic pneumonia and hardening of the udder become apparent.

Infectiousness is expected to be a function of viral load, and although we do not know the timecourse of viral load, it seems most likely that it continues to increase throughout the development of chronic disease. This suggests to us that the infectiousness of an individual is not constant, but is likely to increase as the disease progresses and symptoms emerge.

We are interested in learning how infectiousness changes over the course of an individual’s infection because of the implications at the epidemiological level. Time delays in seroconversion merely make the disease more difficult to detect and control, but the epidemiological significance of a time delay in the development of infectiousness is that it acts to slow the spread of the virus. And if ewes with long-standing infections are the most infectious, they pose the greatest risk to uninfected sheep. This would present an opportunity for the management of exposure to slow the spread of disease. For example, if ewes in their last year of production are the most infectious, then young ewes should be kept away from them when housed – an idea supported by preliminary analysis using individual-based modelling (IBM – see Box 2). Separation of younger animals from older ones may reduce the prevalence of infection to the point where the costs of disease, in terms of lost production and poor welfare, are not overwhelming or at least are less than the cost of attempting to eliminate the disease – we discuss this later in this blog.

So far, there is only very limited and tentative evidence of increasing infectiousness in the literature, and direct experimental evidence would be very hard to come by. But it is plausible that disease severity, viral load and impaired productivity are all related to the extent of inflammatory lesions in the lungs. This suggests that measurably-impaired productivity in infected sheep could be used as a proxy for viral load, and hence infectiousness. And that brings us to our current project.

#### 2            Individual-based modelling

This is a technique to explore the consequences of probabilistic events, such as becoming infected by SRLV. The flock of ewes is modelled as individuals, and their progress through life is followed. Flock replacements are taken from the ewe lambs born to the flock; all other lambs being sold.

The figure shows the mean results (green line) of 1000 iterations of a stochastic simulation of SRLV prevalence in a flock of 400 ewes housed in groups of 100 for one month per year. The probability that an infected ewe will transmit the virus is modelled as rising exponentially with time since infection. The management regime we modelled was to segregate ewes during housing into each of their four age groups (2, 3, 4 and 5 years old) in separate pens, and to sell all the lambs of the oldest ewes, rather than retain any as flock replacements. From an initial starting prevalence of 275 infected ewes, the virus is virtually eliminated from the flock.

#### Using Genstat to model the cost of MV

Eliminating SRLV from an infected flock involves either repeated testing of the whole flock, culling reactors and perhaps also artificially rearing lambs removed at birth, or entirely replacing the flock with accredited disease-free stock. So, the cost of eliminating the virus from a flock can be huge. But what about the costs of living with it? These costs arise from poor welfare leading to lost production: lactation failure, reduced lamb growth and excess lamb and ewe mortality. But under what conditions are they so substantial as to warrant an elimination strategy? That depends, again, on answers at two levels: what are the production losses for each infected ewe, and how prevalent is the disease in the flock?

We have a number of reasons to want to quantify how the costs of being SRLV+ vary over the time-course of the disease. First, it is reasonable to assume that production losses will be related to the emergence of symptoms in affected sheep, but this has never been adequately quantified. Second, if production losses are a function of the duration of infection, and we can regard them as a proxy for viral load, then it would support the idea that infectiousness also increases as the disease progresses. And third, if production losses are only apparent in sheep with long-standing infections, which is therefore restricted to older sheep, then management could focus on early detection of symptoms and culling of older ewes.

We are quantifying these processes using a large dataset from colleagues at the Lublin University of Life Sciences. Their six-year study was designed to assess the response of production parameters to SRLV infection in a flock of breeding sheep kept under standard Polish husbandry conditions. They published results suggesting that infection with SRLV was associated with higher rates of age-specific ewe mortality, especially in older ewes.

The data comprise lambing records for about 800 ewes from three breeds, with over 300 ewes being present each year and a few being present every year. There are also records from about 2800 lambs born during the trial. Ewes were blood-tested in November and June each year, and all SRLV+ ewes were housed together following the November test until the lambs were weaned in April. SRLV- ewes were housed in the same shed, but segregated from the SRLV+ group. We were able to group the ewes on the basis of the series of blood test results as: (1) seroconverted before the trial began, (2) had not seroconverted by the end, and (3) seroconverted during the trial and for whom a time since seroconversion can be estimated to within about six months.

Given the nature of the data - unbalanced design, multiple observations from individual ewes and rams over several years, and different breeds – we used Genstat to fit mixed models to distinguish random and fixed effects. We were given access to Genstat release 22 beta, which adds greater functionality for displaying and saving output, producing predictions and visualising the fit of the model.

The example below addresses pre-weaning lamb mortality (mort, either 0 or 1). We are using a generalized linear mixed model where significant fixed terms were added stepwise. The ewes and rams used to produce these lambs are obvious random terms because they can be regarded as being drawn at random from a large population. There also appears to be a strong ewe.ram interaction, with some combinations faring differently from others. We included ‘year’ as a random term because, over the six years in which data were collected, factors such as flock size and forage quality varied somewhat randomly.

The fixed terms show that the probability of mortality is strongly affected by lamb birthweight (lambbirthwt). A quadratic term (lb2) models the well-known reduction in lamb survival in very large lambs - a consequence of birth difficulties. The age of the ewe, fitted as a factor (eweageF), is the next most significant fixed effect, followed by the SRLV status of the ewe tested in November prior to lambing (ewetestNov). The interaction term of ewe age with SRLV status is highly significant, showing that the way the ageing process in ewes affects the probability of their lambs’ mortality differs according to SRLV status. From the table of back-transformed means, we see that the probability of lamb mortality ranges between about 0.02 to 0.04 in SRLV- ewes aged from 2 to 5 years, perhaps declining in older ewes. SRLV+ ewes show similar lamb mortality in ages 2-4, but a progressive increase as ewes age further, almost doubling each year.

This preliminary analysis provides some evidence that the costs of being infected by SRLV are, indeed, progressive with age. There is some way to go yet to show whether sheep with longer-standing SRLV infection have higher viral loads and are more infectious, but our current research does point to a way to potential better disease control by targeting older animals. Maedi-visna isn’t called a progressive disease for anything, and we should be able to exploit that.

#### Example output from Genstat

We finally submitted our paper for publication in November 2019, just before the Covid 19 pandemic. One might have thought that a paper on the epidemiology and control of a respiratory virus spread by aerosol, especially indoors in close proximity and with a recommendation for social distancing, would have seemed quite topical. Ironically, early 2020 saw the subeditors and reviewers of such work all being urgently re-allocated to analysing information on the burgeoning pandemic. But we got there eventually and by October we were proudly writing a press release recommending social distancing … in sheep.

#### Reflecting on Genstat

Andrew Illius writes, “My experience of Genstat dates back to the early 1980s when I think Release 3 was current. It was hard going, and we queued up to have our fault codes diagnosed at the Genstat Clinic. But having learnt Fortran programming in the punched cards era, I was used to it taking several days to get a job to run. Genstat’s exacting requirements were reassuring and it became indispensable over the following years of agricultural and ecological research. By the 1990s we had learnt that mixed models were required to account separately for random and fixed effects in unbalanced data, and I’d been on a REML course. I was especially proud to use REML as my main analytical procedure thereafter because Robin Thompson invented it just down the corridor where we work in the Zoology building at Edinburgh University, and where he worked with the Animal Breeding group. It’s been a tremendous pleasure to get back to Genstat recently after many years away – like greeting an old friend. In the past, I wrote and submitted batch jobs on a UNIX mainframe before collecting some line-printer output on my way home. Now things have really speeded up, with the menu-driven environment of the Windows version. It’s a fantastic improvement, and a pleasure to use.”

Andrew Illius is Emeritus Prof of Animal Ecology in the Institute of Evolutionary Biology, University of Edinburgh, where he taught animal production and animal ecology from 1978 to 2008 and was latterly Head of the School of Biological Sciences. Most of his work has been on the ecology and management of grazing systems and the ecophysiology and behaviour of grazing animals. He retired in 2008 to spend more time with his sheep, keeping about 400 breeding ewes. Familiarity with sheep diseases led to collaboration with Nick Savill since 2018 on the epidemiology and control of MV.

Nick Savill is a Senior Lecturer at the Institute of Immunology and Infection Research, University of Edinburgh. He teaches a range of quantitative skills to undergraduate biological science students including maths, stats, data analysis and coding. His research interests are in mathematical modelling of infectious disease epidemiology. He has worked on foot and mouth disease, avian influenza, malaria, trypanosomiasis and, most recently, maedi-visna with Andrew Illius.

#### Reference

Illius AW, Lievaart-Peterson K, McNeilly TN, Savill NJ (2020) Epidemiology and control of maedi-visna virus: Curing the flock. PLoS ONE 15 (9): e0238781. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0238781