Virulence of the Pathogen Porphyromonas gingivalis Is Controlled by the CRISPR-Cas Protein Cas3

Porphyromonas gingivalis is a key pathogen of periodontitis, a polymicrobial disease characterized by a chronic inflammation that destroys the tissues supporting the teeth. Thus, understanding the virulence potential of P. gingivalis is essential to maintaining a healthy oral microbiome. In nonoral organisms, CRISPR-Cas systems have been shown to modulate a variety of microbial processes, including protection from exogenous nucleic acids, and, more recently, have been implicated in bacterial virulence. Previously, our clinical findings identified activation of the CRISPR-Cas system in patient samples at the transition to disease; however, the mechanism of contribution to disease remained unknown. The importance of the present study resides in that it is becoming increasingly clear that CRISPR-associated proteins have broader functions than initially thought and that those functions now include their role in the virulence of periodontal pathogens. Studying a P. gingivalis cas3 mutant, we demonstrate that at least one of the CRISPR-Cas systems is involved in the regulation of virulence during infection.

observed an 8-fold increase in expression of the gene after 6 h of incubation in the intracellular bacteria, while no increase in expression was observed in planktonic growth after 2 and 6 h of incubation (Fig. S1b). Next, we assessed whether there were changes in the growth of the wild-type and mutant in planktonic growth and intracellularly in THP-1 cells. Our results showed no significant differences in either planktonic growth (Fig. S2a) or intracellular growth and no changes in the persistence of either P. gingivalis genotype in the two environments (Fig. 1). In our analysis, we distinguished between intracellular P. gingivalis cells obtained after antibiotic treatment and the fraction of P. gingivalis attached to the surface of THP-1 (cell-associated) plus intracellular P. gingivalis obtained after cell washing with no antibiotic treatment. We observed no effect of cas3 deletion on growth rates of P. gingivalis under any of those conditions (Fig. 1).
Deletion of cas3 leads to significant changes in expression profiles of P. gingivalis when intracellular in THP-1 cells and an increase in activities associated with pathogenesis. To begin determining the potential role of the Cas3 protein in the overall metabolism of P. gingivalis, we performed a transcriptome analysis of intracellular bacteria in THP-1 cells and in a human serum-based medium broth. The number of sequences in the THP-1 infection and planktonic growth experiments ranged from 49,592,338 to 82,270,371 sequences.
P. gingivalis growing in the serum-containing medium did not show any significant differences in gene expression profiles between wild-type and mutant strains. Principalcomponent analysis (PCA) of the complete profiles of gene expression showed similar patterns for the wild-type strain and the mutant (Fig. 2). In contrast to the results observed in planktonic growth, the patterns of the two strains growing intracellularly were very different (Fig. 2), suggesting that cas3 controls certain aspects of the metabolism of P. gingivalis when the bacterium is residing inside a host cell but has little to no significant effect when the organism is growing planktonically. More importantly, PCA showed that the expression profiles of the wild-type and mutant strains, growing intracellularly, clustered tightly based on the strain of origin but were very different between strains (Fig. 2) (see Table S1 in the supplemental material). Comparisons of the wild-type transcriptome to the mutant transcriptome in P. gingivalis in the planktonic and THP-1 infection experiments showed differential expression of many genes: 150 differentially expressed genes in the planktonic experiment and 814 in the THP-1 infection experiments (Table S1). Among the most highly upregulated genes in the Δcas3 mutant in the THP-1 infection experiments, we found the genes for rubrerythrin, which plays a role in the oxidative stress response of P. gingivalis (23); for an ATPase; for an XRE family transcriptional regulator; and for a hypothetical protein. Interestingly, several Cas-associated protein genes (cas6, cas1, cas2, and cmr2) were found to be upregulated and a large number of mobile element proteins, 56 in total, were also upregulated in the cas3 mutant compared to the wild-type strain, among them a large number of loci of transposase in intracellular serine protease g1 (ISPg1), whose expression is induced when P. gingivalis is treated with H 2 O 2 (24) (Table S1). Additional genes that have been associated with responses to oxidative stress and iron acquisition in P. gingivalis were also upregulated, including those encoding ferritin (25,26), DNA-binding protein from starved cells (Dps) (25), ferredoxin (25), vimA (virulence modulating gene A) (27), putative universal stress protein UspA (28), hemagglutinin protein HagA (29), heme-binding protein FetB (29), and "upregulated in stationary phase protein A" (30).
Two genes with notably elevated expression were linked to the following proteins involved in modulation of the immune response: fig|431947.7.peg.735, which has been described as a gingipain-sensitive ligand A protein (GslA) (31), and the immunoreactive 46-kDa antigen fig|431947.7.peg.1742. Gingipain RgpA (fig|431947. 7.peg.1936), an essential virulence factor in P. gingivalis, was significantly upregulated in the cas3 mutant. Downregulated genes in the cas3 mutant included the Tra protein genes of the conjugative transposons (traE, traF, traI, traJ, traK, traL, traN, traO, and traQ) and a large number of genes encoding ribosomal proteins both from the large subunit of the ribosome (LSU) and from the small subunit of the ribosome (SSU) ( Table S1).
Gene set enrichment analysis (GSEA) of the differentially expressed genes showed that enriched biological processes were linked to pathogenesis and proteolytic activities in the P. gingivalis Δcas3 mutant under conditions of growth inside the THP-1 cells (Fig. 2b, top panel). The same analysis showed enrichment of molecular functions associated with peptidase activity and binding to metal ions (Fig. 2b, bottom panel).
The cas3 mutant induced changes in gene expression of THP-1 cells associated with immune response. In parallel, we performed transcriptome analysis on the THP-1 cells infected with the two strains of P. gingivalis. PCA showed no apparent clustering of the expression profiles based on whether the THP-1 cells were infected with the wild-type or the mutant strain (Fig. S3a). We found that 2,526 genes were differentially expressed in comparisons of THP-1 cells infected with the wild-type and Δcas3 mutant strains (Table S1). Intriguingly, the most abundant differentially expressed transcripts, all of them upregulated, were either microRNAs (miRNAs) or novel transcripts of unknown function (Table S1). A large number of the microRNAs identified in our study were previously linked to periodontal disease, including microRNA 130b, microRNA 19a, microRNA 20a, microRNA 27a, microRNA 302b, microRNA 30e, microRNA 374a, mi-croRNA 374b, microRNA 593, microRNA 29b-1, and microRNA 21 (32,33).
The upregulated genes were mainly associated with cytoskeleton organization (actin-related proteins and ankyrin), cell death (BCL2 interacting proteins, caspase 9, and CASP8-like apoptosis regulators), and immune response (chemokine ligands and receptors; NF-B inhibitors; interferon alpha and gamma receptors; interleukin 1 beta [IL-1␤] and interleukin 1, 10, and 9 receptors; TEC proteins involved in intracellular signaling mechanisms of cytokine receptors; Toll-like receptor [TLR-2], TLR-4, and TLR-8; tetrahydrofuran [THF] alpha receptor; and tumor necrosis factor alpha [TNF-␣]). Table S1 shows fold changes of the specific differentially expressed genes. Only one gene, the RNA5S9 gene, was found to have been downregulated when the cells were infected with the Δcas3 mutant strain (Table S1). When we performed gene set enrichment analysis (GSEA) on Gene Ontology (GO) terms for those genes, we observed enrichment in genes associated with the immune response. Enriched biological processes were associated with neutrophil and leukocyte chemotaxis and also with gene silencing (Fig. 3), while enriched molecular functions were again found to be associated with gene silencing and chemokine and cytokine activities (Fig. S3b).
⌬cas3 mutant had a modest effect on proinflammatory cytokine and chemokine production in the initial immune response of THP-1 to P. gingivalis. As inflammation is a crucial driver of the bacterium-elicited tissue destruction that char- acterizes periodontal disease, we assessed the effect of mutating cas3 on the cytokine profiles of THP-1 macrophage-like cells. Using Luminex multiplex immunoassay, we measured cell culture supernatant fluid levels of the cytokines TNF-␣, IL-1␤, IL-6, IL-10, RANTES, and IL-8 at 6 h and 24 h following infection. Overall, we observed a slight increase in the levels of IL-1␤, IL-6, and IL-10 in cells infected with the cas3 mutant compared to those infected with the wild-type strain, especially at 6 h, while RANTES levels showed decreases compared to the control, suggesting a potential inhibitory effect (Fig. 4). TNF-␣ expression levels were significantly higher in the mutant-infected cells than in the wild-type-infected cells (Fig. 4). The differences in the levels of the cytokines that we evaluated did not reach statistical significance. Intriguingly, the expression of cas3 intracellularly had its peak at 6 h as shown by quantification of the expression by real-time quantitative PCR (RT-qPCR) (Fig. S1b).
Deletion of cas3 increases bacterial virulence of P. gingivalis ATCC 33277 in the Galleria mellonella killing assay. To directly assess the contribution of cas3 in virulence of P. gingivalis, we challenged groups of 15 G. mellonella larvae by injection with different dilutions (from 10 6 to 10 8 CFU) of the different strains. We found significantly higher mortality (P Ͻ 0.0001) for the groups of worms infected with the Δcas3 mutant, with 100% mortality within the first 24 h when larvae were injected with 2.15 ϫ 10 8 CFU (Fig. 5a) and 90% mortality after 40 h when larvae were injected with 2.15 ϫ 10 7 CFU of the cas3 mutant (Fig. 5c). In contrast, only 20% of the larvae infected with the P. . TNF-␣ levels measured at 6 h showed significant differences in the two-way ANOVA. *, adjusted P ϭ 0.032; **, adjusted P ϭ 0.0228; ***, adjusted P ϭ 0.00061.
gingivalis wild-type strain died in the same period when larvae were injected with 3.85 ϫ 10 8 CFU (Fig. 5a), and larvae died after 40 h when injected with 3.85 ϫ 10 6 CFU of the wild-type strain (Fig. 5c). No larva mortality was observed in the control group until 3 days after inoculation (Fig. 5). Time-series dual-transcriptome results show that the transient response of G. mellonella larvae infected with the P. gingivalis wild-type strain is very different from the response seen with those infected with the ⌬cas3 mutant. To assess the "infection transcriptome" of wild-type P. gingivalis and the Δcas3 mutant in the G. mellonella infection model, we performed a time-series dual-transcriptome analysis, selecting infected G. mellonella larvae at 0, 1, 2, and 6 h. We decided to use G. mellonella as a model because a strong correlation has been observed in microbial pathogenicity in G. mellonella and mammalian systems (34,35) and because such an approach has been widely used as an initial screening model to assess virulence in a large number of different pathogens, including oral pathogens (36)(37)(38). The number of sequences in the G. mellonella infection experiments ranged between 74,281,712 and 86,843,647.
We identified 282 differentially expressed genes when we compared the infection transcriptomes of P. gingivalis wild-type strain and Δcas3 mutant infecting G. mellonella (Table S2). Dirichlet process Gaussian process mixture model (DPGP) analysis identified 6 clusters for the differentially expressed genes in the wild-type strain and 5 in the mutant. The cluster with the larger number of genes included 103 genes differentially expressed in the wild-type strain and mutant (Table S2). Both strains showed increased activity in the first hour, but the wild-type strain showed a steeper upregulation (Fig. 6). Between the first and second hours, both showed a trend of decreased gene expression, and finally, after the second hour, another upregulation which, this time, was more pronounced in the mutant (Fig. 6). Enrichment analysis of these genes identified activities associated with nucleic acid binding, nitrogen metabolism, and translation (Fig. 6). The remaining clusters did not share any significant fraction of genes. However, when we associated those genes with their correspondent GO terms, aminopeptidase activity was found to have decreased in the mutant but not in the wild type and, interestingly, copper-exporting ATPase activity was found to have increased in the wild type whereas there was a peak at 2 h and a steep drop after that in the mutant (Fig. S4).
In the transcriptome analysis of G. mellonella during infection, we found 4,490 differentially expressed genes along the course of the experiment. DPGP analysis  identified 10 clusters for the set of G. mellonella larvae infected with the wild type and 8 clusters for the group infected with the Δcas3 mutant. We performed Fisher's test enrichment analysis on the genes listed in the clusters to assess global activity changes. Only one cluster (cluster 1) in the wild-type strain and the mutant shared a considerable number of differentially expressed genes; 90% (3,477) of the differentially expressed genes were identical (Fig. S5). Cluster 1 was significantly depleted of genes that are associated with nucleic acid binding. In both cases, the GO term "immune system processes" was enriched as well. Genes corresponding to this GO term are associated with activities of the innate immune response, among them those encoding the following: inactivation of cytokines (membrane alanyl aminopeptidases-like, LOC113515582, LOC113515584, and LOC113517732), phagocytic processes in dendritic cells and macrophages (aminopeptidases N-like, LOC113515678, LOC113514207, and LOC113515696), IL-1 receptor signaling pathway (myeloid differentiation primary response protein MyD88, LOC113513443), Toll-like receptor 6 (LOC113512725), immune cell lineage commitment and maintenance (RNA-binding protein lark, transcript variant X2, LOC113513283, and LOC113516103), and antimicrobial peptide (AMP) production (virescein-like, LOC113509608). Interestingly, the differential expression of MyD88 and Toll-like receptor 6, the latter of which is known to heterodimerize with Toll-like receptor 2 (TLR2), may suggest activation of TLR2 in a MyD88-dependent manner. However, the expression profiles were very different; while the wild-type strain showed Normalized % sequences wild-type expected % wild-type observed % Δcas3 observed % Δcas3 expected %

FIG 6
Heat map of GSEA of biological processes in Gene Ontology (GO) of P. gingivalis differentially expressed genes in cluster 1. DPGP analysis identified 6 clusters for the differentially expressed genes in the wild-type strain and 5 in the mutant. The more significant of the clusters (cluster 1) shared 103 differentially expressed genes and are represented by the clusters shown in the figure. Heat map data are clustered vertically based on the normalized frequency of expected or observed sequences and are horizontally arranged based on biological processes. The analysis produces two sets of results, the expected frequencies of sequences (from the reference set) and the observed frequencies of sequences (from the test set). Significance was assessed by Fisher's exact test with an FDR value of Ͻ0.05. Test set, list of differentially expressed genes; reference set, list of all genes in the genome of P. gingivalis. a drop in expression of these genes in the first 2 h (Fig. S5a), the mutant showed an increase in expression during the first hour followed by a drop (Fig. S5b). Interestingly, two genes (LOC113517945 and LOC113515116) involved in the synthesis of phenoloxidase (PO), a key enzyme in melanin biosynthesis during melanization (39,40), were differentially expressed in both cluster 1 of the wild type and cluster 3 of the Δcas3 mutant (Fig. S5).
The remaining clusters did not share a significant number of genes. Two clusters, with very different profiles, one from the group infected with the wild-type strain (Fig. 7, upper right panel) and one from the group infected with the Δcas3 mutant ( Fig. 7, lower right panel), showed similar activities. These activities were associated with apoptosis, autophagy, response to infection, response to stress, hormone metabolism, melanogenesis, and cytoskeleton organization (Fig. 7). Genes in wild-type cluster 3 maintained the same expression level for the first 2 h and did not show increased expression after that time (Fig. 7, upper right panel). In contrast, genes in cluster 3 in the Δcas3 mutant showed increased expression levels immediately (Fig. 7, lower right panel). In both groups of worms, we observed activation of innate immune system genes in the Toll signaling pathway as well as activities associated with responses to stress, such as superoxide metabolic processes and response to hydrogen peroxide (  (Table S2).
In contrast, in the case of the group infected with the Δcas3 mutant, we observed a response meditated by cytokines and phagocytosis (Fig. 7).
Additionally, while the G. mellonella group infected with the wild-type strain showed enrichment in a large number of activities associated with hormone metabolism (serotonin, catecholamine, corticosterone, and dopamine), that was not the case for the group infected with the Δcas3 mutant (Fig. 7). On the other hand, apoptosis and autophagy were enriched in the group infected with the Δcas3 mutant but not in the group infected with the wild-type strain (Fig. 7).
The last activity significantly altered during infection was cytoskeleton organization, with upregulation of genes associated with cytoskeleton organization in the wild-type group during the first hours of infection ( Fig. S6, clusters 2 and 4). In contrast, the increase in activity in the Δcas3 mutant group occurred latterly, and it was associated with repression of cytoskeleton turnover, such as negative regulation of actin filament depolymerization (Fig. S7, clusters 2 and 3).
Integrating dual transcriptomes of host-pathogen expression profiles in response to infection by coexpression network analysis. Using the package 'mixOmics,' we performed coexpression analyses on the dual transcriptomes during infection of G. mellonella. 'mixOmics' allows us to distinguish between positive and negative correlations in the patterns of expression. As can be observed, heat maps corresponding to the correlations of the transcriptomes of the P. gingivalis wild-type strain (Fig. S8a) and the Δcas3 mutant ( Fig. S8b) with the host transcriptome reveal that the associations of gene expression of the two strains with the host are very different.
As a result of our analysis, we obtained two coexpression networks, one for the wild type and one for the mutant. We further separated those correlation networks into positive-coexpression networks and negative-coexpression networks. Using these networks, we assessed both commonalities and differences that could indicate the different ways in which the two strains exerted their pathogenic effects in survival experiments. The first result that we observed was that the coexpression networks of negative associations in the wild-type strain and the mutant were completely different, with no overlap in structure. The analysis of differential networks did not result in any consensus network.
Moreover, although the P. gingivalis genes of those networks were associated with similar biological processes (cell redox homeostasis, lipid A synthesis, acetyl coenzyme A [acetyl-CoA] synthesis, and pathogenesis), the associated genes from the host were enriched in different activities. While the negative associations in the wild-type strain were mainly with membrane invagination and autophagic cell death in the host (Fig. S8c), the negative associations in the Δcas3 mutant were mostly with positive regulation of apoptosis and innate immune response (Fig. S8d). Interestingly, although the GO terms associated with the P. gingivalis genes were similar, the genes were strain specific. Only 37% of the P. gingivalis genes in those coexpression networks were the same in the two strains (Fig. S8e).
The positive associations of coexpression were very different. Employing Diffany analysis, we found a consensus coexpression network in both the wild-type and mutant strains. This consensus network showed that expression profiles of genes associated with cell redox homeostasis and pathogenesis in P. gingivalis were positively correlated with the patterns of genes related to muscle tissue development, pigment synthesis, and cell death regulation in the host (Fig. 8a). Additionally, we identified differential networks, defined as networks with identical edges but with opposite significantly different correlations (41). The two differential networks are shown in Fig. 8. Panel b of Fig. 8 shows the differential network where the connections were weaker in the mutant than in the wild-type connections or were missing in the mutant connections. P. gingivalis genes in this network were associated with response to stress, cell redox homeostasis, acetyl-CoA synthesis, lipid A synthesis, proteolysis, and pathogenesis. These activities were positively correlated with genes in the host enriched in GO terms associated with membrane invagination, myeloid leukocyte activation, and regeneration (Fig. 8b). Conversely, for the differential network where the connections were weaker than or missing in the links in the wild type compared with the links in the cas3 mutant, P. gingivalis genes were similarly associated with response to stress, cell redox homeostasis, acetyl-CoA synthesis, and pathogenesis. However, the enrichment of host coexpressed genes was mainly related to the regulation of apoptosis, vesicle-mediated transport, and small-molecule metabolic processes (Fig. 8c). As in the case of the negative-coexpression networks, despite the similarities in the GO term activities of the P. gingivalis genes, they were somewhat different. Only 36% of the genes in those networks overlapped in the two strains (Fig. 8d).

DISCUSSION
CRISPR-Cas systems are now known to modulate a variety of microbial processes; however, the contributions of these systems to the virulence of periodontal pathogens are thus far unknown. Our prior work focusing on the transition from oral health to disease showed that within individual patients, the majority of genes belonging to the CRISPR-Cas loci in P. gingivalis in sites during the transition to periodontal disease were highly upregulated compared to the sites in the same patients that did not progress, where CRISPR-associated gene expression did not change (20). In this article, we report that an in-frame mutant with a deletion of the cas3 gene, encoding an essential nuclease in type I-C systems, increased the virulence of the mutated strain, demonstrating that at least one of the CRISPR-Cas systems of P. gingivalis (type I-C system) is involved in regulating the virulence of P. gingivalis during infection.
The descriptions of mechanistic use of CRISPR-Cas systems by bacteria have focused mainly on its innate defense nature in protecting microbes from exogenous nucleic acids and phages (2). However, Lum et al. showed that although transcription of oral CRISPR loci is relatively ubiquitous, highly expressed CRISPR spacers do not necessarily target the most abundant oral phages (42). Moreover, recently, new functions, other  than defense, have been attributed to CRISPR systems of the oral microbiome (43). Even though bacteriophages have been isolated from oral bacteria, a bacteriophage for P.
gingivalis has yet to be isolated despite the efforts in this area (44,45). Although the P. gingivalis CRISPR-Cas systems may defend against oral phages, that prior work and our clinical data do not support the notion that oral phages alone are responsible for the changed expression occurring mainly in the context of the transition from periodontal health to disease. CRISPR-Cas genes were highly upregulated at sites with clinical evidence of transition to disease compared to healthy sites in the same individual; it would be expected that under both conditions at both sites, oral phages would be present. An expansion of our understanding of the importance of CRISPR-Cas systems to bacteria beyond protection from nucleic acid insult has come from emerging literature on their role in intracellular persistence of pathogens such as Campylobacter jejuni (6), Francisella novicida (7), Legionella pneumophila (46), Neisseria meningitidis (47), and Listeria monocytogenes (9). P. gingivalis, as well as several other oral microbiome members, is a regular component of the intracellular communities of oral epithelial cells (48,49). Intracellular survival is an essential part of the P. gingivalis lifestyle in the periodontal pocket as well a key element in transmission from cell to cell in the epithelium, enabling avoidance of the humoral immune response (49)(50)(51). CRISP-Cas systems are ideal elements for regulating intracellular survival, given how quickly they can control gene expression in the short time available for cell invasion. Therefore, we hypothesized that CRISPR-Cas systems in P. gingivalis might be involved in intracellular survival. Interestingly, mining a data set from an independent study which examined the genome of P. gingivalis to identify genes essential for colonization of epithelial cells, we noted that all genes in the type I-C system locus, except for cas2 (fig|431947.7.peg.1924), were significantly associated with successful infection in that prior investigation (52). Here, focusing on cas3, our findings confirm this prior observation as we found that cas3 transcription was induced in P. gingivalis growing intracellularly but not planktonically.
Nonetheless, cas3 mutation did not alter the planktonic or intracellular growth of P. gingivalis compared to the wild type. A recent study on the effects of deleting cas3 in Salmonella enterica showed no differences in liquid growth. In contrast to our results, those researchers found a decrease in intracellular survival in the mutant (53). The expression profiles of wild-type P. gingivalis and the cas3 mutant were very similar when the bacteria were placed in human serum-based medium, but when they colonized the human macrophage-like THP-1 cell line, the transcriptome profiles of the wild-type and mutant strains were very different. The only activities that were altered under both sets of growing conditions were activities associated with adhesion, probably because fimbrial genes, including fimA, were downregulated in the mutant. In S. enterica, deletion of cas3 leads to downregulation of safA, safB, safC, and safD genes in the fimbrial operon (53). The mechanisms by which cas3 controls levels of transcription of fimbrial genes remain unknown. However, given the importance of fimbria in attachment of bacteria to host cells and the similarity of our findings of reduced transcription of fimbrial genes due to cas3 deletion to those observed with S. enterica, this mechanism of CRISPR control of infection should be investigated further.
Deletion of cas3 had a significant effect on the pathogenesis of intracellular P. gingivalis, as well as on the immune response of THP-1 cells. We observed an upregulation of genes associated with responses to oxidative stress and iron acquisition in the mutant. These activities are at the core of the mechanisms of pathogenicity in P. gingivalis (23-26, 30, 54, 55), and an increase in their levels increase virulences. Unexpectedly, we did not observe changed secretion of TNF-␣ and IL-1␤ or other proinflammatory cytokines or chemokines between the cells infected with the wild type and those infected with the mutant. However, a modest trend of increased inflammatory molecule production coincides in time with the peak in the number of intracellular transcripts for cas3. These results agreed with the transcriptome analysis, where IL-1␤ and TNF-␣ genes were upregulated in the cells infected with the mutant. In a C57BL/6J mouse model, Li et al. found that secretion of TNF-␣ and IL-1␤ after 2 h of infection was higher when mice were infected with a CRISPR-Cas P. aeruginosa mutant than when they were infected with the wild-type strain (8).
Furthermore, infection with the Δcas3 mutant strain altered neutrophil chemotaxis and gene silencing by microRNA (miRNA). We observed a high number of microRNAs upregulated in the cells infected with the Δcas3 mutant. Previous studies had also shown that during periodontitis, the expression of microRNAs in humans is altered (32,33), a fraction of which we found differentially expressed in our study.
Although the number of studies is still limited, the dampening of the immune response seems to be an element common to the mechanisms by which CRISPR-Cas systems control virulence. F. novicida evades the innate immune response by regulating levels of the endogenous transcript of bacterial lipoprotein (BLP; FTN_1103) via a CRISPR/Cas-associated RNA (small Cajal body-specific RNA [scaRNA]) (7), and P. aeruginosa evades host defenses by targeting mRNA of the quorum-sensing regulator LasR to dampen the recognition by TLR-4 (8). Our results agreed with those two studies (7,8) in that mutation of cas3 increased virulence. The exception to this is S. enterica, where deletion of cas3 weakened bacterial virulence (53). Interestingly, both P. aeruginosa and S. enterica possess type I CRISPR-Cas systems but target quorum-sensing systems, and yet different outcomes are seen (8,53). P. gingivalis employs an AI-2/LuxS-mediated quorum-sensing system. We did not see any differences in the levels of expression of luxS in our experiments performed with THP-1 cells. However, in our time-series analysis of infection in G. mellonella, luxS was differentially expressed, with a peak in the first hour in the wild type and a peak at 6 h in the mutant. Note that P. gingivalis does not possess bacterial lipoprotein (BLP) or LasR homologs (data not shown). Thus, it would appear that neither BPL nor LasR represents common mechanisms of action for P. gingivalis Cas3, thus supporting the notion that another, unknown system is likely involved for at least this organism. Interestingly, 11 of the 1,188 P. gingivalis spacers identified by Watanabe et al. (56) had perfect matches with genes in P. gingivalis ATCC 33277. Most of those genes are hypothetical proteins, but we found that one of the targets upregulated in our analysis was a putative TonB gene, which may suggest a potential role in regulating iron metabolism by the CRISPR-Cas system in P. gingivalis.
The most direct evidence of the increase of virulence in the mutant strain came from our experiments using G. mellonella as an infection model. Previous work with several microbial pathogens demonstrated a positive correlation between this model's results and results from other mammalian disease models (35,57,58). Innate immune responses of G. mellonella are comparable with vertebrate innate immune responses and involve recognition of the bacteria and production of antimicrobial molecules (59). Moreover, there is good evidence that the G. mellonella model is suitable for studying pathogenesis and immune responses in human oral pathogens (60), including P. gingivalis (38,61). Our results showed a significant increase in virulence measured as the death rate of the infected larvae with the cas3 mutant during the 48 h of the experiment.
All previous studies on the role of CRISPR-Cas systems in virulence were performed at a single time point for comparisons of conditions. Time-series analysis of the G. mellonella model of infection gave us a dynamic picture of the process rather than just a snapshot of the differences in gene expression at a given time. The most striking results of those experiments are those demonstrating the great differences in the host's infection transcriptomes. Activities associated with apoptosis, autophagy, response to infection, response to stress, hormone metabolism, melanogenesis, and cytoskeleton organization showed distinct time-dependent behavior. Larvae infected with the wild type maintain the same level for the first 2 h, and it is not until after that time that they increased their expression. In contrast, the Δcas3 mutant increased the expression level of those activities immediately. In human coronary artery endothelial cells (HCAEC) and human aortic endothelial cells (HAEC), P. gingivalis activates cellular autophagy to provide a replicative niche while suppressing apoptosis (62,63). The increase in autophagic activity in the mutant-infected larvae seems to indicate a faster infectious process. Moreover, P. gingivalis additionally induces endoplasmic reticulum (ER) stressmediated apoptosis in human umbilical vein endothelial cells (HUVEC), followed by an autophagic response (64).
The immune response of G. mellonella consists of two tightly interconnected components, i.e., cellular and humoral components, that collaborate to ensure the best protection to the insects (65,66). The cellular response is mediated by hemocytes and involves responses such as phagocytosis, encapsulation, and nodulation (67). The humoral response consists of the synthesis of defense molecules, including antimicrobial peptides and melanin, which are involved in encapsulating the invading pathogen and also stimulate the defense activity of other antimicrobial molecules (67).
We observed two distinct response profiles of the larvae depending on the infecting strain. The first pattern of response was seen with the G. mellonella group infected with the wild type and seemed driven primarily by a humoral response corresponding to the synthesis of AMPs. In periodontitis, neutrophils are recruited to sites of injury and secrete a rich mixture of AMPs. Over 45 AMPs have been identified in saliva and gingival crevicular fluid (GCF). Still, their biological function with respect to disease may be more complex than just the killing of bacteria since the levels of several AMPs are below the MIC for oral pathogens (68). Although they are effective against planktonic bacteria, effectiveness against biofilms is not apparent (69).
As a part of the immune response, G. mellonella synthesizes a wide variety of AMPs (67,(70)(71)(72). Additionally, the G. mellonella group infected with the wild-type strain showed enrichment in a large number of activities associated with hormone synthesis and response. This increase in hormone production is most likely linked to melanogenesis. Melanization is one of the primary defense reactions to the encapsulation of pathogens in G. mellonella (66). Bacterial infections induce bursts of dopamine that can be related to the direct participation of dopamine in the phenoloxidase cascade, a major defense system in many invertebrates, ultimately leading to melanization of pathogens and damaged tissues (73).
In contrast, defense against infection with the Δcas3 mutant was primarily mediated by cellular responses, including synthesis of cytokines and phagocytosis, with enrichment in apoptosis and autophagy activities. The presence of IL-1␣ and TNF-␣-like molecules in the hemocytes of G. mellonella (74) and transcriptome analysis of this organism identified homologs of a cytokine called Spätzle, which functions in the antimicrobial immune response in larvae and adults (75). In Drosophila melanogaster, Spätzle functions as a circulating cytokine and specifically binds with high affinity to Toll, thereby activating the Toll signaling pathway (76). G. mellonella possesses hemocytes that can phagocyte intruders or capture them in multicellular structures called nodules of capsules (67). Granulocytes, representing a type of hemocyte in G. mellonella, showed high levels of expression of the autophagy response through microtubule-associated LC3 when G. mellonella was infected with a high-virulence strain of Actinobacillus pleuropneumoniae, whereas treatment with a low-virulence strain induced lower levels of expression of this protein (36). Moreover, in Helicobacter pylori, the killing of larvae always correlated with the induction of apoptosis in hemocytes (37).
We also observe changes in cytoskeleton metabolism as part of the cellular defense against infection. Infection of P. aeruginosa affects G. mellonella hemocyte cytoskeleton rearrangement, disabling cellular immunity (77). Hemocyte motility relies on the continual restructuring of the cytoskeleton. In the wild-type-infected larvae, we observed upregulation of genes associated with cytoskeleton reorganization in the first hours of infection. Simultaneously, in the Δcas3 mutant group, the increase in activity occurred later, and it was associated with the cytoskeletal organization's inhibitory activities.
Although the two models that we used in our transcriptome infection experiments (human THP-1 cells and G. mellonella) and the experiment designs (cross-sectional data and time-series data) were very different, the transcriptome results suggest that the differential immune responses elicited by P. gingivalis wild-type and cas3 mutant strains in G. mellonella and THP-1 cells are somehow controlled by the CRISPR-Cas system, controlling nucleic acid binding (mRNA binding), migration of immune cells, and cytokine and chemokine responses, which are activities altered in the transcriptomes of both THP-1 cells and G. mellonella.
Most studies of infectious disease pathogenesis focus on either the host or the pathogen, even though their interaction determines the outcome. Moreover, such studies generally analyze transcriptomes at only a single stage of the infection process, missing the changes across different stages of the disease, which is a dynamic process. The presence of coexpression networks does not necessarily imply causality, but they are increasingly employed in bioinformatics applications to explore the system-level functionality of genes in host-pathogen interactions (78,79). Analyses of our study's relevant coexpression networks showed that the levels of expression of P. gingivalis genes associated with G. mellonella were similar with respect to the GO categories in wild-type and Δcas3 mutant infections. However, the specific genes related to those GO categories did not overlap. We identified a consensus network response in G. mellonella that probably acts as a core response against both strains. Not surprisingly, it includes activities essential in the immune response of larvae against pathogens: pigment synthesis and regulation of cell death in the host (36,66,67,80). To identify the specific differences in coexpression between the two strains, we generated differential coexpression networks, defined as networks with shared edges but with opposite significantly different correlations (41). The connections that were weaker or missing in the mutant compared with the connections in the wild-type activities were positively correlated with genes in the host enriched in GO terms associated with membrane invagination, myeloid leukocyte activation, and regeneration. G. mellonella does not possess leukocytes, and the GO terms are probably related to hemocyte activity. Conversely, the differential network where the connections were weaker or missing in the wild type compared with the connections in the mutant, enrichment of host coexpressed genes, was mainly associated with regulation of apoptosis, vesiclemediated transport, and small-molecule metabolic processes.
In conclusion, we have demonstrated a direct link of the presence of cas3 gene with an increase in the virulence of P. gingivalis, as has been described previously for other pathogens (7,8). In addition to revealing the importance of a CRISPR-Cas system in controlling the virulence of P. gingivalis, we found that the lack of cas3 had significant impacts leading to the rearrangement of the host's response, mainly when evaluated in our time-series infection modeling. These findings extend and further support the notion that CRISPR-Cas systems can effectively reshape bacterial virulence and escape the host immune defense.

MATERIALS AND METHODS
Bacterial growth conditions. Porphyromonas gingivalis ATCC 33277 was cultured anaerobically at 37°C. Cells were maintained on brain heart infusion (BHI)-blood agar plates supplemented with 5 g/ml hemin and 1 g/ml menadione (vitamin K). Two different liquid growth media were used in the experiments. For the construction of the mutant and the infection experiments, we used BHI broth supplemented with 1 mg/ml yeast extract, 5 g/ml hemin, and 1 g/ml menadione. For the transcriptome experiments, we grew P. gingivalis in 20% heat-inactivated human serum (Sigma-Aldrich H3667-20ML) supplemented with 5 g/ml hemin and 1 g/ml menadione as described previously by Grenier et al. (81).
Construction of a CRISPR-Cas3 gene knockout strain of Porphyromonas gingivalis. We constructed the fig|431947.7.peg.1929 (PATRIC annotation [82]) (see Fig. S1a in the supplemental material) gene knockout strain by replacing the whole gene with an erythromycin resistance cassette. First, we constructed a plasmid (pUC19) carrying the erythromycin resistance cassette (ermF gene from pVA2198) flanked by the 1-kb region upstream and downstream of fig|431947.7.peg.1929. The fig|431947.7.peg.1929 gene was replaced by the erythromycin cassette and was kept in frame with the rest of the genes in the operon. The Cas3 Controls Virulence in P. gingivalis September/October 2020 Volume 5 Issue 5 e00852-20 msystems.asm.org 17 construct was made in such a way that the complete fig|431947.7.peg.1929 gene was replaced by the erythromycin cassette and was kept in frame with the rest of the genes in the operon. This construct was made using a NEBuilder HiFi DNA assembly kit. The kit was used with 100 ng of P. gingivalis genomic DNA according to the manufacturer's protocol. The NEBuilder Assembly Tool (NEB) was used to design the primers for the NEBuilder HiFi DNA kit protocol. Dh5-alpha chemical competent cells were used for the transformation. The resulting plasmid-transformed Dh5-alpha cells were selected on ampicillin LB plates. The plasmid was extracted using EZNA plasmid DNA minikit II (Omega) and sequenced. After verification by sequencing, the plasmid was named pFLUF001. Primers JS_Cas3KOPCR1-F (5=-GGAAGTGACCGTTATCGAAGAT-3=) and JS_Cas3KOPCR1-R (5=-GCCTTA CGAATAGGCCATAAGA-3=) were used to amplify, from pFLUF001, the 2.7-kb DNA fragment containing the erythromycin cassette and its flanking regions. Pfu polymerase (Fermentas) was used according to the manufacturer's protocol. The amplified fragment was cleaned using an EZNA gel extraction kit (Omega) and used for electroporation of P. gingivalis electrocompetent cells. P. gingivalis ATCC 33277 electrocompetent cells were made by growing the cells on tryptic soy broth (TSB) supplemented with hemin and vitamin K to an optical density at 600 nm (OD 600 ) of 0.6 to 0.7. After centrifugation, the cells were washed twice in ice-cold electroporation buffer (10% glycerol, 1 mM MgCl 2 ) and finally resuspended in a minimal amount of buffer. A 100-l sample of P. gingivalis electrocompetent cells was electroporated with different amounts of the purified DNA fragment. Cells were plated on BHI-blood agar plates supplemented with hemin and vitamin K and 10 g/ml of erythromycin. After anaerobic incubation at 37°C for 7 days, the resulting colonies were streaked on new plates, and single colonies were obtained. The fig|431947.7.peg.1929 (Cas3) gene knockout strain was verified by colony PCR using primers from the erythromycin cassette, and the adjacent region to the flanking region was amplified and cloned in pFLUF001. The amplified product was confirmed by sequencing. The correct gene knockout strain was grown on liquid media, and glycerol and dimethyl sulfoxide (DMSO) stocks were prepared and stored in an Ϫ80°C freezer.
Bacterial infection experiments. P. gingivalis wild-type and Δcas3 mutant cells were harvested from the BHI broth culture by centrifugation, and the bacteria were washed three times in RPMI 1640 medium and adjusted to an OD 660 of 1.0 (approximately 1 ϫ 10 9 CFU/ml) and added to PMA-activated THP-1 cells at a multiplicity of infection of 100, in the antibiotic-free medium, and after 2 h of infection, the THP-1 cells were washed to remove nonadherent bacteria. After 2 h of incubation, the THP-1 cells were treated for 1 h with metronidazole (200 g/ml)/gentamicin (300 g/ml) to kill extracellular bacteria (51), were washed to remove antibiotics, and then were lysed with sterile water, followed by addition of an equal volume of 2ϫ phosphate-buffered saline (PBS). Serial 10-fold dilutions of each lysate were plated on blood agar plates. CFU/ml of P. gingivalis was calculated. In separate experiments, THP-1 cells were cultured with the P. gingivalis ATCC 33277 wild-type strain or the Δcas3 mutant, and both cell culture supernatant fluids and RNA were harvested for measurement of cytokine expression and RNA was extracted for dual transcription. Differences in bacterial counts were assessed by performing a Kruskal-Wallis test corrected for multiple comparisons using the function 'kruskal' from the 'agricolae' package in R, with a false-discovery rate (FDR) value of Ͻ0.05.
RNA extraction. Total RNA was extracted from those samples using a mirVana RNA isolation kit (Life Technologies). Samples were bead-beaten for 1 min at maximum speed with 300 l of 0.1-mm diethyl pyrocarbonate (DEPC)-treated zirconia-silica beads (BioSpec Products) in the mirVana lysis buffer. Bacterial rRNA was depleted using Ribo-Zero gold rRNA removal kits (Illumina) following the manufacturer's protocol. In the THP-1 infection experiments, sequencing of total RNA, including human and bacterial mRNA, was performed. Their identification was performed bioinformatically, as described below.
G. mellonella RNA was extracted using a mirVana RNA extraction kit (Thermo Fisher Scientific) with minor modifications. Individual larvae were washed with autoclaved PBS and were cut at the distal end. The internal content was snap-frozen in liquid nitrogen, homogenized and pulverized, and immediately incubated at room temperature for 1 h in lysis buffer (mirVana kit; Thermo Fisher Scientific). After the lysis steps, the manufacturer's instructions were followed.
In the experiments performed for analysis of G. mellonella infection, we divided the total RNA into two subsamples. One half was used for transcriptome analysis of G. mellonella, for which we used a Dynabeads mRNA purification kit to isolate eukaryotic mRNA for transcriptome analysis. The other half was depleted of eukaryotic mRNA using a MICROBEnrich kit (Thermo Fisher Scientific, catalog no. AM1901).
RT-qPCR quantification of cas3 transcripts. From the initial RNA extraction, possible traces of DNA were removed using Ambion's Turbo DNA-free kit (Ambion) following the manufacturer's instructions. The volume of Turbo DNase I (Ambion's Turbo DNA-free) was increased to 3 l, and the reaction mixture was incubated at 37°C for 60 min. RNA (100 ng) was reverse transcribed with random hexamer primers and with SuperScript II reverse transcriptase (Invitrogen) following the manufacturer's instructions. Reverse transcription was performed at 42°C for 2 h, after an initial incubation step of 10 min at 25°C. The synthesized cDNA was used in a DyNAmo Flash SYBR green qPCR kit (New England Biolabs, Ipswich, MA) according to manufacturer's instructions, using the specific primers for the genes of interest. To compare the relative expression levels of genes, we modified the threshold cycle (2 ϪΔΔCT ) method (84), and we used the formula cDNA mutant /cDNA control ϭ (1 ϩ E cDNAcontrol Only two ) CtcDNAcontrol /(1 ϩ E cDNAmutant ) CtcDNAmutant to take into consideration the different amplification efficiencies in the separate qPCR runs.
RNA-seq library construction. Sequencing was performed at the Interdisciplinary Center for Biotechnology Research (ICBR) at the University of Florida using a HiSeq 2500 machine. First, rRNAs were removed from total RNA by the use of an Illumina Ribo-Zero gold rRNA removal kit following the manufacturer's protocol and eluted into 10 l EB buffer. The transcriptome sequencing (RNA-seq) library was then processed using a NEBNext ultradirectional RNA library prep kit for Illumina (NEB, USA) following the manufacturer's recommendations. A 5-l volume of depleted RNA mix was used together with 5 l of first-strand synthesis reaction mix (NEBNext first-strand synthesis reaction buffer [5ϫ] and NEBNext random primers), fragmented by heating at 94°C for the desired time. This step is followed by first-strand cDNA synthesis performed using reverse transcriptase and oligo(dT) primers. Synthesis of double-stranded (ds) cDNA is performed using the second-strand master mix provided in the kit, followed by end repair and adaptor ligation. Finally, the library is enriched (each library has a unique barcode, and each primer has a common adaptor sequence which was added in the previous adaptor ligation step and a unique overhang index unique to each sample) by a certain number of cycles of amplification and purified by the use of Agencourt AMPure beads (Beckman Coulter, catalog no. A63881). Finally, individual libraries were pooled with equimolar volumes and sequenced by the use of an Illumina HiSeq 3000 system (Illumina Inc., CA, US) and a run of 2 ϫ 100 cycles.
Illumina instrument run. In preparation for sequencing, barcoded libraries were sized on an Agilent 2200 TapeStation system. Quantitation was done by the use of both QUBIT and qPCR (Kapa Biosystems, catalog no. KK4824). Individual samples were pooled in equimolar volumes at 2.5 nM. This "working pool" was used as the input in the HiSeq 3000 instrument sample preparation protocol (Illumina material no. 20015630, document no. 15066496 v04, January 2017). Typically, a 250 pM library concentration was used for clustering on a cBOT amplification system. This resulted in an optimum clustering density at which the proportions of clusters passing the filters ranged between 65% and 75%. Six RNA-seq barcoded libraries were pooled for sequencing in multiplex on a single flow cell lane, using a configuration of 2 ϫ 100 (paired-end) cycles. Such an sequencing configuration was achieved by pooling the reagents from 150-cycle and 50-cycle Illumina HiSeq 3000 SBS kits. A typical sequencing run in the HiSeq 3000 instrument produced Ͼ300 million reads from each end, per lane, with a Q30 value of Ն90%. For RNA-seq, the use of 50 million reads per end per sample provided sufficient depth for transcriptome analysis. The sequencing run was performed at the NextGen research facility of the Interdisciplinary Center for Biotechnology Research (University of Florida).
Assessment of cytokine and chemokine production. Frozen cell culture supernatant fluids were thawed, and the levels of TNF-␣, IL-1␤, IL-6, IL-8, IL-10, and RANTES were determined by the use of Milliplex multiplex assays (EMD, Millipore). Data were acquired on a Luminex 200 system running xPONENT 3.1 software (Luminex) and analyzed using a 5-parameter logistic spline-curve fitting method and Milliplex Analyst V5.1 software (Vigene Tech). Statistical differences were assessed by two-way analysis of variance (ANOVA) using the 'emmeans' package in R, applying an FDR value of Ͻ0.05 for multiple-comparison corrections. Experiments were performed in triplicate. Galleria mellonella infection model. For all of the G. mellonella experiments, insects in the final instar larval stage were purchased from Vanderhorst, Inc. (St. Marys, OH). Upon arrival, any dead larvae present were separated from healthy larvae, which were then weighed and placed randomly into groups and kept at room temperature until the next day, when infection was performed. Seven groups of 15 larvae, ranging in weight from 200 to 300 mg, and with no signs of melanization, were randomly chosen and used for subsequent infection. A 25-l Hamilton syringe was used to inject 5-l aliquots of bacterial inoculum into each larva's hemocoel via the last left proleg. Three groups received wild-type P. gingivalis (3.85 ϫ 10 8 , 7.7 ϫ 10 7 , and 3.85 ϫ 10 6 CFU per larvae), and three groups received the cas3 mutant (2.15 ϫ 10 8 , 4.3 ϫ 10 7 , and 2.15 ϫ 10 6 CFU per larvae). Three control groups included THSB medium alone, tryptic soy broth (TSB) (BD, Becton, Dickinson and Co.) plus P. gingivalis wild-type heat-killed (10 min at 75°C), and THSB plus Δcas3 mutant heat-killed (10 min at 75°C). After injection, larvae were incubated at 37°C, and the appearance of melanization and survival were recorded at 0.5, 1, 1.5, 2, 3, 3.75, 4.25, 6.25, 22, 24, 28, and 42 h. After injection, larvae were incubated in the dark at 37°C, and appearance (signs of melanization) and survival were recorded at selected intervals. Larvae were scored as dead when they displayed no movement in response to touch. Kaplan-Meier killing curves were plotted, and estimations of differences in survival were compared using a log-rank test. A P value of Յ0.05 was considered significant. All data were analyzed with the 'survival' and 'survminer' packages in R. Experiments were repeated three independent times with similar results. We followed the protocol described above for the infection transcriptome experiments, but the initial concentrations of P. gingivalis injected were 7.0 ϫ 10 8 CFU/ml for the wild-type strain and 1.0 ϫ 10 8 CFU/ml for the mutant strain.
In the case of the THP-1 infection experiments, differential expression analysis was performed using NOISeqBio (89). After exploratory analysis with NOISeqBio, we selected reads per kilobase per million (RPKM) normalization for the THP-1 cell data, tmm normalization with length correction for the P. gingivalis intracellular transcriptome, and tmm normalization without length correction for the P. gingivalis planktonic transcriptome. Posterior analyses were used only with significant differential expression and log changes Ͼ2.
The genome of G. mellonella has been sequenced and contains 14,067 protein-coding genes (90). Time-series analysis of the G. mellonella infection transcriptomes was performed in two steps. First, we identified differentially expressed genes using a DESeq2 spline approach, as recommended by Spies et al. for short-time series data (Ͻ8 time points) (91). In the second step, we clustered those genes based on their trajectories during the experimental period. For the following step, we used the Dirichlet process Gaussian process mixture model (DPGP) and DP_GP cluster software with 500 iterations (92).
We used Gene Ontology (GO) terms for gene set enrichment analysis (GSEA). In the case of P. gingivalis and G. mellonella, we generated our own GO annotation using the pipeline in OmicsBox (93). The GO annotation was extracted for the human genome directly from its gff3 file. GSEA was performed using GOrilla (94) in the case of the human sequences and OmicsBox for any other enrichment analysis, including enrichment using Fisher's exact test. In all cases, we consider an FDR value of Ͻ0.05. Summaries of GO terms and GO treemaps were obtained using REVIGO (95). GSEA results corresponding to enriched terms from OmixBox, quantified as percentages of sequences in the reference and test sets, were represented as heat maps using the 'pheatmap' package in R, after normalization with the function 'decostand' from the 'vegan' package.
Coexpression networks of host and P. gingivalis genes. The integration of host-microbe expression profiles was performed using the R package 'mixOmics' (96). We calculated the sparse partial least-square (sPLS) correlations between the differentially expressed genes from eukaryotic cells and the different mutants of P. gingivalis. Relevance networks showing correlations between genes from eukaryotic cells and microorganisms were visualized in Cytoscape (97) with a threshold correlation of 0.90.
We used the Cytoscape plugin Diffany (41) to obtain and analyze consensus and differential coexpression networks. On those networks, we visualized enriched GO terms on the G. mellonella genes at the biological process level using the Cytoscape plugin GOlorize (98), which uses the Cytoscape BiNGO (99) plugin to perform GSEA on the nodes of the network. For GSEA performed on BiNGO, we generated the BiNGO annotation based on the GO annotation described above. The tests were performed using the default settings with an FDR significance value of Ͻ0.05. The selection of the enriched terms to be visualized was performed by summarizing the results with REVIGO and selecting for the representative GO terms. In the case of P. gingivalis genes present in the networks, we associated them with GO terms, removed the singletons, and summarized the results using REVIGO.
Data availability. The data sets used in these analyses were deposited at the Gene Expression Omnibus (GEO) data repository of NCBI. Sequences have been deposited at the GEO database (https:// www.ncbi.nlm.nih.gov/geo/) with submission number GSE154569.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.

ACKNOWLEDGMENTS
This work was supported by the National Institute of Dental and Craniofacial Research of the National Institutes of Health (NIDCR/NIH) under award number DE021553.
We declare that we have no competing interests. J.F.-L. conceived the study; J.S., A.D.-P., and F.G.R. performed the experiments; J.F.-L. and F.C.G. III analyzed the data. All of us wrote and revised the manuscript. All of us read and approved the final manuscript.