Strains of the same microbial species (conspecific strains) often present large phenotypic differences, despite having very similar genotypes (Van Rossum et al. 2020) . Previously published examples for phenotypic differences between strains of the same species include pathogenicity (Pierce and Bernstein 2016) , commensalism (Leimbach, Hacker, and Dobrindt 2013) , drug response (Maier et al. 2018) and susceptibility to infection by phages (Holmfeldt et al. 2007) . In host associated microbial communities, species can stably coexist for years (Faith et al. 2013; Lloyd-Price et al. 2017) , potentially evolving into host-specific strains (S. Zhao et al. 2019) . Occasionally, host specific strains could be transferred to or from other hosts, along familial networks (Yassour et al. 2018) , through the built environment (Brooks et al. 2017) or following fecal microbiota transplantation (Li et al. 2016; Smillie et al. 2018) . The ability to identify and follow conspecific strains is required to understand the mechanisms of between-host strain transmission, within-host evolution, and how these forces interact to shape microbial communities.

Methods to track strains using short-read data currently belong to one of two main classes: de-novo assembly of contigs from metagenomes, and methods relying on alignment of genomic sequences to a reference database (reviewed in (Anyansi et al. 2020) ). Methods in both classes usually rely on detection of single nucleotide polymorphisms (SNPs). In assembly-based methods, high sequencing depth is required to overcome sequencing errors and natural variation in the population, making these tools more suitable for the study of low-complexity microbial communities (Anyansi et al. 2020) . Moreover, identifying SNPs based on metagenomic assembled genomes (MAGs) can introduce errors in low quality MAGs (Van Rossum et al. 2020) . On the other hand, methods based on comparisons to a reference database usually require lower sequencing depth, although SNP detection in these methods can be limited by natural variation in the population and the degree of similarity between the community members and the reference genome (Bush et al. 2020) . Moreover, reference-based methods can only track strains belonging to well-studied species, for which a suitable reference database has been generated.

To complement methods relying on SNP information, we developed SynTracker, an approach to identify and track closely-related strains using genome microsynteny (the local conservation of genetic-marker order in genomic regions). Gene synteny (organization of genes along two chromosomes) has been used to estimate evolutionary distances between genomes (Lemoine, Lespinet, and Labedan 2007; Alexeev and Alekseyev 2017; T. Zhao et al. 2021) and to identify horizontal gene transfer events (Adato et al. 2015) . SynTracker uses pairwise comparisons of homologous genomic regions in either metagenomic, or genomic assemblies, followed by scoring the average synteny per pair of strains. SynTracker is relatively insensitive to SNPs, and requires only a single reference genome per species (either complete or draft and with no regard to its annotation level). Here, we apply SynTracker to compare the within-population synteny of in-silico evolved bacterial populations and to reproduce known within-species phylogenies of E. coli, using a fraction of the entire genome length. Additionally, we define a synteny-score cutoff to identify strains residing in the same hosts’ gut microbiomes over time. Finally, we apply SynTracker to a gut microbiome metagenomic dataset consisting of samples obtained from mothers and their infants ( Bäckhed et al. 2015 ), and describe a high degree of strain sharing between mothers and their infants in species which colonize the infant gut.

After verifying that the synteny method can be used to track host-specific strains in human gut metagenomes with high accuracy, we used this method to identify strain-sharing events between hosts. Our overarching goal was to study the role of vertical strain transmission (i.e., transmission from mother to infant) in the colonization of the human gut at early infanthood. Previously, vertical strain transmission was studied using culture-based techniques (Milani et al. 2015; Makino et al. 2013) , which are inherently limited to specific taxa. More recently, vertical strain transmission was studied by identifying SNP profiles in metagenomic data (Nayfach et al. 2016; Yassour et al. 2018) .

To study mother-infant strain transmission using strain synteny, we analyzed the dataset collected by Bäckhed and colleagues ( Bäckhed et al. 2015 ). These data contain stool-derived metagenomes obtained from 98 mothers and their infants, sampled at ages of 4 days, 4 months and 1 year post-birth. We assembled the metagenomic samples de-novo (see methods) and calculated the APSS scores, for a collection of 38 bacterial species (Figure 1). In order to maximize the number of samples included in per-species analyses, we used 30 regions per pairwise comparison.

We expected the gut microbiome of 4-day old infants to contain strains vertically transmitted from their mothers. Therefore, we predicted that true mother-infant-pairs (MIP) will show significantly higher APSS values compared to unrelated MIP, in comparisons of mothers and newborn infants. Therefore, for each analyzed species, we grouped mother-infant comparisons by the relatedness of the pair and compared the APSS values of the two groups (i.e., true MIP and unrelated MIP). We only considered species with at least 6 pairwise comparisons as suitable for statistical hypothesis testing (Wilcoxon rank test, single tailed, Benjamini-Hochberg multiple testing corrected (Benjamini and Hochberg 1995) ) . Only a small subset of 9 species passed our criterion in the newborn age group, which could be explained by the low complexity of the newborn gut microbiome and is in agreement with previous findings regarding the maturation of the infant gut microbiome (Koenig et al. 2011) . As expected, most species in this age group (7/9) had significantly higher APSS values in true MIP (q-value<0.05, figure 5). Additionally, 80/84 of the strain comparisons in true MIP had APSS value > 0.94, and therefore were considered as a vertical strain transmission, based on the APSS threshold described above.

We observed that in the 4- and 12- month age groups, strains of species identified in the newborn group remained similar to those of the mothers, in true MIP (adjusted p-value < 5*10-4 in the 12-months group, Figure 5). In contrast, late colonising species, identified only in later samples, did not have significantly higher APSS values in true MIP compared to unrelated MIP. Overall, in the 12-months old group, we compared 240 strains in true MIP pairs, out of which 126 could be considered as resulting from a strain-sharing event.

In this report we introduce SynTracker, a method for tracking closely related microbial strains using genome synteny. SynTracker requires as input a collection of genomes or per-sample assembled metagenomic contigs, a reference genome file, and a metadata file. The SynTracker code can be used as a standalone tool or as a part of a custom pipeline.

We designed SynTracker to complement other existing methods that track closely related strains, which mostly rely on SNP profiles or on analysis of specific sets of genes. In studies tracking strains across individuals, reliance on SNP information could be potentially limiting: environmental changes can spur the emergence of hypermutator strains (Travis and Travis 2002; Swings et al. 2017) , increasing the point-mutation rate by a factor of up to 150-fold (Wielgoss et al. 2013) . To avoid the limitations of these tools, we intentionally set the default pairwise-alignment parameters in our pipeline to have low sensitivity to SNPs, which may be erroneously identified due to sequencing and amplification errors. When examining SynTracker’s performance using in-silico evolved bacterial populations, we observed that, as expected, populations that evolve exclusively through introduction of SNPs had a marginal reduction in the synteny scores compared to populations that evolved through introduction of insertions and deletions at a lower mutation frequency. This characteristic of our approach makes it a good candidate to complement existing SNP-based tools and also makes it ideal to track closely related strains in data produced using long-read sequencing methods, as the error rate of these methods is higher than in methods based on short-reads (Amarasinghe et al. 2020) .

While some popular tools for conspecific strain analysis require a pre-existing gene database, our approach requires only a single reference genome per species, either fully assembled or as a collection of contigs. This feature is advantageous, as it allows for tracking strains of relatively understudied species. As state-of-the-art methods for assembly of genomes from metagenomes yield ever larger collections of MAGs, we propose a potential workflow, in which the MAG collection is clustered to create “species representative genomes”, which could be used as the reference genome in our pipeline. This approach, which is also utilized in the InStrain program (Olm et al. 2021) , can expand our ability to study strains of novel species.

One of the most important assets of our approach is its ability to track closely related strains using only a small fraction of the full length of the genome. We were able to reconstruct the phylogeny of 140 E.coli strains using <20% of the length of the reference genome and to identify within-host strains with an average sensitivity of 97% and specificity of 99%, using the same accumulative length of the compared regions. The ability to track strains using a fraction of the full genome length is especially important when analyzing MAGs with low completeness values or less abundant taxa in metagenomic assemblies.

We examined the performance of our method by identifying strains residing in the same human hosts over time periods of a few weeks to two years. We observed a decrease in the performance of the approach with the reduction in the number of regions used per pairwise comparison. On the other hand, reducing the number of regions could increase the number of samples included in the final analysis. The SynTracker pipeline provides a number of average-synteny score tables, prepared using 20-200 regions/pairwise-comparison. It is up to the user to select the relevant table, based on their specific needs and dataset.

When investigating low abundance taxa, metagenome assembly might yield relatively short contigs. In those instances, the likelihood of identifying a sufficient number of overlapping 5 kbp regions in any given two metagenomic assemblies is reduced. In such cases, it could be beneficial for the user to perform the synteny-based analysis on shorter genomic regions. This could be easily achieved by reducing both the length of the flanking-regions and the spacing between the “central-regions” (Figure 1A, methods).

To demonstrate the use of SynTracker we analyzed the metagenomic dataset collected by Bäckhed et al (Bäckhed et al. 2015 ), who followed a cohort of mothers and their infants from birth to one year of age. Since SynTracker uses pairwise comparisons of homologous genomic regions, the number of pairwise comparisons increases exponentially with increasing numbers of samples. To reduce the overall running time of our pipeline without losing relevant information (i.e., comparisons of true MIP and comparisons of longitudinal samples), we divided the dataset into 20 bins, while keeping all same-family samples in the same bin. Using this strategy, we were able to reduce the number of pairwise comparisons by a factor of ~21 at the cost of losing some between-family comparisons, which were only used in our analysis as a control group, relative to the true MIP group. We strongly recommend this strategy to researchers analyzing larger datasets consisting of hundreds of metagenomic samples and dozens of reference genomes.

Our analysis of the Backhed et al. dataset showed that early colonizing species, that inhabit the guts of both the mothers and the newborn infants, had higher APSS values in comparisons of MIPs, compared to unrelated MIPs. Moreover, a striking majority of the strains tracked in true MIPs (newborn age group) had APSS scores high enough to be considered within-host strains. In the 12-month-old group, early colonizing species maintained the higher APSS values in true MIPs, compared to unrelated MIPs, while no significant difference was observed for late colonizing species. These results suggest that vertical strain transmission plays a role in the acquisition of early colonizing species, while late colonizers could be obtained from additional sources as well.

We have introduced SynTracker, a tool for tracking conspecific strains and to evaluate their relatedness using genome synteny, in both genomes and metagenomic assemblies. To our knowledge, this is the first tool which is entirely based on this level of biological organization. SynTracker’s attractive features include that it does not require pre existing databases, and has a minimal sensitivity to sequencing errors and natural variation in microbial populations. SynTracker performs well when classifying isolate genomes and when tracking strains in longitudinal metagenomes. SynTracker could be used as a standalone tool or combined with existing tools in a multi-tool pipeline setup. SynTracker is available at: https://github.com/leylabmpi/SynTracker Acknowledgements: We thank Nick Youngblut for providing comments on a previous version of this work. phylogroups, based on Average Pairwise Synteny Scores (APSS). A. A tree based on 200 regions/pair (accumulative length of ~1 mbp) correctly classified 139 genomes. B. A tree based on 40 regions/pair shows correct classification of 136 genomes.

The SynTracker pipeline consists of three main parts. In the first part, SynTracker accepts a collection of reference genomes (a single genome per species), either fully assembled or as a collection of contigs. Each per-species reference is fragmented into a collection of 1kbp central-regions, which are binned and stored together. In the second part SynTracker creates a blast Database, based on a user-provided collection of metagenomic assemblies or genomes. Next, it performs a blast search, for each of the central regioins, against the newly created blast database with a minimal identity of 97% and a minimal query coverage of 70%, i.e., 700bp. In the final step of this part, hits for each blast search are retrieved, using the blastcmddb command, in addition to a 2kbp region on each side of the blast hit. Hits with <2kbp both downstream and downstream to the hit are excluded from further analysis. Each retrieved sequence is denoted by its sample of origin and matching region in the reference genome. In the third part of the pipeline genomic fragments are grouped by their matching region in the reference genome, and pairwise alignment is conducted to identify synteny blocks in each pair of sequences. The identification of synteny blocks in each pairwise alignment is performed using the “FindSynteny” function, in the “DECIPHER” R package (Wright 2016) , with parameters “maxGap” and “maxSep” both set to 15. Additionally, only pairwise comparisons with a minimal overlap 4800 bp are considered for downstream analysis. Next, per each pairwise alignment, a synteny score is calculated, as described in equation 1: Eq. 1

SynScore=1+ log10〖(Ov/len)/B〗 Where Ov stands for the accumulative length of the overlapping synteny blocks identified in the pairwise alignment, len denotes the length of the shorter sequence in each pair and B stands for the number of synteny blocks identified in each pairwise alignment.

In the final step of the third part of the pipeline, for each reference genome n genomic regions are randomly selected, per pair of metagenomic samples or genomes. APSS (average pairwise synteny scores) are calculated by averaging the individual pairwise synteny scores. Pairs of samples/genomes with <n regions are excluded from downstream analysis.

Calculation of the synteny scores per group of sampled cells was performed as described above, however, as the length of the genomic fragments used in the simulation was limited to 20kbp, synteny scores were based on a single alignment of the ~20kbp region, per pair of simulated genomes.

Calculation of APSS values for E.Coli strain pairs was performed as described above and using the E.coli str. K-12 substr. MG1655 genome as a reference (NCBI Reference Sequence: NC_000913.3).

Phylogenetic trees were generated by conversion of APSS values to synteny distances, which equal to 1-APSS. All pairwise synteny distances were placed in a symmetric matrix which was used to calculate UPGMA phylogenetic trees, by employing the “phangorn” R-package (Schliep 2011) .

Longitudinal metagenomes were obtained from the NCBI-SRA database, and were quality filtered as described previously (Youngblut et al. 2020) . Metagenomic samples were de-novo assembled using metaSPades (Nurk et al. 2017) , with a maximal number of 20M reads/sample. ROC curves and matching APSS thresholds for each combination of species and sampling depth, in the testing set, were calculated using the R-programing language “pROC” package (Robin et al. 2011) .

Metagenomic samples were downloaded from the NCBI-SRA database, and were quality filtered and assembled as described above, however, as for some samples only one of the two matching read files passed our quality filtration, we performed the metagenomic assembly using single-end reads. genomic regions per pairwise comparison.

Specificity (%) 15 20 30 40 60 80 100 200

NC_016077.1 NZ_CP021420.1 NZ_CP021420.1 NC_018011.1 NZ_AP019734.1 NZ_DS499581.1 NC_021030.1 NZ_CP012801.1 NZ_UFSX01000001.1 NC_003228.3 NZ_KB905475.1 NZ_SPFU01000010.1 NZ_JH724307.1 NZ_CP012937.1 NZ_CZAF01000001.1 NC_009614.1 NZ_RCXZ01000001.1 NZ_JH815206.1 NZ_CP028341.1 NZ_AKCA01000001.1 NZ_CP026999.1 NZ_CYZN01000001.1 NZ_CP048433.1