Transcriptional profile of Trichomonas vaginalis in response to metronidazole
BMC Genomics volume 24, Article number: 318 (2023)
Trichomoniasis caused by Trichomonas vaginalis, combined with its complications, has long frequently damaged millions of human health. Metronidazole (MTZ) is the first choice for therapy. Therefore, a better understanding of its trichomonacidal process to ultimately reveal the global mechanism of action is indispensable. To take a step toward this goal, electron microscopy and RNA sequencing were performed to fully reveal the early changes in T. vaginalis at the cellular and transcriptome levels after treatment with MTZ in vitro.
The results showed that the morphology and subcellular structures of T. vaginalis underwent prominent alterations, characterized by a rough surface with bubbly protrusions, broken holes and deformed nuclei with decreased nuclear membranes, chromatin and organelles. The RNA-seq data revealed a total of 10,937 differentially expressed genes (DEGs), consisting of 4,978 upregulated and 5,959 downregulated genes. Most DEGs for the known MTZ activators, such as pyruvate:ferredoxin oxidoreductase (PFOR) and iron-sulfur binding domain, were significantly downregulated. However, genes for other possible alternative MTZ activators such as thioredoxin reductase, nitroreductase family proteins and flavodoxin-like fold family proteins, were dramatically stimulated. GO and KEGG analyses revealed that genes for basic vital activities, proteostasis, replication and repair were stimulated under MTZ stress, but those for DNA synthesis, more complicated life activities such as the cell cycle, motility, signaling and even virulence were significantly inhibited in T. vaginalis. Meanwhile, increased single nucleotide polymorphism (SNP) and insertions - deletions (indels) were stimulated by MTZ.
The current study reveals evident nuclear and cytomembrane damage and multiple variations in T. vaginalis at the transcriptional level. These data will offer a meaningful foundation for a deeper understanding of the MTZ trichomonacidal process and the transcriptional response of T. vaginalis to MTZ-induced stress or even cell death.
Trichomoniasis caused by Trichomonas vaginalis is the most prevalent nonviral sexually transmitted disease worldwide. A recent report by the World Health Organization (WHO) estimated approximately 156 million new cases of trichononasis globally . Usually, T. vaginalis infection can lead to variable clinical manifestations, ranging from asymptomatic carriers to patients with severe inflammation of the genital system in women or men . Rather than these, it may also increase the risk of HIV  and cervical neoplasia  acquisition in women and probably result in several adverse pregnancy consequences . Although infected males are widely asymptomatic and self-limited, the increased risk of nongonococcal urethritis , impaired infertility , prostate hyperplasia and even prostate cancer  associated with male trichomoniasis should not be ignored.
Metronidazole (MTZ), as one of the mainstay drugs for killing anaerobic bacteria, microaerophilic bacteria, and protozoa such as Giardia lamblia, Entamoeba histolytica, and T. vaginalis, has been widely used for the treatment of related infectious diseases in the clinic for several decades . It is also considered the first choice for the systemic treatment of trichmononasis with an excellent cure rate, and there are no replacement drugs with superior advantages to date . Therefore, a deeper understanding of its trichomonacidal process to ultimately fully reveal the global mechanism of action is indispensable and meaningful for the better and more accurate use of this drug against trichmononasis and even other pathogens. After decades of unremitting efforts by scientists, research on the mechanism by which MTZ kills T. vaginalis has made tremendous progress, however, it has not yet been fully illustrated .
As a prodrug, MTZ is inactive until it enters T. vaginalis and becomes activated after its nitro group is reduced to cytotoxic nitro radical anions [12, 13]. Previous studies revealed several MTZ activation-related molecules, such as pyruvate:ferredoxin oxidoreductase (PFOR), ferredoxin (Fd) [12, 14, 15] and the flavin enzyme thioredoxin reductase (TrxR) . Studies also indicate that rapid responses of T. vaginalis exposure to MTZ occur in vitro, such as DNA synthesis interruption at approximately 30 min, cell division ceasing one hour later, motility stopping within one to two hours, and ultimately parasite death occurring within 5 to 8 h [17,18,19]. However, most studies about the mechanism of MTZ against T. vaginalis seemed to be focused on the activation of this drug and the final ending of the parasites. Regarding the influence of MTZ on T. vaginalis after MTZ treatment, but before T. vaginalis widely dies, especially how MTZ affects the gene expression of T. vaginalis post coincubation, there are limited documents.
Furthermore, there are frequently refractory or recurrent trichomoniasis cases reported in the clinic in recent years . Investigate the reasons, downregulated genes or reduced activity of PFOR, Fd, flavin reductase and nitroreductase (NTR), in some aerobic or anaerobic resistant strains are probably closely related to it [14, 21,22,23]. However, changes in these genes cannot explain all cases of MTZ resistance . The deeper and overall mechanism of MTZ resistance in trichomoniasis needs further investigation. Certainly, more explicit recognition of the MTZ action is beneficial for a better understanding of MTZ resistance.
Recently, transcriptome sequencing has provided an excellent tool to dissect the differentially expressed genes (DEGs) between different research conditions, which could supply global and important information reflecting the influence of one or more inducements on the target cells or tissues at the transcriptional level . Therefore, the present study aims to reveal the changes in gene expression of T. vaginalis in response to MTZ at an early stage after MTZ treatment but before T. vaginalis widely dies depending on RNA-seq. The results will help us gain a more comprehensive understanding of the drug action model and the transcriptional response of T. vaginalis, or even other pathogens, to MTZ-induced stress and even cell death, which could also offer a clue to help understanding the resistant mechanism and develop new drugs.
Structural changes in T. vaginalis after treatment with MTZ
According to Fig. S1A, the minimum lethal concentration (MLC) of MTZ on TV-THS1 was 12.5 µg/ml, lower than the resistance threshold (50 µg/ml) . Therefore, TV-THS1 was considered to be a sensitive strain. The drug concentration at which approximately half of the parasites compared to the control (MTZ = 0 µg/ml) remained after treatment with MTZ for 3 h was 25 µg/ml (Fig. S1B). Therefore, T. vaginalis treated with MTZ at 25 µg/ml for 3 h was chosen for morphology and subcellular structure analysis and RNA-seq.
Observation by scanning electron microscopy showed that the trophozoites of normal T. vaginalis were oval or pear-shaped. There were folds of different depths on the cell surface, and the membrane was clear (Fig. 1A). After MTZ treatment for 3 h, T. vaginalis was still ovoid, but the surface was rough, and there were many bubbly protrusions and small depressions of varying sizes. There were also obvious broken holes on the surface of the cell membrane with off-white granular substances attached to the surface of the trophozoites (Fig. 1B). For the severely damaged parasites, the cell membrane was completely broken, leaving the cauliflower-like cytoplasmic contents (Fig. 1C).
Transmission electron microscopy presented that normal T. vaginalis had an intact cell membrane. The nuclear membrane was also intact with abundant nuclear chromatin inside. The endoplasmic reticulum, Golgi apparatus, and axis column were all clearly visible (Fig. 1D). After treatment with MTZ for 3 h, T. vaginalis underwent prominent alterations. The nucleus was deformed in morphology. Part of the nuclear membrane disappeared, and the boundary between the nucleus and cytoplasm was blurred (Fig. 1E). The Golgi and endoplasmic reticulum were also blurred and even disappeared. For some trophozoites, the chromatin in the nucleus decreased, or the electron density in the entire cell was significantly reduced (Fig. 1F). Furthermore, some trophozoites were obviously disrupted, leaving only a network structure. However, compared with the control (untreated group), the morphology and number of hydrogenases and vacuoles did not change significantly.
General property and quality assessment of transcriptome sequencing
First, the quality of sequencing data was assessed, which consisted of the data output statistics, base content across all bases, quality scores across all bases, and sequence GC content. As shown in Table 1, the sequencing finally produced an average of 20,945,600 and 19,660,581 clean read pairs with extremely low error rates in the MTZ-treated group (Q20: 96.9%) and untreated group (Q20: 97.1%), respectively. The highly matched base pairs (A-T, C-G), good quality for each base, and the approximately 45% GC content all demonstrated the extremely high quality and accuracy of the current sequencing data (Fig. S2). Of these data sets (Table S1), 36,469,263 (87.07%, MTZ treated) and 35,284,989 (89.73%, untreated) clean reads were mapped to the reference genome of T. vaginalis. A total of 22,810 (MTZ-treated) and 23,208 (untreated) expressed genes were clustered by the fragments per kilobase per million bases (FPKM) values for five intervals. Of these, gene numbers within the FPKM > 1 intervals all dramatically deceased after MTZ treatment (Table S2).
High coverage depth and homogenization in data sets revealed the high quality of the RNA-seq in the untreated group, but RNA degradation appeared in the MTZ-treated group (Fig. S3). The saturability differentiation of an expression level less than 15% guaranteed quantitative accuracy (Fig. S3). The high repeatability and correlation of gene expression within the three repeats in both the untreated and MTZ-treated groups further manifested the experimental reliability and reasonability of sample selection (Fig. 2A and B). Additionally, the density distribution of genes all disclosed high consistency within the three repeats in each group and total downregulated expression in MTZ-treated samples (Fig. 2 C and 2D).
Analysis of differentially expressed genes
To globally assess the quantitative changes in the MTZ-treated transcriptome, the DEGs (red dots) were screened out by edgeR, and the selective condition was defined as false discovery rate (FDR) < 0.05 and Log2 Fold Change (MTZ-treated/MTZ-untreated) (Log2 FC) > 1 or Log2 FC <-1. The MA plot (Fig. 3A) and volcano plot (Fig. 3B) visually show that more genes were inhibited than stimulated in T. vaginalis after MTZ treatment. Only a relatively small portion of the genes were not differentially expressed (black dots). Statistically, a total of 10,937 genes were differentially expressed, including 4,978 upregulated and 5,959 downregulated genes (Fig. 3C). The heatmap and clustering tree revealed a distinct gene expression pattern between the MTZ treated and untreated groups, but similar within a group (Fig. 3D). The hierarchical clustering revealed three clusters of genes.
Differential expression of known MTZ activators
Of these genes (Table 2), three of the five annotated genes for PFOR were differentially expressed, with one upregulated approximately 4.0-fold and two downregulated approximately 3.0-fold and 2.4-fold, respectively. Four genes for the 4Fe-4 S binding domain-containing protein showed differential expression, with only one upregulated approximately 2.2-fold and the remaining three downregulated approximately 3.0- to 10.3-fold. However, six annotated genes for Fd did not change significantly. Additionally, one of three genes mapped to thioredoxin reductase, an alternative MTZ activator, was significantly upregulated by approximately 3.8-fold. Five of six DEGs encoding nitroreductase family proteins were significantly upregulated by approximately 9.2- to even 739.3-fold. Meanwhile, 24 genes annotated to the flavodoxin-like fold family proteins were differentially expressed, including two downregulated approximately 3.0-fold and 22 upregulated approximately 2.1- to even 685.0-fold (Table 3).
GO and KEGG enrichment analysis of the DEGs
To further explore the changes in DEGs of T. vaginalis exposed to MTZ, GO and KEGG analysis were performed. The significantly enriched (Q value < 0.05) upregulated DEGs consisted of 74 GO terms in biological processes, 25 terms in cellular components and 31 terms in molecular functions (Table S3). Accordingly, the downregulated DEGs were significantly enriched in 78 biological process terms, 31 cellular component terms and 43 molecular function terms (Table S4). As shown in Fig. 4, the top 10 significantly enriched terms in each primary displayed that genes for metabolism and biosynthetic processes, intracellular components (especially for organelles and cytoplasm), structural molecules, and oxidoreductase activity were specifically upregulated together with genes involved in transmembrane-related transport, such as drug transmembrane transporters. In contrast, a relatively smaller number of genes involved in metabolic processes and organelles/organellar lumen were significantly downregulated. Instead, large numbers of genes contributing to macromolecule/protein modification, biological regulation, the nucleus, and bandings of various intracellular molecules, especially nucleotide binding, were significantly inhibited.
However, KEGG analysis revealed many more significantly enriched downregulated pathways than upregulated pathways. As shown in Fig. 5, among the 15 significantly upregulated pathways, genes were mainly mapped to ribosome (337 out of 421), aminoacyl-tRNA biosynthesis (15 out of 33) and proteasome (30 out of 46). In accordance with the upregulated GO terms, a total of 173 genes were significantly enriched in metabolism-related pathways involving amino acid carbon, carbohydrate, lipid, cofactors and vitamin. Interestingly, 16 out of the 31 genes involved in homologous recombination during replication and repair and 34 out of the 88 genes involved in ABC transporters were significantly upregulated. For the top 20 significantly downregulated pathways (Fig. 5), genes were widely mapped to the following various levels: translation (e.g., mRNA surveillance pathway), nucleotide metabolism (e.g., purine metabolism), cell growth and death (e.g., cell cycle, meiosis), cell motility (e.g., regulation of actin cytoskeleton), cellular community-eukaryotes (e.g., adherens junction, focal adhesion, gap junction), transport and catabolism (e.g., mitophagy, autophagy), signal transduction in environmental information processing sets and organismal systems (e.g., circulatory system, endocrine system, environmental adaptation, excretory system, immune system, nervous system and sensory system). Further analysis of the other 37 significantly enriched KEGG pathways showed that many genes involved in system regulation, transcription, translation, signal transduction, transport, catabolism, community, cell growth and death, as well as those involved in various signaling pathways, cellular junctions and the cell cycle, were dramatically downregulated (Table S5).
Single nucleotide polymorphism (SNP) and insertion - deletions (indel) analysis
After MTZ treatment, the total SNPs increased but were not significant (Table 4). However, a more interesting finding in our study showed that dramatically increased indels (from 568 to 2161, P = 0.0317), especially heterozygous indels (from 289 to 1939, P = 0.0275), appeared in the transcripts after MTZ treatment. Further analysis showed that the main types of these indels were deletions, and most of them were G/C deletions. For both the untreated (62%) and MTZ treated groups (56%), most SNPs appeared in exonic regions. However, SNPs were significantly increased in other regions including the 3’UTR, introns and [upstream; downstream] (3 kb region of upstream from transcription start site; 3 kb region of downstream from the transcription termination site) after MTZ treatment. Similarly, the exonic region contained the most indels for both groups. By MTZ treatment, indels in regions of downstream (3 kb region of downstream from the transcription termination site), exonic, intronic and upstream (3 kb region of upstream from transcription start site) were all significantly increased (Fig. S4).
Confirmation of gene expression by RT‒qPCR
To validate the RNA-seq results, expression level of 14 significant DEGs closely related to the known MTZ activators were relatively quantified by RT‒qPCR (Fig. 6A), including PFOR E (EAY02588.1), 4Fe-4 S binding domain containing protein (EAY19262.1 and EAY22303.1), thioredoxin reductase (EAY04700.1) and ten flavodoxin-like fold family proteins (EAY20871.1, EAX97338.1, EAY00354.1, EAY15747.1, XP_001307810.1, XP_001322702.1, XP_001317528.1, XP_001317526.1, XP_001306204.1 and XP_001319426.1). All the primers for target genes and reference gene  are shown in Table S6. The results of RT‒qPCR showed that the selected genes were significantly differentially expressed, and their upregulated or downregulated trends were consistent with those obtained by RNA-Seq. Additionally, a high correlation (R2 = 0.8078) measured by Log2 FC between RNA-Seq and RT‒qPCR was observed (Fig. 6B), indicating the reliability of DEGs obtained from transcriptome sequencing.
Revealing the transcriptional changes in T. vaginalis after MTZ treatment, but before T. vaginalis widely dies, can assist in a deeper understanding of the MTZ trichomonacidal process and the transcriptional response of T. vaginalis to MTZ-induced stress or even cell death. After MTZ treatment for 3 h, variations in the morphology and subcellular structures of T. vaginalis were observed. RNA-seq revealed more downregulated DEGs than the upregulated ones. Of these, most DEGs referring to the well-known MTZ activators, such as PFOR and iron-sulfur binding domain, were significantly inhibited, but gene sets annotated to thioredoxin reductase, nitroreductase family proteins and flavodoxin -like fold family proteins, the other possible alternative MTZ activators, were dramatically stimulated. GO and KEGG analyses revealed that genes for basic vital activities, proteostasis, replication and repair were stimulated under MTZ stress, but those for DNA synthesis, more complicated life activities such as the cell cycle, motility, signaling and even virulence were significantly inhibited in T. vaginalis.
Of the 10,937 DEGs, 4,978 upregulated and 5,959 downregulated DEGs were discovered in T. vaginalis, suggesting that more genes were significantly downregulated in T. vaginalis in response to MTZ-induced stress. This is different from a previous study in which 14,072 DEGs were found, with 8432 upregulated and 6270 downregulated . A possible cause of the difference should be the diversity in experimental materials and conditions between the two studies, such as T. vaginalis isolates, MTZ concentrations and the time for treatment. This, accompanied by the inherent variability of the transcriptome, indicated that the results of one transcriptome-based study would be theoretically difficult to resolve the target problem in one study. In other words, the shared findings among different studies probably have greater significance.
The DEGs encoding PFOR and 4Fe-4 S binding domain-containing proteins were almost all downregulated by 2.4- to even 10.3-fold, except for two that were upregulated. PFORs from T. vaginalis are known as homodimers of 240 to 280 kDa containing 2[4Fe-4 S] clusters . A previous study concluded that a strong correlation existed between the presence of PFOR activity and MTZ among different microorganisms . Previous studies also reported that genes for PFOR were significantly downregulated in resistant isolates of T. vaginalis  and Giardia lamblia . Therefore, as the key proteins for MTZ activation, the wide downregulation of PFOR was likely to be a self-rescue response at the early stage. However, one of three genes mapped to thioredoxin reductase and five of six DEGs encoding nitroreductase family proteins were significantly upregulated by approximately 3.8-to to even 739.3-fold. Thioredoxin reductase (with nitroreductase activity) and nitroreductase have been reported to play roles in activating MTZ against T. vagianlis depending on an alternative theory . A “covalent adduct” between MTZ with thioredoxin reductase and other enzymes can be formed, leading to the failure of thioredoxin activation and then blocking thioredoxin peroxidase from reducing H2O2, which is cytotoxic. This, perhaps, is an important mechanism of action of MTZ other than DNA disruption [16, 20]. Interestingly, studies have confirmed that nitroreductase gene is also an important sign of MTZ resistance, which mainly refers to SNPs in two nitroreductase genes (ntr4Tv and ntr6Tv). These SNPs introduced nonsense mutations, leading to stop codons in ntr4Tv and ntr6Tv associated with MTZ resistance [20, 21]. Therefore, their upregulation in the present study may further demonstrate the important alternative roles of these enzymes in MTZ action. Additionally, another family including 24 DEGs was annotated to flavodoxin-like fold sets, with 22 upregulated approximately 2.1- to even 685.0-fold. Flavodoxin - dependent MTZ reduction has been reported in Helicobacter pylori , indicating a possible functional prediction of these flavodoxin-like fold genes in T. vaginalis. Actually, a study proved that flavodoxin-like proteins in T. vaginalis, as functionally redundant electron donors with Fd (no differential expression in the present study), were indeed capable of activating MTZ . Its wide and high upregulation presents an interesting target involved in MTZ activation that deserves further concern, especially in MTZ-resistant T. vaginalis isolates.
From the significantly enriched GO terms, genes responsible for metabolism, intracellular components, structural molecules and oxidoreductase activity were the most upregulated, indicating that responses supporting basic vital activity were stimulated under MTZ stress. Instead, large numbers of genes contributing to macromolecule/protein modification and biological regulation, especially cell nuclear and nucleotide binding, were inhibited. This was in accordance with findings in a previous study showing that most genes involved in nucleotide metabolism were dramatically suppressed in MTZ-sensitive parasites after MTZ treatment . These results, combined with the broken nuclear membrane and decreased chromatin observed by TEM in the present study, further verify the DNA damage-dependent mechanism of MTZ.
In KEGG analysis, the significantly enriched upregulated genes were mainly involved in the basic life process, similar to the GO analysis. However, three special pathways were also stimulated under MTZ stress: proteasome for folding, sorting degradation, homologous recombination for replication and repair, and ABC transporters. Of these, all 14 standard proteasome subunits of core particles for the 20 S proteasome were significantly upregulated after MTZ treatment by approximately 2- to 4-fold. The 20 S proteasome, as a large protease complex consisting of seven α subunits (α1 to α7) and seven β subunits (β1 to β7), can widely degrade intracellular proteins in all eukaryotes . This is meaningful for proteostasis because the accumulation of unfolded or misfolded proteins in the endoplasmic reticulum lumen might be harmful to cells . Short-term incubation of T. vaginalis with a proteasome inhibitor leads to autophagy of parasites, implying its important function in T. vaginalis survival [35, 36]. Upregulation of the replication and repair gene pathway is also important for T. vaginalis survival under emergency conditions. For membrane transport, genes in ABC transporter pathways with a rich factor of 38.6% (34/88) were upregulated. Of these, a total of 18/34 genes were enriched in the ATP-binding cassette, subfamily B (MDR/TAP) and subfamily A (ABC1), and member 3 pathways with high Log2 FC even reached up to 8.0. Interestingly, these two ATP-binding cassette subfamilies were all reported to be multidrug resistance-associated proteins in humans [37, 38]. Therefore, the upregulation of ABC transporters may be a self-saving route for T. vaginalis to resist damage by MTZ.
Purines and pyrimidines are essential precursors for nucleotide and nucleic acid synthesis. However, T. vaginalis lacks the enzymes to synthesize de novo nucleotides, making it an obligate parasite reliant on salvage pathways. Purine salvage is mediated by some nucleoside phosphorylases and kinases . KEGG data in the present study revealed that 90 genes in the purine metabolism pathway were significantly suppressed. This indicated that the normal nucleic acid metabolic pathway was sharply disturbed, which may be another mechanism by which MTZ functions. Widely downregulated transcripts of mRNA surveillance pathways (85/222) indicate the collapse in fidelity and quality of mRNA. This is consistent with an interesting finding in our study that a dramatic increase in indels (from 568 to 2161), especially the G/C deletion, was observed in the transcripts after MTZ treatment. Previous studies revealed that the nitro radical of MTZ is partial to targeting sections of DNA rich in the thymine and adenine residues (A + T) [13, 19]. In genome of T. vaginalis, approximately 71% of the A + T content explained the specificity and sensitivity of MTZ for this parasite . The significantly increased deletion of G/C, may imply another mechanism of MTZ actions. In which, a higher proportion of A + T residues were left and enriched after G/C deletion at the gene level, resulting in a more severe impairment in these regions. Of course, this speculation needs the direct and reliable evidence. Additionally, the increased SNPs and indels in specific regions may introduce more coding variations, such as missense mutation or nonsense mutation. For example, the formation of premature stop codons in genes due to SNPs or indels, especially for those closely related to the drug’s action, perhaps leads to the inactivation of key enzymes and then drug resistance appears [24, 40, 41]. As reported in another study, 72 SNPs in T. vaginalis were proven to be associated with MTZ resistance .
The downregulation of cell meiosis and cell cycle genes indicated that the cell cycle of T. vaginalis was rapidly disrupted and may even arrest. This is in accordance with our observation of the structure of T. vaginalis by SEM/TEM and the previous finding that DNA synthesis could be inhibited within 30 min and T. vaginalis died within 5 h after incubation with MTZ in vitro . It is also worth mentioning that the wide downregulation of genes for actin cytoskeleton regulation should be responsible for the rapid decrease in cell motility. At 3 h postincubation with MTZ in the present study, almost all T. vaginalis had stopped moving. The downregulated DEGs in adhesion and junctions suggested that the adhesion-based virulence of T. vaginalis was reduced. After all, motility and adhesion to host cells are pivotal steps in the pathogenesis of T. vaginalis . Remarkably, genes for many signaling pathways were downregulated. In brief, the genes for purine metabolism, mRNA surveillance, cell cycle, various signaling pathways, motility and adhesion were inhibited, indicating that widely transcriptional disorders have appeared from DNA synthesis to more complicated life activities and even virulence in T. vaginalis.
Although our study revealed early changes in T. vaginalis at the cellular and transcriptome levels after treatment with MTZ for 3 h in vitro, there were also limitations in the study. For example, dynamic observation at more time points, including the early, middle, and universal dying stages, and the involvement of more different strains of T. vaginalis should have more advantages for better understanding the trichomonacidal process of MTZ.
The present study focuses on the total changes in the morphology, subcellular structures and transcriptome of T. vaginalis after MTZ treatment in vitro (Fig. 7A). The results revealed evident nuclear and cytomembrane damage beyond other changes and multiple changes in T. vaginalis at the transcriptional level. Specifically, key genes for MTZ activation were highly concerned (Fig. 7B). Other genes for basic vital activities, proteostasis and replication and repair were stimulated under MTZ stress, but those for DNA synthesis, more complicated life activities and even virulence were significantly inhibited in T. vaginalis (Fig. 7C). Meanwhile, SNPs and indels in specific regions of genes were increased after MTZ treatment (Fig. 7C). These data will offer a meaningful foundation for a deeper understanding of the MTZ trichomonacidal process and the transcriptional response of T. vaginalis to MTZ-induced stress or even cell death, which potentially takes a step toward ultimately revealing the global mechanism of MTZ action and then helps understanding the resistant mechanism.
Materials and methods
The T. vaginalis strain TV-THS1 used in this work was isolated from a female patient undergoing routine examination in Taihe Hospital (Shiyan, China). The leucorrhea specimen was obtained with cotton swabs, kept in normal saline and passaged for 48 generations before being used in vitro. Ethical approval for this study was obtained from the Medical Ethics Committee of the Hubei University of Medicine (2018-TH-003). Informed consent was supplied by the participating individual.
Parasite culture and treatment by MTZ
T. vaginalis were cultured in modified hepar-peptone-glucose medium supplemented with 10% heat-inactivated fetal calf serum, 100 U/ml penicillin, and 100 µg/ml streptomycin and maintained at 37 °C. Only the logarithmic-phase parasites were used in this work. Before RNA-seq, the MTZ susceptibility of T. vaginalis was determined by microscopy. T. vaginalis, approximately 2 × 105/ml in 200 µl, was incubated with different concentrations (from 200 to 0 µg/ml) of MTZ in 96-well plates for 24 h. The drug concentration at which no viable parasites were observed after trypan blue staining was defined as the MLC. The MTZ IC50 of T. vaginalis was determined after treatment for 3 h in vitro. The drug concentration at which approximately half of the parasites compared to the control (MTZ = 0 µg/ml) was chosen for electron microscopic observation and RNA-seq.
Preparation of scanning electron microscopy and transmission electron microscopy
Approximately 5 × 107T. vaginalis trophozoites at the logarithmic phase were collected and fixed in 2.5% glutaraldehyde. The parasites were postfixed for 2 h at room temperature in 1% OsO4. Then, they were dehydrated in graded ethanol and isoamyl acetate and dried with a critical point dryer. The specimens were attached to metallic stubs using carbon stickers and sputter-coated with gold for 30 s. Images were collected with a scanning electron microscope (Regulus 8100, HITACHI).
For transmission electron microscopy, T. vaginalis fixed in 2.5% glutaraldehyde were first preembedded in agarose and then postfixed in 1% OsO4 for 2 h at room temperature. After being dehydrated in graded ethanol and acetone, they were embedded in EMBed 812. The resin blocks were cut to 60–80 nm, and the tissues were fished onto 150-mesh cuprum grids with formvar film. Then, the sections were stained with 2% uranium acetate and 2.6% lead citrate. Images were finally collected with a transmission electron microscope (HT7700, HITACHI).
RNA extraction and cDNA library preparation for illumina sequencing
The parasite materials were maintained with TRIzol at -80 °C until RNA extraction. Sequencing was performed by the Frasergen Corporation (Wuhan, China). After RNA contamination and degradation were monitored on 1.5% agarose gels, RNA purity and concentration were measured with a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 3.0 Fluorometer (Life Technologies, CA, USA), respectively. Then, RNA integrity was assessed by using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
Qualified RNA was submitted to generate sequencing libraries using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s instructions and index codes. Briefly, mRNA was first gathered from the total RNA with poly-T oligo-attached magnetic beads. Then, mRNA was fragmented with divalent cations under elevated temperature in NEB Next First Strand Synthesis Reaction Buffer (5×). By using random hexamer primers and M-MuLV Reverse Transcriptase, first-strand cDNA was synthesized, which was followed by the synthesis of second-strand cDNA using DNA Polymerase I and RNase H. After double-stranded cDNA was purified by AMPure XP beads, USER enzyme was used to degrade the second strand cDNA containing U. The remaining purified double-stranded cDNA was processed with an end-repair reaction, the addition of a single ‘A’ base and adapter ligation. To select and purify cDNA fragments with the right length, PCR and AMPure XP beads (Beckman Coulter, Beverly, USA) were applied to obtain the final cDNA library. The qualified library assessed by the Agilent Bioanalyzer 2100 System was finally sequenced using the Illumina HiSeq2000.
RNA-seq data analysis
For quality control of the sequencing data, the software FraserQC (v1.2), an in-house software developed by Frasergen Co. Ltd. (Wuhan, China), was used. Clean data (clean reads) were selected from the raw data after a series of routine processes. For read alignment, the T. vaginalis genome obtained by de novo sequencing on the PacBio Sequel platform (Pacific Biosciences of California, Menlo Park, CA, USA) conducted by our laboratory was chosen as the reference genome by using Tophat2 (v2.1.1)  and Bowtie2 (v2.2.2)  with default parameters. The expression levels of genes were quantified using the software package RSEM (RNASeq by Expectation Maximization v1.3.0) . After FPKM transformation , screening of DEGs with the EdgeR (v3.6.8) package  was performed based on the following criteria: fold change (FC = condition 2/condition 1 for a gene) ≥ 2 and FDR < 0.05. The upregulated DEGs met the requirements of FDR < 0.05 and Log2 FC > 1, and the downregulated DEGs met the criteria of FDR < 0.05 and Log2 FC< -1. To visualize the distribution of FC and FDR values of all genes between the two groups, MA plot and volcano plot were drawn. Then, the numbers of upregulated and downregulated DEGs were sent for statistical analysis and heatmap drawing of the known gene expression pattern clusters.
GO and KEGG enrichment analysis
After screening out target DEGs, Gene Ontology (GO) enrichment analysis was accomplished using GOseq (v1.22) based on Wallenius noncentral hypergeometric distribution, which contains three ontologies: molecular function, biological process, and cellular component. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs was performed in KOBAS (v3.0) [48,49,50]. For both GO enrichment and KEGG pathway analysis, FDR ≤ 0.05 was considered significantly enriched.
Quantification and region-based annotation of SNPs and indels
The Picard tools version 1.96 were used to sort the bam files obtained by comparing them with the reference genomes. After duplication removal, SNPs and indels were called using SAMtools version 0.1.19 and Genome Analysis Toolkit (GATK) version 3.7. Finally, the obtained SNPs and indels with high quality were sent for region-based annotation using ANNOVAR version 2.1.
Quantitative real-time PCR (RT‒qPCR)
The parasite materials and RNA extraction were the same as those used for RNA-Seq. Reverse transcription and qPCR were conducted using the PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) and TB Green® Premix Ex Taq™ II (Tli RNaseH Plus) kit according to the manufacturer’s protocols (Takara, Dalian, China). Briefly, after the genomic DNA removal reaction, total RNA was reverse transcribed to cDNA, which was used as the template for the following qPCRs. The reaction systems and conditions of RT‒qPCR were all in accordance with the manufacturer’s suggestion on a CFX96 Real-Time PCR Detection System (Bio-Rad, USA). Melting curve analysis was performed to validate specific amplification. The primers used for RT‒qPCR were designed with primer-BLAST. The actin gene was used as the reference gene for normalization, and the relative expression was calculated with the 2−∆∆Ct method. Three technical replicates and two negative controls without template were set up for each biological replicate.
The plots and statistical analysis of the MLC and IC50 of MTZ against T. vaginalis, the statistics for distribution of SNP and indel, the relative quantification by RT-qPCR and the correlation between RNA-Seq with RT-qPCR were all performed using GraphPad Prism 8.0.2. The determinations were expressed as the means ± standard deviations (SDs). Correlation analysis was performed with linear regression. SNP and indels with significant differences between groups were analyzed by T-test (F test, P>0.05), or Welch’s t test (F test, P < 0.05) in GraphPad Prism 8.0.2.
The clean data sets for sequencing presented in this study were uploaded into the Sequence Read Archive (SRA) in online repositories. The BioProject accession code PRJNA849385 can be found in NCBI through the website https://www.ncbi.nlm.nih.gov/sra/PRJNA849385.
Differentially expressed gene
World Health Organization
Single nucleotide polymorphism
Minimum lethal concentration
Fragments per kilobase per million bases
False discovery rate
Kyoto Encyclopedia of Genes and Genomes
Sequence read archive
Nicotinamide adenine dinucleotide phosphate
WHO. Report on global sexually transmitted infection surveillance., 2018. Geneva: World Health Organization;. 2018; Licence: CC BY-NC-SA 3.0 IGO.
Bouchemal K, Bories C, Loiseau PM. Strategies for Prevention and Treatment of Trichomonas vaginalis infections. Clin Microbiol Rev. 2017;30(3):811–25.
Masha SC, Cools P, Sanders EJ, Vaneechoutte M, Crucitti T. Trichomonas vaginalis and HIV infection acquisition: a systematic review and meta-analysis. Sex Transm Infect. 2019;95(1):36–42.
Ghosh I, Muwonge R, Mittal S, Banerjee D, Kundu P, Mandal R, et al. Association between high risk human papillomavirus infection and co-infection with Candida spp. and Trichomonas vaginalis in women with cervical premalignant and malignant lesions. J Clin Virol. 2017;87:43–8.
Workowski KA, Bachmann LH, Chan PA, Johnston CM, Muzny CA, Park I, et al. Sexually transmitted Infections Treatment Guidelines, 2021. MMWR recommendations and reports: morbidity and mortality weekly report. Recommendations and reports. 2021;70(4):1–187.
Ziaei Hezarjaribi H, Saberi R, Fakhar M, Sadeghian N. Is There Any Relationship between Trichomonas vaginalis Infection and Male Urethritis Risk? A Systematic Review and Meta-Analysis. Interdiscip Perspect Infect Dis. 2022; 2022:8359859.
Zhang Z, Li Y, Lu H, Li D, Zhang R, Xie X, et al. A systematic review of the correlation between Trichomonas vaginalis infection and infertility. Acta Trop. 2022;236:106693.
Yang HY, Su RY, Chung CH, Huang KY, Lin HA, Wang JY, et al. Association between trichomoniasis and prostate and bladder diseases: a population-based case-control study. Sci Rep. 2022;12(1):15358.
Hernández Ceruelos A, Romero-Quezada LC, Ruvalcaba Ledezma JC, López Contreras L. Therapeutic uses of metronidazole and its side effects: an update. Eur Rev Med Pharmacol Sci. 2019;23(1):397–401.
Weir CB, Le JK. Metronidazole. StatPearls Treasure Island (FL): StatPearls Publishing Copyright © 2022,StatPearls Publishing LLC.; 2022.
Dingsdag SA, Hunter N. Metronidazole: an update on metabolism, structure-cytotoxicity and resistance mechanisms. J Antimicrob Chemother. 2018;73(2):265–79.
Edwards DI. Nitroimidazole drugs–action and resistance mechanisms. I. Mechanisms of action. J Antimicrob Chemother. 1993;31(1):9–20.
Cudmore SL, Delgaty KL, Hayward-McClelland SF, Petrin DP, Garber GE. Treatment of infections caused by metronidazole-resistant Trichomonas vaginalis. Clin Microbiol Rev. 2004;17(4):783–93. table of contents.
Yarlett N, Yarlett NC, Lloyd D. Metronidazole-resistant clinical isolates of Trichomonas vaginalis have lowered oxygen affinities. Mol Biochem Parasitol. 1986;19(2):111–6.
Vidakovic M, Crossnoe CR, Neidre C, Kim K, Krause KL, Germanas JP. Reactivity of reduced [2Fe-2S] ferredoxins parallels host susceptibility to nitroimidazoles. Antimicrob Agents Chemother. 2003;47(1):302–8.
Leitsch D, Kolarich D, Binder M, Stadlmann J, Altmann F, Duchêne M. Trichomonas vaginalis: metronidazole and other nitroimidazole drugs are reduced by the flavin enzyme thioredoxin reductase and disrupt the cellular redox system. Implications for nitroimidazole toxicity and resistance. Mol Microbiol. 2009;72(2):518–36.
Ings RM, McFadzean JA, Ormerod WE. The mode of action of metronidazole in Trichomonas vaginalis and other micro-organisms. Biochem Pharmacol. 1974;23(9):1421–9.
Nielsen MH. In vitro effect of metronidazole on the ultrastructure of Trichomonas vaginalis donné. Acta Pathol Microbiol Scand B. 1976;84(2):93–100.
Edwards DI. Mechanisms of selective toxicity of metronidazole and other nitroimidazole drugs. Br J Vener Dis. 1980;56(5):285–90.
Graves KJ, Novak J, Secor WE, Kissinger PJ, Schwebke JR, Muzny CA. A systematic review of the literature on mechanisms of 5-nitroimidazole resistance in Trichomonas vaginalis. Parasitology. 2020;147(13):1383–91.
Paulish-Miller TE, Augostini P, Schuyler JA, Smith WL, Mordechai E, Adelson ME, et al. Trichomonas vaginalis metronidazole resistance is associated with single nucleotide polymorphisms in the nitroreductase genes ntr4Tv and ntr6Tv. Antimicrob Agents Chemother. 2014;58(5):2938–43.
Kulda J. Trichomonads, hydrogenosomes and drug resistance. Int J Parasitol. 1999;29(2):199–212.
Leitsch D, Janssen BD, Kolarich D, Johnson PJ, Duchêne M. Trichomonas vaginalis flavin reductase 1 and its role in metronidazole resistance. Mol Microbiol. 2014;91(1):198–208.
Bradic M, Warring SD, Tooley GE, Scheid P, Secor WE, Land KM, et al. Genetic indicators of Drug Resistance in the highly repetitive genome of Trichomonas vaginalis. Genome Biol Evol. 2017;9(6):1658–72.
Yu Y, Li J, Wang W, Wang T, Qi W, Zheng X, et al. Transcriptome analysis uncovers the key pathways and candidate genes related to the treatment of Echinococcus granulosus protoscoleces with the repurposed drug pyronaridine. BMC Genomics. 2021;22(1):534.
Lossick JG, Müller M, Gorrell TE. In vitro drug susceptibility and doses of metronidazole required for cure in cases of refractory vaginal trichomoniasis. J Infect Dis. 1986;153(5):948–55.
dos Santos O, de Vargas Rigo G, Frasson AP, Macedo AJ, Tasca T. Optimal reference genes for gene expression normalization in Trichomonas vaginalis. PLoS ONE. 2015;10(9):e0138331.
Huang P-J, Huang C-Y, Li Y-X, Liu Y-C, Chu L-J, Yeh Y-M, et al. Dissecting the Transcriptomes of multiple metronidazole-resistant and sensitive trichomonas vaginalis strains identified distinct genes and Pathways Associated with Drug Resistance and Cell Death. Biomedicines. 2021;9(12):1817.
Ali V, Nozaki T. Current therapeutics, their problems, and sulfur-containing-amino-acid metabolism as a novel target against infections by “amitochondriate” protozoan parasites. Clin Microbiol Rev. 2007;20(1):164–87.
Leitsch D, Burgess AG, Dunn LA, Krauer KG, Tan K, Duchêne M, et al. Pyruvate:ferredoxin oxidoreductase and thioredoxin reductase are involved in 5-nitroimidazole activation while flavin metabolism is linked to 5-nitroimidazole resistance in Giardia lamblia. J Antimicrob Chemother. 2011;66(8):1756–65.
Kaihovaara P, Höök-Nikanne J, Uusi-Oukari M, Kosunen TU, Salaspuro M. Flavodoxin-dependent pyruvate oxidation, acetate production and metronidazole reduction by Helicobacter pylori. J Antimicrob Chemother. 1998;41(2):171–7.
Land KM, Delgadillo-Correa MG, Tachezy J, Vanacova S, Hsieh CL, Sutak R, et al. Targeted gene replacement of a ferredoxin gene in Trichomonas vaginalis does not lead to metronidazole resistance. Mol Microbiol. 2004;51(1):115–22.
Huber EM, Basler M, Schwab R, Heinemeyer W, Kirk CJ, Groettrup M, et al. Immuno- and constitutive proteasome crystal structures reveal differences in substrate and inhibitor specificity. Cell. 2012;148(4):727–38.
Kocaturk NM, Gozuacik D. Crosstalk between mammalian autophagy and the ubiquitin-proteasome system. Front cell Dev biology. 2018;6:128.
Huang KY, Chen RM, Lin HC, Cheng WH, Lin HA, Lin WN, et al. Potential role of autophagy in proteolysis in Trichomonas vaginalis. J Microbiol Immunol Infect. 2019;52(2):336–44.
O’Donoghue AJ, Bibo-Verdugo B, Miyamoto Y, Wang SC, Yang JZ, Zuill DE et al. 20S Proteasome as a Drug Target in Trichomonas vaginalis. Antimicrob Agents Chemother 2019; 63(11).
Klugbauer N, Hofmann F. Primary structure of a novel ABC transporter with a chromosomal localization on the band encoding the multidrug resistance-associated protein. FEBS Lett. 1996;391(1–2):61–5.
Chen CJ, Chin JE, Ueda K, Clark DP, Pastan I, Gottesman MM, et al. Internal duplication and homology with bacterial transport proteins in the mdr1 (P-glycoprotein) gene from multidrug-resistant human cells. Cell. 1986;47(3):381–9.
Miller RL, Lindstead D. Purine and pyrimidine metabolizing activities in Trichomonas vaginalis extracts. Mol Biochem Parasitol. 1983;7(1):41–51.
Palmieri N, de Jesus Ramires M, Hess M, Bilic I. Complete genomes of the eukaryotic poultry parasite Histomonas meleagridis: linking sequence analysis with virulence / attenuation. BMC Genomics. 2021;22(1):753.
Liu Y, Wang S, Yang F, Chi W, Ding L, Liu T, et al. Antimicrobial resistance patterns and genetic elements associated with the antibiotic resistance of Helicobacter pylori strains from Shanghai. Gut Pathog. 2022;14(1):14.
Petrin D, Delgaty K, Bhatt R, Garber G. Clinical and microbiological aspects of Trichomonas vaginalis. Clin Microbiol Rev. 1998;11(2):300–17.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.
Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics 2010; Chap. 11:Unit 11.7.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–d92.
This work was partially supported by the Principle Investigator Program of Hubei University of Medicine (HBMUPI202101) and foundation from the Research Project of Hubei Provincial Department of Education (Q20192102) and Post-Graduates of Hubei University of Medicine (2018QDJZR17). The funders had no role in the study design and data analysis, the decision to publish, or reparation of the manuscript.
Ethics approval and consent to participate
The collection and use of the parasite strain in this study were approved by the Medical Ethics Committee of the Hubei University of Medicine (2018-TH-003) and methods were carried out in compliance with relevant guidelines and regulations. Informed consent was supplied by the participating individual for strain collection.
Consent for publication
The authors declare that they have no competing interests.
The authors would like to thank the Clinical Lab and Department of Gynecology and Obstetrics, Affiliated Taihe Hospital of Hubei University of Medicine, and the participants who contributed their leucorrhea samples.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
About this article
Cite this article
Xie, Y., Zhong, P., Guan, W. et al. Transcriptional profile of Trichomonas vaginalis in response to metronidazole. BMC Genomics 24, 318 (2023). https://doi.org/10.1186/s12864-023-09339-9