Molecular Ecology Resources (2010) 10, 393–396

doi: 10.1111/j.1755-0998.2009.02764.x

COMPUTER PROGRAM NOTE

SOFSOG: a suite of programs to avoid inbreeding in

plantation designs ´ NDEZ* and S . C . G O N Z A ´ L E Z - M A R T I´ N E Z † J. FERNA *Departamento de Mejora Gene´tica Animal, I. N. I. A. Ctra. Corun˜a Km 7,5, 28040 Madrid, Spain, †Departamento de Sistemas y Recursos Forestales, Centro de Investigacio´n Forestal, I. N. I. A. Ctra. Corun˜a Km 7,5, 28040 Madrid, Spain

Abstract Cost-effective ways of controlling inbreeding in conservation or productive plantations imply the allocation of individuals reducing the possibility of close relatives’ mating and, consequently, limiting inbreeding. SOFSOG is a suite of programs, which helps to design plantation sites. First, if the plantation scheme involves several plots, it allows distribution of individuals available among different sites minimizing within-site global coancestry. Then, it yields a plantation design for each site, either following the classical permutated neighbourhood strategy or the recently developed method by Ferna´ndez and Gonza´lez-Martı´nez. This new method allows the implementation of different pollen dispersion kernels, and to include in the designing strategy any available information on individual relationships, reproductive success, differences in phenology, etc., via weighting or penalization matrices. Additionally, the package includes a tool for calculating the molecular coancestry (Identity By State) from codominant marker data. Keywords: ex situ plantations, inbreeding, kinship, plant conservation, seed orchards Received 14 May 2009; revision received 7 July 2009; accepted 3 August 2009

Introduction Plant conservation programs often focus on species or populations with a low effective population size, and, normally, with severely threatened native environments. Avoiding mating among relatives in plant population management, in particular in plantations established as ex situ conservation units, is important to reduce inbreeding and its deleterious consequences (i.e. inbreeding depression) in subsequent generations. Usually, managers of open-pollinated plantations use, to avoid mating among relatives, physical isolation (e.g. implementing ‘barriers’ of other species) or, more commonly, plantation designs that separate like-individuals in space. Indeed, the larger the distance that separates two individuals, the less likely it is they mate (see reviews in Adams & Burczyk 2000; Smouse & Sork 2004; Sork & Smouse 2006). However, better designs for ex situ conservation plantations can be obtained by considering pairwise relatedness and specific dispersal kernels, and by the use of optimization algorithms.

Correspondence: Jesu´s Ferna´ndez Martı´n, Fax: 34 913478743; E-mail: [email protected]

 2009 Blackwell Publishing Ltd

Traditionally, plantations designed to avoid mating among relatives have used permutated neighbourhood and related methods (e.g. La Bastide 1967; Bell & Fletcher 1978; Chakravarty & Bagchi 1993; Nester 1994). However, these methods are limited in several ways. First, they only consider pollination among the closest individuals, not taking into account the possibility of long distance pollination, which is well documented (Oddou-Muratorio et al. 2005; Robledo-Arnuncio & Gil 2005; de Lucas et al. 2008). Second, they do not use any other available information, such as known relationships among individuals, pollen dispersal kernels or phenology. Finally, userfriendly implementations of permutated neighbourhood methods are not currently available. A recently developed method (Ferna´ndez & Gonza´lez-Martı´nez 2009), which can be of general use in conservation programs, provides a solution to all the limitations aforementioned: it considers all individuals in a plantation, it includes different pollen dispersal kernels, and it takes into account other sources of information too. To provide an example of the utility of this new method in conservation genetics, we will consider the case of the endemic variety Pinus sylvestris var. nevadensis, which is severely threatened within its native range by substantial

394 C O M P U T E R P R O G R A M N O T E levels of pollen introgression (Robledo-Arnuncio et al. 2009). Establishment of ex situ plantations is recommended for this species. Ideally, seeds produced in these plantations could be used to reinforce the natural populations if, through control crosses and ⁄ or management measures, inbreeding were avoided. Additionally, pollen dispersal kernels are available for this species (RobledoArnuncio & Gil 2005). A user-friendly plantation designing tool able to consider dispersal kernels to reduce offspring inbreeding, such as SOFSOG, would be a valuable tool for managers of this relict variety. SOFSOG is a suite of programs, which helps to design plantation sites with the ultimate goal of minimizing the probability of generating inbred offspring in open-pollination populations for conservation purposes. The programs can also be used to design seed orchards within breeding programs. The package includes three programs that can be helpful for managers facing the task of designing such plots. The executable files for all of them (compiled both for Windows and Macintosh), as well as some example files and a user manual, can be freely downloaded from http://www.uvigo.es/webs/c03/ webc03/XENETICA/XB2/Jesus/Fernandez.htm.

Description of the software is a suite of programs that helps to design conservation or breeding plantations. The package contains three programs: Coancestry, to calculate molecular coancestry between individuals; Divide, to distribute candidates among different available sites (if more than one); and Sofsog, to determine the plantation scheme within plots to avoid inbred offspring, following either the classical neighbourhood design or the new method by Ferna´ndez & Gonza´lez-Martı´nez (2009).

SOFSOG

Coancestry In natural plant populations, there is usually a lack of genealogies (at least of paternal lineages), and, therefore, it is not possible to construct a pedigree-based coancestry matrix. However, the recent development of DNA screening techniques has provided an increasingly large number of molecular markers that can be used to determine the genetic relationships between individuals. One of the simplest ways to quantify the relationship between individuals from markers is molecular coancestry (Toro et al. 2002), which measures the probability of Identity By State (or Alike In State) between pairs of individuals for one or several loci. This program calculates the molecular coancestry (Identity By State) between each pair of individuals (or genets) from the genotypes for codominant markers (for example microsatellites or single nucleotide

polymorphisms). Multilocus values are the average of single locus measures. Alleles should be coded for each locus separately and correlatively from 1 to the number of alleles; for example, if the locus has 10 different alleles these must be coded 1, 2, 3,…, 9, 10. Notwithstanding, the program can automatically recode from other kinds of numeric values (e.g. number of repeats, fragment length, etc.). Output file contains the upper triangle of the molecular coancestry matrix between all pairs of individuals, in the format the other programs of the suite will need as input of relationship information.

Divide Sometimes managers deal with the task of planting a number of individuals in more than a single plot. This is particularly common when dealing with shrubs and trees, as they occupy substantial amounts of space. Consequently, the first decision to be taken is which individuals are to be located in each site. The most sensible strategy is to distribute the individuals so that the mean global coancestry among these individuals selected for the same plot is minimized. By acting this way, the probability of generating inbred offspring is reduced even if the plantation scheme within plots is not optimized. A direct interpretation of this strategy would lead to the following objective function: n X xik xjk fij min ; Nk2 k¼1 where xik is a variable which takes the value 1 if the individual i is to be planted in the plot k and 0 otherwise, fij is the coancestry between individuals i and j (genealogical or marker based estimate), Nk is the number of plants to be allocated in plot k, and n is the number of considered plots. Notice that the minimization function considers also self-ancestry. There are three reasons for this. First, the function reflects the inbreeding arising from self-fertilization (in compatible plants). Second, self-coancestry is a measure of the inbreeding levels of the individuals themselves. And, finally, it also affects the amount of genetic diversity stored (even with no self-fertilization). This is important to determine the fate of the population beyond the first generation, especially when dealing with threatened populations. When sizes of the different plots are very different (and, therefore, also the number of individuals to be planted in each), the above strategy may lead to some sites with high global coancestry, even if the average for all the plots is low. To avoid this, the program actually minimizes the coancestry of the site with the highest global coancestry. Consequently, in the optimal configuration, the coancestry levels of every plot are quite similar.

 2009 Blackwell Publishing Ltd

C O M P U T E R P R O G R A M N O T E 395 An extra advantage of this procedure is that different sites become ‘genetic copies’ of the others (they have similar constitutions) protecting against the loss of genetic information if any particular plot became eventually extinct due to fire, disease breakthrough, etc.

Sofsog This is the main program performing the optimization of a single plantation site. The software assumes that individuals of a monoecious species (but other mating systems can also be simulated in the F&G method through the use of a weighting matrix, see below) are to be allocated in a plot with an orthogonal design comprising a certain number of rows and columns with the same distance between adjacent rows (columns). The method is also designed in a way that can take the particularities of the terrain into account, such as places (rocks, bad soil, etc.) not suitable for the plantation. This program has two sections. The first one performs the classical neighbourhood method, where the only information used is the genets each plant belongs to. This section is suitable, for instance, for plantations based on grafting. This tool allows the user to define the number of neighbourhoods to be accounted for (i.e. number of ‘rings’ where individuals from the same genet are avoided, or exclusion area). When it is impossible to fit completely the above restriction (i.e. some individuals from the same genet are to be placed in the ‘influence’ area of each other due to space constraints), the program searches for a solution that minimizes a function inversely proportional to the distance to the target plant. Consequently, solutions in which clone plants are located in the ‘corners’ of the exclusion area are preferred over others with like plants closer. The second section performs the allocation of individuals following the F&G method (Ferna´ndez & Gonza´lezMartı´nez 2009). In an open-pollinated plantation, the probability of a mother plant of generating inbred offspring depends on its coancestry with any other plant and the probability of being pollinated by each of the candidate fathers. Therefore, the optimum strategy to allocate N plants in a conservation or breeding plantation is to minimize the global probability across the whole population of generating inbred offspring, expressed by the function: PN N X y¼1 pðx; yÞf ðx; yÞ wðx; yÞ ; PN x¼1 y¼1 pðx; yÞ wðx; yÞ where p(x, y) is the probability for a plant x of being pollinated by another plant y, f(x, y) is the coancestry between these plants, and w(x, y) is a weight assigned to that particular mating. Note that probabilities have to be stan-

 2009 Blackwell Publishing Ltd

dardized for each plant to ensure that the sum of probabilities equals one. Both the probability associated to the couple (x, y) and the reciprocal combination (y, x) are considered as monoecious species are assumed. The selffertilization rate can be determined by the user depending on the problem species. For two different plants (i.e. two different locations), p(x, y) is inversely related to the distance between them [d(x, y)], because the larger the distance, the lower the probability of pollen reaching the target position. The program assumes orthogonal distribution of plants with rows i = 1,…, R and columns j = 1,…, C. The distance (in Euclidean terms) between plants located at positions ij and kl is calculated as dðij; klÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðm  ði  kÞÞ2 þðm  ðj  lÞÞ2 ;

where m is the distance between adjacent rows (columns). The particular relationship between p(x, y) and d(x, y) depends on the pollen dispersal probability function of the species under study. In the present software, the assumed pollen dispersal kernel is always power exponential, because this function seems to fit most plant species well (see Tufto et al. 1997; Austerlitz et al. 2004). It can also easily accommodate the normal (b = 2) and exponential (b = 1) kernels. The expression for calculating pollination probability in this case is   ! b dðx; yÞ b pðx; yÞ ¼ exp  ; 2pa2 Cð2=bÞ a where a and b are the scale and shape parameters of the distribution, respectively, and C is the gamma function (Tufto et al. 1997; Clark et al. 1998). Under this probability distribution, the mean pollination distance is d ¼ a½Cð3=bÞ=Cð2=bÞ: The parameters a and b can be defined by the user to conform to the particularities of the species he ⁄ she is dealing with and, thus, different functions can be modelled fitting almost every possible situation. Relationships between individuals are defined via the coancestry matrix, either obtained from genealogies or inferred from marker information. Only coancestries between genets are to be provided to the software and they will be applied to all plants obtained from those genets. The F&G method can also account for other factors affecting the probability of mating, such as differences in male fecundity (for instance, due to quantity and ⁄ or quality of pollen; Burczyk et al. 1996; Oddou-Muratorio et al. 2005) or phenological overlapping (RobledoArnuncio et al. 2004 and references therein) through the matrix of weights or penalizations [w(x, y)]. For example,

396 C O M P U T E R P R O G R A M N O T E the weights in the probability of mating may be of the P form: ki uij = k2N kk ukj ; where ki is the male fecundity of the ith individual, uij the phenological overlapping of the ith pollen donor with the jth mother plant and N is the total number of plants in the plantation. Nevertheless, any particular mating can be avoided (promoted) by including a high (low) weight in the appropriate position of the penalization matrix. All procedures implying an optimization (Divide and both methods in the Sofsog program) are solved using a simulated annealing method (Kirkpatrick et al. 1983) to search across the feasible space of solutions. The performance and accuracy of a simulated annealing algorithm depend on a series of parameters (initial ‘temperature’, number of steps, number of solutions per step, etc.), which can be modified by the user to adapt to the particular situation. See Kirkpatrick et al. (1983) for a thorough explanation of the annealing method and for considerations about the definition of the different parameters and Ferna´ndez & Gonza´lez-Martı´nez (2009) for deeper insights in the implementation to the plantation design. The software we present is a powerful tool for plantation design in the ex situ conservation and plant breeding context, in particular for forest trees and shrubs, replacing other applications advantageously. It can also facilitate the integration of ecological, life-history and pedigree information in ex situ conservation strategies, and encourage managers to adopt objective-oriented plantation designs.

References Adams WT, Burczyk J (2000) Magnitude and implications of gene flow in gene conservation reserves. In: Forest Conservation Genetics. Principles and Practice (eds Young A, Boshier D, Boyle T), pp. 215–224. CSIRO-CABI Publishing, Collingwood. Austerlitz F, Dick CW, Dutech C et al. (2004) Using genetic markers to estimate the pollen dispersal curve. Molecular Ecology, 13, 937–954. Bell G, Fletcher AM (1978) Computer organised orchard layouts (COOL) based on the permutated neighbourhood design concept. Silvae Genetica, 27, 223–225.

Burczyk J, Adams WT, Shimizu JY (1996) Mating patterns and pollen dispersal in a natural knobcone pine (Pinus attenuata Lemmon.) stand. Heredity, 77, 251–260. Chakravarty N, Bagchi SK (1993) A computer program for permutated neighbourhood seed orchard design. Silvae Genetica, 42, 1–5. Clark JS, Macklin E, Wood L (1998) Stages and spatial scales of recruitment limitation in southern Appalachian forests. Ecological Monographs, 68, 213–235. Ferna´ndez J, Gonza´lez-Martı´nez SC (2009) Allocating individuals to avoid inbreeding in ex situ conservation plantations: so far, so good. Conservation Genetics, 10, 45–57. Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optimization by simulated annealing. Science, 220, 671. La Bastide JGA (1967) A computer program for the layout of seed orchards. Euphytica, 16, 321–323. de Lucas AI, Robledo-Arnuncio JJ, Hidalgo E, Gonza´lez-Martı´nez SC (2008) Mating system and pollen gene flow in Mediterranean maritime pine. Heredity, 100, 390–399. Nester R (1994) Modulo tile construction for systematic seed orchard designs. Silvae Genetica, 43, 312–321. Oddou-Muratorio S, Klein EK, Austerlitz F (2005) Real-time patterns of pollen flow in the wildservice tree, Sorbus torminalis (L.) Crantz. II. Spatial patterns of pollen flow and variance in male reproductive success inferred from parent-offspring. Molecular Ecology, 14, 4441–4452. Robledo-Arnuncio JJ, Gil L (2005) Patterns of pollen dispersal in a small population of Pinus sylvestris L. revealed by totalexclusion paternity analysis. Heredity, 94, 13–22. Robledo-Arnuncio JJ, Alı´a R, Gil L (2004) Increased selfing and correlated paternity in a small population of a predominantly outcrossing conifer, Pinus sylvestris. Molecular Ecology, 13, 2567–2577. Robledo-Arnuncio JJ, Navascue´s M, Gonza´lez-Martı´nez SC, Gil L (2009) Estimating gametic introgression rates in a risk assessment context: a case study with Scots pine relicts. Heredity, DOI: 10.1038/hdy.2009.78 Smouse P, Sork VL (2004) Measuring pollen flow in forest trees: an exposition of alternative approaches. Forest Ecology and Management, 197, 21–38. Sork VL, Smouse PE (2006) Genetic analysis of landscape connectivity in tree populations. Landscape Ecology, 21, 821–836. ´ vilo C, Rodriga´n˜ez J, Rodrı´guez C, Silio´ Toro MA, Barraga´n C, O L (2002) Estimation of coancestry in Iberian pigs using molecular markers. Conservation Genetics, 3, 309–320. Tufto J, Engen S, Hindar K (1997) Stochastic dispersal processes in plant populations. Theoretical Population Biology, 52, 16–26.

 2009 Blackwell Publishing Ltd