Comparative genomics of the parasite Trichomonas vaginalis reveals genes involved in spillover from birds to humans

Abstract

Trichomonas vaginalis, the causative agent of the venereal disease trichomoniasis, infects men and women globally and is associated with serious outcomes during pregnancy, increased risk of HIV-1 infection, and cancers of the human reproductive tract. Species of trichomonad parasitize a range of hosts in addition to humans, including birds, livestock, and pets. Genetic analysis of trichomonads recovered from columbid birds has provided evidence that they undergo frequent host-switching, and that a spillover event from columbids likely gave rise to T. vaginalis in humans. Here we describe a comparative genomics study of seven trichomonad species, generating chromosome-scale reference genomes for T. vaginalis and its avian sister species Trichomonas stableri, and assemblies of five other species that infect birds and mammals. Human-infecting trichomonad lineages have undergone recent and convergent genome size expansions compared to their avian sister species, a result of extensive repeat expansions specifically of multicopy gene families and transposable elements, with genetic drift likely a driver due to relaxed selection. Trichomonads are thought to have independently host-switched twice from birds to mammals/humans. We identify gene functions implicated in the transition, including host tissue adherence and phagocytosis, extracellular vesicle formation, and CAZyme virulence factors, which are all associated with pathogenesis phenotypes.

Fig. 1. Architecture and genome features of T. vaginalis G3 across its six chromosomes.

The concentric rings, from innermost to outermost, represent: (1) chromosome size in Mb; (2) gene density (green plot) shown in 20 Kb windows; vertical blue lines represent 11 rRNA cassettes, and the vertical black line represents the 47.5 Kb block from an LGT event of the bacterium Peptoniphilus harei; (3) TE density (pink plot) shown in 20 Kb windows; (4) transcript abundance (brown plot) of all genes shown as transcripts per million (TPM) in 100 Kb windows; (5) TE transcript abundance (orange plot) of annotated TE genes shown as TPM in 100 Kb windows; and (6) dN/dS values (grey dots). The axes are shown next to chromosome I.

Fig. 2. Synteny plot of human parasite T. vaginalis and its closest relative in birds T. stableri.

Each of the six T. vaginalis chromosomes I–VI are colored uniquely, and synteny blocks between the two species are indicated by ribbons connecting the chromosomes. Chromosomes are not shown in numbered order for visualization purposes. Blue graph plots show normalized TE density (genomic sequence classified as containing TEs) in 100 Kb windows on the top and bottom of each species’ chromosomes. Hi-C interaction maps are shown as red triangles above (T. vaginalis) and below (T. stableri) the TE density plots, with lengths of contigs from the assemblies shown as white squares for T. vaginalis and green squares for T. stableri.

Fig. 3. Genome content distribution.

A An Upset plot displaying the intersection of orthogroups identified in OrthoFinder across seven trichomonad species. Each vertical bar represents the number of orthogroups shared at each species intersection, the set size indicates the number of orthogroups found in each species, and the connected dots represent the species in the intersection. B Ultrametric tree from 6226 concatenated single-copy genes. Black dots at terminal nodes are proportional to estimated genome size, and hosts are denoted by cartoons. Estimated gene family expansions/contractions from 12,345 genes are denoted as + or − values on the tree. The heat map shows the log10 transformed count of TE family members for each tree branch.

Fig. 4. Analysis of orthologs across seven trichomonad species.

A Graphs showing count of all single-copy orthologs (SCOs; left panel) and BUSCO genes (right panel) identified by RELAX as being under relaxed, neutral, or intensified selection in the species with expanded genomes (T. vaginalis and T. tenax) for a range of P-values. B Left panel: mean dN/dS values (plotted from 0.0 to 1.0) for SCOs under relaxed purifying selection for bird-infecting species (n = 141) and mammal-infecting species (n = 31), and right panel: dN/dS values (plotted from 1.0 to 100.0) for SCOs under relaxed positive selection for bird-infecting species (n = 401) and mammal-infecting species (n = 104). Boxplots represent the interquartile range (IQR) with the median as a horizontal line and whiskers extending to 1.5 × IQR. For relaxed purifying selection, the mammal group had a Q1 of 0.38, median of 0.52, Q3 of 0.68, and whiskers from 0.09 to 0.98; the bird group had a Q1 of 0.38, median of 0.56, Q3 of 0.83, and whiskers from 0.08 to 1.00. For relaxed positive selection, the mammal group had a Q1 of 2.36, median of 5.79, Q3 of 33.92, and whiskers from 1.02 to 80.48; the bird group had a Q1 of 2.40, median of 7.32, Q3 of 31.15, and whiskers from 1.03 to 73.13.

Fig. 5. GO enrichment of 140 expanded gene families in T. vaginalis.

Dot size represents the number of genes with a specific GO term. Biological process (BP), cellular component (CC), and molecular function (MF) are plotted. Only significant GO enrichments after FDR correction (0.05 < ) are reported.

Fig. 6. T. vaginalis shares the largest number of expanded gene families with human/mammal-infecting T. tenax A.

Heatmap showing number of shared expanded orthogroups between seven trichomonad species. T. vaginalis shares its largest number of expanded orthogroups with T. tenax (n = 33) and T. sp. 2a (n = 35), two non-sister species that infect humans and birds, respectively. B GO enrichment of convergently expanded multicopy gene families in T. vaginalis and T. tenax. C GO enrichment of convergently expanded multicopy gene families in T. vaginalis and T. sp. 2a. For both (B and C): size of dot represents the number of genes with a specific GO term; biological process (BP), cellular component (CC) and molecular function (MF) are plotted; only significant GO enrichments after FDR correction (0.05 < ) are reported.