Epigenomic Landscape of Lyme Disease Spirochetes Reveals Novel Motifs

Abstract

Borrelia burgdorferi, the etiological agent of Lyme disease, persists in nature through an enzootic cycle consisting of a vertebrate host and an Ixodes tick vector. The sequence motifs modified by two well-characterized restriction/modification (R/M) loci of B. burgdorferi type strain B31 were recently described, but the methylation profiles of other Lyme disease Borrelia bacteria have not been characterized. Here, the methylomes of B. burgdorferi type strain B31 and 7 clonal derivatives, along with B. burgdorferi N40, B. burgdorferi 297, B. burgdorferi CA-11, B. afzelii PKo, B. afzelii BO23, and B. garinii PBr, were defined through PacBio single-molecule real-time (SMRT) sequencing. This analysis revealed 9 novel sequence motifs methylated by the plasmid-encoded restriction/modification enzymes of these Borrelia strains. Furthermore, while a previous analysis of B. burgdorferi B31 revealed an epigenetic impact of methylation on the global transcriptome, the current data contradict those findings; our analyses of wild-type B. burgdorferi B31 revealed no consistent differences in gene expression among isogenic derivatives lacking one or more restriction/modification enzymes. IMPORTANCE The principal causative agent of Lyme disease in humans in the United States is Borrelia burgdorferi, while B. burgdorferi, B. afzelii, and B. garinii, collectively members of the Borrelia burgdorferi sensu lato species complex, cause Lyme disease in Europe and Asia. Two plasmid-encoded restriction/modification systems have been shown to limit the genetic transformation of B. burgdorferi type strain B31 with foreign DNA, but little is known about the restriction/modification systems of other Lyme disease Borrelia bacteria. This paper describes the methylation motifs present on genomic DNAs of multiple B. burgdorferi, B. afzelii, and B. garinii strains. Contrary to a previous report, we did not find evidence for an epigenetic impact on gene expression by methylation. Knowledge of the motifs recognized and methylated by the restriction/modification enzymes of Lyme disease Borrelia will facilitate molecular genetic investigations of these important human pathogens. Additionally, the similar motifs methylated by orthologous restriction/modification systems of Lyme disease Borrelia bacteria and the presence of these motifs within recombinogenic loci suggest a biological role for these ubiquitous restriction/modification systems in horizontal gene transfer.

Keywords: Borrelia; Lyme disease; epigenetic regulation; methylation; motifs.

FIG 1

Phylogenetic trees depicting homologs of the BBE02 (Pfam1) (A) and BBQ67 (Pfam2) (B) type II R/M systems of Borrelia burgdorferi B31. The trees were rooted by the outgroup (GenBank accession number YP_002344108). The CAT approximations of the Shimodaira-Hasegawa (SH) test were used to determine the local support values (SH support values) with 1,000 resamplings. The B. burgdorferi sensu stricto (s. s.) single nucleotide polymorphism (SNP) groups are indicated by the boxes next to the R/M protein GeneIDs. Genes belonging to B. burgdorferi B31 are indicated by stars, while additional Borrelia strains that underwent SMRT sequencing are indicated with triangles. Any truncations identified in GenBank and any indels identified through sequencing are indicated next to the strain name. The identified motifs are indicated in the rightmost column next to the R/M enzyme responsible for its methylation; however, if it is not known which protein is responsible for motif recognition, all motifs identified within that strain are shown in parentheses. The colored circles and similarly colored lines indicate orthologs whose clade matches the expected strain phylogeny. The sequences used to create the tree are located in Table S1 in the supplemental material (R/M homologs).

FIG 2

Heat map depicting the methylated motifs per genome and their interpulse duration (IPD) ratios. The horizontal heat map (red) indicates the detected motifs and the pairwise motif differences (based on the Pearson correlation coefficient), with the m6A-methylated bases underlined. The vertical bar plot demonstrates the number of R/M enzymes encoded within each Borrelia strain and derivative, identified to the right of the heat map. Ambiguity codes for nucleotides are as follows: H is A, C, or T; K is G or T; M is A or C; N is any base; R is A or G; V is A, C, or G; and Y is C or T. Ba, B. afzelii; Bg, B. garinii; Bb, B. burgdorferi.

FIG 3

Pearson correlation of the number of motifs per gene (per 1,000 bp) and the average number of motifs per 1,000 synonymously shuffled genes. An R value close to 1 with a significant P value (P < 0.01) indicates a correlation between the number of motifs and gene length. The regression line and confidence interval along with the R and P values for each motif of B. burgdorferi derivatives (A) and B. afzelii derivatives and B. garinii (B) are shown. Other than the HMAAG motif of B. afzelii PKo and the CAGC and CMAAYC motifs of B. garinii PBr, the majority of genes contain fewer motifs than their synonymously shuffled counterparts.

FIG 4

Heat map of methylated motifs per 1,000 bp within select genes of sequenced Lyme disease Borrelia spirochetes. The numbers of methylated motifs within each gene are displayed as the ratio of methylated motifs to the gene length per 1,000 bp so that each column (strain) totals 1. Each strain and derivative is displayed as a different color below the motif.

FIG 5

Pearson correlation of the number of motifs to the length of the replicon in B. burgdorferi strains B31, N40, CA-11, and 297 (A) and B. afzelii strains PKo and BO23 and B. garinii strain PBr (B). An R value close to 1 with a significant P value (P < 0.01) indicates a correlation between the number of motifs and the length of the replicon. The regression line and confidence interval along with the R and P values for each motif are shown. Only replicons that lie outside the confidence interval are labeled. The number of sites for each motif within their respective genomes significantly correlates with the length of the replicon.

FIG 6

bbe02 and bbq67 transcript levels of B. burgdorferi B31-A3 grown in vitro or within I. scapularis. cDNA of total RNA from B. burgdorferi B31-A3 in vitro or within I. scapularis at different stages underwent qPCR. bbe02 and bbq67 were expressed at lower levels within I. scapularis than in in vitro-grown organisms, with bbq67 being present at undetectable levels in larvae and feeding nymphs. The exception is the expression of bbe02 at nymphal repletion. Significance was determined by Dunn’s multiple comparison of the Kruskal-Wallis test (n = 3 biological and 3 technical replicates each). **, P value of <0.002; ***, P value of <0.0002.

FIG 7

Volcano plots and qRT-PCR graphs of differentially expressed genes. (A to C) Volcano plots of B. burgdorferi B31-A3 Δbbe02 (A) and B. burgdorferi B31-A3 Δbbq67 (B) gene expression compared to B31-A3 gene expression and B. burgdorferi B31-A3 gene expression compared to B31-A3 Δbbe02 Δbbq67 gene expression (C). Differentially expressed genes with a false discovery rate (FDR) of <0.01 as determined by DESeq2, those with a log fold change (|logFC|) of >2 as determined by DESeq2, and genes that have an FDR of <0.01 and a log fold change of >2 as determined by DESeq2 are shown. Genes that have an FDR of <0.01 and a log fold change of >2 in both DESeq2 and edgeR are labeled and shown in blue. (D and E) qPCR of B31-A3, B31-A3 Δbbe02, B31-A3 Δbbq67, and B31-A3 Δbbe02 Δbbq67 to verify decreases in bba25 gene expression within B31-A3 Δbbe02 (D) and bbh26 gene expression within B31-A3Δbbq67 (E) compared to B31-A3 from RNA-seq data. Decreased bba25 expression in B31-A3 Δbbe02 was confirmed by qPCR, but B31-A3 Δbbe02 Δbbq67 cells do not exhibit significantly lower expression levels of bba25. Decreased bbh26 expression in B31-A3 Δbbq67 was not confirmed by qPCR. Significance was determined by Dunn’s multiple comparison of the Kruskal-Wallis test (n = 3 biological and 3 technical replicates each). *, P value of <0.05; **, P value of <0.002; ***, P value of <0.0002.

FIG 8

Volcano plots of B. burgdorferi B31-A3 gene expression compared to B31-A3-68 (lp56⁻) gene expression (A) and B. burgdorferi B31-A3 gene expression compared to B31-A3-68 Δbbe02 (Δbbe02/lp56⁻) gene expression (B). Differentially expressed genes with an FDR of <0.01 as determined by DESeq2, those with a log fold change of >2 as determined by DESeq2, and genes that have an FDR of <0.01 and a log fold change of >2 in DESeq2 are shown. Genes that have an FDR of <0.01 and a log fold change of >2 by both DESeq2 and edgeR are labeled and shown in blue.

FIG 9

bbq67 isogenic deletion. To create the B31-A3 Δbbq67 and B31-A3 Δbbe02 Δbbq67 strains, the 3,261-bp bbq67 gene was replaced with the 919-bp flgB::aphA1 antibiotic resistance cassette in the same orientation as bbq67.