Skyline‐plot methods for estimating demographic history from nucleotide sequences

Accounting for reporting delays in real-time phylodynamic analyses with preferential sampling

The COVID-19 pandemic demonstrated that fast and accurate analysis of continually collected infectious disease surveillance data is crucial for situational awareness and policy making. Coalescent-based phylodynamic analysis can use genetic sequences of a pathogen to estimate changes in its effective population size, a measure of genetic diversity. These changes in effective population size can be connected to the changes in the number of infections in the population of interest under certain conditions. Phylodynamics is an important set of tools because its methods are often resilient to the ascertainment biases present in traditional surveillance data (e.g., preferentially testing symptomatic individuals). Unfortunately, it takes weeks or months to sequence and deposit the sampled pathogen genetic sequences into a database, making them available for such analyses. These reporting delays severely decrease precision of phylodynamic methods closer to present time, and for some models can lead to extreme biases. Here we present a method that affords reliable estimation of the effective population size trajectory closer to the time of data collection, allowing for policy decisions to be based on more recent data. Our work uses readily available historic times between sampling and reporting of sequenced samples for a population of interest, and incorporates this information into the sampling model to mitigate the effects of reporting delay in real-time analyses. We illustrate our methodology on simulated data and on SARS-CoV-2 sequences collected in the state of Washington in 2021.

Ancestral process for infectious disease outbreaks with superspreading

When an infectious disease outbreak is of a relatively small size, describing the ancestry of a sample of infected individuals is difficult because most ancestral models assume large population sizes. Given a set of infected individuals, we show that it is possible to express exactly the probability that they have the same infector, either inclusively (so that other individuals may have the same infector too) or exclusively (so that they may not). To compute these probabilities requires knowledge of the offspring distribution, which determines how many infections each infected individual causes. We consider transmission both without and with superspreading, in the form of a Poisson and a Negative-Binomial offspring distribution, respectively. We show how our results can be incorporated into a new Lambda-coalescent model which allows multiple lineages to coalesce together. We call this new model the Omega-coalescent, we compare it with previously proposed alternatives, and advocate its use in future studies of infectious disease outbreaks.

Legacy parasite collections reveal species-specific population genetic patterns among three species of zoonotic schistosomes

Studies estimating genetic diversity and population structure in multi-host parasites are often constrained by temporally and spatially limited sampling. This study addresses these limitations by analyzing globally distributed samples of three congeneric avian schistosomes (Trematoda: Schistosomatidae: Trichobilharzia), including collections spanning 20 years archived at The Museum of Southwestern Biology, Parasites Division. The three species exhibited significant differences in population genetic parameters across one nuclear and two mitochondrial loci. Trichobilharzia querquedulae (TQ) maintained a well-connected, globally diverse metapopulation, with an effective population size approximately three times larger than that of the other two species, T. physellae (TP) and Trichobilharzia sp. A (TA). TP and TA had lower overall genetic diversity and greater population structure. These differences are likely shaped by the ecologies of the duck definitive hosts that disperse these parasites. This study highlights the value of natural history collections, particularly since Trichobilharzia is a key agent of zoonotic cercarial dermatitis, a disease whose etiology and epidemiology remain poorly understood. Within a comparative congeneric framework, population genetic data can provide insights into host-parasite natural history and its influence on microevolutionary patterns, including contributions to zoonotic disease.

Multiple merger coalescent inference of effective population size

Variation in a sample of molecular sequence data informs about the past evolutionary history of the sample's population. Traditionally, Bayesian modelling coupled with the standard coalescent is used to infer the sample's bifurcating genealogy and demographic and evolutionary parameters such as effective population size and mutation rates. However, there are many situations where binary coalescent models do not accurately reflect the true underlying ancestral processes. Here, we propose a Bayesian non-parametric method for inferring effective population size trajectories from a multifurcating genealogy under the [Formula: see text]-coalescent. In particular, we jointly estimate the effective population size and the model parameter for the Beta-coalescent model, a special type of [Formula: see text]-coalescent. Finally, we test our methods on simulations and apply them to study various viral dynamics as well as Japanese sardine population size changes over time. The code and vignettes can be found in the phylodyn package.This article is part of the theme issue '"A mathematical theory of evolution": phylogenetic models dating back 100 years'.

High genetic variation, low differentiation, and Pleistocene expansions of the migratory and endangered long-nosed tequila bat, Leptonycteris nivalis, inferred using both maternal and paternal genetic markers

Tequila bats (genus Leptonycteris) have gained attention for their critical role in pollinating different plant species, especially Agave spp. and columnar cacti. Leptonycteris nivalis is the largest nectar-feeding bat in the Americas, and the females exhibit migratory behavior during the breeding season. Due to its relatively small and seemingly declining population sizes, this species is protected by government agencies in the United States and Mexico. We conducted population genetics and phylogeographic analyses to elucidate the genetic structure and demographic history of the species using two mitochondrial markers and a Y chromosome-associated gene, to describe both maternal and paternal lineages. We estimated high haplotypic diversity measures for the different markers (Dloop-Hd = 0.775; Cyt-b-Hd = 0.937; DBY -Hd = 0.946). We found that geographic genetic differentiation is very low, and there is high connectivity among localities. The estimated divergence time between L. nivalis and L. yerbabuenae, the other species in the genus found in Mexico, aligns with previous estimates for the genus (6.91-9.43 mya). A demographic expansion was detected approximately at 600 ka-700 ka (thousands of years ago). The historical demographic changes observed in L. nivalis appear to be associated with environmental shifts during the Pleistocene, which likely impacted the distribution range of the plants that these bats feed on, such as Agave species.

Novel B2 mitogenomes from Continental southern Patagonia's Late Holocene: New insights into the peopling of the Southern Cone

OBJECTIVES

The main aim of this study is to discuss the migratory processes and peopling dynamics that shaped the genetic variability of populations during the settlement of the Southern Cone, through the analysis of complete mitogenomes of individuals from southern Patagonia.

MATERIALS AND METHODS

Complete mitogenomes were sequenced through massively parallel sequencing from two late Holocene individuals (SAC 1-1-3 and SAC 1-1-4) buried in the same chenque at Salitroso Lake Basin (Santa Cruz province, Argentina). To evaluate matrilineal phylogenetic affinities with other haplotypes, maximum likelihood and Bayesian phylogenetic reconstructions were performed, as well as a haplotype median-joining network.

RESULTS

The mitogenomes were assigned to haplogroups B2 and B2b, exhibiting an average depth of 54X and 89X (≥1X coverage of 98.6% and 100%), and a high number of nucleotide differences among them. The phylogenetic analyses showed a relatively close relationship between the haplotype found in SAC 1-1-4 and those retrieved from a Middle Holocene individual from Laguna Chica (Buenos Aires province), and from a group of individuals from the Peruvian coast. For the SAC 1-1-3, no clear affiliations to any other haplotype were established.

DISCUSSION

The large divergence between the haplotypes presented in this study suggests either a highly variable founder gene pool, or a later enrichment by frequent biological contact with other populations. Our results underline the persistence of genetic signals related to the first waves of peopling in South America, suggesting that the regional settlement of the southern end of the continent has been much more complex than initially thought.

Modeling the velocity of evolving lineages and predicting dispersal patterns

Accurate estimation of the dispersal velocity or speed of evolving organisms is no mean feat. In fact, existing probabilistic models in phylogeography or spatial population genetics generally do not provide an adequate framework to define velocity in a relevant manner. For instance, the very concept of instantaneous speed simply does not exist under one of the most popular approaches that models the evolution of spatial coordinates as Brownian trajectories running along a phylogeny. Here, we introduce a family of models-the so-called Phylogenetic Integrated Velocity (PIV) models-that use Gaussian processes to explicitly model the velocity of evolving lineages instead of focusing on the fluctuation of spatial coordinates over time. We describe the properties of these models and show an increased accuracy of velocity estimates compared to previous approaches. Analyses of West Nile virus data in the United States indicate that PIV models provide sensible predictions of the dispersal of evolving pathogens at a one-year time horizon. These results demonstrate the feasibility and relevance of predictive phylogeography in monitoring epidemics in time and space.

Assessing the effect of model specification and prior sensitivity on Bayesian tests of temporal signal

Our understanding of the evolution of many microbes has been revolutionised by the molecular clock, a statistical tool to infer evolutionary rates and timescales from analyses of biomolecular sequences. In all molecular clock models, evolutionary rates and times are jointly unidentifiable and 'calibration' information must therefore be used. For many organisms, sequences sampled at different time points can be employed for such calibration. Before attempting to do so, it is recommended to verify that the data carry sufficient information for molecular dating, a practice referred to as evaluation of temporal signal. Recently, a fully Bayesian approach, BETS (Bayesian Evaluation of Temporal Signal), was proposed to overcome known limitations of other commonly used techniques such as root-to-tip regression or date randomisation tests. BETS requires the specification of a full Bayesian phylogenetic model, posing several considerations for untangling the impact of model choice on the detection of temporal signal. Here, we aimed to (i) explore the effect of molecular clock model and tree prior specification on the results of BETS and (ii) provide guidelines for improving our confidence in molecular clock estimates. Using microbial molecular sequence data sets and simulation experiments, we assess the impact of the tree prior and its hyperparameters on the accuracy of temporal signal detection. In particular, highly informative priors that are inconsistent with the data can result in the incorrect detection of temporal signal. In consequence, we recommend: (i) using prior predictive simulations to determine whether the prior generates a reasonable expectation of parameters of interest, such as the evolutionary rate and age of the root node, (ii) conducting prior sensitivity analyses to assess the robustness of the posterior to the choice of prior, and (iii) selecting a molecular clock model that reasonably describes the evolutionary process.

Modeling the velocity of evolving lineages and predicting dispersal patterns

Accurate estimation of the dispersal velocity or speed of evolving organisms is no mean feat. In fact, existing probabilistic models in phylogeography or spatial population genetics generally do not provide an adequate framework to define velocity in a relevant manner. For instance, the very concept of instantaneous speed simply does not exist under one of the most popular approaches that models the evolution of spatial coordinates as Brownian trajectories running along a phylogeny (Lemey et al., 2010). Here, we introduce a new family of models - the so-called "Phylogenetic Integrated Velocity" (PIV) models - that use Gaussian processes to explicitly model the velocity of evolving lineages instead of focusing on the fluctuation of spatial coordinates over time. We describe the properties of these models and show an increased accuracy of velocity estimates compared to previous approaches. Analyses of West Nile virus data in the U.S.A. indicate that PIV models provide sensible predictions of the dispersal of evolving pathogens at a one-year time horizon. These results demonstrate the feasibility and relevance of predictive phylogeography in monitoring epidemics in time and space.

Genetic Insights into the Giant Keyhole Limpet (Megathura crenulata), an Eastern Pacific Coastal Endemic: Complete Mitogenome, Phylogenetics, Phylogeography, and Historical Demography

BACKGROUND

The giant keyhole limpet Megathura crenulata is a gastropod mollusk (Fissurella superfamily) that is endemic to the eastern Pacific coast from southern California, USA, to Baja California Sur, Mexico. M. crenulata is socioeconomically important as it produces a potent immune-stimulating protein, called Keyhole Limpet Hemocyanin, which is extracted in vivo and utilized for vaccine development. However, ecological studies are scarce and genetic knowledge of the species needs to be improved. Our objectives were to assemble and annotate the mitogenome of M. crenulata, and to assess its phylogenetic relationships with other marine gastropods and to evaluate its population genetic diversity and structure.

METHODS

Samples were collected for mitogenome assembly (n = 3) spanning its geographic range, Puerto Canoas (PCA) and Punta Eugenia (PEU), Mexico, and California (CAL), USA. Total DNA was extracted from gills sequenced using Illumina paired-end 150-bp-read sequencing. Reads were cleaned, trimmed, assembled de novo, and annotated. In addition, 125 samples from eight locations were analyzed for genetic diversity and structure analysis at the 16s rRNA and COX1 genes.

RESULTS

The M. crenulata mitogenomes had lengths of 16,788 bp (PCA) and 16,787 bp (PEU) and were composed of 13 protein-coding regions, 22 tRNAs, two rRNAs, and the D-Loop region. In terms of phylogeographic diversity and structure, we found a panmictic population that has experienced recent demographic expansion with low nucleotide diversity (0.002), high haplotypic diversity (0.915), and low φST (0.047).

CONCLUSIONS

Genetic insights into the giant keyhole limpet provides tools for its management and conservation by delimiting fishing regions with low genetic diversity and/or genetically discrete units.

A comparison of genomic diversity and demographic history of the North Atlantic and Southwest Atlantic southern right whales

Right whales (genus Eubalaena) were among the first, and most extensively pursued, targets of commercial whaling. However, understanding the impacts of this persecution requires knowledge of the demographic histories of these species prior to exploitation. We used deep whole genome sequencing (~40×) of 12 North Atlantic (E. glacialis) and 10 Southwest Atlantic southern (E. australis) right whales to quantify contemporary levels of genetic diversity and infer their demographic histories over time. Using coalescent- and identity-by-descent-based modelling to estimate ancestral effective population sizes from genomic data, we demonstrate that North Atlantic right whales have lived with smaller effective population sizes (Ne) than southern right whales in the Southwest Atlantic since their divergence and describe the decline in both populations around the time of whaling. North Atlantic right whales exhibit reduced genetic diversity and longer runs of homozygosity leading to higher inbreeding coefficients compared to the sampled population of southern right whales. This study represents the first comprehensive assessment of genome-wide diversity of right whales in the western Atlantic and underscores the benefits of high coverage, genome-wide datasets to help resolve long-standing questions about how historical changes in effective population size over different time scales shape contemporary diversity estimates. This knowledge is crucial to improve our understanding of the right whales' history and inform our approaches to address contemporary conservation issues. Understanding and quantifying the cumulative impact of long-term small Ne, low levels of diversity and recent inbreeding on North Atlantic right whale recovery will be important next steps.

Impact of pneumococcal conjugate vaccines on invasive pneumococcal disease-causing lineages among South African children

Invasive pneumococcal disease (IPD) due to non-vaccine serotypes after the introduction of pneumococcal conjugate vaccines (PCV) remains a global concern. This study used pathogen genomics to evaluate changes in invasive pneumococcal lineages before, during and after vaccine introduction in South Africa. We included genomes (N = 3104) of IPD isolates from individuals aged <18 years (2005-20), spanning four periods: pre-PCV, PCV7, early-PCV13, and late-PCV13. Significant incidence reductions occurred among vaccine-type lineages in the late-PCV13 period compared to the pre-PCV period. However, some vaccine-type lineages continued to cause invasive disease and showed increasing effective population size trends in the post-PCV era. A significant increase in lineage diversity was observed from the PCV7 period to the early-PCV13 period (Simpson's diversity index: 0.954, 95% confidence interval 0.948-0.961 vs 0.965, 0.962-0.969) supporting intervention-driven population structure perturbation. Increases in the prevalence of penicillin, erythromycin, and multidrug resistance were observed among non-vaccine serotypes in the late-PCV13 period compared to the pre-PCV period. In this work we highlight the importance of continued genomic surveillance to monitor disease-causing lineages post vaccination to support policy-making and future vaccine designs and considerations.

Reconstruction of Avian Reovirus History and Dispersal Patterns: A Phylodynamic Study

Avian reovirus (ARV) infection can cause significant losses to the poultry industry. Disease control has traditionally been attempted mainly through vaccination. However, the increase in clinical outbreaks in the last decades demonstrated the poor effectiveness of current vaccination approaches. The present study reconstructs the evolution and molecular epidemiology of different ARV genotypes using a phylodynamic approach, benefiting from a collection of more than one thousand sigma C (σC) sequences sampled over time at a worldwide level. ARVs' origin was estimated to occur several centuries ago, largely predating the first clinical reports. The origins of all genotypes were inferred at least one century ago, and their emergence and rise reflect the intensification of the poultry industry. The introduction of vaccinations had only limited and transitory effects on viral circulation and further expansion was observed, particularly after the 1990s, likely because of the limited immunity and the suboptimal and patchy vaccination application. In parallel, strong selective pressures acted with different strengths and directionalities among genotypes, leading to the emergence of new variants. While preventing the spread of new variants with different phenotypic features would be pivotal, a phylogeographic analysis revealed an intricate network of viral migrations occurring even over long distances and reflecting well-established socio-economic relationships.

Comparison of Bayesian Coalescent Skyline Plot Models for Inferring Demographic Histories

Bayesian coalescent skyline plot models are widely used to infer demographic histories. The first (non-Bayesian) coalescent skyline plot model assumed a known genealogy as data, while subsequent models and implementations jointly inferred the genealogy and demographic history from sequence data, including heterochronous samples. Overall, there exist multiple different Bayesian coalescent skyline plot models which mainly differ in two key aspects: (i) how changes in population size are modeled through independent or autocorrelated prior distributions, and (ii) how many change-points in the demographic history are used, where they occur and if the number is pre-specified or inferred. The specific impact of each of these choices on the inferred demographic history is not known because of two reasons: first, not all models are implemented in the same software, and second, each model implementation makes specific choices that the biologist cannot influence. To facilitate a detailed evaluation of Bayesian coalescent skyline plot models, we implemented all currently described models in a flexible design into the software RevBayes. Furthermore, we evaluated models and choices on an empirical dataset of horses supplemented by a small simulation study. We find that estimated demographic histories can be grouped broadly into two groups depending on how change-points in the demographic history are specified (either independent of or at coalescent events). Our simulations suggest that models using change-points at coalescent events produce spurious variation near the present, while most models using independent change-points tend to over-smooth the inferred demographic history.

Ancient mitogenomes from Pre-Pottery Neolithic Central Anatolia and the effects of a Late Neolithic bottleneck in sheep (Ovis aries)

Occupied between ~10,300 and 9300 years ago, the Pre-Pottery Neolithic site of Aşıklı Höyük in Central Anatolia went through early phases of sheep domestication. Analysis of 629 mitochondrial genomes from this and numerous sites in Anatolia, southwest Asia, Europe, and Africa produced a phylogenetic tree with excessive coalescences (nodes) around the Neolithic, a potential signature of a domestication bottleneck. This is consistent with archeological evidence of sheep management at Aşıklı Höyük which transitioned from residential stabling to open pasturing over a millennium of site occupation. However, unexpectedly, we detected high genetic diversity throughout Aşıklı Höyük's occupation rather than a bottleneck. Instead, we detected a tenfold demographic bottleneck later in the Neolithic, which caused the fixation of mitochondrial haplogroup B in southwestern Anatolia. The mitochondrial genetic makeup that emerged was carried from the core region of early Neolithic sheep management into Europe and dominates the matrilineal diversity of both its ancient and the billion-strong modern sheep populations.

Deep-pelagic fishes: Demographic instability in a stable environment

Demographic histories are frequently a product of the environment, as populations expand or contract in response to major environmental changes, often driven by changes in climate. Meso- and bathy-pelagic fishes inhabit some of the most temporally and spatially stable habitats on the planet. The stability of the deep-pelagic could make deep-pelagic fishes resistant to the demographic instability commonly reported in fish species inhabiting other marine habitats, however the demographic histories of deep-pelagic fishes are unknown. We reconstructed the historical demography of 11 species of deep-pelagic fishes using mitochondrial and nuclear DNA sequence data. We uncovered widespread evidence of population expansions in our study species, a counterintuitive result based on the nature of deep-pelagic ecosystems. Frequency-based methods detected potential demographic changes in nine species of fishes, while extended Bayesian skyline plots identified population expansions in four species. These results suggest that despite the relatively stable nature of the deep-pelagic environment, the fishes that reside here have likely been impacted by past changes in climate. Further investigation is necessary to better understand how deep-pelagic fishes, by far Earth's most abundant vertebrates, will respond to future climatic changes.

The effect of missing data on evolutionary analysis of sequence capture bycatch, with application to an agricultural pest

Sequence capture is a genomic technique that selectively enriches target sequences before high throughput next-generation sequencing, to generate specific sequences of interest. Off-target or 'bycatch' data are often discarded from capture experiments, but can be leveraged to address evolutionary questions under some circumstances. Here, we investigated the effects of missing data on a variety of evolutionary analyses using bycatch from an exon capture experiment on the global pest moth, Helicoverpa armigera. We added > 200 new samples from across Australia in the form of mitogenomes obtained as bycatch from targeted sequence capture, and combined these into an additional larger dataset to total > 1000 mitochondrial cytochrome c oxidase subunit I (COI) sequences across the species' global distribution. Using discriminant analysis of principal components and Bayesian coalescent analyses, we showed that mitogenomes assembled from bycatch with up to 75% missing data were able to return evolutionary inferences consistent with higher coverage datasets and the broader literature surrounding H. armigera. For example, low-coverage sequences broadly supported the delineation of two H. armigera subspecies and also provided new insights into the potential for geographic turnover among these subspecies. However, we also identified key effects of dataset coverage and composition on our results. Thus, low-coverage bycatch data can offer valuable information for population genetic and phylodynamic analyses, but caution is required to ensure the reduced information does not introduce confounding factors, such as sampling biases, that drive inference. We encourage more researchers to consider maximizing the potential of the targeted sequence approach by examining evolutionary questions with their off-target bycatch where possible-especially in cases where no previous mitochondrial data exists-but recommend stratifying data at different genome coverage thresholds to separate sampling effects from genuine genomic signals, and to understand their implications for evolutionary research.

Past volcanic activity predisposes an endemic threatened seabird to negative anthropogenic impacts

Humans are regularly cited as the main driver of current biodiversity extinction, but the impact of historic volcanic activity is often overlooked. Pre-human evidence of wildlife abundance and diversity are essential for disentangling anthropogenic impacts from natural events. Réunion Island, with its intense and well-documented volcanic activity, endemic biodiversity, long history of isolation and recent human colonization, provides an opportunity to disentangle these processes. We track past demographic changes of a critically endangered seabird, the Mascarene petrel Pseudobulweria aterrima, using genome-wide SNPs. Coalescent modeling suggested that a large ancestral population underwent a substantial population decline in two distinct phases, ca. 125,000 and 37,000 years ago, coinciding with periods of major eruptions of Piton des Neiges. Subsequently, the ancestral population was fragmented into the two known colonies, ca. 1500 years ago, following eruptions of Piton de la Fournaise. In the last century, both colonies declined significantly due to anthropogenic activities, and although the species was initially considered extinct, it was rediscovered in the 1970s. Our findings suggest that the current conservation status of wildlife on volcanic islands should be firstly assessed as a legacy of historic volcanic activity, and thereafter by the increasing anthropogenic impacts, which may ultimately drive species towards extinction.

Inference of Infectious Disease Transmission through a Relaxed Bottleneck Using Multiple Genomes Per Host

In recent times, pathogen genome sequencing has become increasingly used to investigate infectious disease outbreaks. When genomic data is sampled densely enough amongst infected individuals, it can help resolve who infected whom. However, transmission analysis cannot rely solely on a phylogeny of the genomes but must account for the within-host evolution of the pathogen, which blurs the relationship between phylogenetic and transmission trees. When only a single genome is sampled for each host, the uncertainty about who infected whom can be quite high. Consequently, transmission analysis based on multiple genomes of the same pathogen per host has a clear potential for delivering more precise results, even though it is more laborious to achieve. Here, we present a new methodology that can use any number of genomes sampled from a set of individuals to reconstruct their transmission network. Furthermore, we remove the need for the assumption of a complete transmission bottleneck. We use simulated data to show that our method becomes more accurate as more genomes per host are provided, and that it can infer key infectious disease parameters such as the size of the transmission bottleneck, within-host growth rate, basic reproduction number, and sampling fraction. We demonstrate the usefulness of our method in applications to real datasets from an outbreak of Pseudomonas aeruginosa amongst cystic fibrosis patients and a nosocomial outbreak of Klebsiella pneumoniae.

Parallel signatures of Mycobacterium tuberculosis and human Y-chromosome phylogeography support the Two Layer model of East Asian population history

The Two Layer hypothesis is fast becoming the favoured narrative describing East Asian population history. Under this model, hunter-gatherer groups who initially peopled East Asia via a route south of the Himalayas were assimilated by agriculturalist migrants who arrived via a northern route across Eurasia. A lack of ancient samples from tropical East Asia limits the resolution of this model. We consider insight afforded by patterns of variation within the human pathogen Mycobacterium tuberculosis (Mtb) by analysing its phylogeographic signatures jointly with the human Y-chromosome. We demonstrate the Y-chromosome lineages enriched in the traditionally hunter-gatherer groups associated with East Asia's first layer of peopling to display deep roots, low long-term effective population size, and diversity patterns consistent with a southern entry route. These characteristics mirror those of the evolutionarily ancient Mtb lineage 1. The remaining East Asian Y-chromosome lineage is almost entirely absent from traditionally hunter-gatherer groups and displays spatial and temporal characteristics which are incompatible with a southern entry route, and which link it to the development of agriculture in modern-day China. These characteristics mirror those of the evolutionarily modern Mtb lineage 2. This model paves the way for novel host-pathogen coevolutionary research hypotheses in East Asia.

The raw material of evolution

Estimates of whale mutation rates contribute to understanding evolutionary processes.

Phylogeography of the desert scorpion illuminates a route out of Central Asia

A comprehensive understanding of phylogeography requires the integration of knowledge across different organisms, ecosystems, and geographic regions. However, a critical knowledge gap exists in the arid biota of the vast Asian drylands. To narrow this gap, here we test an "out-of-Central Asia" hypothesis for the desert scorpion Mesobuthus mongolicus by combining Bayesian phylogeographic reconstruction and ecological niche modeling. Phylogenetic analyses of one mitochondrial and three nuclear loci and molecular dating revealed that M. mongolicus represents a coherent lineage that diverged from its most closely related lineage in Central Asia about 1.36 Ma and underwent radiation ever since. Bayesian phylogeographic reconstruction indicated that the ancestral population dispersed from Central Asia gradually eastward to the Gobi region via the Junggar Basin, suggesting that the Junggar Basin has severed as a corridor for Quaternary faunal exchange between Central Asia and East Asia. Two major dispersal events occurred probably during interglacial periods (around 0.8 and 0.4 Ma, respectively) when climatic conditions were analogous to present-day status, under which the scorpion achieved its maximum distributional range. M. mongolicus underwent demographic expansion during the Last Glacial Maximum, although the predicted distributional areas were smaller than those at present and during the Last Interglacial. Development of desert ecosystems in northwest China incurred by intensified aridification might have opened up empty habitats that sustained population expansion. Our results extend the spatiotemporal dimensions of trans-Eurasia faunal exchange and suggest that species' adaptation is an important determinant of their phylogeographic and demographic responses to climate changes.

Testing the efficacy of different molecular tools for parasite conservation genetics: a case study using horsehair worms (Phylum: Nematomorpha)

In recent years, parasite conservation has become a globally significant issue. Because of this, there is a need for standardized methods for inferring population status and possible cryptic diversity. However, given the lack of molecular data for some groups, it is challenging to establish procedures for genetic diversity estimation. Therefore, universal tools, such as double-digest restriction-site-associated DNA sequencing (ddRADseq), could be useful when conducting conservation genetic studies on rarely studied parasites. Here, we generated a ddRADseq dataset that includes all 3 described Taiwanese horsehair worms (Phylum: Nematomorpha), possibly one of the most understudied animal groups. Additionally, we produced data for a fragment of the cytochrome c oxidase subunit I (COXI) for the said species. We used the COXI dataset in combination with previously published sequences of the same locus for inferring the effective population size (Ne) trends and possible population genetic structure.We found that a larger and geographically broader sample size combined with more sequenced loci resulted in a better estimation of changes in Ne. We were able to detect demographic changes associated with Pleistocene events in all the species. Furthermore, the ddRADseq dataset for Chordodes formosanus did not reveal a genetic structure based on geography, implying a great dispersal ability, possibly due to its hosts. We showed that different molecular tools can be used to reveal genetic structure and demographic history at different historical times and geographical scales, which can help with conservation genetic studies in rarely studied parasites.

Inferring Evolutionary Timescale of Omsk Hemorrhagic Fever Virus

Until 2020, there were only three original complete genome (CG) nucleotide sequences of Omsk hemorrhagic fever virus (OHFV) in GenBank. For this reason, the evolutionary rate and divergence time assessments reported in the literature were based on the E gene sequences, but notably without temporal signal evaluation, such that their reliability is unclear. As of July 2022, 47 OHFV CG sequences have been published, which enables testing of temporal signal in the data and inferring unbiased and reliable substitution rate and divergence time values. Regression analysis in the TempEst software demonstrated a stronger clocklike behavior in OHFV samples for the complete open reading frame (ORF) data set (R² = 0.42) than for the E gene data set (R² = 0.11). Bayesian evaluation of temporal signal indicated very strong evidence, with a log Bayes factor of more than 5, in favor of temporal signal in all data sets. Our results based on the complete ORF sequences showed a more precise OHFV substitution rate (95% highest posterior density (HPD) interval, 9.1 × 10^-5-1.8 × 10^-4 substitutions per site per year) and tree root height (416-896 years ago) compared with previous assessments. The rate obtained is significantly higher than tick-borne encephalitis virus by at least 3.8-fold. The phylogenetic analysis and past population dynamics reconstruction revealed the declining trend of OHFV genetic diversity, but there was phylogenomic evidence that implicit virus subpopulations evolved locally and underwent an exponential growth phase.

Comparative phylogenetic analysis and transcriptomic profiling of Dengue (DENV-3 genotype I) outbreak in 2021 in Bangladesh

Background

The next-generation sequencing (NGS) technology facilitates in-depth study of host-pathogen metatranscriptome. We, therefore, implicated phylodynamic and transcriptomic approaches through NGS technology to know/understand the dengue virus (DENV) origin and host response with dengue fever.

Methods

In this study, blood serum RNA was extracted from 21 dengue patients and 3 healthy individuals. Total transcriptomic data were analyzed for phylogenetic, phylodynamic, differential express gene (DEG), and gene ontology (GO) using respective bioinformatics tools.

Results

The viral genome sequence revealed dengue viral genome size ranges 10647 to 10707 nucleotide. Phylogenetic and phylodynamic analysis showed that the 2021 epidemic isolates were DENV-3 genotype-I and maintained as a new clade in compared to 2019 epidemic. Transcriptome analysis showed a total of 2686 genes were DEG in dengue patients compared to control with a q-value < 0.05. DESeq2 plot counts function of the top 24 genes with the smallest q-values of differential gene expression of RNA-seq data showed that 11 genes were upregulated, whereas 13 genes were downregulated. GO analysis showed a significant upregulation (p = < 0.001) in a process of multicellular organismal, nervous system, sensory perception of chemical stimulus, and G protein-coupled receptor signaling pathways in the dengue patients. However, there were a significant downregulation (p = < 0.001) of intracellular component, cellular anatomical entity, and protein-containing complex in dengue patients. Most importantly, there was a significant increase of a class of immunoregulatory proteins in dengue patients in compared to the controls, with increased GO of immune system process. In addition, upregulation of toll receptor (TLR) signaling pathways were found in dengue patients. These TLR pathways were particularly involved for the activation of innate system coupled with adaptive immune system that probably involved the rapid elimination of dengue virus infected cells. These differentially expressed genes could be further investigated for target based prophylactic interventions for dengue.

Conclusion

This is a first report describing DENV complete genomic features and differentially expressed genes in patients in Bangladesh. These genes may have diagnostic and therapeutic values for dengue infection. Continual genomic surveillance is required to further investigate the shift in dominant genotypes in relation to viral pathogenesis.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12985-023-02030-1.

Estimating the fitness cost and benefit of antimicrobial resistance from pathogen genomic data

Increasing levels of antibiotic resistance in many bacterial pathogen populations are a major threat to public health. Resistance to an antibiotic provides a fitness benefit when the bacteria are exposed to this antibiotic, but resistance also often comes at a cost to the resistant pathogen relative to susceptible counterparts. We lack a good understanding of these benefits and costs of resistance for many bacterial pathogens and antibiotics, but estimating them could lead to better use of antibiotics in a way that reduces or prevents the spread of resistance. Here, we propose a new model for the joint epidemiology of susceptible and resistant variants, which includes explicit parameters for the cost and benefit of resistance. We show how Bayesian inference can be performed under this model using phylogenetic data from susceptible and resistant lineages and that by combining data from both we are able to disentangle and estimate the resistance cost and benefit parameters separately. We applied our inferential methodology to several simulated datasets to demonstrate good scalability and accuracy. We analysed a dataset of Neisseria gonorrhoeae genomes collected between 2000 and 2013 in the USA. We found that two unrelated lineages resistant to fluoroquinolones shared similar epidemic dynamics and resistance parameters. Fluoroquinolones were abandoned for the treatment of gonorrhoea due to increasing levels of resistance, but our results suggest that they could be used to treat a minority of around 10% of cases without causing resistance to grow again.

Molecular epidemiology of Canine circovirus in domestic dogs and wildlife in Namibia, Africa

Canine circovirus (CanineCV) is a DNA virus affecting domestic dogs and other wild carnivore species. Despite the potential implications for dogs' health and wildlife conservation, data on CanineCV presence, epidemiology and genetic features from Africa is still poor. In the present study, biological specimens collected between 2020 and 2022 from a total of 32 jackals and 575 domestic dogs were tested for the presence of CanineCV DNA to evaluate its frequency. Furthermore, sequencing was conducted on positive samples to characterize the strains and compare them with publicly available sequences through phylogenetic analysis. A high CanineCV prevalence was observed both in jackals (43.75%; 95 CI: 28.17% - 60.67%) and domestic dogs (27.13%; 95 CI: 23.66% - 30.91%). All aside from one Namibian strain formed an independent clade, suggestive of extremely rare introduction events, followed by local persistence, circulation, and evolution. Remarkably, different recombination events were observed involving strains from both jackals and domestic dogs, which testify to the likely strain exchange between these populations. Distinctive amino acid residues were also observed in jackals. The limitations of the considered host populations however prevent a definitive conclusion on host adaptation, biological, and clinical features. Further studies should be performed to expand our current knowledge of the CanineCV disease scenario in Namibia, other African regions, and associated host species in Africa.

Model design for nonparametric phylodynamic inference and applications to pathogen surveillance

Inference of effective population size from genomic data can provide unique information about demographic history and, when applied to pathogen genetic data, can also provide insights into epidemiological dynamics. The combination of nonparametric models for population dynamics with molecular clock models which relate genetic data to time has enabled phylodynamic inference based on large sets of time-stamped genetic sequence data. The methodology for nonparametric inference of effective population size is well-developed in the Bayesian setting, but here we develop a frequentist approach based on nonparametric latent process models of population size dynamics. We appeal to statistical principles based on out-of-sample prediction accuracy in order to optimize parameters that control shape and smoothness of the population size over time. Our methodology is implemented in a new R package entitled mlesky. We demonstrate the flexibility and speed of this approach in a series of simulation experiments and apply the methodology to a dataset of HIV-1 in the USA. We also estimate the impact of non-pharmaceutical interventions for COVID-19 in England using thousands of SARS-CoV-2 sequences. By incorporating a measure of the strength of these interventions over time within the phylodynamic model, we estimate the impact of the first national lockdown in the UK on the epidemic reproduction number.

Evolution and epidemic success of Mycobacterium tuberculosis in eastern China: evidence from a prospective study

BACKGROUND

Lineage distribution of Mycobacterium tuberculosis (Mtb) isolates is strongly associated with geographically distinct human populations, and its transmission can be further impacted by the bacterial genome. However, the epidemic success of Mtb isolates at an individual level was unknown in eastern China. Knowledge regarding the emergence and transmission of Mtb isolates as well as relevant factors may offer a new solution to curb the spread of the disease. Thus, this study aims to reveal the evolution and epidemic success of Mtb isolates in eastern China.

RESULTS

Of initial 1040 isolates, 997 were retained after removing duplicates and those with insufficient sequencing depth. Of the final samples, 733 (73.52%) were from Zhejiang Province, and 264 (26.48%) were from Shanghai City. Lineage 2 and lineage 4 accounted for 80.44% and 19.56%, with common ancestors dating around 7017 years ago and 6882 years ago, respectively. Sub-lineage L2.2 (80.34%) contributed the majority of total isolates, followed by L4.4 (8.93%) and L4.5 (8.43%). Additionally, 51 (5.12%) isolates were identified to be multidrug-resistant (MDR), of which 21 (29.17%) were pre-extensively drug-resistant (pre-XDR). One clade harboring katG S315T mutation may date back to 65 years ago and subsequently acquired mutations conferring resistance to another five antibiotic drugs. The prevalence of compensatory mutation was the highest in pre-XDR isolates (76.19%), followed by MDR isolates (47.06%) and other drug-resistant isolates (20.60%). Time-scaled haplotypic density analyses suggested comparable success indices between lineage 2 and lineage 4 (P = 0.306), and drug resistance did not significantly promote the transmission of Mtb isolates (P = 0.340). But for pre-XDR isolates, we found a higher success index in those with compensatory mutations (P = 0.025). Mutations under positive selection were found in genes associated with resistance to second-line injectables (whiB6) and drug tolerance (prpR) in both lineage 2 and lineage 4.

CONCLUSIONS

Our study demonstrates the population expansion of lineage 2 and lineage 4 in eastern China, with comparable transmission capacity, while accumulation of resistance mutations does not necessarily facilitate the success of Mtb isolates. Compensatory mutations usually accompany drug resistance and significantly contribute to the epidemiological transmission of pre-XDR strains. Prospective molecular surveillance is required to further monitor the emergence and spread of pre-XDR/XDR strains in eastern China.

adaPop: Bayesian inference of dependent population dynamics in coalescent models

The coalescent is a powerful statistical framework that allows us to infer past population dynamics leveraging the ancestral relationships reconstructed from sampled molecular sequence data. In many biomedical applications, such as in the study of infectious diseases, cell development, and tumorgenesis, several distinct populations share evolutionary history and therefore become dependent. The inference of such dependence is a highly important, yet a challenging problem. With advances in sequencing technologies, we are well positioned to exploit the wealth of high-resolution biological data for tackling this problem. Here, we present adaPop, a probabilistic model to estimate past population dynamics of dependent populations and to quantify their degree of dependence. An essential feature of our approach is the ability to track the time-varying association between the populations while making minimal assumptions on their functional shapes via Markov random field priors. We provide nonparametric estimators, extensions of our base model that integrate multiple data sources, and fast scalable inference algorithms. We test our method using simulated data under various dependent population histories and demonstrate the utility of our model in shedding light on evolutionary histories of different variants of SARS-CoV-2.

Dates and Rates of Tick-Borne Encephalitis Virus-The Slowest Changing Tick-Borne Flavivirus

We evaluated the temporal signal and substitution rate of tick-borne encephalitis virus (TBEV) using 276 complete open reading frame (ORF) sequences with known collection dates. According to a permutation test, the TBEV Siberian subtype (TBEV-S) data set has no temporal structure and cannot be applied for substitution rate estimation without other TBEV subtypes. The substitution rate obtained suggests that the common clade of TBEV (TBEV-common), including all TBEV subtypes and louping-ill virus (LIV), is characterized by the lowest rate (1.87 × 10-5 substitutions per site per year (s/s/y) or 1 nucleotide substitution per ORF per 4.9 years; 95% highest posterior density (HPD) interval, 1.3-2.4 × 10-5 s/s/y) among all tick-borne flaviviruses previously assessed. Within TBEV-common, the TBEV European subtype (TBEV-E) has the lowest substitution rate (1.3 × 10-5 s/s/y or 1 nucleotide substitution per ORF per 7.5 years; 95% HPD, 1.0-1.8 × 10-5 s/s/y) as compared with TBEV Far-Eastern subtype (3.0 × 10-5 s/s/y or 1 nucleotide substitution per ORF per 3.2 years; 95% HPD, 1.6-4.5 × 10-5 s/s/y). TBEV-common representing the species tick-borne encephalitis virus diverged 9623 years ago (95% HPD interval, 6373-13,208 years). The TBEV Baikalian subtype is the youngest one (489 years; 95% HPD, 291-697 years) which differs significantly by age from TBEV-E (848 years; 95% HPD, 596-1112 years), LIV (2424 years; 95% HPD, 1572-3400 years), TBEV-FE (1936 years, 95% HPD, 1344-2598 years), and the joint clade of TBEV-S (2505 years, 95% HPD, 1700-3421 years) comprising Vasilchenko, Zausaev, and Baltic lineages.

Genetic evidence for widespread population size expansion in North American boreal birds prior to the Last Glacial Maximum

Pleistocene climate cycles are well documented to have shaped contemporary species distributions and genetic diversity. Northward range expansions in response to deglaciation following the Last Glacial Maximum (LGM; approximately 21 000 years ago) are surmised to have led to population size expansions in terrestrial taxa and changes in seasonal migratory behaviour. Recent findings, however, suggest that some northern temperate populations may have been more stable than expected through the LGM. We modelled the demographic history of 19 co-distributed boreal-breeding North American bird species from full mitochondrial gene sets and species-specific molecular rates. We used these demographic reconstructions to test how species with different migratory strategies were affected by glacial cycles. Our results suggest that effective population sizes increased in response to Pleistocene deglaciation earlier than the LGM, whereas genetic diversity was maintained throughout the LGM despite shifts in geographical range. We conclude that glacial cycles prior to the LGM have most strongly shaped contemporary genetic diversity in these species. We did not find a relationship between historic population dynamics and migratory strategy, contributing to growing evidence that major switches in migratory strategy during the LGM are unnecessary to explain contemporary migratory patterns.

Distinguishing imported cases from locally acquired cases within a geographically limited genomic sample of an infectious disease

MOTIVATION

The ability to distinguish imported cases from locally acquired cases has important consequences for the selection of public health control strategies. Genomic data can be useful for this, for example, using a phylogeographic analysis in which genomic data from multiple locations are compared to determine likely migration events between locations. However, these methods typically require good samples of genomes from all locations, which is rarely available.

RESULTS

Here, we propose an alternative approach that only uses genomic data from a location of interest. By comparing each new case with previous cases from the same location, we are able to detect imported cases, as they have a different genealogical distribution than that of locally acquired cases. We show that, when variations in the size of the local population are accounted for, our method has good sensitivity and excellent specificity for the detection of imports. We applied our method to data simulated under the structured coalescent model and demonstrate relatively good performance even when the local population has the same size as the external population. Finally, we applied our method to several recent genomic datasets from both bacterial and viral pathogens, and show that it can, in a matter of seconds or minutes, deliver important insights on the number of imports to a geographically limited sample of a pathogen population.

AVAILABILITY AND IMPLEMENTATION

The R package DetectImports is freely available from https://github.com/xavierdidelot/DetectImports.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

eSMC: a statistical model to infer admixture events from individual genomics data

Background

Inferring historical population admixture events yield essential insights in understanding a species demographic history. Methods are available to infer admixture events in demographic history with extant genetic data from multiple sources. Due to the deficiency in ancient population genetic data, there lacks a method for admixture inference from a single source. Pairwise Sequentially Markovian Coalescent (PSMC) estimates the historical effective population size from lineage genomes of a single individual, based on the distribution of the most recent common ancestor between the diploid’s alleles. However, PSMC does not infer the admixture event.

Results

Here, we proposed eSMC, an extended PSMC model for admixture inference from a single source. We evaluated our model’s performance on both in silico data and real data. We simulated population admixture events at an admixture time range from 5 kya to 100 kya (5 years/generation) with population admix ratio at 1:1, 2:1, 3:1, and 4:1, respectively. The root means the square error is $$\pm 7.61$$ ± 7.61 kya for all experiments. Then we implemented our method to infer the historical admixture events in human, donkey and goat populations. The estimated admixture time for both Han and Tibetan individuals range from 60 kya to 80 kya (25 years/generation), while the estimated admixture time for the domesticated donkeys and the goats ranged from 40 kya to 60 kya (8 years/generation) and 40 kya to 100 kya (6 years/generation), respectively. The estimated admixture times were concordance to the time that domestication occurred in human history.

Conclusion

Our eSMC effectively infers the time of the most recent admixture event in history from a single individual’s genomics data. The source code of eSMC is hosted at https://github.com/zachary-zzc/eSMC.

Consideration of genetic variation and evolutionary history in future conservation of Indian one-horned rhinoceros (Rhinoceros unicornis)

Background

The extant members of the Asian rhinos have experienced severe population and range declines since Pleistocene through a combination of natural and anthropogenic factors. The one-horned rhino is the only Asian species recovered from such conditions but most of the extant populations are reaching carrying capacity. India currently harbours ~ 83% of the global wild one-horned rhino populations distributed across seven protected areas. Recent assessments recommend reintroduction-based conservation approaches for the species, and implementation of such efforts would greatly benefit from detailed genetic assessments and evolutionary history of these populations. Using mitochondrial data, we investigated the phylogeography, divergence and demographic history of one-horned rhinos across its Indian range.

Results

We report the first complete mitogenome from all the extant Indian wild one-horned rhino populations (n = 16 individuals). Further, we identified all polymorphic sites and assessed rhino phylogeography (2531 bp mtDNA, n = 111 individuals) across India. Results showed 30 haplotypes distributed as three distinct genetic clades (F_st value 0.68–1) corresponding to the states of Assam (n = 28 haplotypes), West Bengal and Uttar Pradesh (both monomorphic). The reintroduced population of Uttar Pradesh showed maternal signatures of Chitwan National Park, Nepal. Mitochondrial phylogenomics suggests one-horned rhino diverged from its recent common ancestors ~ 950 Kya and different populations (Assam, West Bengal and Uttar Pradesh/Nepal) coalesce at ~ 190–50 Kya, corroborating with the paleobiogeography history of the Indian subcontinent. Further, the demography analyses indicated historical decline in female effective population size ~ 300–200 Kya followed by increasing trends during ~ 110–60 Kya.

Conclusion

The phylogeography and phylogenomic outcomes suggest recognition of three ‘Evolutionary Significant Units (ESUs)’ in Indian rhino. With ongoing genetic isolation of the current populations, future management efforts should focus on identifying genetically variable founder animals and consider periodic supplementation events while planning future rhino reintroduction programs in India. Such well-informed, multidisciplinary approach will be the only way to ensure evolutionary, ecological and demographic stability of the species across its range.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12862-022-02045-2.

A scalable analytical approach from bacterial genomes to epidemiology

Recent years have seen a remarkable increase in the practicality of sequencing whole genomes from large numbers of bacterial isolates. The availability of this data has huge potential to deliver new insights into the evolution and epidemiology of bacterial pathogens, but the scalability of the analytical methodology has been lagging behind that of the sequencing technology. Here we present a step-by-step approach for such large-scale genomic epidemiology analyses, from bacterial genomes to epidemiological interpretations. A central component of this approach is the dated phylogeny, which is a phylogenetic tree with branch lengths measured in units of time. The construction of dated phylogenies from bacterial genomic data needs to account for the disruptive effect of recombination on phylogenetic relationships, and we describe how this can be achieved. Dated phylogenies can then be used to perform fine-scale or large-scale epidemiological analyses, depending on the proportion of cases for which genomes are available. A key feature of this approach is computational scalability and in particular the ability to process hundreds or thousands of genomes within a matter of hours. This is a clear advantage of the step-by-step approach described here. We discuss other advantages and disadvantages of the approach, as well as potential improvements and avenues for future research. This article is part of a discussion meeting issue 'Genomic population structures of microbial pathogens'.

Demographic Expansions and the Emergence of Host Specialization in Genetically Distinct Ecotypes of the Tick-Transmitted Bacterium Anaplasma phagocytophilum

In Europe, genetically distinct ecotypes of the tick-vectored bacterium Anaplasma phagocytophilum circulate among mammals in three discrete enzootic cycles. To date, potential ecological factors that contributed to the emergence of these divergent ecotypes have been poorly studied. Here, we show that the ecotype that predominantly infects roe deer (Capreolus capreolus) is evolutionarily derived. Its divergence from a host generalist ancestor occurred after the last glacial maximum as mammal populations, including roe deer, recolonized the European mainland from southern refugia. We also provide evidence that this host specialist ecotype's effective population size (N_e) has tracked changes in the population of its roe deer host. Specifically, both host and bacterium have undergone substantial increases in N_e over the past 1,500 years. In contrast, we show that while it appears to have undergone a major population expansion starting ~3,500 years ago, in the past 500 years, the contemporary host generalist ecotype has experienced a substantial reduction in genetic diversity levels, possibly as a result of reduced opportunities for transmission between competent hosts. IMPORTANCE The findings of this study reveal specific events important for the evolution of host specialization in a naturally occurring, obligately intracellular bacterial pathogen. Specifically, they show that host range shifts and the emergence of host specialization may occur during periods of population growth in a generalist ancestor. Our results also demonstrate the close correlation between demographic patterns in host and pathogen for a specialist system. These findings have important relevance for understanding the evolution of host range diversity. They may inform future work on host range dynamics, and they provide insights for understanding the emergence of pathogens that have human and veterinary health implications.

Assessing model adequacy for Bayesian Skyline plots using posterior predictive simulation

Bayesian skyline plots (BSPs) are a useful tool for making inferences about demographic history. For example, researchers typically apply BSPs to test hypotheses regarding how climate changes have influenced intraspecific genetic diversity over time. Like any method, BSP has assumptions that may be violated in some empirical systems (e.g., the absence of population genetic structure), and the naïve analysis of data collected from these systems may lead to spurious results. To address these issues, we introduce P2C2M.Skyline, an R package designed to assess model adequacy for BSPs using posterior predictive simulation. P2C2M.Skyline uses a phylogenetic tree and the log file output from Bayesian Skyline analyses to simulate posterior predictive datasets and then compares this null distribution to statistics calculated from the empirical data to check for model violations. P2C2M.Skyline was able to correctly identify model violations when simulated datasets were generated assuming genetic structure, which is a clear violation of BSP model assumptions. Conversely, P2C2M.Skyline showed low rates of false positives when models were simulated under the BSP model. We also evaluate the P2C2M.Skyline performance in empirical systems, where we detected model violations when DNA sequences from multiple populations were lumped together. P2C2M.Skyline represents a user-friendly and computationally efficient resource for researchers aiming to make inferences from BSP.

Genomic epidemiology of Salmonella Typhi in Central Division, Fiji, 2012 to 2016

BACKGROUND

Typhoid fever is endemic in some Pacific Island Countries including Fiji and Samoa yet genomic surveillance is not routine in such settings. Previous studies suggested imports of the global H58 clade of Salmonella enterica var Typhi (Salmonella Typhi) contribute to disease in these countries which, given the MDR potential of H58, does not auger well for treatment. The objective of the study was to define the genomic epidemiology of Salmonella Typhi in Fiji.

METHODS

Genomic sequencing approaches were implemented to study the distribution of 255 Salmonella Typhi isolates from the Central Division of Fiji. We augmented epidemiological surveillance and Bayesian phylogenomic approaches with a multi-year typhoid case-control study to define geospatial patterns among typhoid cases.

FINDINGS

Genomic analyses showed Salmonella Typhi from Fiji resolved into 2 non-H58 genotypes with isolates from the two dominant ethnic groups, the Indigenous (iTaukei) and non-iTaukei genetically indistinguishable. Low rates of international importation of clones was observed and overall, there were very low levels an antibiotic resistance within the endemic Fijian typhoid genotypes. Genomic epidemiological investigations were able to identify previously unlinked case clusters. Bayesian phylodynamic analyses suggested that genomic variation within the larger endemic Salmonella Typhi genotype expanded at discreet times, then contracted.

INTERPRETATION

Cyclones and flooding drove 'waves' of typhoid outbreaks in Fiji which, through population aggregation, poor sanitation and water safety, and then mobility of the population, spread clones more widely. Minimal international importations of new typhoid clones suggest that targeted local intervention strategies may be useful in controlling endemic typhoid infection. These findings add to our understanding of typhoid transmission networks in an endemic island country with broad implications, particularly across Pacific Island Countries.

FUNDING

This work was supported by the Coalition Against Typhoid through the Bill and Melinda Gates Foundation [grant number OPP1017518], the Victorian Government, the National Health and Medical Research Council Australia, the Australian Research Council, and the Fiji Ministry of Health and Medical Services.

Measles Epidemiology and Viral Nucleoprotein Gene Evolution in Shandong Province, China

Measles, caused by measles virus (MeV), has not been eradicated in many regions and countries, threatening human health. Thus, it is beneficial for measles elimination to understand measles epidemiology and molecular evolution of key viral genes, such as nucleoprotein (N) gene. Based on public data, measles epidemiological information and MeV N gene sequences reported in Shandong Province, China were comprehensively collected and systematically analyzed. The results showed a positive correlation between population density and measles incidence (r = + 0.31), while negative correlations were found between measles incidence and healthcare condition (r = - 0.21) as well as average routine vaccination rate (r = - 0.11). Additionally, the predominant lineage of MeV in Shandong was formed by genotype H1 strains, and the time of the most recent common ancestor of the N gene of MeV genotype H1 in Shandong traced back to 1987 (95% highest posterior density, 1984-1990) with relatively rapid evolution (mean rate, 1.267×10^-3 substitutions/site/year). The genetic diversity of MeV N gene increased with the substantial emergence of major divergent clades of genotype H1 before 2005 and then remained relatively stable. In summary, these findings provided a significant insight into the measles elimination. This article is protected by copyright. All rights reserved.

The bounded coalescent model: Conditioning a genealogy on a minimum root date

The coalescent model represents how individuals sampled from a population may have originated from a last common ancestor. The bounded coalescent model is obtained by conditioning the coalescent model such that the last common ancestor must have existed after a certain date. This conditioned model arises in a variety of applications, such as speciation, horizontal gene transfer or transmission analysis, and yet the bounded coalescent model has not been previously analysed in detail. Here we describe a new algorithm to simulate from this model directly, without resorting to rejection sampling. We show that this direct simulation algorithm is more computationally efficient than the rejection sampling approach. We also show how to calculate the probability of the last common ancestor occurring after a given date, which is required to compute the probability density of realisations under the bounded coalescent model. Our results are applicable in both the isochronous (when all samples have the same date) and heterochronous (where samples can have different dates) settings. We explore the effect of setting a bound on the date of the last common ancestor, and show that it affects a number of properties of the resulting phylogenies. All our methods are implemented in a new R package called BoundedCoalescent which is freely available online.

Epidemiological inference from pathogen genomes: A review of phylodynamic models and applications

Phylodynamics requires an interdisciplinary understanding of phylogenetics, epidemiology, and statistical inference. It has also experienced more intense application than ever before amid the SARS-CoV-2 pandemic. In light of this, we present a review of phylodynamic models beginning with foundational models and assumptions. Our target audience is public health researchers, epidemiologists, and biologists seeking a working knowledge of the links between epidemiology, evolutionary models, and resulting epidemiological inference. We discuss the assumptions linking evolutionary models of pathogen population size to epidemiological models of the infected population size. We then describe statistical inference for phylodynamic models and list how output parameters can be rearranged for epidemiological interpretation. We go on to cover more sophisticated models and finish by highlighting future directions.

Cross-reactive immunity potentially drives global oscillation and opposed alternation patterns of seasonal influenza A viruses

Several human pathogens exhibit distinct patterns of seasonality and circulate as pairs. For instance, influenza A virus subtypes oscillate and peak during winter seasons of the world’s temperate climate zones. Alternation of dominant strains in successive influenza seasons makes epidemic forecasting a major challenge. From the start of the 2009 influenza pandemic we enrolled influenza A virus infected patients (n = 2980) in a global prospective clinical study. Complete hemagglutinin sequences were obtained from 1078 A/H1N1 and 1033 A/H3N2 viruses. We used phylodynamics to construct high resolution spatio-temporal phylogenetic hemagglutinin trees and estimated global influenza A effective reproductive numbers (R) over time (2009–2013). We demonstrate that R oscillates around R = 1 with a clear opposed alternation pattern between phases of the A/H1N1 and A/H3N2 subtypes. Moreover, we find a similar alternation pattern for the number of global viral spread between the sampled geographical locations. Both observations suggest a between-strain competition for susceptible hosts on a global level. Extrinsic factors that affect person-to-person transmission are a major driver of influenza seasonality. The data presented here indicate that cross-reactive host immunity is also a key intrinsic driver of influenza seasonality, which determines the influenza A virus strain at the onset of each epidemic season.

What is N e, anyway?

Few doubt that effective population size (Ne) is one of the most important parameters in evolutionary biology, but how many can say they really understand the concept? Ne is the evolutionary analogue of the number of individuals (or adults) in the population, N. Whereas ecological consequences of population size depend on N, evolutionary consequences (rates of loss of genetic diversity and increase in inbreeding; relative effectiveness of selection) depend on Ne. Formal definitions typically relate effective size to a key population genetic parameter, such as loss of heterozygosity or variance in allele frequency. However, for practical application to real populations, it is more useful to define Ne in terms of three demographic parameters: number of potential parents (adult N), and mean and variance in offspring number. Defined this way, Ne determines the rate of random genetic drift across the entire genome in the offspring generation. Other evolutionary forces (mutation, migration, selection)—together with factors such as variation in recombination rate—can also affect genetic variation, and this leads to heterogeneity across the genome in observed rates of genetic change. For some, it has been convenient to interpret this heterogeneity in terms of heterogeneity in Ne, but unfortunately this has muddled the concepts of genetic drift and effective population size. A commonly-repeated misconception is that Ne is the number of parents that actually contribute genes to the next generation (NP). In reality, NP can be smaller or larger than Ne, and the NP/Ne ratio depends on the sex ratio, the mean and variance in offspring number, and whether inbreeding or variance Ne is of interest.

Phylogeography of the parasitic mite Laelaps agilis in Western Palearctic shows lineages lacking host specificity but possessing different demographic histories

Background

Laelaps agilis C.L. Koch, 1836 is one the most abundant and widespread parasitic mite species in the Western Palearctic. It is a permanent ectoparasite associated with the Apodemus genus, which transmits Hepatozoon species via the host’s blood. Phylogenetic relationships, genealogy and host specificity of the mite are uncertain in the Western Palearctic. Here, we investigated the population genetic structure of 132 individual mites across Europe from their Apodemus and Clethrionomys hosts. Phylogenetic relationships and genetic variation of the populations were analyzed using cytochrome c oxidase subunit I (COI) gene sequences.

Results

We recovered three main mtDNA lineages within L. agilis in the Western Palearctic, which differentiated between 1.02 and 1.79 million years ago during the Pleistocene period: (i) Lineage A, including structured populations from Western Europe and the Czech Republic, (ii) Lineage B, which included only a few individuals from Greece and the Czech Republic; and (iii) Lineage C, which comprised admixed populations from Western and Eastern Europe. Contrary to their population genetic differentiation, the lineages did not show signs of specificity to different hosts. Finally, we confirmed that the sympatric congener L. clethrionomydis is represented by a separated monophyletic lineage.

Conclusion

Differences in the depth of population structure between L. agilis Lineages A and C, corroborated by the neutrality tests and demographic history analyses, suggested a stable population size in the structured Lineage A and a rapid range expansion for the geographically admixed Lineage C. We hypothesized that the two lineages were associated with hosts experiencing different glaciation histories. The lack of host specificity in L. agilis lineages was in contrast to the co-occurring highly host-specific lineages of Polyplax serrata lice, sharing Apodemus hosts. The incongruence was attributed to the differences in mobility between the parasites, allowing mites to switch hosts more often.

Variation in global distribution, population structures, and demographic history for four Trichiurus cutlassfishes

Background

Species-specific information on distribution and demographic patterns provides important implications for conservation and fisheries management. However, such information is often lacking for morphologically-similar species, which may lead to biases in the assessments of these species and even decrease effort towards sustainable management. Here, we aimed to uncover the distribution range, population structure and demographic history for four exploited Trichiurus cutlassfishes using genetics. These cutlassfishes contribute substantial global fisheries catch, with a high proportion of catch harvested from the NW Pacific.

Methods

We chose the widely available mitochondrial 16S ribosomal RNA (16S) as the genetic marker for cutlassfishes. We compiled the 16S sequence data from both the GenBank and a survey of trawler catch samples along the NW Pacific coasts 22-39°N. Genealogical relationships within each species was visualized with haplotype networks and potential population differentiations were further evaluated with AMOVA. Demographic histories were estimated using neutrality test, mismatch analysis, and the Bayesian skyline plot. The reconstructed phylogenetic trees were used to delimit and estimate the divergence time of species and included populations.

Results

In each of two cosmopolitan species, T. lepturus and T. nanhaiensis, we observed distinct populations along the coasts of warm oceans; such population differentiation might result from historical geographic barriers in the Pleistocene. In the NW Pacific, four Trichiurus species vary in their distribution habitats, which reflect differential ecological niches among these species. The small-sized T. brevis are primarily found in nearshore habitats; the warm-affiliated T. nanhaiensis are present along the path of the Kuroshio Current; the cold-affiliated T. japonicus spatially diverged from the widely-distributed T. lepturus, with the latter mainly occupy in warmer regions. Despite these differences, a single well-mixing fish stock, thus one management unit, was identified in each of the four species, presumably due to expansion of their population sizes predated the Last Glacial Maximum and a lack of distribution barrier. The most dominant T. japonicus, which have at least one magnitude higher effective population size than the others, show a unique abrupt size expansion event at 75 to 50-kilo years ago when the low sea level occurred during the ice age.

Main conclusions

The demographic history revealed by our genetic analyses advances understanding of the current distribution and population structure for these congeneric species. Moreover, the uncovered population structure provides insight into the assessment and management of these species. Such information complements contemporary knowledge about these species and enables us to forecast their ability to resist future environmental and anthropogenic disturbances.

Bayesian inference of clonal expansions in a dated phylogeny

Microbial population genetics models often assume that all lineages are constrained by the same population size dynamics over time. However, many neutral and selective events can invalidate this assumption and can contribute to the clonal expansion of a specific lineage relative to the rest of the population. Such differential phylodynamic properties between lineages result in asymmetries and imbalances in phylogenetic trees that are sometimes described informally but which are difficult to analyze formally. To this end, we developed a model of how clonal expansions occur and affect the branching patterns of a phylogeny. We show how the parameters of this model can be inferred from a given dated phylogeny using Bayesian statistics, which allows us to assess the probability that one or more clonal expansion events occurred. For each putative clonal expansion event, we estimate its date of emergence and subsequent phylodynamic trajectory, including its long-term evolutionary potential which is important to determine how much effort should be placed on specific control measures. We demonstrate the applicability of our methodology on simulated and real data sets. Inference under our clonal expansion model can reveal important features in the evolution and epidemiology of infectious disease pathogens. [Clonal expansion; genomic epidemiology; microbial population genomics; phylodynamics.]

Reductive evolution of virulence repertoire to drive the divergence between community- and hospital-associated methicillin-resistant Staphylococcus aureus of the ST1 lineage

Methicillin-resistant Staphylococcus aureus (MRSA) of the ST1-SCCmecIV lineage has been associated with community-acquired (CA) infections in North America and Australia. In Brazil, multi-drug resistant ST1-SCCmecIV MRSA has emerged in hospital-associated (HA) diseases in Rio de Janeiro. To understand these epidemiological differences, genomic and phylogenetic analyses were performed. In addition, virulence assays were done for representative CA - and HA-MRSA strains. Despite the conservation of the virulence repertoire, some genes were missing in Brazilian ST1-SCCmecIV including lukSF-PV, fnbB, and several superantigen-encoded genes. Additionally, CA-MRSA lost the splDE while HA-MRSA strains conserved the complete operon. Most of these variable genes were located in mobile genetic elements (MGE). However, conservation and maintenance of MGEs were often observed despite the absence of their associated virulence markers. A Bayesian phylogenetic tree revealed the occurrence of more than one entrance of ST1 strains in Rio de Janeiro. The tree shape and chronology allowed us to infer that the hospital-associated ST1-SCCmecIV from Brazil and the community-acquired USA400 from North America are not closely related and that they might have originated from different MSSA strains that independently acquired SCCmecIV cassettes. As expected, representatives of ST1 strains from Brazil showed lower cytotoxicity and a greater ability to survive inside human host cells. We suggest that Brazilian ST1-SCCmecIV strains have adapted to the hospital setting by reducing virulence and gaining the ability to persist and survive inside host cells. Possibly, these evolutionary strategies may balance the biologic cost of retaining multiple antibiotic resistance genes.

Biogeographic and demographic history of the Mediterranean snakes Malpolon monspessulanus and Hemorrhois hippocrepis across the Strait of Gibraltar

BACKGROUND

The contribution of North Africa to the assembly of biodiversity within the Western Palaearctic is still poorly documented. Since the Miocene, multiple biotic exchanges occurred across the Strait of Gibraltar, underlying the high biogeographic affinity between the western European and African sides of the Mediterranean basin. We investigated the biogeographic and demographic dynamics of two large Mediterranean-adapted snakes across the Strait and assess their relevance to the origin and diversity patterns of current European and North African populations.

RESULTS

We inferred phylogeographic patterns and demographic history of M. monspessulanus and H. hippocrepis, based on range-wide multilocus data, combined with fossil data and species distribution modelling, under present and past bioclimatic envelopes. For both species we identified endemic lineages in the High Atlas Mountains (Morocco) and in eastern Iberia, suggesting their persistence in Europe during the Pleistocene. One lineage is shared between North Africa and southern Iberia and likely spread from the former to the latter during the sea-level low stand of the last glacial stage. During this period M. monspessulanus shows a sudden demographic expansion, associated with increased habitat suitability in North Africa. Lower habitat suitability is predicted for both species during interglacial stages, with suitable areas restricted to coastal and mountain ranges of Iberia and Morocco. Compiled fossil data for M. monspessulanus show a continuous fossil record in Iberia at least since the Pliocene and throughout the Pleistocene.

CONCLUSIONS

The previously proposed hypothesis of Pleistocene glacial extinction of both species in Europe is not supported based on genetic data, bioclimatic envelopes models, and the available fossil record. A model of range retraction to mountain refugia during arid periods and of glacial expansion (demographic and spatial) associated to an increase of Mediterranean habitats during glacial epochs emerges as a general pattern for mesic vertebrates in North Africa. Moreover, the phylogeographic pattern of H. hippocrepis conforms to a well-established biogeographic partition between western and eastern Maghreb.