Clues on the dynamics of DNA replication in Giardia lamblia

Eukaryotic DNA replication is a high-fidelity process that undergoes precise regulation and follows specific steps that vary slightly among the diverse taxonomic groups (Kelly and Callegari, 2019; Machida et al., 2005). In general, DNA replication begins with the firing (activation) of replication origins which initiates the S phase of the cell cycle.

Each fired origin generates two replication forks that move in opposite directions (bidirectional movement) and are responsible for synthesizing the DNA at an average rate (velocity) that varies from cell to cell (Stanojcic et al., 2016; Coster and Diffley, 2017; da Silva et al., 2020). The time required for chromosome duplication determines the duration of the cell cycle S phase, which is a strategy to regulate cell cycle progression (Günesdogan et al., 2014; Turrero-García et al., 2016). However, the entry into the S phase, as well as progression throughout the cell cycle are mostly dependent on the expression and activation of cyclin-dependent kinases (CDKs), which are turned into holoenzymes in the presence of their coenzymes termed cyclins. Briefly, at the end of the G1 phase, after the cell crosses the restriction point (Johnson and Skotheim, 2013), the activated cyclin-E–CDK2 complex promotes the G1/S transition. The cyclin-A–CDK1 and cyclin-A–CDK2 complexes control the progression throughout the S phase (DNA replication), whereas cyclin-B–CDK1 promotes the orderly development of mitosis, ensuring the correct segregation of the chromosomes (Hochegger et al., 2008).

The pervasive view is that, during the S phase, all eukaryotes replicate their chromosomes from multiple origins (Fragkos et al., 2015; Kelly and Callegari, 2019). Each eukaryotic chromosome has a minimum number of replication origins (MO) that need to be fired to complete replication within the S phase (da Silva et al., 2019, 2020; da Silva, 2020). The MO varies according to the organism and cell type but essentially depends on three factors: chromosome size, S phase duration and the average replication rate. The MO also contributes to the establishment of a standard threshold where the estimate of a smaller number of fired origins relative to the MO is implausible, whereas the opposite might suggest the presence of replication stress (da Silva et al., 2020).

Most eukaryotes tightly regulate origin licensing and firing (origin usage) during the cell cycle, through several conserved mechanisms, to avoid aberrant origin firing and DNA re-replication (Parker et al., 2017; Sclafani and Holzen, 2007). These mechanisms are predominantly regulated by differential expression of cell-cycle regulatory proteins, such as cyclins and CDKs (Gaggioli et al., 2020; Tanaka et al., 2007), and licensing and firing factors (Boos and Ferreira, 2019; Köhler et al., 2016; Truong and Wu, 2011). Furthermore, impairments during DNA replication, such as head-on replication–transcription conflicts and secondary structure formation might lead to intra-S checkpoint activation. This process primarily regulates the replication forks by slowing down and/or arresting forks together with their stabilization. The intra-S checkpoint can also block the firing of new origins while promoting the activation of dormant origins close to the region that presents replication stress (Hu et al., 2012; García-Muse and Aguilera, 2016; Iyer and Rhind, 2017; Prorok et al., 2019).

Most DNA metabolism processes have been determined using experiments with cells containing a single nucleus (mononucleated). Some of these processes might be imprecise when analyzing organisms containing two nuclei (binucleated), such as Giardia spp. This single-celled eukaryote is the etiological agent of giardiasis, a diarrheic illness that affects millions of people, predominantly from low-income regions (Choy et al., 2014; Coelho et al., 2017; Horton et al., 2019; Vivancos et al., 2018). The two nuclei in Giardia are morphologically indistinguishable and are transcriptionally active (Kabnick and Peattie, 1990). Moreover, each nucleus has two complete copies of the haploid genome arranged in five chromosomes, although constitutive aneuploidy is reported (Tůmová et al., 2016, 2019). Nuclei are partitioned equationally during the cytokinesis phase of the cell division cycle, with each daughter cell inheriting one nucleus from each of the parental nuclei (Sagolla et al., 2006; Yu et al., 2002). The unusual feature of a cell containing two nuclei raises questions: is the replication process for both nuclei synchronized, as previously suggested (Wiesehahn et al., 1984)? If this is the case, does the presence of varying quantities of DNA damage in each nucleus during the G1 phase cause the S phase to begin at different times? How frequently do endogenous replication issues arise that cause different nucleus replication rates?

Here, we thoroughly investigated DNA replication in Giardia lamblia trophozoites, and analyzed the spatial and temporal organization of replication units and fired origins to better understand the dynamics of DNA synthesis that might ensure synchronous whole-genome duplication between nuclei.

Our findings demonstrated a high replication rate, long inter-origin distance and the utilization of a small number of origins. This suggests that Giardia trophozoites are subjected to modest levels of DNA replication stress during the cell cycle. Also, an examination of DNA replication patterns using a double pulse of two different thymidine analogs showed that ∼80% of Giardia cells synchronously replicate their nuclei. Furthermore, single-molecule analysis of moving replication forks in long nanopore reads (DNAscent) was performed to better understand the regulation of Giardia DNA replication. This approach allowed the identification of possible head-on replication–transcription (HoRT) collisions and their target genes. DNAscent forksense analysis detected 5561 genes (49.65% of analyzed sequencing reads) containing potential HoRT collision sites.

These combined analyses provide the first comprehensive atlas of whole-genome replication dynamics in Giardia. Furthermore, our findings provide a body of specific knowledge that might help answer questions about DNA replication in binucleated organisms.

DNA combing was used to investigate replication features in G. lamblia. This technique allows the visualization of replication origins in replicated DNA molecules stretched onto coverslips through subsequent short-pulse incorporation of the halogenated thymidine analogs IdU (red) and CldU (green). The replication rate was measured by determining the ratio between the length of specific green tracks (in kb) and the period of CldU incorporation (20 min). The specific green tracks continued from red tracks on one side, and blue tracks – non-incorporated DNA on the other side. In other words, the length of the green tracks was measured from blue-green-red (BGR), red-green-blue (RGB), and blue-green-red-green-blue (BGRGB) patterns (Fig. 1A). A histogram showing the replication rate frequency of the analyzed molecules is presented (Fig. 1B). The mean±s.d. replication rate was 3.64±1.32 kb min⁻¹ (Fig. 1C) and the median rate among the molecules analyzed was 3.7 kb min⁻¹. Interestingly, this value is higher relative to those in other single-celled protozoa including Trypanosoma cruzi (2.05 kb min⁻¹) (de Araujo et al., 2019), Leishmania donovani (2.37 kb min⁻¹) (Stanojcic et al., 2016), Leishmania major (2.45 kb min⁻¹) (Stanojcic et al., 2016), Leishmania mexicana (2.48 kb min⁻¹) (Stanojcic et al., 2016), and different strains of Trypanosoma brucei [1.84 kb min⁻¹ for Lister 427 (Stanojcic et al., 2016) and 3.06 kb min⁻¹ for TREU927 (da Silva et al., 2019)].

The origin distribution along the G. lamblia chromosomes was characterized by measuring the inter-origin distance (IOD). The IOD is defined as the length between two adjacent replication domains. Each replication domain is determined by having two adjacent replication forks in progression (Stanojcic et al., 2016).

A representative image of two replication domains with white arrows representing origins is shown in Fig. 1D. DNA fibers ranging from 86 to 975 kb were combed; there was more than one origin per molecule, allowing an easy estimate of the IOD. The median IOD for G. lamblia was 294.5 kb (mean±s.d. of 335.7±36.49 kb) (Fig. 1E,F). The IOD from G. lamblia is slightly higher relative to those in other single-celled eukaryotes [T. cruzi IOD=171 kb (de Araujo et al., 2019), different species of Leishmania (L. donovani IOD=179.7 kb, L. major IOD=192.7 kb, L. mexicana IOD=203.3 kb), and T. brucei Lister 427 IOD=148.8 kb], and substantially different from that of budding yeast (IOD=46 kb) (Stanojcic et al., 2016).

Although further assays are necessary to determine the mechanism associated with this genomic feature, the high IOD value of G. lamblia infers that it replicates its genome using fewer origins than do other single-celled eukaryotes. The high IOD value suggests that few origins are fired (activated) during a standard G. lamblia cell cycle, which indicates low levels of endogenous replication stress in this organism.

To obtain more evidence supporting the use of a few replication origins by Giardia, we followed up on the above experiment by characterizing the S phase length and estimating the minimum number of origins needed to complete genome duplication during S phase. DNA replication was initially monitored to further characterize and update the S phase duration. A recent study in protozoa has revealed that significant differences are seen when using standard approaches with BrdU or EdU to monitor DNA replication (da Silva et al., 2017a,b). In summary, EdU is more sensitive than BrdU for monitoring DNA replication. Therefore, we used EdU to monitor DNA replication to improve the accuracy of our analyses.

Analysis of growth curves determined a doubling time of 6 h (r²⁼0.984) for G. lamblia trophozoites (Fig. 2A,B). This was incorporated into William's equation, and Stanner's and Till's equation (Stanners and Till, 1960; Williams, 1971) to estimate the cell cycle phase lengths (see Material and Methods). The percentage of cells performing cytokinesis (C) was estimated at 8.77±0.93% 4N trophozoites (mean±s.d.; n=308) (Fig. 2C), and the duration of the C phase was estimated at 0.73 h or 0.122 ccu (cell cycle unit) through assessing the morphology of nuclei stained with DAPI (4N configuration) combined with the morphology of the parasite using differential interference contrast (DIC).

Duration of the G2+M phase was estimated at 1.25 h, indicating that cells at the end of S phase required 1.25 h to proceed through G2 and M phases (Fig. 2D). The results were consistent amongst the three biological replicates; all cells contained four EdU-labeled nuclei at the same time.

Edu-labeling of trophozoites (1 h pulse) showed that 35.7±3.5% (mean±s.d.) of parasites can replicate DNA (Fig. 2E). This value was used to calculate the S phase duration and the estimated duration of the G2+M+C phases (see Materials and Methods). The S phase was estimated at 1.16 h or 0.194 ccu (Fig. 2F), which is required to estimate the minimum number of origins used by G. lamblia.

Our group recently developed a formula to estimate the MO required to duplicate an entire chromosome within the S phase (da Silva et al., 2019). This formula can be applied to any organism (da Silva, 2020) and is based on the bidirectional movement of the replication forks, replication rate, S phase duration and the chromosome size in question (see Materials and Methods for more details).

The previously obtained parameters (Figs 1 and 2) were used to estimate the MO for each chromosome from G. lamblia (Fig. 3A). There was a positive correlation between the number of origins and the size of chromosomes, which correlates with the results of other organisms (da Silva, 2020) (Fig. 3A,B). The MO was compared with the number of origins estimated by DNA combing (Fig. 3B). G. lamblia uses slightly more origins than the minimum number needed (compare the red and black dots in Fig. 3B). The MO is estimated from relatively constant parameters in a wild-type population; this allows the establishment of a threshold (dashed black line in Fig. 3B) that serves as a parameter to validate the number of origins estimated by other techniques, such as DNA combing (da Silva, 2020).

The ratio between the angular coefficient (a value) of linear equations (y=ax+b) of origins was estimated by DNA combing and MO (red and black dashed lines, respectively in Fig. 3B). This indicated that G. lamblia uses 1.55 times more origins than the MO (green bar in Fig. 3C). Compared with other single-celled eukaryotes, G. lamblia fits between the yeasts Schizosaccharomyces pombe and Saccharomyces cerevisiae, and with a quite different value to that of other protozoa parasites, such as trypanosomatids (compare the green bar with the other gray bars in Fig. 3C).

These results corroborate the previously obtained DNA combing data: G. lamblia uses a few origins above the MO to complete replication during the S phase. This suggests that Giardia activates only a few backup origins during DNA replication, most likely due to the lack of significant replication stresses. Therefore, an intriguing question can be raised: given that DNA replication in Giardia appears to be subjected to low levels of stress, is the replication behavior in both nuclei similar?

Previous autoradiographic analysis of thymidine [³H]-labeled trophozoites has indicated that replication occurs simultaneously in both nuclei (Wiesehahn et al., 1984). We revisited this data using a double-label strategy to monitor DNA replication. This approach consisted of sequentially labeling the DNA that was being replicated using the two different halogenated thymidine analogs IdU and CldU. Giardia presents a unique ventral disc that is responsible for its adhesion to different surfaces (Schwartz et al., 2012). A basic positioning pattern for the nuclei was established (RN and LN for the right and left nucleus, respectively) given that they always hold their position relative to the ventral disc.

A column containing illustrations of all possible patterns (left side) from the representative images of different G. lamblia trophozoite models found in our analyses is shown in Fig. 4A. The DIC image shows the morphological integrity of the parasites. DAPI (blue) was used to stain the DNA from both nuclei. The specific antibodies recognize their respective thymidine analogs (red for IdU and green for CldU) in different patterns according to the nuclei position. The measurement of all obtained nuclei patterns (n=388) is shown on a bar graph (Fig. 4B). The highest proportion of trophozoites (78.49%) exhibited stunning replication synchronicity (replicating both nuclei concomitantly), with 62.47±4% (mean±s.d.) showing both nuclei as yellow, 7±1.5% as red and 9±1.5% as green, while the remaining cells (21.51%) exhibited an unusual variable pattern (asynchronous replication). A more precise grouping of the patterns was undertaken according to when the beginning and end of DNA replication occurred in each nucleus as shown in Fig. 4C: nuclei replicate together (yellow, 78.49±0.7%), undetermined patterns (gray, 7.47±3.4%), one nucleus (left or right) completes replication first (blue, 11.54%±3.6), and one nucleus (left or right) initiates replication first (green, 2.51% ±0.5).

Curiously, all the undetermined parasites had at least one nucleus without replication activity, even after 2 h in the presence of thymidine analogs. DAPI fluorescence is inaccurate for quantitative analyses of DNA content. Therefore, further assays are required to determine whether these few apparently inactive nuclei are arrested before the S phase due to DNA damage (G1/S checkpoint), exhibit slower replication due to a possible intra-S checkpoint, or are quiescent relative to its twin nucleus. Also, the more than 20% of Giardia trophozoites exhibiting asynchronous replication makes us wonder about the origin of this asynchrony given that the single-molecule approach did not find evidence of the presence of significant replication stress.

Synchronized nuclei replication in most parts of the population and a high replication rate relative to other single-celled eukaryotes suggests that DNA replication in G. lamblia occurs smoothly. Therefore, it appears that few replication stresses occur, allowing both nuclei (from ∼80% of the population) to start and conclude DNA replication in the same time frame. One of the sources of replication stress is the head-on collision between replication and transcription machineries (HoRT collisions), which occurs when these machineries are traveling in opposite and convergent directions. We previously showed that such HoRT collisions occur in the trypanosomatid T. brucei, due to its genome organization (da Silva et al., 2019). Therefore, we asked whether the low occurrence of replication stress in G. lamblia is due to genome organization. We hypothesized that HoRT collisions are minimized due to favorable guidance of replication and transcription machinery, in contrast to what is seen in T. brucei (da Silva et al., 2019), but similar to what is seen in bacteria (Wu et al., 2020). However, the possible absence of HoRT collisions does not explain the presence of more than 20% of Giardia trophozoites exhibiting asynchronous replication. To further investigate the dynamics of Giardia genomic organization, we applied the cutting-edge technique DNAscent.

DNAscent methodology allows us to visualize the direction of replication forks in the genome by the MinION sequencing (long reads) of BrdU-incorporated DNA fragments. Briefly, asynchronous trophozoites are pulsed with BrdU thymidine analogs to label nascent DNA. The DNA is extracted, sequenced and analyzed by DNAscent software considering the BrdU gradient, which provides replication fork direction (Müller et al., 2019; Boemo, 2021). The transcription direction was obtained from GiardiaDB (https://giardiadb.org/), and the direction of replication and transcription events were compared within the same DNA molecule. A potential HoRT collision was considered when transcription and replication presented convergent directions in the same molecule (see scheme of potential HoRT collisions in Fig. 5A).

Two independent experiments totalizing 545,880 MinION sequencing reads were performed with an average read length of 30 kb (Fig. S1). The reads were aligned with the G. lamblia genome using Minimap2 software with an alignment rate of 94.9%. BrdU incorporation of the MinION raw, base called and mapped reads were then submitted for DNAscent analysis.

Common base calling and alignment software cannot identify the BrdU signals, which are interpreted as thymidine nucleotides. DNAscent software combines the electrical signal characteristic of BrdU with the thymidine positions in the read and determines the probability that the base is, in fact, BrdU. After the program identifies the thymidine analog regions, it calculates the probability of replication fork direction through the decay of BrdU on the read (Boemo, 2021). We selected all molecules with a probability ≥50% to determine the fork direction. This probability was chosen by verification of false positives on the Integrative Genomics Viewer (IGV), which provides a large set of tools for inspection, validation and interpretation of next-generation sequencing (NGS) datasets, as well as other types of genomic data. A sample of DNAscent output reads with a ≥50% fork probability was visualized in IGV to compare BrdU decay and fork sense (Fig. S2). False positives were assigned to reads with divergent directions between BrdU decay and fork direction. IGV analysis demonstrated only 19% false positives considering Nieduszynski parameters (Dr Conrad Nieduszynski, Earlham Institute, UK, personal communication) and a ≥50% probability. The number of false positives did not decrease with increasing fork probability. As a result, 3663 BrdU-positive DNA molecules (0.67% of total reads) were obtained that corresponded to regions containing fork directions that were mapped with the annotated genome of G. lamblia. Those that corresponded to gene sequences (3622) are plotted in Table S1; a single-cell methodology was used, and the same gene can appear more than once. DNAscent only detected 0.03% of total reads as BrdU-positive in the negative control (performed in the absence of BrdU). BrdU-positive DNA molecules with a determined fork direction were analyzed concerning potential HoRT collisions. The establishment of replication fork direction allows us to infer the transcription direction in the same DNA molecule and identify probable conflict sites between the replication and transcription machinery. Replication forks were detected in ∼2000 reads from nanopore sequencing ranging from 6664 to 100,800 bp with an average size of 23,406 bp, and forks were observed in both directions (Fig. S3). Annotation of fork direction coordinates in the G. lamblia genome identified a total of 11,400 genes: 5739 (50.3%) were in genomic regions without potential HoRT collision, while 5661 (49.7%) could be a source of replicative stress by potential HoRT collisions. This indicates that G. lamblia genome organization does not show a co-oriented bias to avoid replication–transcription conflicts, which might explain the presence of more than 20% of G. lamblia trophozoites exhibiting asynchronous replication.

Head-on collisions damage DNA, leading to genetic variability through repair mechanisms (Sankar et al., 2016). As it was expected, we found that genes in regions of potential HoRT collisions were enriched in GO categories linking to chromatin (dis)assembly and organization, DNA repair, microtubule-based processes, regulation of signal transduction, cell cycle and glycerolipid biosynthetic processes (Fig. 5C,D; Table S2). This is due to the fact that these GO terms are related to the only genes that are actively transcribed during the S phase, as well as those required for mitosis. Genes located in genomic regions without HoRT collisions were represented by GO terms associated with the maintenance of protein location within the cell, metabolic processes, proteasome processing and regulation of G1/S transition in mitotic cells. Giardia undergoes antigenic variation to overcome the host immune system (Prucca and Lujan, 2009). Therefore, we expected that the variant-specific surface proteins (VSP) genes would be positioned in regions prone to conflict as a way to increase variability. However, there were no significant differences in the distribution of VSP genes to regions that were or were not affected by HoRT conflicts, despite its heterogeneity (data not shown). This finding supports the proposal that antigenic variation in Giardia mainly occurs at the post-transcriptional level (Prucca and Lujan, 2009). To the best of our knowledge, there is no transcriptome data for Giardia during the S phase available. If there was, it would be interesting to assess whether the HoRT collision sites correlate with loci containing highly transcribed genes. Although our findings suggest that genes related to chromatin dynamics and DNA repair might be altered due to replication-transcription conflicts, further assays are needed to show an increase or decrease in variability in these genes. However, the enrichment of these genes in HoRT conflict regions might explain the presence of asynchronous replication between nuclei observed in ∼20% of the population.

This study systematically investigated the dynamics of DNA replication in G. lamblia trophozoites by single-molecule techniques, including DNA combing and nanopore-based sequencing, followed by DNAscent analysis to determine the spatial and temporal activation of replication origins and fork progression. G. lamblia replication parameters (including mean IODs, replication rate and MO) were substantially different relative to other single-celled eukaryotes (da Silva et al., 2019, 2020). Remarkably, DNA replication occurred faster in Giardia (Fig. 1) than in T. cruzi (de Araujo et al., 2019), T. brucei (da Silva et al., 2019), L. major (Lombraña et al., 2016; Stanojcic et al., 2016) and S. cerevisiae (Sekedat et al., 2010), which is consistent with low levels of replication stress (Zhong et al., 2013).

Giardia uses on average 1.55 more origins than the minimum required. This indicates that Giardia probably does not activate backup origins in response to endogenous replication stress, unlike T. brucei (da Silva et al., 2019). A possible explanation for this result is the reduction of genes transcribed polycistronically by RNA pol II (Clayton, 2019). Most individual genes in Giardia have individual promoters (Holberton and Marshall, 1995; Teodorovic et al., 2007). Giardia has bidirectional transcription, which is an inherent feature of its promoters and contributes to an abundance of sterile antisense transcripts throughout the genome (Elmendorf et al., 2001; Teodorovic et al., 2007). In other words, there is abundant evidence that transcription in Giardia is a well-regulated process (Holberton and Marshall, 1995; Knodler et al., 1999; Yee et al., 2000; Davis-Hayman et al., 2003; Teodorovic et al., 2007) and, therefore, transcriptional control might contribute to avoiding conflicts with replication machinery, at least in part of the genome, thereby resulting in a low rate of dormant origin activation.

Curiously, our precise analysis of DNA replication in single trophozoite cells confirmed that both nuclei replicate synchronously in 78.49% of cells, with a variation in stochasticity of replication timing in the other 21.51% of the population. Synchronous DNA replication between nuclei suggests the simultaneous firing of different origins during the S phase. However, the sequence composition of these origins is still unknown. The asynchronous DNA replication for the other 21.51% of the Giardia population might be due to (1) early entry of one of the nuclei into the S phase, and/or (2) replication arrest on one nucleus but not the other.

The first situation could occur if endogenous DNA damage occurs in the G1 phase in only one nucleus, which would activate a G1/S checkpoint that blocks the beginning of the S phase for that nucleus. We believe that this is theoretically possible because the G1/S checkpoint is predominantly activated in the nuclear space when Cdc25A phosphorylation renders it inactive by nuclear exclusion and degradation. In the absence of Cdc25A, the S phase promoting cyclin-E–Cdk2 complex is inactivated (Sancar et al., 2004). Another possibility that would help explain this first situation is a possible unequal partitioning of some unidentified essential DNA replication factor, such as ORC, Cdc6, MCM_2-7, Cdt1, Cdc45 or other factors not yet characterized in Giardia (Morrison et al., 2007; Tiengwe et al., 2012). There is evidence that the inhibition and/or removal of some of these factors impairs the progression of DNA replication in mononucleated protozoa (Dang and Li, 2011; Benmerzouga et al., 2013; Deshmukh et al., 2015; Kim, 2019; da Silva et al., 2017a,b). We can hypothesize that an analogous situation could be occurring in just one of the Giardia nuclei. Further cutting-edge assays, such as a high-resolution single-cell type approach to map only one Giardia nucleus, could bring evidence to support this hypothesis.

The second situation would occur during replication stress in only one nucleus. In this case, the intra-S checkpoint (involving only nuclear players) (Iyer and Rhind, 2017) is activated, which arrests or delays DNA replication in only one nucleus. Furthermore, ∼1% of cells (4/388) contained three nuclei (data not shown). This might represent cells where mitosis occurred in one nucleus but not in the other. The asynchrony of nuclei division is not expected because mitosis is controlled by a cytoplasmic cyclin–Cdk (the equivalent of mammalian CyclinB–Cdk1) that is translocated to the nucleus after its activation at the beginning of mitosis (late prophase) (Takizawa and Morgan, 2000). Therefore, non-finalized replication in even one nucleus would activate the G2/M checkpoint and block the entrance of mitosis activator Cyclin–Cdk to both nuclei. The finding that only 1% of cells exhibit non-synchronized mitosis corroborates this scenario. Alternatively, non-finalized replication in just one nucleus might not be enough to activate the G2/M checkpoint. Further experiments are needed to explain this question.

Our findings allowed us to improve our understanding of DNA replication timing during the S phase by analyzing replication fork movement and potential head-on collisions between replication and transcription in G. lamblia. We provide evidence that the organization of Giardia genome imparts a co-orientation of replication and transcription, probably to preserve the expression of conserved genes involved in metabolic and biosynthetic processes like those observed in bacteria and humans (Srivatsan et al., 2010; Liu et al., 2017). Interestingly, genes enriched in locations of potential HoRT collisions are related to microtubule-based processes, cell cycle regulation and DNA replication mechanisms (Fig. 5C,D). Most of these genes are essential for Giardia pathogenesis, including cell division, transport of cellular cargo, organelle positioning and cell migration. For example, the dynein light intermediate chain is an essential subunit that contributes to spindle checkpoint inactivation (Mische et al., 2008) and the assembly and activation of the dynein–dynactin complex (Schroeder and Vale, 2016). Most interestingly, our GO analysis of genes located in sites with potential HoRT collisions showed enrichment for ATP binding and adenyl nucleotide binding. Previous studies have shown that transcription of these genes peaked in the G0/G1 phase in human cells (Liu et al., 2017). In addition, genes highly transcribed in the M phase and enriched for the GO term ‘cell cycle’ in humans are also localized in genomic regions with HoRT collisions in Giardia. This suggests that these genes are not transcribed (or present at a reduced transcription rate) during the S phase, which might contribute to minimizing replication stress. On the other hand, the genes related to chromatin dynamics, cell cycle regulation, and DNA replication/repair mechanisms located in regions prone to potential HoRT collisions might suffer DNA damage due to these conflicts, such as double-strand breaks and R-loop formation (da Silva et al., 2019), impairing the cell cycle of the affected nucleus and, consequently, the synchronicity of the replication in part of the population, which might explain the observed asynchronous replication of part of the population. The most intriguing question we still cannot answer is why these HoRT collisions do not seem to activate dormant origins to keep the robustness of the S phase duration. Further in-depth single-cell genomics approaches are necessary to isolate the Giardia cells exhibiting asynchronous replication and investigate their genome.

In summary, our findings provide important clues regarding DNA replication dynamics in Giardia, which paves the way for a better understanding of the mechanisms regulating DNA replication and gene transcription in this peculiar parasite that affect millions of people worldwide.

G. lamblia trophozoites strain WB (American Type Culture Collection, catalog number 50582) were cultivated at 37°C in sterilized TYI-S-33 medium (Keister, 1983) supplemented with 10% (v/v) adult bovine serum, 0.5 mg ml⁻¹ bovine bile (Sigma-Aldrich) and 1× antibiotic-antimycotic solution (Thermo Fisher Scientific). Daily growth curves were generated by culturing trophozoites at 10⁴ cells ml⁻¹ and collecting samples up to the 6th day. Hourly curves were created by culturing trophozoites at 10⁴ cells ml⁻¹ and collecting samples up to 7 h.

Exponentially growing trophozoites were sequentially labeled with two halogenated thymidine analogs dissolved in fresh culture medium: 100 μM of 5-iodo-2′-deoxyuridine (IdU; Sigma-Aldrich) for 20 min, and 100 μM of 5-chloro-2′ deoxyuridine (CldU; Sigma-Aldrich) for an additional 20 min. There was no intermediate washing; instead, the IdU-containing medium was removed by centrifugation (500 g, 5 min, at 21°C) before adding the CldU-containing medium. After labeling, the cells were immediately harvested by centrifugation at 500 g for 10 min at 4°C, washed once with cold PBS (137 mM NaCl, 2.7 mM KCl, 10 mM Na₂HPO₄ and 2 mM KH₂PO₄, pH 7.4), and resuspended in 100 μl of warmed PBS containing 1% low-melting agarose to embed the cells in agarose plugs. The plugs were suspended in 500 μl of lysis solution [0.5 M EDTA pH 8.0, 1% (v/v) N-lauryl-sarcosyl, and 4 μg ml⁻¹ proteinase K] and incubated at 50°C for 24 h. Fresh lysis solution was added, and the plugs were incubated for an additional 24 h. The plugs were carefully washed several times using 0.5 M EDTA at pH 8.0, to enable the complete removal of digested proteins and other degradation products. Protein-free DNA plugs were stored in 0.5 M EDTA pH 8.0 at 4°C or immediately used. Plug samples were washed in TE buffer (10 mM Tris-HCl, 1 mM EDTA, pH 8.0), resuspended in 1 ml of 0.5 M MES buffer [2-(N-morpholino) ethanesulfonic acid pH 5.5] and melted at 68°C for 20 min. The solution was maintained at 42°C for 10 min and digested overnight with 2 U of β-agarase (New England Biolabs, Ipswich, MA, USA); 1 ml of 0.5 M MES was added, and DNA fibers were regularly stretched (2 kb μm⁻¹) on silanized coverslips using a molecular combing system (Genomic Vision) as previously described (Michalet et al., 1997).

Combed DNA was fixed onto coverslips at 65°C for at least 2 h, denatured in 1 M NaOH for 20 min, and washed several times in PBS. Coverslips containing the DNA fibers were blocked with 1% (w/v) BSA in PBS. Immunodetection was performed using primary antibodies mouse anti-BrdU antibody clone B44 (AB_400443) (Becton Dickinson, Franklin Lakes, NJ, USA); and rat anti-BrdU antibody clone BU1/75 (ICR1) (AB_1523225, Abcam, Cambridge, UK), with both antibodies diluted 1:20 in 1% (w/v) BSA, and incubated at 37°C in a humid chamber for 1 h. It is worth highlighting that mouse anti-BrdU reacts with IdU and BrdU (Bakker et al., 1989; Gratzner, 1982), and rat anti-BrdU antibody reacts with CldU and BrdU but does not cross-react with thymidine or IdU (Bakker et al., 1989). The coverslips were incubated with secondary goat anti-rat-IgG conjugated to Alexa Fluor 488 diluted 1:20 (Life Technologies, Carlsbad, CA, USA), and goat anti-mouse-IgG conjugated to Alexa Fluor 568 (1:20 dilution, Life Technologies) antibodies. Each antibody incubation step was followed by extensive washes with PBS. DNA immunodetection was performed using anti-ssDNA antibody (1:50 dilution, Chemicon, São Paulo, Brazil) and goat anti-mouse-IgG conjugated to Alexa Fluor 350 (MAB3868, Merck Millipore; 1:10 dilution, Life Technologies). Coverslips were mounted with 20 μl of Prolong Gold antifade mountant (Thermo Fisher Scientific), dried at 21°C for at least 2 h, and processed for image acquisition using an Olympus BX51 fluorescent microscope (Japan) with a 100× oil objective attached to an EXFO Xcite series 120Q lamp and an Olympus XM10 digital camera (Japan). Images were superimposed using the Cell F software (Olympus, Tokyo, JP) when required. The observation of longer DNA fibers required the capture of adjacent fields. Fibers <100 kb were excluded from the analysis. The percentage of origins activated during the thymidine pulses among the patterns of randomly collected DNA fibers was manually measured using ImageJ software (National Institutes of Health, Bethesda, MD, USA). The formula/time of the CldU pulse was used to estimate the average replication rate in the analyzed population. Statistical analyses of the replication rate and inter-origin distance (IOD) were performed using Graph Pad Prism (version 5.0, GraphPad Software, La Jolla, CA, USA). Three independent combing experiments were performed.

Exponentially growing trophozoites were incubated with 100 µM 5-ethynyl-2′-deoxyuridine (EdU; Thermo Fisher Scientific) for different periods (30–105 min) at 37°C. The parasites (∼10⁵ total cells) were harvested by centrifugation at 500 g for 10 min at 4°C, washed three times in PBS, and the pellet was resuspended in 150 µl of cell fixative solution [1% (v/v) Triton X-100 , 40 mM citric acid, 20 mM dibasic sodium phosphate, 0.2 M sucrose, pH 3.0] for 10 min. Diluent buffer (300 µl) (125 mM MgCl₂ in PBS) was added, and the mixture was incubated for 20 min. The samples were centrifuged at 500 g for 10 min at 4°C, washed and resuspended in 100 µl of PBS. Cell suspension (20 µl) was loaded onto poly-L-lysine pre-treated microscope slides (Tekdon, Miakka City, FL, USA) and washed three times with 3% (w/v) BSA (Sigma-Aldrich, St. Louis, MO, USA) diluted in PBS. Incorporated EdU was detected using Click-iT EdU detection solution for 45 min protected from light. This solution is composed of 5 µl of 100 mM copper sulfate (CuSO₄), 2.5 µl of Alexa Fluor 488 azide (Thermo Fisher Scientific), 25 µL of 500 mM ascorbic acid (C₆H₈O₆) and 467.5 µl of distilled water (for details regarding the EdU procedure, see Salic and Mitchison, 2008). Finally, the cells were washed five times with PBS. Vectashield mounting medium (Vectorlabs, Burlingame, CA, USA) containing DAPI was used as an anti-fade mounting solution and to stain nuclei DNA. Images were acquired with the Olympus BX51 fluorescent microscope with a 100× oil objective attached to an EXFO Xcite series 120Q lamp and a digital camera (Olympus XM10, Olympus, Tokyo, JP). Images were further analyzed using ImageJ software (National Institute of Health) to count the numbers of EdU-positive cells. The percentage of proliferating parasites was calculated for each sample relative to the total number of DAPI-positive parasites. Images were superimposed using the Cell F software (Olympus, Tokyo, JP) when necessary.

Finally, the duration of the G1 phase was estimated by calculating the difference between the doubling time (dt) and the sum of the other phases (S+G2+M+C).

The MO analyses used the estimated S phase duration and replication rate determined in this study; the size of each chromosome was taken from the GiardiaDB database (www.giardiadb.org).

The initial steps of this protocol are similar to those described for DNA combing. Exponentially growing trophozoites (∼5×10⁴ cells ml⁻¹) were sequentially incubated with 100 μM IdU for 60 min, centrifuged at 500 g for 10 min at 4°C, then treated with 100 μM CldU for another 60 min. The parasites (∼10⁵ total cells) were harvested by centrifugation at 500 g for 10 min, washed three times in PBS, and the pellet was resuspended in 150 µl of fixation solution [1% (v/v) Triton X-100 , 40 mM citric acid, 20 mM dibasic sodium phosphate, 0.2 M sucrose, pH 3.0] for 10 min. Diluent buffer (300 µl of 125 mM MgCl₂ in PBS) was added and the mixture was incubated for 20 min at room temperature. The samples were centrifuged at 500 g for 10 min, washed and resuspended in 100 µl of PBS. The cell suspension (20 µl) was loaded onto poly-L-lysine pre-treated microscope slides (Tekdon, Miakka City, FL, USA) and washed three times with PBS. The samples were treated with 2.5 M HCl for 20 min at 21°C to expose the incorporated thymidine analogs. Samples were washed and neutralized by incubating with 0.1 M borate buffer (100 mM H₃BO₃, 75 mM NaCl, 25 mM Na₂B₄O₇·10H₂O, pH=7.4) for 15 min at 21°C.

Immunodetection was performed using primary antibodies [mouse anti-BrdU antibody clone B44, (Becton Dickinson; Franklin Lakes, NJ, USA) and rat anti-BrdU antibody clone BU1/75 (Abcam, Cambridge, UK)] diluted 1:250 in 4% (w/v) BSA and incubated at 37°C for 2 h. The samples were incubated with secondary goat anti-rat-IgG conjugated to Alexa Fluor 488 (Life Technologies) and goat anti-mouse-IgG conjugated to Alexa Fluor 568 (Life Technologies, Carlsbad, CA, USA) antibodies, both diluted 1:500 in 4% (w/v) BSA, and incubated at 37°C for 2 h. Each antibody incubation step was followed by extensive washes with PBS. Vectashield mounting medium (Vectorlabs, Burlingame, CA, USA) containing DAPI was used as an anti-fade mounting solution, and to stain nuclei.

Images were acquired using the Olympus BX51 fluorescent microscope with a 100× oil objective attached to an EXFO Xcite series 120Q lamp and an Olympus XM10 digital camera (Olympus, Tokyo, JP). Images were further analyzed using ImageJ software (National Institute of Health, Bethesda, Maryland, USA) to count and establish the labeled patterns shown by trophozoites. Images were superimposed using the Cell F software (Olympus, Tokyo, JP) when necessary.

Exponentially growing trophozoites (10⁷ cells ml⁻¹) were harvested by centrifugation at 500 g for 10 min at 4°C, and the pellet was resuspended in the same volume of PBS with or without 300 µM of 5-bromo-2′-deoxyuridine (BrdU). The culture was incubated for 30 min at 28°C and harvested by centrifugation at 500 g for 10 min at 4°C. The BrdU-incorporated pellet was resuspended in 20 ml of PBS and harvested by centrifugation at 2000 g for 10 min. The pellet was resuspended in 4 ml of lysis buffer [200 mM Tris-HCl pH 8.0, 10 mM EDTA, 50 µg ml⁻¹ proteinase K and 0.5% (w/v) SDS] and incubated for 3 h at 50°C with periodic homogenization. DNA extraction was performed by adding one volume of phenol:chloroform:isoamyl alcohol (25:24:1) and the sample was gently homogenized by inversion. The solution was centrifuged at 1000 g for 10 min at 21°C and the aqueous phase was transferred to a new tube. The extraction was repeated using one volume of phenol:chloroform, followed by centrifugation, and the aqueous phase was transferred to a new tube. Sodium acetate (300 mM) was added, and the solution was homogenized by inversion. Two volumes of 100% ice-cold ethanol were added, and the sample was carefully homogenized by inversion until the formation of a ‘cloud’ of DNA. The DNA was collected with a glass hook and dipped three times into ice-cold 70% ethanol. The ethanol was air dried and the DNA was added to a tube containing 200 µl of 10 mM Tris-HCl pH 7.5 and 0.2 mg ml⁻¹ RNase A and incubated for 30 min at 37°C. All DNA quantifications were performed using the Qubit fluorometric double-stranded DNA High Sensitivity method (Life Technologies).

Unsheared high-molecular-mass DNA samples were treated with the Short Read Eliminator kit XL (Pacific Biosciences, Menlo Park, CA, USA) and prepared for nanopore sequencing using the 1D ligation-based library kit (SQK-LSK109) and the NEBNext^® Companion Module for Oxford Nanopore Technologies^® ligation sequencing (New England Biolabs, Ipswich, MA, USA) according to recommendations by Oxford Nanopore Technologies. Nanopore libraries were sequenced on a MinION Mk1B sequencer with Spot-ON flow cells (FLO-MIN106D R9 version; Oxford Nanopore Technologies, Oxford, UK).

Nanopore fast5 reads were base-called using Guppy base-calling software (version 4.2.2; Oxford Nanopore Technologies). Base-called reads were aligned to the G. lamblia WB genome assembly (Xu et al., 2020) using minimap2, version 2.17 (Li, 2018) and converted into bam format using samtools, version 1.9 (Li et al., 2009). Raw-, base-called- and mapped-nanopore reads were then used as input data for BrdU detection by DNAscent software (version 2.0.2), followed by analysis using forkSense to determine the direction of moving forks (considering 50–100% moving fork probability) (Boemo, 2021). Only reads with a mapping length ≥1000 bp and mapping quality ≥20 were used.

Mapped reads were fitted using a BrdU minimum incorporation rate of 5% across the whole read and a 350 bp window to determine the regions of potential conflicts between the replication and transcription machineries. ForkSense coordinates were then clustered and analyzed using scripts written in house available at https://github.com/davidsilvapires/potentialHoRTCollisions and the G. lamblia assemblage A genome annotation (release-28, GiardiaDB) because of chromosome correspondence (Warrenfeltz et al., 2018).

HoRT collisions occur when a replication fork encounters the transcription machinery head-on in the lagging strand at specific gene loci. Replication forks can move in opposite directions from the origin; therefore, genes might have more than one positive (lagging strand) or negative position (leading strand) for replication–transcription collisions. Genes identified with at least one positive position were considered to occur in a genomic region with a potential HoRT collision.

Gene Ontology analysis was performed on the GiardiaDB platform using gene IDs identified with potential collisions. Repeated gene IDs were removed, and Gene Ontology enrichment was performed through the GiardiaDB database (www.giardiadb.org) for all three ontologies: biological process, cellular component and molecular function. The same strategy was used for genes containing either the term ‘variant’ or ‘vsp’ in the product description field to evaluate the impact of HoRT on the expression of VSPs. Curated and computed evidence was considered, and enrichment analysis was not limited based on the GO Slim generic subset. The P-value cutoff was 0.05.

We thank Dr Conrad Nieduszynski (Earlham Institute, UK) for his assistance with DNAscent analysis.

The peer review history is available online at https://journals.biologists.com/jcs/lookup/doi/10.1242/jcs.260828.reviewer-comments.pdf.