Prokaryotic genomes are usually densely packed with intact and functional genes. However, in certain contexts, such as after recent ecological shifts or extreme population bottlenecks, broken and nonfunctional gene fragments can quickly accumulate and form a substantial fraction of the genome. Identification of these broken genes, called pseudogenes, is a critical step for understanding the evolutionary forces acting upon, and the functional potential encoded within, prokaryotic genomes. Here, we present Pseudofinder, an open-source software dedicated to pseudogene identification and analysis in bacterial and archaeal genomes. We demonstrate that Pseudofinder's multi-pronged, reference-based approach can detect a wide variety of pseudogenes, including those that are highly degraded and typically missed by gene-calling pipelines, as well newly formed pseudogenes containing only one or a few inactivating mutations. Additionally, Pseudofinder can detect genes that lack inactivating substitutions but experiencing relaxed selection. Implementation of Pseudofinder in annotation pipelines will allow more precise estimations of the functional potential of sequenced microbes, while also generating new hypotheses related to the evolutionary dynamics of bacterial and archaeal genomes.

The selection of the best-fitting substitution model of molecular evolution is a traditional step for phylogenetic inferences, including ancestral sequence reconstruction (ASR). However, a few recent studies suggested that applying this procedure does not affect the accuracy of phylogenetic tree reconstruction. Here, we revisited this debate topic by analyzing the influence of selection among substitution models of protein evolution, with focus on exchangeability matrices, on the accuracy of ASR using simulated and real data. We found that the selected best-fitting substitution model produces the most accurate ancestral sequences, especially if the data present large genetic diversity. Indeed, ancestral sequences reconstructed under substitution models with similar exchangeability matrices were similar, suggesting that if the selected best-fitting model cannot be used for the reconstruction, applying a model similar to the selected one is preferred. We conclude that selecting among substitution models of protein evolution is recommended for reconstructing accurate ancestral sequences.

Alternative splicing underpins functional diversity in proteins and the complexity and diversity of eukaryotes. An example is the doublesex gene, the key transcriptional factor in arthropod sexual differentiation. doublesex is controlled by sex-specific splicing and promotes both male and female differentiation in holometabolan insects, whereas in hemimetabolan species, doublesex has sex-specific isoforms but is not required for female differentiation. How doublesex evolved to be essential for female development remains largely unknown. Here, we investigate ancestral states of doublesex using Thermobia domestica belonging to Zygentoma, the sister group of Pterygota, that is, winged insects. We find that, in T. domestica, doublesex expresses sex-specific isoforms but is only necessary for male differentiation of sexual morphology. This result supports the hypothesis that doublesex initially promoted male differentiation during insect evolution. However, T. domestica doublesex has a short female-specific region and upregulates the expression of vitellogenin homologs in females, suggesting that doublesex may already play some role in female morphogenesis of the common ancestor of Pterygota. Reconstruction of the ancestral sequence and prediction of protein structures show that the female-specific isoform of doublesex has an extended C-terminal disordered region in holometabolan insects but not in nonholometabolan species. We propose that doublesex acquired its function in female morphogenesis through a change in the protein motif structure rather than the emergence of the female-specific exon.

The divergence of regulatory regions and gene regulatory network (GRN) rewiring is a key driver of cichlid phenotypic diversity. However, the contribution of miRNA-binding site turnover has yet to be linked to GRN evolution across cichlids. Here, we extend our previous studies by analyzing the selective constraints driving evolution of miRNA and transcription factor (TF)-binding sites of target genes, to infer instances of cichlid GRN rewiring associated with regulatory binding site turnover. Comparative analyses identified increased species-specific networks that are functionally associated to traits of cichlid phenotypic diversity. The evolutionary rewiring is associated with differential models of miRNA- and TF-binding site turnover, driven by a high proportion of fast-evolving polymorphic sites in adaptive trait genes compared with subsets of random genes. Positive selection acting upon discrete mutations in these regulatory regions is likely to be an important mechanism in rewiring GRNs in rapidly radiating cichlids. Regulatory variants of functionally associated miRNA- and TF-binding sites of visual opsin genes differentially segregate according to phylogeny and ecology of Lake Malawi species, identifying both rewired, for example, clade-specific and conserved network motifs of adaptive trait associated GRNs. Our approach revealed several novel candidate regulators, regulatory regions, and three-node motifs across cichlid genomes with previously reported associations to known adaptive evolutionary traits.

Transcription regulatory sequences (TRSs), which occur upstream of structural and accessory genes as well as the 5' end of a coronavirus genome, play a critical role in discontinuous transcription in coronaviruses. We introduce two problems collectively aimed at identifying these regulatory sequences as well as their associated genes. First, we formulate the TRS Identification problem of identifying TRS sites in a coronavirus genome sequence with prescribed gene locations. We introduce CORSID-A, an algorithm that solves this problem to optimality in polynomial time. We demonstrate that CORSID-A outperforms existing motif-based methods in identifying TRS sites in coronaviruses. Second, we demonstrate for the first time how TRS sites can be leveraged to identify gene locations in the coronavirus genome. To that end, we formulate the TRS and Gene Identification problem of simultaneously identifying TRS sites and gene locations in unannotated coronavirus genomes. We introduce CORSID to solve this problem, which includes a web-based visualization tool to explore the space of near-optimal solutions. We show that CORSID outperforms state-of-the-art gene finding methods in coronavirus genomes. Furthermore, we demonstrate that CORSID enables de novo identification of TRS sites and genes in previously unannotated coronavirus genomes. CORSID is the first method to perform accurate and simultaneous identification of TRS sites and genes in coronavirus genomes without the use of any prior information.

Primate herpes simplex viruses are species-specific and relatively harmless to their natural hosts. However, cross-species transmission is often associated with severe disease, as exemplified by the virulence of macacine herpesvirus 1 (B virus) in humans. We performed a genome-wide scan for signals of adaptation of simplexviruses to their hominin hosts. Among core genes, we found evidence of episodic positive selection in three glycoproteins, with several selected sites located in antigenic determinants. Positively selected noncore genes were found to be involved in different immune-escape mechanisms. The herpes simplex virus (HSV)-1/HSV-2 encoded product (ICP47) of one of these genes is known to down-modulate major histocompatibility complex class I expression. This feature is not shared with B virus, which instead up-regulates Human Leukocyte Antigen (HLA)-G, an immunomodulatory molecule. By in vitro expression of different ICP47 mutants, we functionally characterized the selection signals. Results indicated that the selected sites do not represent the sole determinants of binding to the transporter associated with antigen processing (TAP). Conversely, the amino acid status at these sites was sufficient to determine HLA-G up-regulation. In fact, both HSV-1 and HSV-2 ICP47 induced HLA-G when mutated to recapitulate residues in B virus, whereas the mutated version of B virus ICP47 failed to determine HLA-G expression. These differences might contribute to the severity of B virus infection in humans. Importantly, they indicate that the evolution of ICP47 in HSV-1/HSV-2 led to the loss of an immunosuppressive effect. Thus, related simplexviruses finely tune the balance between immunosuppressive and immunostimulatory pathways to promote successful co-existence with their primate hosts.

The popular trap model assumes that the invasions of transposable elements (TEs) in mammals and invertebrates are stopped by piRNAs that emerge after insertion of the TE into a piRNA cluster. It remains, however, still unclear which factors influence the dynamics of TE invasions. The activity of the TE (i.e., transposition rate) is one frequently discussed key factor. Here we take advantage of the temperature-dependent activity of the P-element, a widely studied eukaryotic TE, to test how TE activity affects the dynamics of a TE invasion. We monitored P-element invasion dynamics in experimental Drosophila simulans populations at hot and cold culture conditions. Despite marked differences in transposition rates, the P-element reached very similar copy numbers at both temperatures. The reduction of the insertion rate upon approaching the copy number plateau was accompanied by similar amounts of piRNAs against the P-element at both temperatures. Nevertheless, we also observed fewer P-element insertions in piRNA clusters than expected, which is not compatible with a simple trap model. The ping-pong cycle, which degrades TE transcripts, becomes typically active after the copy number plateaued. We generated a model, with few parameters, that largely captures the observed invasion dynamics. We conclude that the transposition rate has at the most only a minor influence on TE abundance, but other factors, such as paramutations or selection against TE insertions are shaping the TE composition.

In Drosophila melanogaster, a key germline stem cell (GSC) differentiation factor, bag of marbles (bam) shows rapid bursts of amino acid fixations between sibling species D. melanogaster and Drosophila simulans, but not in the outgroup species Drosophila ananassae. Here, we test the null hypothesis that bam's differentiation function is conserved between D. melanogaster and four additional Drosophila species in the melanogaster species group spanning approximately 30 million years of divergence. Surprisingly, we demonstrate that bam is not necessary for oogenesis or spermatogenesis in Drosophila teissieri nor is bam necessary for spermatogenesis in D. ananassae. Remarkably bam function may change on a relatively short time scale. We further report tests of neutral sequence evolution at bam in additional species of Drosophila and find a positive, but not perfect, correlation between evidence for positive selection at bam and its essential role in GSC regulation and fertility for both males and females. Further characterization of bam function in more divergent lineages will be necessary to distinguish between bam's critical gametogenesis role being newly derived in D. melanogaster, D. simulans, Drosophila yakuba, and D. ananassae females or it being basal to the genus and subsequently lost in numerous lineages.

African wild pigs have a contentious evolutionary and biogeographic history. Until recently, desert warthog (Phacochoerus aethiopicus) and common warthog (P. africanus) were considered a single species. Molecular evidence surprisingly suggested they diverged at least 4.4 million years ago, and possibly outside of Africa. We sequenced the first whole-genomes of four desert warthogs and 35 common warthogs from throughout their range. We show that these two species diverged much later than previously estimated, 400,000-1,700,000 years ago depending on assumptions of gene flow. This brings it into agreement with the paleontological record. We found that the common warthog originated in western Africa and subsequently colonized eastern and southern Africa. During this range expansion, the common warthog interbred with the desert warthog, presumably in eastern Africa, underlining this region's importance in African biogeography. We found that immune system-related genes may have adaptively introgressed into common warthogs, indicating that resistance to novel diseases was one of the most potent drivers of evolution as common warthogs expanded their range. Hence, we solve some of the key controversies surrounding warthog evolution and reveal a complex evolutionary history involving range expansion, introgression, and adaptation to new diseases.

The p53 tumor suppressor is a transcription factor with roles in cell development, apoptosis, oncogenesis, aging, and homeostasis in response to stresses and infections. p53 is tightly regulated by the MDM2 E3 ubiquitin ligase. The p53-MDM2 pathway has coevolved, with MDM2 remaining largely conserved, whereas the TP53 gene morphed into various isoforms. Studies on prevertebrate ancestral homologs revealed the transition from an environmentally induced mechanism activating p53 to a tightly regulated system involving cell signaling. The evolution of this mechanism depends on structural changes in the interacting protein motifs. Elephants such as Loxodonta africana constitute ideal models to investigate this coevolution as they are large and long-living as well as having 20 copies of TP53 isoformic sequences expressing a variety of BOX-I MDM2-binding motifs. Collectively, these isoforms would enhance sensitivity to cellular stresses, such as DNA damage, presumably accounting for strong cancer defenses and other adaptations favoring healthy aging. Here we investigate the molecular evolution of the p53-MDM2 system by combining in silico modeling and in vitro assays to explore structural and functional aspects of p53 isoforms retaining the MDM2 interaction, whereas forming distinct pools of cell signaling. The methodology used demonstrates, for the first time that in silico docking simulations can be used to explore functional aspects of elephant p53 isoforms. Our observations elucidate structural and mechanistic aspects of p53 regulation, facilitate understanding of complex cell signaling, and suggest testable hypotheses of p53 evolution referencing Peto's Paradox.

Hybridization between lineages that have not reached complete reproductive isolation appears more and more like a common phenomenon. Indeed, speciation genomic studies have now extensively shown that many species' genomes have hybrid ancestry. However, genomic patterns of introgression are often heterogeneous across the genome. In many organisms, a positive correlation between introgression levels and recombination rate has been observed. It is usually explained by the purging of deleterious introgressed material due to incompatibilities. However, the opposite relationship was observed in a North American population of Drosophila melanogaster with admixed European and African ancestry. In order to explore how directional and epistatic selection can impact the relationship between introgression and recombination, we performed forward simulations of whole D. melanogaster genomes reflecting the North American population's history. Our results revealed that the simplest models of positive selection often yield negative correlations between introgression and recombination such as the one observed in D. melanogaster. We also confirmed that incompatibilities tend to produce positive introgression-recombination correlations. And yet, we identify parameter space under each model where the predicted correlation is reversed. These findings deepen our understanding of the evolutionary forces that may shape patterns of ancestry across genomes, and they strengthen the foundation for future studies aimed at estimating genome-wide parameters of selection in admixed populations.

Multiheme cytochromes play key roles in diverse biogeochemical cycles, but understanding the origin and evolution of these proteins is a challenge due to their ancient origin and complex structure. Up until now, the evolution of multiheme cytochromes composed by multiple redox modules in a single polypeptide chain was proposed to occur by gene fusion events. In this context, the pentaheme nitrite reductase NrfA and the tetraheme cytochrome c554 were previously proposed to be at the origin of the extant octa- and nonaheme cytochrome c involved in metabolic pathways that contribute to the nitrogen, sulfur, and iron biogeochemical cycles by a gene fusion event. Here, we combine structural and character-based phylogenetic analysis with an unbiased root placement method to refine the evolutionary relationships between these multiheme cytochromes. The evidence show that NrfA and cytochrome c554 belong to different clades, which suggests that these two multiheme cytochromes are products of truncation of ancestral octaheme cytochromes related to extant octaheme nitrite reductase and MccA, respectively. From our phylogenetic analysis, the last common ancestor is predicted to be an octaheme cytochrome with nitrite reduction ability. Evolution from this octaheme framework led to the great diversity of extant multiheme cytochromes analyzed here by pruning and grafting of protein modules and hemes. By shedding light into the evolution of multiheme cytochromes that intervene in different biogeochemical cycles, this work contributes to our understanding about the interplay between biology and geochemistry across large time scales in the history of Earth.

Sublineages (SLs) within microbial species can differ widely in their ecology and pathogenicity, and their precise definition is important in basic research and for industrial or public health applications. Widely accepted strategies to define SLs are currently missing, which confuses communication in population biology and epidemiological surveillance. Here, we propose a broadly applicable genomic classification and nomenclature approach for bacterial strains, using the prominent public health threat Klebsiella pneumoniae as a model. Based on a 629-gene core genome multilocus sequence typing (cgMLST) scheme, we devised a dual barcoding system that combines multilevel single linkage (MLSL) clustering and life identification numbers (LINs). Phylogenetic and clustering analyses of >7,000 genome sequences captured population structure discontinuities, which were used to guide the definition of 10 infraspecific genetic dissimilarity thresholds. The widely used 7-gene multilocus sequence typing (MLST) nomenclature was mapped onto MLSL SLs (threshold: 190 allelic mismatches) and clonal group (threshold: 43) identifiers for backwards nomenclature compatibility. The taxonomy is publicly accessible through a community-curated platform (https://bigsdb.pasteur.fr/klebsiella), which also enables external users' genomic sequences identification. The proposed strain taxonomy combines two phylogenetically informative barcode systems that provide full stability (LIN codes) and nomenclatural continuity with previous nomenclature (MLSL). This species-specific dual barcoding strategy for the genomic taxonomy of microbial strains is broadly applicable and should contribute to unify global and cross-sector collaborative knowledge on the emergence and microevolution of bacterial pathogens.

Due to extensive pleiotropy, trans-acting elements are often thought to be evolutionarily constrained. While the impact of trans-acting elements on gene expression evolution has been extensively studied, relatively little is understood about the contribution of a single trans regulator to interspecific expression and phenotypic divergence. Here, we disentangle the effects of genomic context and miR-983, an adaptively evolving young microRNA, on expression divergence between Drosophila melanogaster and D. simulans. We show miR-983 effects promote interspecific expression divergence in testis despite its antagonism with the often-predominant context effects. Single-cyst RNA-seq reveals that distinct sets of genes gain and lose miR-983 influence under disruptive or diversifying selection at different stages of spermatogenesis, potentially helping minimize antagonistic pleiotropy. At the round spermatid stage, the effects of miR-983 are weak and distributed, coincident with the transcriptome undergoing drastic expression changes. Knocking out miR-983 causes reduced sperm length with increased within-individual variation in D. melanogaster but not in D. simulans, and the D. melanogaster knockout also exhibits compromised sperm defense ability. Our results provide empirical evidence for the resolution of antagonistic pleiotropy and also have broad implications for the function and evolution of new trans regulators.

Opossums in the tribe Didelphini are resistant to pit viper venoms and are hypothesized to be coevolving with venomous snakes. Specifically, a protein involved in blood clotting (von Willebrand factor [vWF] which is targeted by snake venom C-type lectins [CTLs]) has been found to undergo rapid adaptive evolution in Didelphini. Several unique amino acid changes in vWF could explain their resistance; however, experimental evidence that these changes disrupt binding to venom CTLs was lacking. Furthermore, without explicit testing of ancestral phenotypes to reveal the mode of evolution, the assertion that this system represents an example of coevolution rather than noncoevolutionary adaptation remains unsupported. Using expressed vWF proteins and purified venom CTLs, we quantified binding affinity for vWF proteins from all resistant taxa, their venom-sensitive relatives, and their ancestors. We show that CTL-resistant vWF is present in opossums outside clade Didelphini and likely across a wider swath of opossums (family Didelphidae) than previously thought. Ancestral reconstruction and in vitro testing of vWF phenotypes in a clade of rapidly evolving opossums reveal a pattern consistent with trench warfare coevolution between opossums and their venomous snake prey.

Several types of shell matrix proteins (SMPs) have been identified in molluskan shells. Their diversity is the consequence of various molecular processes, including domain shuffling and gene duplication. However, the evolutionary origin of most SMPs remains unclear. In this study, we investigated the evolutionary process EGF-like and zona pellucida (ZP) domains containing SMPs. Two types of the proteins (EGF-like protein (EGFL) and EGF-like and ZP domains containing protein (EGFZP)) were found in the pearl oyster, Pinctada fucata. In contrast, only EGFZP was identified in the gastropods. Phylogenetic analysis and genomic arrangement studies showed that EGFL and EGFZP formed a clade in bivalves, and their encoding genes were localized in tandem repeats on the same scaffold. In P. fucata, EGFL genes were expressed in the outer part of mantle epithelial cells are related to the calcitic shell formation. However, in both P. fucata and the limpet Nipponacmea fuscoviridis, EGFZP genes were expressed in the inner part of the mantle epithelial cells are related to aragonitic shell formation. Furthermore, our analysis showed that in P. fucata, the ZP domain interacts with eight SMPs that have various functions in the nacreous shell mineralization. The data suggest that the ZP domain can interact with other SMPs, and EGFL evolution in pterimorph bivalves represents an example of neo-functionalization that involves the acquisition of a novel protein through gene duplication.

The mutation rate is a fundamental evolutionary parameter with direct and appreciable effects on the health and function of individuals. Here, we examine this important parameter in the domestic cat, a beloved companion animal as well as a valuable biomedical model. We estimate a mutation rate of 0.86 × 10-8 per bp per generation for the domestic cat (at an average parental age of 3.8 years). We find evidence for a significant paternal age effect, with more mutations transmitted by older sires. Our analyses suggest that the cat and the human have accrued similar numbers of mutations in the germline before reaching sexual maturity. The per-generation mutation rate in the cat is 28% lower than what has been observed in humans, but is consistent with the shorter generation time in the cat. Using a model of reproductive longevity, which takes into account differences in the reproductive age and time to sexual maturity, we are able to explain much of the difference in per-generation rates between species. We further apply our reproductive longevity model in a novel analysis of mutation spectra and find that the spectrum for the cat resembles the human mutation spectrum at a younger age of reproduction. Together, these results implicate changes in life-history as a driver of mutation rate evolution between species. As the first direct observation of the paternal age effect outside of rodents and primates, our results also suggest a phenomenon that may be universal among mammals.

Most vertebrate sex-determining genes (SDGs) emerge as neofunctionalized genes through duplication and/or mutation of ancestral genes that are involved with sexual differentiation. We previously demonstrated dm-W to be the SDG in the African clawed frog Xenopus laevis and found that a portion of this gene emerged from the masculinization gene dmrt1 after allotetraploidization by interspecific hybridization between two ancestral species around 17-18 Ma. dm-W has four exons consisting of a noncoding exon 1, dmrt1-derived exons 2 and 3, and an orphan exon 4 (Ex4) of unknown origin that includes coding sequence (CDS). In this study, we searched for the origin of Ex4 and investigated the function of the CDS of this exon. We found that the Ex4-CDS is derived from a noncoding portion of the hAT-10 family of DNA transposon. Evolutionary analysis of transposons and determination of the Ex4 sequences from three other species indicated that Ex4 was generated before the diversification of most or all extant allotetraploid species in subgenus Xenopus, during which time we hypothesize that transposase activity of this hAT superfamily was active. Using DNA-protein binding and transfection assays, we further demonstrate that the Ex4-encoded amino acid sequence increases the DNA-binding ability and transrepression activity of DM-W. These findings suggest that the conversion of the noncoding transposon sequence to the CDS of dm-W contributed to neofunctionalization of a new chimeric SDG in the ancestor of the allotetraploid Xenopus species, offering new insights into de novo origin and functional evolution of chimerical genes.

Understanding intratumor heterogeneity is critical for studying tumorigenesis and designing personalized treatments. To decompose the mixed cell population in a tumor, subclones are inferred computationally based on variant allele frequency (VAF) from bulk sequencing data. In this study, we showed that sequencing depth, mean VAF, and variance of VAF of a subclone are confounded. Without considering this effect, current methods require deep-sequencing data (>300× depth) to reliably infer subclones. Here, we present a novel algorithm that incorporates depth-variance and mean-variance dependencies in a clustering error model and successfully identifies subclones in tumors sequenced at depths of as low as 30×. We implemented the algorithm as a model-based adaptive grouping of subclones (MAGOS) method. Analyses of computer simulated data and empirical sequencing data showed that MAGOS outperformed existing methods on minimum sequencing depth, decomposition accuracy, and computation efficiency. The most prominent improvements were observed in analyzing tumors sequenced at depths between 30× and 200×, whereas the performance was comparable between MAGOS and existing methods on deeply sequenced tumors. MAGOS supports analysis of single-nucleotide variants and copy number variants from a single sample or multiple samples of a tumor. We applied MAGOS to whole-exome data of late-stage liver cancers and discovered that high subclone count in a tumor was a significant risk factor of poor prognosis. Lastly, our analysis suggested that sequencing multiple samples of the same tumor at standard depth is more cost-effective and robust for subclone characterization than deep sequencing a single sample. MAGOS is available at github (https://github.com/liliulab/magos).

Genetic recombination plays a critical role in the emergence of pathogens with phenotypes such as drug resistance, virulence, and host adaptation. Here, we tested the hypothesis that recombination between sympatric ancestral populations leads to the emergence of divergent variants of the zoonotic parasite Cryptosporidium parvum with modified host ranges. Comparative genomic analyses of 101 isolates have identified seven subpopulations isolated by distance. They appear to be descendants of two ancestral populations, IIa in northwestern Europe and IId from southwestern Asia. Sympatric recombination in areas with both ancestral subtypes and subsequent selective sweeps have led to the emergence of new subpopulations with mosaic genomes and modified host preference. Subtelomeric genes could be involved in the adaptive selection of subpopulations, while copy number variations of genes encoding invasion-associated proteins are potentially associated with modified host ranges. These observations reveal ancestral origins of zoonotic C. parvum and suggest that pathogen import through modern animal farming might promote the emergence of divergent subpopulations of C. parvum with modified host preference.