Whole genome resequencing and phenotyping of MAGIC population for high resolution mapping of drought tolerance in chickpea

Terminal drought is one of the major constraints to crop production in chickpea (Cicer arietinum L.). In order to map drought tolerance related traits at high resolution, we sequenced multi-parent advanced generation intercross (MAGIC) population using whole genome resequencing approach and phenotyped it under drought stress environments for two consecutive years (2013-14 and 2014-15). A total of 52.02 billion clean reads containing 4.67 TB clean data were generated on the 1136 MAGIC lines and eight parental lines. Alignment of clean data on to the reference genome enabled identification of a total, 932,172 of SNPs, 35,973 insertions, and 35,726 deletions among the parental lines. A high-density genetic map was constructed using 57,180 SNPs spanning a map distance of 1606.69 cM. Using compressed mixed linear model, genome-wide association study (GWAS) enabled us to identify 737 markers significantly associated with days to 50% flowering, days to maturity, plant height, 100 seed weight, biomass, and harvest index. In addition to the GWAS approach, an identity-by-descent (IBD)-based mixed model approach was used to map quantitative trait loci (QTLs). The IBD-based mixed model approach detected major QTLs that were comparable to those from the GWAS analysis as well as some exclusive QTLs with smaller effects. The candidate genes like FRIGIDA and CaTIFY4b can be used for enhancing drought tolerance in chickpea. The genomic resources, genetic map, marker-trait associations, and QTLs identified in the study are valuable resources for the chickpea community for developing climate resilient chickpeas.

A conceptual framework for the dynamic modeling of time-resolved phenotypes for sets of genotype-environment-management combinations: a model library

Introduction

Dynamic crop growth models are an important tool to predict complex traits, like crop yield, for modern and future genotypes in their current and evolving environments, as those occurring under climate change. Phenotypic traits are the result of interactions between genetic, environmental, and management factors, and dynamic models are designed to generate the interactions producing phenotypic changes over the growing season. Crop phenotype data are becoming increasingly available at various levels of granularity, both spatially (landscape) and temporally (longitudinal, time-series) from proximal and remote sensing technologies.

Methods

Here we propose four phenomenological process models of limited complexity based on differential equations for a coarse description of focal crop traits and environmental conditions during the growing season. Each of these models defines interactions between environmental drivers and crop growth (logistic growth, with implicit growth restriction, or explicit restriction by irradiance, temperature, or water availability) as a minimal set of constraints without resorting to strongly mechanistic interpretations of the parameters. Differences between individual genotypes are conceptualized as differences in crop growth parameter values.

Results

We demonstrate the utility of such low-complexity models with few parameters by fitting them to longitudinal datasets from the simulation platform APSIM-Wheat involving in silico biomass development of 199 genotypes and data of environmental variables over the course of the growing season at four Australian locations over 31 years. While each of the four models fits well to particular combinations of genotype and trial, none of them provides the best fit across the full set of genotypes by trials because different environmental drivers will limit crop growth in different trials and genotypes in any specific trial will not necessarily experience the same environmental limitation.

Discussion

A combination of low-complexity phenomenological models covering a small set of major limiting environmental factors may be a useful forecasting tool for crop growth under genotypic and environmental variation.

Genetic Mapping of Genotype-by-Ploidy Effects in Arabidopsis thaliana

Plants can express different phenotypic responses following polyploidization, but ploidy-dependent phenotypic variation has so far not been assigned to specific genetic factors. To map such effects, segregating populations at different ploidy levels are required. The availability of an efficient haploid inducer line in Arabidopsis thaliana allows for the rapid development of large populations of segregating haploid offspring. Because Arabidopsis haploids can be self-fertilised to give rise to homozygous doubled haploids, the same genotypes can be phenotyped at both the haploid and diploid ploidy level. Here, we compared the phenotypes of recombinant haploid and diploid offspring derived from a cross between two late flowering accessions to map genotype × ploidy (G × P) interactions. Ploidy-specific quantitative trait loci (QTLs) were detected at both ploidy levels. This implies that mapping power will increase when phenotypic measurements of monoploids are included in QTL analyses. A multi-trait analysis further revealed pleiotropic effects for a number of the ploidy-specific QTLs as well as opposite effects at different ploidy levels for general QTLs. Taken together, we provide evidence of genetic variation between different Arabidopsis accessions being causal for dissimilarities in phenotypic responses to altered ploidy levels, revealing a G × P effect. Additionally, by investigating a population derived from late flowering accessions, we revealed a major vernalisation-specific QTL for variation in flowering time, countering the historical bias of research in early flowering accessions.

statgenMPP: an R package implementing an IBD-based mixed model approach for QTL mapping in a wide range of multi-parent populations

MOTIVATION

Multi-parent populations (MPPs) are popular for QTL mapping because they combine wide genetic diversity in parents with easy control of population structure, but a limited number of software tools for QTL mapping are specifically developed for general MPP designs.

RESULTS

We developed an R package called statgenMPP, adopting a unified identity-by-descent (IBD)-based mixed model approach for QTL analysis in MPPs. The package offers easy-to-use functionalities of IBD calculations, mixed model solutions and visualizations for QTL mapping in a wide range of MPP designs, including diallele, nested-association mapping populations, multi-parent advanced genetic inter-cross populations and other complicated MPPs with known crossing schemes.

AVAILABILITY AND IMPLEMENTATION

The R package statgenMPP is open-source and freely available on CRAN at https://CRAN.R-project.org/package=statgenMPP.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Identification of environment types and adaptation zones with self-organizing maps; applications to sunflower multi-environment data in Europe

We evaluate self-organizing maps (SOM) to identify adaptation zones and visualize multi-environment genotypic responses. We apply SOM to multiple traits and crop growth model output of large-scale European sunflower data. Genotype-by-environment interactions (G × E) complicate the selection of well-adapted varieties. A possible solution is to group trial locations into adaptation zones with G × E occurring mainly between zones. By selecting for good performance inside those zones, response to selection is increased. In this paper, we present a two-step procedure to identify adaptation zones that starts from a self-organizing map (SOM). In the SOM, trials across locations and years are assigned to groups, called units, that are organized on a two-dimensional grid. Units that are further apart contain more distinct trials. In an iterative process of reweighting trial contributions to units, the grid configuration is learnt simultaneously with the trial assignment to units. An aggregation of the units in the SOM by hierarchical clustering then produces environment types, i.e. trials with similar growing conditions. Adaptation zones can subsequently be identified by grouping trial locations with similar distributions of environment types across years. For the construction of SOMs, multiple data types can be combined. We compared environment types and adaptation zones obtained for European sunflower from quantitative traits like yield, oil content, phenology and disease scores with those obtained from environmental indices calculated with the crop growth model Sunflo. We also show how results are affected by input data organization and user-defined weights for genotypes and traits. Adaptation zones for European sunflower as identified by our SOM-based strategy captured substantial genotype-by-location interaction and pointed to trials in Spain, Turkey and South Bulgaria as inducing different genotypic responses.

Two-dimensional P-spline smoothing for spatial analysis of plant breeding trials

Large agricultural field trials may display irregular spatial trends that cannot be fully captured by a purely randomization-based analysis. For this reason, paralleling the development of analysis-of-variance procedures for randomized field trials, there is a long history of spatial modeling for field trials, starting with the early work of Papadakis on nearest neighbor analysis, which can be cast in terms of first or second differences among neighboring plot values. This kind of spatial modeling is amenable to a natural extension using splines, as has been demonstrated in recent publications in the field. Here, we consider the P-spline framework, focusing on model options that are easy to implement in linear mixed model packages. Two examples serve to illustrate and evaluate the methods. A key conclusion is that first differences are rather competitive with second differences. A further key observation is that second differences require special attention regarding the representation of the null space of the smooth terms for spatial interaction, and that an unstructured variance-covariance structure is required to ensure invariance to translation and rotation of eigenvectors associated with that null space. We develop a strategy that permits fitting this model with ease, but the approach is more demanding than that needed for fitting models using first differences. Hence, even though in other areas, second differences are very commonly used in the application of P-splines, our conclusion is that with field trials, first differences have advantages for routine use.

A two-stage approach for the spatio-temporal analysis of high-throughput phenotyping data

High throughput phenotyping (HTP) platforms and devices are increasingly used for the characterization of growth and developmental processes for large sets of plant genotypes. Such HTP data require challenging statistical analyses in which longitudinal genetic signals need to be estimated against a background of spatio-temporal noise processes. We propose a two-stage approach for the analysis of such longitudinal HTP data. In a first stage, we correct for design features and spatial trends per time point. In a second stage, we focus on the longitudinal modelling of the spatially corrected data, thereby taking advantage of shared longitudinal features between genotypes and plants within genotypes. We propose a flexible hierarchical three-level P-spline growth curve model, with plants/plots nested in genotypes, and genotypes nested in populations. For selection of genotypes in a plant breeding context, we show how to extract new phenotypes, like growth rates, from the estimated genotypic growth curves and their first-order derivatives. We illustrate our approach on HTP data from the PhenoArch greenhouse platform at INRAE Montpellier and the outdoor Field Phenotyping platform at ETH Zürich.

An IBD-based mixed model approach for QTL mapping in multiparental populations

The identity-by-descent (IBD)-based mixed model approach introduced in this study can detect quantitative trait loci (QTLs) referring to the parental origin and simultaneously account for multilevel relatedness of individuals within and across families. This unified approach is proved to be a powerful approach for all kinds of multiparental population (MPP) designs. Multiparental populations (MPPs) have become popular for quantitative trait loci (QTL) detection. Tools for QTL mapping in MPPs are mostly developed for specific MPPs and do not generalize well to other MPPs. We present an IBD-based mixed model approach for QTL mapping in all kinds of MPP designs, e.g., diallel, Nested Association Mapping (NAM), and Multiparental Advanced Generation Intercross (MAGIC) designs. The first step is to compute identity-by-descent (IBD) probabilities using a general Hidden Markov model framework, called reconstructing ancestry blocks bit by bit (RABBIT). Next, functions of IBD information are used as design matrices, or genetic predictors, in a mixed model approach to estimate variance components for multiallelic genetic effects associated with parents. Family-specific residual genetic effects are added, and a polygenic effect is structured by kinship relations between individuals. Case studies of simulated diallel, NAM, and MAGIC designs proved that the advanced IBD-based multi-QTL mixed model approach incorporating both kinship relations and family-specific residual variances (IBD.MQMkin_F) is robust across a variety of MPP designs and allele segregation patterns in comparison to a widely used benchmark association mapping method, and in most cases, outperformed or behaved at least as well as other tools developed for specific MPP designs in terms of mapping power and resolution. Successful analyses of real data cases confirmed the wide applicability of our IBD-based mixed model methodology.

Temperature response of wheat affects final height and the timing of stem elongation under field conditions

In wheat, temperature affects the timing and intensity of stem elongation. Genetic variation for this process is therefore important for adaptation. This study investigates the genetic response to temperature fluctuations during stem elongation and its relationship to phenology and height. Canopy height of 315 wheat genotypes (GABI wheat panel) was scanned twice weekly in the field phenotyping platform (FIP) of ETH Zurich using a LIDAR. Temperature response was modelled using linear regressions between stem elongation and mean temperature in each measurement interval. This led to a temperature-responsive (slope) and a temperature-irresponsive (intercept) component. The temperature response was highly heritable (H2=0.81) and positively related to a later start and end of stem elongation as well as final height. Genome-wide association mapping revealed three temperature-responsive and four temperature-irresponsive quantitative trait loci (QTLs). Furthermore, putative candidate genes for temperature-responsive QTLs were frequently related to the flowering pathway in Arabidopsis thaliana, whereas temperature-irresponsive QTLs corresponded to growth and reduced height genes. In combination with Rht and Ppd alleles, these loci, together with the loci for the timing of stem elongation, accounted for 71% of the variability in height. This demonstrates how high-throughput field phenotyping combined with environmental covariates can contribute to a smarter selection of climate-resilient crops.

Temperature response of wheat affects final height and the timing of stem elongation under field conditions

In wheat, temperature affects the timing and intensity of stem elongation. Genetic variation for this process is therefore important for adaptation. This study investigates the genetic response to temperature fluctuations during stem elongation and its relationship to phenology and height. Canopy height of 315 wheat genotypes (GABI wheat panel) was scanned twice weekly in the field phenotyping platform (FIP) of ETH Zurich using a LIDAR. Temperature response was modelled using linear regressions between stem elongation and mean temperature in each measurement interval. This led to a temperature-responsive (slope) and a temperature-irresponsive (intercept) component. The temperature response was highly heritable (H2=0.81) and positively related to a later start and end of stem elongation as well as final height. Genome-wide association mapping revealed three temperature-responsive and four temperature-irresponsive quantitative trait loci (QTLs). Furthermore, putative candidate genes for temperature-responsive QTLs were frequently related to the flowering pathway in Arabidopsis thaliana, whereas temperature-irresponsive QTLs corresponded to growth and reduced height genes. In combination with Rht and Ppd alleles, these loci, together with the loci for the timing of stem elongation, accounted for 71% of the variability in height. This demonstrates how high-throughput field phenotyping combined with environmental covariates can contribute to a smarter selection of climate-resilient crops.

Assessment of Multi-Image Unmanned Aerial Vehicle Based High-Throughput Field Phenotyping of Canopy Temperature

Canopy temperature (CT) has been related to water-use and yield formation in crops. However, constantly (e.g., sun illumination angle, ambient temperature) as well as rapidly (e.g., clouds) changing environmental conditions make it difficult to compare measurements taken even at short time intervals. This poses a great challenge for high-throughput field phenotyping (HTFP). The aim of this study was to i) set up a workflow for unmanned aerial vehicles (UAV) based HTFP of CT, ii) investigate different data processing procedures to combine information from multiple images into orthomosaics, iii) investigate the repeatability of the resulting CT by means of heritability, and iv) investigate the optimal timing for thermography measurements. Additionally, the approach was v) compared with other methods for HTFP of CT. The study was carried out in a winter wheat field trial with 354 genotypes planted in two replications in a temperate climate, where a UAV captured CT in a time series of 24 flights during 6 weeks of the grain-filling phase. Custom-made thermal ground control points enabled accurate georeferencing of the data. The generated thermal orthomosaics had a high spatial accuracy (mean ground sampling distance of 5.03 cm/pixel) and position accuracy [mean root-mean-square deviation (RMSE) = 4.79 cm] over all time points. An analysis on the impact of the measurement geometry revealed a gradient of apparent CT in parallel to the principle plane of the sun and a hotspot around nadir. Averaging information from all available images (and all measurement geometries) for an area of interest provided the best results by means of heritability. Correcting for spatial in-field heterogeneity as well as slight environmental changes during the measurements were performed with the R package SpATS. CT heritability ranged from 0.36 to 0.74. Highest heritability values were found in the early afternoon. Since senescence was found to influence the results, it is recommended to measure CT in wheat after flowering and before the onset of senescence. Overall, low-altitude and high-resolution remote sensing proved suitable to assess the CT of crop genotypes in a large number of small field plots as is required in crop breeding and variety testing experiments.

From QTLs to Adaptation Landscapes: Using Genotype-To-Phenotype Models to Characterize G×E Over Time

Genotype by environment interaction (G×E) for the target trait, e.g. yield, is an emerging property of agricultural systems and results from the interplay between a hierarchy of secondary traits involving the capture and allocation of environmental resources during the growing season. This hierarchy of secondary traits ranges from basic traits that correspond to response mechanisms/sensitivities, to intermediate traits that integrate a larger number of processes over time and therefore show a larger amount of G×E. Traits underlying yield differ in their contribution to adaptation across environmental conditions and have different levels of G×E. Here, we provide a framework to study the performance of genotype to phenotype (G2P) modeling approaches. We generate and analyze response surfaces, or adaptation landscapes, for yield and yield related traits, emphasizing the organization of the traits in a hierarchy and their development and interactions over time. We use the crop growth model APSIM-wheat with genotype-dependent parameters as a tool to simulate non-linear trait responses over time with complex trait dependencies and apply it to wheat crops in Australia. For biological realism, APSIM parameters were given a genetic basis of 300 QTLs sampled from a gamma distribution whose shape and rate parameters were estimated from real wheat data. In the simulations, the hierarchical organization of the traits and their interactions over time cause G×E for yield even when underlying traits do not show G×E. Insight into how G×E arises during growth and development helps to improve the accuracy of phenotype predictions within and across environments and to optimize trial networks. We produced a tangible simulated adaptation landscape for yield that we first investigated for its biological credibility by statistical models for G×E that incorporate genotypic and environmental covariables. Subsequently, the simulated trait data were used to evaluate statistical genotype-to-phenotype models for multiple traits and environments and to characterize relationships between traits over time and across environments, as a way to identify traits that could be useful to select for specific adaptation. Designed appropriately, these types of simulated landscapes might also serve as a basis to train other, more deep learning methodologies in order to transfer such network models to real-world situations.

Combining Crop Growth Modeling and Statistical Genetic Modeling to Evaluate Phenotyping Strategies

Genomic prediction of complex traits, say yield, benefits from including information on correlated component traits. Statistical criteria to decide which yield components to consider in the prediction model include the heritability of the component traits and their genetic correlation with yield. Not all component traits are easy to measure. Therefore, it may be attractive to include proxies to yield components, where these proxies are measured in (high-throughput) phenotyping platforms during the growing season. Using the Agricultural Production Systems Simulator (APSIM)-wheat cropping systems model, we simulated phenotypes for a wheat diversity panel segregating for a set of physiological parameters regulating phenology, biomass partitioning, and the ability to capture environmental resources. The distribution of the additive quantitative trait locus effects regulating the APSIM physiological parameters approximated the same distribution of quantitative trait locus effects on real phenotypic data for yield and heading date. We use the crop growth model APSIM-wheat to simulate phenotypes in three Australian environments with contrasting water deficit patterns. The APSIM output contained the dynamics of biomass and canopy cover, plus yield at the end of the growing season. Each water deficit pattern triggered different adaptive mechanisms and the impact of component traits differed between drought scenarios. We evaluated multiple phenotyping schedules by adding plot and measurement error to the dynamics of biomass and canopy cover. We used these trait dynamics to fit parametric models and P-splines to extract parameters with a larger heritability than the phenotypes at individual time points. We used those parameters in multi-trait prediction models for final yield. The combined use of crop growth models and multi-trait genomic prediction models provides a procedure to assess the efficiency of phenotyping strategies and compare methods to model trait dynamics. It also allows us to quantify the impact of yield components on yield prediction accuracy even in different environment types. In scenarios with mild or no water stress, yield prediction accuracy benefitted from including biomass and green canopy cover parameters. The advantage of the multi-trait model was smaller for the early-drought scenario, due to the reduced correlation between the secondary and the target trait. Therefore, multi-trait genomic prediction models for yield require scenario-specific correlated traits.

Construction of Genetic Linkage Maps in Multiparental Populations

Construction of genetic linkage maps has become a routine step for mapping quantitative trait loci (QTL), particularly in animal and plant breeding populations. Many multiparental populations have recently been produced to increase genetic diversity and QTL mapping resolution. However, few software packages are available for map construction in these populations. In this paper, we build a general framework for the construction of genetic linkage maps from genotypic data in diploid populations, including bi- and multiparental populations, cross-pollinated (CP) populations, and breeding pedigrees. The framework is implemented as an automatic pipeline called magicMap, where the maximum multilocus likelihood approach utilizes genotypic information efficiently. We evaluate magicMap by extensive simulations and eight real datasets: one biparental, one CP, four multiparent advanced generation intercross (MAGIC), and two nested association mapping (NAM) populations, the number of markers ranging from a few hundred to tens of thousands. Not only is magicMap the only software capable of accommodating all of these designs, it is more accurate and robust to missing genotypes and genotyping errors than commonly used packages.

Modelling strategies for assessing and increasing the effectiveness of new phenotyping techniques in plant breeding

New types of phenotyping tools generate large amounts of data on many aspects of plant physiology and morphology with high spatial and temporal resolution. These new phenotyping data are potentially useful to improve understanding and prediction of complex traits, like yield, that are characterized by strong environmental context dependencies, i.e., genotype by environment interactions. For an evaluation of the utility of new phenotyping information, we will look at how this information can be incorporated in different classes of genotype-to-phenotype (G2P) models. G2P models predict phenotypic traits as functions of genotypic and environmental inputs. In the last decade, access to high-density single nucleotide polymorphism markers (SNPs) and sequence information has boosted the development of a class of G2P models called genomic prediction models that predict phenotypes from genome wide marker profiles. The challenge now is to build G2P models that incorporate simultaneously extensive genomic information alongside with new phenotypic information. Beyond the modification of existing G2P models, new G2P paradigms are required. We present candidate G2P models for the integration of genomic and new phenotyping information and illustrate their use in examples. Special attention will be given to the modelling of genotype by environment interactions. The G2P models provide a framework for model based phenotyping and the evaluation of the utility of phenotyping information in the context of breeding programs.

Accurate Genotype Imputation in Multiparental Populations from Low-Coverage Sequence

Many different types of multiparental populations have recently been produced to increase genetic diversity and resolution in QTL mapping. Low-coverage, genotyping-by-sequencing (GBS) technology has become a cost-effective tool in these populations, despite large amounts of missing data in offspring and founders. In this work, we present a general statistical framework for genotype imputation in such experimental crosses from low-coverage GBS data. Generalizing a previously developed hidden Markov model for calculating ancestral origins of offspring DNA, we present an imputation algorithm that does not require parental data and that is applicable to bi- and multiparental populations. Our imputation algorithm allows heterozygosity of parents and offspring as well as error correction in observed genotypes. Further, our approach can combine imputation and genotype calling from sequencing reads, and it also applies to called genotypes from SNP array data. We evaluate our imputation algorithm by simulated and real data sets in four different types of populations: the F2, the advanced intercross recombinant inbred lines, the multiparent advanced generation intercross, and the cross-pollinated population. Because our approach uses marker data and population design information efficiently, the comparisons with previous approaches show that our imputation is accurate at even very low ([Formula: see text]) sequencing depth, in addition to having accurate genotype phasing and error detection.

Modelling spatial trends in sorghum breeding field trials using a two-dimensional P-spline mixed model

A flexible and user-friendly spatial method called SpATS performed comparably to more elaborate and trial-specific spatial models in a series of sorghum breeding trials. Adjustment for spatial trends in plant breeding field trials is essential for efficient evaluation and selection of genotypes. Current mixed model methods of spatial analysis are based on a multi-step modelling process where global and local trends are fitted after trying several candidate spatial models. This paper reports the application of a novel spatial method that accounts for all types of continuous field variation in a single modelling step by fitting a smooth surface. The method uses two-dimensional P-splines with anisotropic smoothing formulated in the mixed model framework, referred to as SpATS model. We applied this methodology to a series of large and partially replicated sorghum breeding trials. The new model was assessed in comparison with the more elaborate standard spatial models that use autoregressive correlation of residuals. The improvements in precision and the predictions of genotypic values produced by the SpATS model were equivalent to those obtained using the best fitting standard spatial models for each trial. One advantage of the approach with SpATS is that all patterns of spatial trend and genetic effects were modelled simultaneously by fitting a single model. Furthermore, we used a flexible model to adequately adjust for field trends. This strategy reduces potential parameter identification problems and simplifies the model selection process. Therefore, the new method should be considered as an efficient and easy-to-use alternative for routine analyses of plant breeding trials.

Fine mapping of a major QTL for awn length in barley using a multiparent mapping population

Awn length was mapped using a multiparent population derived from cv. Morex and four wild accessions. One QTL was fine mapped and candidate genes were identified in NILs by RNA-seq. Barley awns are photosynthetically active and contribute to grain yield. Awn length is variable among both wild and cultivated barley genotypes and many mutants with alterations in awn length have been identified. Here, we used a multiparent mapping population derived from cv. Morex and four genetically diverse wild barley lines to detect quantitative trait loci (QTLs) for awn length. Twelve QTLs, distributed over the barley genome, were identified with the most significant one located on chromosome arm 7HL (QTL AL7.1). The effect of AL7.1 was confirmed using near isogenic lines (NILs) and fine-mapped in two independent heterogeneous inbred families to a < 0.9 cM interval. With exception of a small effect on grain width, no other traits such as plant height or flowering time were affected by AL7.1. Variant calling on transcripts obtained from RNA sequencing reads in NILs was used to narrow down the list of candidate genes located in the interval. This data may be used for further characterization and unravelling of the mechanisms underlying natural variation in awn length.

First step in using molecular data for microbial food safety risk assessment; hazard identification of Escherichia coli O157:H7 by coupling genomic data with in vitro adherence to human epithelial cells

The potential for using whole genome sequencing (WGS) data in microbiological risk assessment (MRA) has been discussed on several occasions since the beginning of this century. Still, the proposed heuristic approaches have never been applied in a practical framework. This is due to the non-trivial problem of mapping microbial information consisting of thousands of loci onto a probabilistic scale for risks. The paradigm change for MRA involves translation of multidimensional microbial genotypic information to much reduced (integrated) phenotypic information and onwards to a single measure of human risk (i.e. probability of illness).

In this paper a first approach in methodology development is described for the application of WGS data in MRA; this is supported by a practical example. That is, combining genetic data (single nucleotide polymorphisms; SNPs) for Shiga toxin-producing Escherichia coli (STEC) O157 with phenotypic data (in vitro adherence to epithelial cells as a proxy for virulence) leads to hazard identification in a Genome Wide Association Study (GWAS).

This application revealed practical implications when using SNP data for MRA. These can be summarized by considering the following main issues: optimum sample size for valid inference on population level, correction for population structure, quantification and calibration of results, reproducibility of the analysis, links with epidemiological data, anchoring and integration of results into a systems biology approach for the translation of molecular studies to human health risk.

Future developments in genetic data analysis for MRA should aim at resolving the mapping problem of processing genetic sequences to come to a quantitative description of risk. The development of a clustering scheme focusing on biologically relevant information of the microbe involved would be a useful approach in molecular data reduction for risk assessment.

Genotype-environment interactions affecting preflowering physiological and morphological traits of Brassica rapa grown in two watering regimes

Plant growth and productivity are greatly affected by drought, which is likely to become more threatening with the predicted global temperature increase. Understanding the genetic architecture of complex quantitative traits and their interaction with water availability may lead to improved crop adaptation to a wide range of environments. Here, the genetic basis of 20 physiological and morphological traits is explored by describing plant performance and growth in a Brassica rapa recombinant inbred line (RIL) population grown on a sandy substrate supplemented with nutrient solution, under control and drought conditions. Altogether, 54 quantitative trait loci (QTL) were identified, of which many colocated in 11 QTL clusters. Seventeen QTL showed significant QTL-environment interaction (Q×E), indicating genetic variation for phenotypic plasticity. Of the measured traits, only hypocotyl length did not show significant genotype-environment interaction (G×E) in both environments in all experiments. Correlation analysis showed that, in the control environment, stomatal conductance was positively correlated with total leaf dry weight (DW) and aboveground DW, whereas in the drought environment, stomatal conductance showed a significant negative correlation with total leaf DW and aboveground DW. This correlation was explained by antagonistic fitness effects in the drought environment, controlled by a QTL cluster on chromosome A7. These results demonstrate that Q×E is an important component of the genetic variance and can play a great role in improving drought tolerance in future breeding programmes.

Studying the genetic basis of drought tolerance in sorghum by managed stress trials and adjustments for phenological and plant height differences

Managed environments in the form of well watered and water stressed trials were performed to study the genetic basis of grain yield and stay green in sorghum with the objective of validating previously detected QTL. As variations in phenology and plant height may influence QTL detection for the target traits, QTL for flowering time and plant height were introduced as cofactors in QTL analyses for yield and stay green. All but one of the flowering time QTL were detected near yield and stay green QTL. Similar co-localization was observed for two plant height QTL. QTL analysis for yield, using flowering time/plant height cofactors, led to yield QTL on chromosomes 2, 3, 6, 8 and 10. For stay green, QTL on chromosomes 3, 4, 8 and 10 were not related to differences in flowering time/plant height. The physical positions for markers in QTL regions projected on the sorghum genome suggest that the previously detected plant height QTL, Sb-HT9-1, and Dw2, in addition to the maturity gene, Ma5, had a major confounding impact on the expression of yield and stay green QTL. Co-localization between an apparently novel stay green QTL and a yield QTL on chromosome 3 suggests there is potential for indirect selection based on stay green to improve drought tolerance in sorghum. Our QTL study was carried out with a moderately sized population and spanned a limited geographic range, but still the results strongly emphasize the necessity of corrections for phenology in QTL mapping for drought tolerance traits in sorghum.

QTL linkage analysis of connected populations using ancestral marker and pedigree information

The common assumption in quantitative trait locus (QTL) linkage mapping studies that parents of multiple connected populations are unrelated is unrealistic for many plant breeding programs. We remove this assumption and propose a Bayesian approach that clusters the alleles of the parents of the current mapping populations from locus-specific identity by descent (IBD) matrices that capture ancestral marker and pedigree information. Moreover, we demonstrate how the parental IBD data can be incorporated into a QTL linkage analysis framework by using two approaches: a Threshold IBD model (TIBD) and a Latent Ancestral Allele Model (LAAM). The TIBD and LAAM models are empirically tested via numerical simulation based on the structure of a commercial maize breeding program. The simulations included a pilot dataset with closely linked QTL on a single linkage group and 100 replicated datasets with five linkage groups harboring four unlinked QTL. The simulation results show that including parental IBD data (similarly for TIBD and LAAM) significantly improves the power and particularly accuracy of QTL mapping, e.g., position, effect size and individuals' genotype probability without significantly increasing computational demand.

Gene and QTL detection in a three-way barley cross under selection by a mixed model with kinship information using SNPs

Quantitative trait locus (QTL) detection is commonly performed by analysis of designed segregating populations derived from two inbred parental lines, where absence of selection, mutation and genetic drift is assumed. Even for designed populations, selection cannot always be avoided, with as consequence varying correlation between genotypes instead of uniform correlation. Akin to linkage disequilibrium mapping, ignoring this type of genetic relatedness will increase the rate of false-positives. In this paper, we advocate using mixed models including genetic relatedness, or 'kinship' information for QTL detection in populations where selection forces operated. We demonstrate our case with a three-way barley cross, designed to segregate for dwarfing, vernalization and spike morphology genes, in which selection occurred. The population of 161 inbred lines was screened with 1,536 single nucleotide polymorphisms (SNPs), and used for gene and QTL detection. The coefficient of coancestry matrix was estimated based on the SNPs and imposed to structure the distribution of random genotypic effects. The model incorporating kinship, coancestry, information was consistently superior to the one without kinship (according to the Akaike information criterion). We show, for three traits, that ignoring the coancestry information results in an unrealistically high number of marker-trait associations, without providing clear conclusions about QTL locations. We used a number of widely recognized dwarfing and vernalization genes known to segregate in the studied population as landmarks or references to assess the agreement of the mapping results with a priori candidate gene expectations. Additional QTLs to the major genes were detected for all traits as well.

A mixed-model quantitative trait loci (QTL) analysis for multiple-environment trial data using environmental covariables for QTL-by-environment interactions, with an example in maize

Complex quantitative traits of plants as measured on collections of genotypes across multiple environments are the outcome of processes that depend in intricate ways on genotype and environment simultaneously. For a better understanding of the genetic architecture of such traits as observed across environments, genotype-by-environment interaction should be modeled with statistical models that use explicit information on genotypes and environments. The modeling approach we propose explains genotype-by-environment interaction by differential quantitative trait locus (QTL) expression in relation to environmental variables. We analyzed grain yield and grain moisture for an experimental data set composed of 976 F(5) maize testcross progenies evaluated across 12 environments in the U.S. corn belt during 1994 and 1995. The strategy we used was based on mixed models and started with a phenotypic analysis of multi-environment data, modeling genotype-by-environment interactions and associated genetic correlations between environments, while taking into account intraenvironmental error structures. The phenotypic mixed models were then extended to QTL models via the incorporation of marker information as genotypic covariables. A majority of the detected QTL showed significant QTL-by-environment interactions (QEI). The QEI were further analyzed by including environmental covariates into the mixed model. Most QEI could be understood as differential QTL expression conditional on longitude or year, both consequences of temperature differences during critical stages of the growth.

Heteroclinic orbits indicate overexploitation in predator–prey systems with a strong Allee effect

Species establishment in a model system in a homogeneous environment can be dependent not only on the parameter setting, but also on the initial conditions of the system. For instance, predator invasion into an established prey population can fail and lead to system collapse, an event referred to as overexploitation. This phenomenon occurs in models with bistability properties, such as strong Allee effects. The Allee effect then prevents easy re-establishment of the prey species. In this paper, we deal with the bifurcation analyses of two previously published predator-prey models with strong Allee effects. We expand the analyses to include not only local, but also global bifurcations. We show the existence of a point-to-point heteroclinic cycle in these models, and discuss numerical techniques for continuation in parameter space. The continuation of such a cycle in two-parameter space forms the boundary of a region in parameter space where the system collapses after predator invasion, i.e. where overexploitation occurs. We argue that the detection and continuation of global bifurcations in these models are of vital importance for the understanding of the model dynamics.

Extending Xu's Bayesian model for estimating polygenic effects using markers of the entire genome

In 2003, Xu obtained remarkably precise estimates of QTL positions despite the many markers simultaneously in his Bayesian model. We extend his model and Gibbs algorithm to ensure a valid posterior distribution and convergence to it, without changing the attractiveness of the method.

A penalized likelihood method for mapping epistatic quantitative trait Loci with one-dimensional genome searches

Epistasis is a common and important phenomenon, as indicated by results from a number of recent experiments. Unfortunately, the discovery of epistatic quantitative trait loci (QTL) is difficult since one must search for multiple QTL simultaneously in two or more dimensions. Such a multidimensional search necessitates many statistical tests, and a high statistical threshold must be adopted to avoid false positives. Furthermore, the large number of (interaction) parameters in comparison with the number of observations results in a serious danger of overfitting and overinterpretation of the data. In this article we present a new statistical framework for mapping epistasis in inbred line crosses. It is based on reducing the high dimensionality of the problem in two ways. First, epistatic QTL are mapped in a one-dimensional genome scan for high interactions between QTL and the genetic background. Second, the dimension of the search is bounded by penalized likelihood methods. We use simulated backcross data to illustrate the new approach.

Numerical bifurcation analysis of a tri-trophic food web with omnivory

We study the consequences of omnivory on the dynamic behaviour of a three species food web under chemostat conditions. The food web consists of a prey consuming a nutrient, a predator consuming a prey and an omnivore which preys on the predator and the prey. For each trophic level an ordinary differential equation describes the biomass density in the reactor. The hyperbolic functional response for single and multi prey species figures in the description of the trophic interactions. There are two limiting cases where the omnivore is a specialist; a food chain where the omnivore does not consume the prey and competition where the omnivore does not prey on the predator. We use bifurcation analysis to study the long-term dynamic behaviour for various degrees of omnivory. Attractors can be equilibria, limit cycles or chaotic behaviour depending on the control parameters of the chemostat. Often multiple attractor occur. In this paper we will discuss community assembly. That is, we analyze how the trophic structure of the food web evolves following invasion where a new invader is introduced one at the time. Generally, with an invasion, the invader settles itself and persists with all other species, however, the invader may also replace another species. We will show that the food web model has a global bifurcation, being a heteroclinic connection from a saddle equilibrium to a limit cycle of saddle type. This global bifurcation separates regions in the bifurcation diagram with different attractors to which the system evolves after invasion. To investigate the consequences of omnivory we will focus on invasion of the omnivore. This simplifies the analysis considerably, for the end-point of the assembly sequence is then unique. A weak interaction of the omnivore with the prey combined with a stronger interaction with the predator seems advantageous.

Multiple attractors and boundary crises in a tri-trophic food chain

The asymptotic behaviour of a model of a tri-trophic food chain in the chemostat is analysed in detail. The Monod growth model is used for all trophic levels, yielding a non-linear dynamical system of four ordinary differential equations. Mass conservation makes it possible to reduce the dimension by 1 for the study of the asymptotic dynamic behaviour. The intersections of the orbits with a Poincaré plane, after the transient has died out, yield a two-dimensional Poincaré next-return map. When chaotic behaviour occurs, all image points of this next-return map appear to lie close to a single curve in the intersection plane. This motivated the study of a one-dimensional bi-modal, non-invertible map of which the graph resembles this curve. We will show that the bifurcation structure of the food chain model can be understood in terms of the local and global bifurcations of this one-dimensional map. Homoclinic and heteroclinic connecting orbits and their global bifurcations are discussed also by relating them to their counterparts for a two-dimensional map which is invertible like the next-return map. In the global bifurcations two homoclinic or two heteroclinic orbits collide and disappear. In the food chain model two attractors coexist; a stable limit cycle where the top-predator is absent and an interior attractor. In addition there is a saddle cycle. The stable manifold of this limit cycle forms the basin boundary of the interior attractor. We will show that this boundary has a complicated structure when there are heteroclinic orbits from a saddle equilibrium to this saddle limit cycle. A homoclinic bifurcation to a saddle limit cycle will be associated with a boundary crisis where the chaotic attractor disappears suddenly when a bifurcation parameter is varied. Thus, similar to a tangent local bifurcation for equilibria or limit cycles, this homoclinic global bifurcation marks a region in the parameter space where the top-predator goes extinct. The 'Paradox of Enrichment' says that increasing the concentration of nutrient input can cause destabilization of the otherwise stable interior equilibrium of a bi-trophic food chain. For a tri-trophic food chain enrichment of the environment can even lead to extinction of the highest trophic level.

Resistance of a food chain to invasion by a top predator

We study the invasion of a top predator into a food chain in a chemostat. For each trophic level, a bioenergetic model is used in which maintenance and energy reserves are taken into account. Bifurcation analysis is performed on the set of nonlinear ordinary differential equations which describe the dynamic behaviour of the food chain. In this paper, we analyse how the ability of a top predator to invade the food chain depends on the values of two control parameters: the dilution rate and the concentration of the substrate in the input. We investigate invasion by studying the long-term behaviour after introduction of a small amount of top predator. To that end we look at the stability of the boundary attractors; equilibria, limit cycles as well as chaotic attractors using bifurcation analysis. It will be shown that the invasibility criterion is the positiveness of the Lyapunov exponent associated with the change of the biomass of the top predator. It appears that the region in the control parameter space where a predator can invade increases with its growth rate. The resulting system becomes more resistant to further invasion when the top predator grows faster. This implies that short food chains with moderate growth rate of the top predator are liable to be invaded by fast growing invaders which consume the top predator. There may be, however, biological constraints on the top predator's growth rate. Predators are generally larger than prey while larger organisms commonly grow slower. As a result, the growth rate generally decreases with the trophic level. This may enable short food chains to be resistant to invaders. We will relate these results to ecological community assembly and the debate on the length of food chains in nature.

Consequences of population models for the dynamics of food chains

A class of bioenergetic ecological models is studied for the dynamics of food chains with a nutrient at the base. A constant influx rate of the nutrient and a constant efflux rate for all trophic levels is assumed. Starting point is a simple model where prey is converted into predator with a fixed efficiency. This model is extended by the introduction of maintenance and energy reserves at all trophic levels, with two state variables for each trophic level, biomass and reserve energy. Then the dynamics of each population are described by two ordinary differential equations. For all models the bifurcation diagram for the bi-trophic food chain is simple. There are three important regions; a region where the predator goes to extinction, a region where there is a stable equilibrium and a region where a stable limit cycle exists. Bifurcation diagrams for tritrophic food chains are more complicated. Flip bifurcation curves mark regions where complex dynamic behaviour (higher periodic limit cycles as well as chaotic attractors) can occur. We show numerically that Shil'nikov homoclinic orbits to saddle-focus equilibria exists. The codimension 1 continuations of these orbits form a 'skeleton' for a cascade of flip and tangent bifurcations. The bifurcation analysis facilitates the study of the consequences of the population model for the dynamic behaviour of a food chain. Although the predicted transient dynamics of a food chain may depend sensitively on the underlying model for the populations, the global picture of the bifurcation diagram for the different models is about the same.

Food chain dynamics in the chemostat

The asymptotic behavior of a tri-trophic food chain model in the chemostat is studied. The Monod-Herbert growth model is used for all trophic levels. The analysis is carried out numerically, by finding both local and global bifurcations of equilibria and of limit cycles with respect to two chemostat control parameters: the dilution rate of the chemostat and the concentration of input substrate. It is shown that the bifurcation structure of the food chain model has much in common with the bifurcation structure of a one-dimensional map with two turning points. This map is used to explain how attractors are created and destroyed under variation of the bifurcation parameters. It is shown that low as well as high concentration of input substrate can lead to extinction of the highest trophic level.