The strains analyzed in this study are listed in Table 1 . The paired-end mHVR amplicon sequences were obtained by Rusman et al. (2019). Reads are available at the Sequence Read Archive (SRA) under the BioProject ID: PRJNA514922. Raw sequence reads for all samples were pre-processed. First, reads were quality filtered using the pair-end mode of Trimmomatic v0.36 (Bolger et al., 2014). Posteriorly, the preserved paired reads were merged into consensus sequences with their associated quality score and chimeras were removed using LeeHom software with default parameters (Renaud et al., 2014).

Predicted edited mRNA sequences for strains Sylvio (mitochondrial clade A, TcI), CL Brener (mitochondrial clade B, TcVI) and Esmeraldo (mitochondrial clade C, TcII) (Ruvalcaba-Trejo and Sturm, 2011) were generated by manually inserting or deleting uridines (U’s) into the unedited sequences of their genes as proposed by (Greif et al., 2015). Briefly, unedited sequences corresponding to Sylvio (FJ203996.1), CL Brener (DQ343645.1) and Esmeraldo (DQ343646.1) obtained from the NCBI database were aligned using MEGA v7 software (Kumar et al., 2016). Then, the alignment was manually manipulated by inserting or deleting uridines (U’s) at specific sites into the unedited sequences following the editing patterns of the corresponding T. brucei and T. vivax edited mRNAs, meanwhile the resulting amino acid sequence was preserved. Predicted sequences were generated for mitochondrial unidentified reading frame 2 (MURF2), NADH dehydrogenase subunit 3 (ND3) and C-rich region 4 (CR4) following the known edited mRNAs for T. vivax and T. brucei ( Supplementary File 1 ). On the other hand, sequence data of predicted edited mRNAs for the following genes were obtained from Ruvalcaba-Trejo et al. (Ruvalcaba-Trejo and Sturm, 2011) for the same three T. cruzi strains: ATPase6, C-rich region 3 (CR3), cytochrome b (CyB), ribosomal protein S12 (RPS12), NADH dehydrogenase subunit 7 (ND7), NADH dehydrogenase subunit 8 (ND8) and NADH dehydrogenase subunit 9 (ND9). CL Brener cytochrome oxidase subunit 3 (COIII) mRNA was downloaded from GenBank (accession number: EF058194.1). In addition, Sylvio and Esmeraldo COIII edited mRNA sequences were manually generated ( Supplementary File 1 ). Additional maxicircle sequences described by Reis-Cunha et al. 2018) were downloaded from GenBank, manually edited, and also used for gRNA search. Predicting gRNAs based on such sequences (shown in Table S1 ) had similar results to the obtained by using Sylvio, Esmeraldo or CL-Brener edited mRNAs.

An algorithm using c++ and the SeqAn library (Reinert et al., 2017) was built to detect potential gRNA in the mHVR reads. First, a hash-table was built from the reads by storing position of all the possible k-mers (seeds) of 15 bases. However, due to non-canonical base pair a hash-table was used considering just purines (R) and pyrimidines (Y) instead the four bases. Although C::A pairing is allowed with this approach, it is considered a mismatch in the following steps. An algorithm scheme is provided in Figure S1 . The hash table only allows to reduce the search time in the next steps. The mRNA sequence for each edited gene was used to identify gRNAs. The estimated edited mRNA of the corresponding mitochondrial clade was used to infer gRNA for each strain. First, a window of 15 bases was selected at the 5’ end in the mRNA and the reverse complement of the purine/pyrimidine pattern was used to determine which mHVRs contains a potential gRNA by using the hash-table. Then, the potential site was evaluated considering canonical and non-canonical base pairing (G::U) with the mRNA and allowing only one mismatch. Finally, the pairing is extended to both sides using an x-drop extension allowing one mismatch. Mismatches were allowed considering some flexibility in editing and considering that mRNAs are estimated and minimal errors in estimations are possible (Thomas et al., 2007). A filter was then applied to discard the predicted gRNA if it has fewer bases than the parameter f. This parameter is the minimum number of bases that a gRNA candidate must have to be accepted. This parameter is required because short sequences with complementarity to mRNAs are expected by simple random. Then, the window is moved over the mRNA one base at a time and the process is repeated until the entire mRNA is scanned. Different values for f were evaluated in order to discard sequences that are not gRNAs. In a first approach, gRNAs with <40 bases were discarded (f = 40). This filter was too stringent but allowed to identify the most probable regions of the mHVR that contain the gRNA, that is the gRNA cassette. Then, the filter was lowered looking for hits outside the cassette and the minimum filter that had as maximum as 5% of the gRNA hits out of the cassette was selected (f = 30) ( Figure S3 ). Alternatively, other approaches for gRNA detection were evaluated. The first one was based on (Koslowsky et al., 2014). Briefly, after looking for 15 bp seeds, the algorithm scored canonical (score = 2) and non-canonical (score = 1) base pairings and gRNAs with scores <45 were discarded. The second approach was based on (Cooper et al., 2019) by looking for gRNAs with 6 bp-anchors (canonical base pairing) in the 5’ end and then discarding gRNAs with less than 25 bp (canonical plus non-canonical base pairing). Both approaches were too stringent, they recovered reduced gRNA coverage in ATPase 6 mRNA and a smaller number of gRNAs than using the above approach f = 40 ( Figure S1 ). Consequently, f = 30 was established given its higher sensitivity than f = 40, and a reasonable rate of gRNA hits out of the cassette. The c++ source code for gRNA detection in mHVRs and example datasets are stored in https://osf.io/kn34z/.

Identified gRNAs were clustered according to sequence similarity and overlapping in the edited-mRNA target region. A SWARM-like algorithm (Mahe et al., 2014) was implemented to cluster gRNAs. SWARM is a fast unsupervised (de novo) single-linkage-clustering algorithm. An algorithm scheme is provided in Figure S2 . The modified algorithm also included an overlapping threshold (o) in addition to a threshold of the maximum number of base mismatches between sequences (d). Basically, the algorithm selects one gRNA as a cluster seed and then it iteratively looks for gRNAs with a number of mismatches lower than d and an overlapping > o and adds them to the cluster. Then, it iteratively looks for unclustered gRNAs with less than d mismatches and overlapping > o to every gRNA previously included in the cluster. When no new gRNAs can be added to the cluster a new cluster is seeded using a gRNA that was not previously clustered. The process is iterated until the all gRNAs were joined to a cluster. o = 0.7 and d = 4 were implemented. This is a relaxed threshold to cluster moderately different gRNAs in the same cluster considering divergency in such kind of sequences. Considering the overlapping region of two gRNAs of ~30-50bp, d = 4 implies 85-92% pairwise similarity. Similar percentages were used to cluster mHVRs in (Rusman et al., 2019). The second phase of breakage in the original SWARM algorithm was not required (Mahe et al., 2014) because of the high distances and reduced overlapping among clusters. Identified gRNA clusters based on SWARM were termed as gRNA classes here. Finally, a consensus gRNA was built for every gRNA class to represent editing cascades. In order to fast calculation, the consensus gRNA was built using a base-by-base majority-rule criterium (i.e. the most abundant base in each position of gRNA alignment was used to build the consensus gRNA). Because artifactual gRNA may be generated by this method the algorithm also calculated the number of gRNAs sharing such sequence and informed it in the output file to manual inspection. The consensus gRNAs were graphed aligned to the edited mRNA and the cluster abundance was represented with a color scale.

In order to assess linkage disequilibrium as a measure of genetic structuring among strains, the index of association (Ia) and the proportion of compatible pairs of loci with strict clonality (PrCompat) were calculated. Briefly, gRNA classes were coded as present (>20 reads) or absent (≥20 reads) emulating RFLP pattern data for each strain and the coded matrix was analyzed as haploid data in Mulilocus v1.3 (Agapow and Burt, 2001) with 10,000 iterations to calculate statistical significance.

A subset of gRNA classes was selected to predict secondary structure. Briefly, gRNA classes with more than 200 reads (high abundance), with more than 90% of the reads being identical to the consensus sequence of the class and with at least 40 bases were selected. The selected class representative gRNAs were analyzed using RNAstructure software (Reuter and Mathews, 2010). The maximum accuracy structures were calculated using default parameters but considering a temperature of 27°C (optimal to temperature for epimastigotes and previously used in T. brucei analysis of gRNA secondary structure) (Schmid et al., 1995).

mHVR reads were processed and clustered according to (Rusman et al., 2019) using a UCLUST algorithm with a threshold of 85% of identity for cluster definition. Clustering at 85% was previously used to define mHVR clusters (Rusman et al., 2019) and such clusters were used to determine relationships between mHVR clusters and gRNA classes. Briefly, a c++ algorithm was designed to add a string with the mHVR cluster ID to each read name in the raw fastq file. Then, this fastq file was used to determine gRNAs using the above-described algorithm. This way, each identified gRNA could be assigned to an mHVR cluster. After clustering gRNAs in classes using the SWARM algorithm, it was evaluated whether each mHVR cluster coded none, one or more than one gRNA. Around 9.5% of the clusters corresponded to two gRNA classes in the same mHVR cluster. Although it is expected none or just one gRNA in an mHVR, such possible method inaccuracy was allowed in order to gain sensitivity to obtain the maximum coverage on mRNA editing cascades.

Guide RNAs coded in mHVRs were inferred based on the available sequences or predicted edited mRNAs of T. cruzi. A total of 7,476,003 gRNAs were detected for nine T. cruzi strains belonging to the six DTUs. The percentage of detected gRNAs from the total number of mHVR reads for each studied strain varied from 38% to 60% ( Figure 1A ), suggesting that some mHVRs would not encode gRNAs. Posteriorly, the gRNAs were clustered in gRNA classes according to the editing region and their sequence similarity ( Supplementary File 2 ). A total of 1334 gRNA classes were identified (every class with >20 reads). The number of gRNA classes per strain was variable ranging from 104 (LL014R1-TcV) to 401 (Esmeraldo-TcII). Particularly, TcV strains had up to three or four times fewer gRNA classes than other DTUs ( Figure 1 ). The number of gRNA classes was also variable within DTUs ( Figure 1 and Figure 2 ) suggesting fast changes in gRNA repertoires and probably gRNA redundancy (i.e., different gRNA classes editing the same mRNA region). Rarefactions of each dataset showed that these differences between strains are not the effect of different sequencing depths ( Figure 1B ).

Just 22.8% of the whole identified gRNA classes in all strains were shared among two or more of them and only four gRNA classes were shared among all strains (sequences provided in Supplementary File 2 ) revealing high divergence. However, the percentage of shared gRNA classes between strains of the same DTU was variable ( Figure 2 ). TcI strains shared a relatively low number of gRNA classes (19.2%) despite that the studied strains are very close phylogenetically (according to nuclear genes) (Schmid et al., 1995) and they were isolated in the same geographical area and time (Diosque et al., 2003). In addition, TcII strains shared even few gRNA classes (9.1%) ( Figure 2 ). Instead, TcV strains (LL014R1 and MNcl2) shared a high proportion of gRNA classes (73.8%) despite they were isolated in geographically distant places and with a difference of around 30 years between the isolation dates. All these results suggest significant changes in gRNA repertoires in relatively short evolutionary times. In addition, linkage disequilibrium was analyzed to address the structuring of such diversity. The index of association (Ia) was 49.3 and showed high significance (p < 0.001). Furthermore, the ratio of pairs of gRNA classes compatible (PrCompat) with unmixed gRNA classes was 0.98 suggesting structuring of diversity among strains. Finally, it was addressed whether probable predicted secondary structures are common in such gRNA sequence diversity. A subset of 89 gRNA classes defined by a stringent criterion (see Material and methods) was studied. Although variable structures were observed ( Supplementary File 4 ), 70.7% of the gRNAs have a stem region intercalated by one or two internal loops (or less frequently bulges) and ending in a hairpin loop ( Figure S4 ). Only 26% of the structures have a double hairpin ( Figure S4 ).

In a previous work, entire mHVR sequences were clustered according to sequence identity (Rusman et al., 2019). Here, the gRNA classes were assigned to mHVR clusters. First, most (97.1%) mHVR clusters of all strains had one of two states: (i) high proportion (> 75%) of the sequences in the mHVR cluster coding one gRNA class; or (ii) low proportion (< 10%) of sequences in the mHVR cluster coding a gRNA class ( Table S2 ). Most exceptions to the rule were observed in TcV strains ( Table S1 ). The term “silent” mHVR clusters was used for those in which less than 20 gRNAs were identified. Silent mHVR clusters were from 25% to 54% of the total mHVR clusters in the studied strains ( Table 2 ). Curiously, strains with the lowest number of gRNA classes (TcV and TcVI strains) also had the lowest percentage of silent mHVR clusters ( Table 2 ).

Furthermore, the analyses showed there is no clear relationship between mHVR cluster size and whether it encodes a gRNA or not since silent mHVR clusters were of quite different abundances ( Figure 3 ).

The gRNA classes identified for each strain were aligned with the edited mRNAs of the corresponding mitochondrial clade. The number of gRNA classes editing each gene for all strains is represented in Figure 4 and Table S3 . Most gRNA classes (more than 50%) edit ATPase 6 and COIII mRNAs in all strains ( Figure 4A ). Interestingly, the percentage of gRNA classes editing each mRNA was nearly conserved among strains ( Figure 4B ). This implies correlated variations in the number of gRNAs classes editing each mRNA (see correlation coefficients in Table S4 ). These results show that observed variations in the number of gRNA classes are balanced, i.e., if one strain decreased gRNA classes for one mRNA also decreased the number of gRNAs classes for editing other mRNAs. This would be an expected result whether each minicircle codes for different gRNAs, each one editing different mRNAs.

The algorithms identified a complete or almost complete set of gRNA classes necessary to edit the mRNAs for the ATPase6 (complex V) and COIII (complex IV) in every strain despite the highly variable number of gRNA classes editing such mRNAs. This clearly shows redundant gRNA classes for some strains. The editing mRNA cascades for ATPase6 and COIII pan-edited genes are shown in Figure 5 and Figure S5 , respectively. Aligned gRNAs-mRNAs are provided in the Supplementary File 2 . Strains PalDa20cl3, Esmeraldo, X109/2 and CANIIIcl1 showed a greater redundancy of gRNA classes than other strains which implied strong variations among and within some DTUs, mainly TcII strains ( Figure 5 and Figure S5 ). On the other side, TcV strains - MNcl2 and LL014R1 - had the lowest redundancy of gRNA classes ( Figure 5 and Figure S5 ).

The percentage of coverage by gRNAs for each mRNA is shown in Figure 4B and Table S5 . Interestingly, the coverage was almost complete for ATPase6 and COIII mRNAs but highly variable for the remaining mRNAs, especially in the complex I members ( Figure 4B ). The editing cascades of the mRNAs coding complex I members —ND9, CR3, CR4, ND3, ND7, ND8— of the respiratory chain are shown in Figure 6 and Figures S6–S10 . Even strains showing high mRNA class redundancy in COIII and ATPase6, such as Esmeraldo and X109/1, have incomplete cascades for complex I members. The exception was ND9 mRNA, in which both strains (X109/2 and Esmeraldo) presented almost a complete coverage ( Table S5 ), although most of the gRNA classes (82.1% and 74.8%, respectively) had less than 20 reads (see also Table S6 ). Breaking points in ND3, ND7 and ND8 editing cascades in these two strains were also observed. Interestingly, TcV strains —LL014R1 and MNcl2— presented relatively few gRNA classes that edit the complex I components messengers ( Table S3 , Figure 6 and Figures S6–S10 ) and relatively low coverage of the mRNAs ( Table S5 and Table S6 ).

Biparental inheritance of minicircles in hybrid DTUs was previously proposed (Rusman et al., 2019). Here, it was analyzed whether gRNAs from one strain may potentially edit mRNAs from a strain from another DTU, particularly another mitochondrial clade. This ability for “cross-editing” is a requirement in hybrid DTUs, where minicircles come from both parentals but maxicircles come from just one parent. Editing cascades were again generated for Esmeraldo gRNAs and X109/2 gRNAs but based on mRNAs of different mitochondrial clades. mRNA coverage by gRNAs was evaluated and both gRNA repertoires have similar coverages on mRNAs from different mitochondrial clades ( Figure 7 ). Cross-editing was also evaluated in more stringent conditions by lowering the allowed number of mismatches in cluster definition (d = 3, d = 2 and d = 1) and it was not detected a significant difference in mRNA coverage percentage. Minor differences were observed in the mRNA coverage percentages for genes CyB and ND9 in d = 1 settings ( Figure S11 ). These results suggest that biparental inheritance of the kDNA would not generate mRNA editing problems even among phylogenetically distant DTUs.