
==== Front
PLoS Pathog
PLoS Pathog
plos
PLOS Pathogens
1553-7366
1553-7374
Public Library of Science San Francisco, CA USA

39178326
10.1371/journal.ppat.1012454
PPATHOGENS-D-24-01051
Research Article
HPV induced R-loop formation represses innate immune gene expression while activating DNA damage repair pathways
HPV induced R-loops alter cellular gene expression
Templeton Conor W. Conceptualization Data curation Formal analysis Funding acquisition Investigation Methodology Resources Software Supervision Validation Visualization Writing – original draft Writing – review & editing
https://orcid.org/0000-0002-6314-623X
Laimins Laimonis A. Conceptualization Formal analysis Funding acquisition Investigation Methodology Project administration Resources Supervision Validation Visualization Writing – original draft Writing – review & editing *
Department of Microbiology-Immunology, Northwestern University Feinberg School of Medicine, Chicago, Illinois, United States of America
Lambert Paul F. Editor
University of Wisconsin Madison School of Medicine and Public Health, UNITED STATES OF AMERICA
The authors have declared that no competing interests exist.

* E-mail: l-laimins@northwestern.edu
23 8 2024
8 2024
20 8 e101245417 5 2024
28 7 2024
© 2024 Templeton, Laimins
2024
Templeton, Laimins
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

R-loops are trimeric nucleic acid structures that form when an RNA molecule hybridizes with its complementary DNA strand, displacing the opposite strand. These structures regulate transcription as well as replication, but aberrant R-loops can form, leading to DNA breaks and genomic instability if unresolved. R-loop levels are elevated in many cancers as well as cells that maintain high-risk human papillomaviruses. We investigated how the distribution as well as function of R-loops changed between normal keratinocytes and HPV positive cells derived from a precancerous lesion of the cervix (CIN I). The levels of R-loops associated with cellular genes were found to be up to 10-fold higher in HPV positive cells than in normal keratinocytes while increases at ALU1 elements increased by up to 500-fold. The presence of enhanced R-loops resulted in altered levels of gene transcription, with equal numbers increased as decreased. While no uniform global effects on transcription due to the enhanced levels of R-loops were detected, genes in several pathways were coordinately increased or decreased in expression only in the HPV positive cells. This included the downregulation of genes in the innate immune pathway, such as DDX58, IL-6, STAT1, IFN-β, and NLRP3. All differentially expressed innate immune genes dependent on R-loops were also associated with H3K36me3 modified histones. Genes that were upregulated by the presence of R-loops in HPV positive cells included those in the DNA damage repair such as ATM, ATRX, and members of the Fanconi Anemia pathway. These genes exhibited a linkage between R-loops and H3K36me3 as well as γH2AX histone marks only in HPV positive cells. These studies identify a potential link in HPV positive cells between DNA damage repair as well as innate immune regulatory pathways with R-loops and γH2AX/H3K36me3 histone marks that may contribute to regulating important functions for HPV pathogenesis.

Author summary

R-loops are trimeric RNA: DNA hybrids that regulate transcription but can also lead to DNA breaks. R-loop levels are elevated in many cancers as well in precancerous cells containing high-risk HPVs. How the distribution and function of R-loops changes due to the presence of viral genomes was examined by comparing cells derived from an HPV 31 positive precancer (CIN 612) to normal keratinocytes. Up to 10-fold higher R-loop levels were detected on cellular genes in CIN 612 cells, with over 500-fold increases at ALU elements. Analysis of the changes in R-loop landscape between normal and HPV positive cells identified a functional association between R-loop formation, gene expression, and chromatin states. R-loop formation in HPV postive cells was equally associated with negatively and positively regulated gene expression; however, coordinated expression of genes in specific pathways by R-loops was observed. This included repression of genes in the innate immune surveillance pathway, while those regulating DNA damage repair and metabolism were activated.

http://dx.doi.org/10.13039/100007316 Division of Cancer Prevention, National Cancer Institute RO1CA142861 https://orcid.org/0000-0002-6314-623X
Laimins Laimonis A. http://dx.doi.org/10.13039/100007316 Division of Cancer Prevention, National Cancer Institute RO1 CA59655 https://orcid.org/0000-0002-6314-623X
Laimins Laimonis A. http://dx.doi.org/10.13039/100015691 Division of Microbiology and Infectious Diseases, National Institute of Allergy and Infectious Diseases R21AI180285 https://orcid.org/0000-0002-6314-623X
Laimins Laimonis A. http://dx.doi.org/10.13039/100000069 National Institute of Arthritis and Musculoskeletal and Skin Diseases 5T32AR060710-10 Templeton Conor W. This work was supported by grants to L.A.L. from the National Cancer Institute RO1CA142861 and RO1 CA59655 along with R21AI180285 from National Institute of Allergy and Infectious Diseases. L.A.L. received salary support from all three grants. C.W.T. salary was supported by the RO1CA142861, R21AI180285, and 5T32AR060710-10 from the National Institute of Arthritis and Musculoskeletal and Skin Diseases. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. PLOS Publication Stagevor-update-to-uncorrected-proof
Publication Update2024-09-05
Data AvailabilityAll relevant data are in the manuscript and its Supporting information files. All raw datasets generated in this report have been submitted to the NCBI GEO database and are public. These datasets can be found at GSE240391, GSE274109, GSE274120, GSE274122, GSE274119, and GSE274005.
Data Availability

All relevant data are in the manuscript and its Supporting information files. All raw datasets generated in this report have been submitted to the NCBI GEO database and are public. These datasets can be found at GSE240391, GSE274109, GSE274120, GSE274122, GSE274119, and GSE274005.
==== Body
pmcIntroduction

R-loops are trimeric nucleic acid structures that are formed when an RNA strand hybridizes with its complementary DNA and displaces the opposite strand [1–5]. These structures are long-lived and regulate normal transcription as well as replication. Aberrant R-loops can form, and failure to efficiently resolve these structures leads to transcription/replication conflicts, resulting in DNA break formation [6–11]. High levels of R-loops have been detected in cell lines derived from precancerous lesions that maintain high risk human papillomaviruses (HPVs) [12–16]. Furthermore, human cancers themselves contain high levels of R-loops, which suggests they contribute to progression [17–20]. Few studies have, however, examined how R-loop distributions and functions differ between cells, such as those that maintain human papillomaviruses (HPVs) and normal cells. Our studies investigated how the landscape of R-loop distributions and functions change between cells that maintain high risk HPV-31 and normal keratinocytes.

HPVs are the etiological agents of cervical cancer and are responsible for ~5% of all human cancers. Cancers and precancers induced by infection with high-risk HPVs provide an excellent model for studying factors influencing progression [21–23]. Cervical lesions caused by high-risk HPVs are characterized as cervical intraepithelial neoplasia grades I to III (CINI-CINIII), and these precede the development of frank cervical cancer [24–26]. Characterization of lesions as CIN I to CINII is made according to the degrees to which epithelia are altered. In precancerous CIN I lesions, HPV genomes are maintained as extrachromosomal elements or episomes that replicate coordinately with cellular replication, while productive viral replication or amplification is restricted to differentiated suprabasal cells [27–29]. CIN 612 is an immortal cell line that was derived from a CINI cervical biopsy and stably maintains high-risk HPV 31 genomes as episomes [30]. Transfection of normal human keratinocytes with cloned HPV sequences leads to their immortalization and stable maintenance of viral episomes [31]. These cell lines are similar to those derived from CIN I lesions and demonstrate that viral genomes are responsible for changes indicative of precancerous lesions. Previous studies demonstrated the presence of high levels of R-loops in CIN612 cells in comparison to normal human keratinocytes (HFKs) [20,32]. Furthermore, R-loops formed on HPV genomes as well as cellular sites, and these high levels were found to be critical for viral transcription as well as replication. Elevated levels of R-loops have also been detected by immunofluorescence and immunohistochemistry analyses of biopsies from HPV positive cervical cancers [20,33]. In this study, we examined how the landscape of R-loops on cellular sites changes between HPV positive cells and normal keratinocytes as well as whether these alterations have functional consequences on cellular gene expression and HPV pathogenesis.

Results

To investigate how the distributions and functions of R-loops change due to the presence of high-risk HPV genomes, we examined cells derived from a biopsy of an HPV 31 positive precancerous cervical lesion (CIN 612) and compared effects in normal keratinocytes (HFKs). Included in this initial analysis was the HFK-31 cell line that was generated by transfection of HFKs with cloned HPV 31 sequences and maintains viral sequences as episomes. Both HPV positive cell lines have been shown to exhibit similar histological changes in organotypic raft cultures consistent with CIN I lesions in vivo [31,34]. The levels of total R-loops in these cells were measured by dot blot assays that utilize the S9.6 antibody, which is specific for R-loops (Fig 1A). This analysis demonstrated that total R-loop levels in both HPV positive cell lines were significantly increased compared to normal keratinocytes. RNase H treatment abrogated these signals, demonstrating the assay was specific for R-loop formation. R-loops often form at promoter as well as transcription termination sequences, and the examination of these regions on a series of representative genes by DRIP-qPCR demonstrated substantially increased levels in HFK-31 and CIN 612 cells (Fig 1B). In this analysis, the levels of R-loops at previously reported sites in MYADM, RPL13a, SLC35B2, and LGAL2 were found to be increased on average by 5-to-10-fold relative to levels detected at the same sites in normal keratinocytes [12,16,35–37]. In addition, sites with low or negligible levels of R-loops associated with genes such as EGR1 and SNRPN in normal keratinocytes showed minimal increases in CIN 612 cells (Fig 1B). Increases of over 500-fold in R-loop levels were also detected in association with ALU elements, which may account for the higher total levels seen by dot blot analysis (Fig 1C). High levels of R-loops were also detected on HPV genomes at the early promoter and termination sites but not at coding sequences or the late poly A site (Fig 1C). These observations indicate there are substantial increases in R-loop levels in cells derived from HPV positive precancers or generated by transfection in comparison to normal keratinocytes and are consistent with previous reports [20].

10.1371/journal.ppat.1012454.g001 Fig 1 Distribution of R-loops in HPV positive cells compared to normal keratinocytes.

(A) S9.6 dot blot analysis from whole cell nucleic acid extracts of normal keratinocytes (HFKs), transfected, HPV positive keratinocytes (HFK-31), and HPV positive keratinocytes derived from a cervical CIN 1 lesion (CIN 612). Total nucleic acid levels were measured via methylene blue staining (top), and the specificity of the S9.6 monoclonal antibody was assessed through RNase H treatment. (B) DNA:RNA immunoprecipitation assays (DRIP) were performed on HFKs, HFK 31, and CIN 612 cells for six representative cellular sites, and immunoprecipitated chromatin was analyzed by quantitative PCR (qPCR). Primers mapped to EGR1, RPL13a, SLC35B2, and LGALS2 were used as positive controls for regions previously characterized to contain R-loops; SNRPN was used as a negative control, while MYADM has variable reports of its association with R-loops [12,16,35–37]. Fold enrichment for each primer set over HFKS9.6 is plotted: (S9.6x/IgGx)/ (S9.6HFK/IgGHFK) where x is Ct values from either HFK 31 or CIN 612 cells. The error bars represent the standard error of the mean (n = 3, ns, not significant; p<0.05, *; p<0.001, ***; p<0.0001, ****). (C) DRIP-qPCR of three regions on the HPV 31 genome in comparison to ALU repetitive cellular elements. DRIP-qPCR was performed on HFKs, HFK 31, and CIN 612 cells using primers mapping to ALU elements and viral genomic elements (early polyA site, upstream regulatory region (URR), and the late polyA site). Percentage input was plotted: Input % = 100/2(ΔCt [normalized to input control]). The error bars represent the standard error of the mean (n = 3, p<0.001, ***; p<0.00001, ****). (D) Metaplot distribution of S9.6 signal (IP—input) through genic regions, including 2kb flanking upstream or downstream (n = 2, top). Heat map of S9.6 intensity through genic and 2kb flanking regions (n = 2, bottom). (E) Depth graphs of input normalized S9.6 reads through MYADM, LGALS2, RPL13a, and ZNF554 in normal keratinocytes (black) and CIN 612 cells (red).

DRIP-sequencing (DRIP-seq) was next performed to investigate how the distributions of R-loops varied between HPV positive cells and normal keratinocytes. This method allows for an unbiased approach to identify where R-loops are present within cells utilizing immunoprecipitations with the S9.6 monoclonal antibody followed by NEXTGen sequencing [38]. We focused this analysis on CIN612 cells in comparison to HFKs. Metaplot analysis of R-loop distribution of 2kb upstream and downstream of coding sequences in normal keratinocytes (HFK) and CIN 612 cells demonstrated that R-loop reads in both cell types peak near the transcription start site (TSS), at the transcription end site (TES), and about 1–1.5kb downstream of the TES. This distribution is similar to the profile published by Promonet et al. [12]. For our downstream analyses, R-loops within these regions were associated with a gene’s coding region and referred to as genic R-loops. Importantly, the overall R-loop distributions at TSS and TES sites are similar in CIN 612 cells (Fig 1D, left) as well as in normal keratinocytes (Fig 1D, right). The primary difference between the two cell types was that the signal in the CIN 612 cells was significantly higher than in the normal keratinocytes. Heatmap analysis demonstrated that genes with high levels of R-loops at the TSS also had high levels of R-loops at the TES (Fig 1D, bottom). Only a minority of genes exhibited distinct patterns of R-loop formation in CIN 612 compared to normal keratinocytes. Examples of R-loops distributions showing the IP/input enrichment on four different genes are shown in Fig 1E. RPL13a, MYADM, and LGALS2 all contained significantly higher R-loops levels over the input background control in the precancerous CIN 612 cells than in normal keratinocytes, consistent with our DRIP-qPCR analysis (Fig 1E to 1B). ZNF554 is a gene not associated with R-loop formation, and it exhibited minimal R-loop levels over the input background in both cell lines examined. S1 Fig shows similar analyses for 2 additional genes (DNA lig IV and CALML5). This overview indicates that enhanced levels of R-loops are present in CIN 612 cells compared to normal keratinocytes, and these form at similar, though not necessarily identical, regions in the proximal 2kb upstream and downstream regions.

Identification of genomic regions enriched with R-loops

Further analysis of the DRIP-seq data was then used to provide an overall picture of which sites were associated with enhanced R-loop formation in CIN 612 cells as compared to normal keratinocytes. Peak calling analysis demonstrated that R-loops were significantly enriched over background at over 90,000 sites in HPV positive cells and at approximately 40,000 sites in normal keratinocytes. About 30,000 sites were shared between both cell lines, leaving over 60,000 unique R-loop sites in CIN 612 cells and approximately 9,500 unique sites in normal keratinocytes (Fig 2A). MA plot analysis of the ~30,000 common sites indicated significantly higher R-loop levels at those sites in CIN 612 cells compared to normal keratinocytes, with some R-loop sites being enriched by almost ~1000 fold in the HPV positive cells (Fig 2B). Although there are significant differences in the numbers of R-loop reads between CIN 612 cells and normal keratinocytes, peak distribution over promoter, exonic, and downstream sequences were very similar, with approximately 60% of the R-loop peaks detected in genic regions in both cell types (Fig 2C).

10.1371/journal.ppat.1012454.g002 Fig 2 R-loops form preferentially on genes in pathways responsible for the cancer progression and viral pathogenesis in CIN 612 cells.

(A) Venn diagram of the genomic regions containing R-loop peaks (MACS) overlapping between normal keratinocytes and HPV positive cells (CIN612) (n = 2, a representative image is shown). Total R-loops in CIN 612 (pink) and HFK (red) with common sites (brown). (B) MA plot analysis of the common R-loop containing genes between CIN 612 cells and normal keratinocytes. Log2 enrichment of R-loop levels over the matched input control is plotted on the x-axis, and Log2 enrichment of R-loop levels in genes present within CIN 612 cells over HFKs is plotted on the y-axis. (C) Distribution of R-loop peaks relative to genomic locations in HFKs and CIN 612 cells. CHIPSEEKER was used to analyze the location of R-loop reads within each sample. (D) HOMER was used to identify the location of where R-loop peaks occurred within HFK and CIN 612 cells. Intergenic R-loops were filtered out, leaving R-loops that fell within introns, exons, TES, TSS, 3’UTR, and 5’UTR. Common genes found in both HFKs and CIN 612 cells were also filtered out. Pathway analysis was then performed on the genes to which these R-loops were assigned to either the CIN 612 cells or the HFKs using Shiny GO 0.80. KEGG pathways or molecular function analyses are shown.

Since the levels of R-loops were substantially increased in HPV positive cells relative to normal keratinocytes, it was possible they were linked to specific genes or pathways. For this analysis, we first examined genes associated with R-loops in both normal keratinocytes and CIN 612 cells, focusing on proximal promoter, gene body, and terminator regions. Pathway analysis on the R-loop containing genes was then performed using ShinyGO 0.80 [39]. The KEGG pathways associated with high levels of R-loops in CIN 612 cells included those involving cancer progression (pathways in cancer, proteoglycans in cancer, and transcriptional misregulation in cancer) and DNA virus infection, many of which were not found to be enhanced in normal keratinocytes (Fig 2D, top). In contrast, pathway analysis of genes containing R-loops in normal keratinocytes identified prominent pathways involved in clathrin binding, lipid binding, and kinase activity (Fig 2D, bottom). These analyses indicated that the increased formation of R-loops in CIN 612 cells is localized to genes in pathways that are distinct from those seen in normal keratinocytes. It was next important to determine if this increased R-loop association correlated with increased transcription.

HPV positive cells have similar numbers of genes upregulated and downregulated, despite high R-loop levels

In order to determine if there was a correlation between high levels of R-loops and increased transcription of the associated genes, RNA sequencing was performed on CIN 612 [20] and normal keratinocytes. This analysis demonstrated that approximately 20% (~4,500) of the genes analyzed were differentially expressed compared to normal cells. Interestingly, these genes were divided almost equally between those upregulated (2,207) and those downregulated (2,280) (Fig 3A). A similar distribution in the fold changes of differentially expressed genes was seen by MA plot analysis (Fig 3B). We then performed pathway analysis of the differentially expressed genes in CIN 612 cells compared to normal keratinocytes to determine which pathways were altered (Fig 3C). Consistent with previous findings, genes in pathways involved in DNA repair, DNA replication, and cell cycle were upregulated, while those in pathways important for epithelial differentiation and epidermal development were downregulated [40,41]. In addition, pathways involved in the immune response and interferon signaling were downregulated in CIN 612 cells.

10.1371/journal.ppat.1012454.g003 Fig 3 Formation of R-loops at unique sites in HPV positive, CIN 612 cells correlates with differential gene expression.

(A) Venn diagram of the differentially expressed genes in precancerous CIN 612 cells compared to normal keratinocytes (n = 2, RPKM > 0) Genes upregulated in CIN 612 (green), upregulated in HFKs (red), and those with no difference (<1 Log2 FC, brown). (B) MA plot analysis of the differentially expressed genes exhibiting a similar distribution of downregulation and upregulation. (C) The differentially expressed genes were divided between those upregulated and downregulated in CIN 612 cells compared to HFKs. Pathway analysis of the biological processes of these genes was performed using Shiny GO 0.80. (D) mRNA levels of genes R-loop positive or negative in HFKs and CIN 612 cells. The line represents the mean (p<0.0001; ****). (E) Around 25% of all differentially expressed genes are associated with R-loops only in precancerous CIN 612 cells. R-loop peaks in the genic or 2kb flanking regions of genes in HFKs or those that were common to both HFKs and CIN 612 cells were filtered out. The remaining genes were screened against the differentially expressed genes in CIN 612 cells compared to HFKs. Fold enrichment of S9.6 reads over input is plotted on the X-axis, while fold change of mRNA levels in CIN 612 cells compared to HFKs is plotted on the Y-axis (left). (F) Pathway analysis of the R-loop containing genes upregulated (left) and downregulated (right) in CIN 612 cells.

DRIP-seq analyses were then used to investigate if altered expression levels were linked to increased R-loop levels. This analysis demonstrated an approximately 2-fold higher level of total transcripts associated with genes that were linked to R-loops compared to those without R-loops (Fig 3D). Importantly, approximately 30% of the differentially expressed genes were found to be associated with R-loops only in CIN 612 cells and absent in normal keratinocytes (Fig 3E). These differentially regulated genes were similarly distributed between those downregulated (672) and those upregulated (594). While there was not a consistent global increase or decrease in the expression of genes that were uniquely associated with R-loops in CIN 612 cells, KEGG analysis identified genes in specific pathways that were coordinately regulated (Fig 3F). The most prominent pathways associated with enhanced R-loop levels were involved in replication and DNA metabolism. At the same time, genes associated with the innate immune surveillance pathway, including IL1B, STAT1, and MYD88, were linked to enhanced R-loop levels but exhibited decreased expression in CIN 612 cells relative to normal keratinocytes. Furthermore, no correlation was found between the enrichment of R-loops and their corresponding mRNA levels with respect to whether these structures formed at either TSS, TES, or intronic locations (S2 Fig).

The above studies indicated there was a correlation between the presence of unique R-loops in CIN 612 cells and altered expression of genes in specific pathways. It was next important to determine if their expression was functionally dependent upon enhanced levels of R-loops. For this analysis, we utilized CIN 612 cells that were generated to overexpress RNase H1, an R-loop processing enzyme, by transfection of CMV-directed tagged expression vector followed by selecting stable cell lines [42]. Overexpression of RNAse H1 has been characterized as the “gold standard method” for reducing R-loop levels [43]. When RNase H1 was overexpressed in CIN 612 cells, viral early gene expression and episome levels were reduced by ~70% and 50%, respectively. As HPV 31 E6 was shown to induce R-loop formation in HFKs, this could further enhance these reductions [20]. RNA-seq analysis was previously performed on these cells and compared to that seen with parental CIN 612 cells [20]. The cells overexpressing RNase H1 exhibited substantially reduced levels of R-loops compared to the parental line. Around 12% of all genes (FKPM > 0) were differentially expressed within CIN 612 cells overexpressing RNase H1 compared to the scramble control cells (S3A Fig). Furthermore, 833 of the genes were associated with R-loops only in CIN 612 cells, further supporting a direct functional relationship between R-loop formation and gene expression (S3B Fig).

Pathway analysis of genes whose expression was dependent upon R-loops present only in CIN 612 cells were linked to innate immune surveillance, including interferon-alpha and interferon-gamma responses, complement signaling, and inflammation (Fig 4A). Genes in these innate immune surveillance pathways that are repressed by R-loops in CIN 612 cells include STAT1 (8-fold), NLRP3 (14-fold), JAK2 (4-fold), AIM2 (34-fold), RIG-I (22-fold), IFNB (greater than 50-fold), and IL6 (190-fold). In contrast, TRIM 14 and STING are only modestly repressed (Fig 4B). The expression of some genes, such as IFN β and RIG-I, increased in response to reductions in levels of R-loops despite not being physically linked with these structures. This likely indicates that R-loops target their upstream regulators, so it is likely that the effects on IFN β and RIG-I are indirect.

10.1371/journal.ppat.1012454.g004 Fig 4 Reduction of R-loops through RNase H1 overexpression in CIN 612 cells identifies genes that are functionally dependent on their formation.

(A) Major pathways whose expression is dependent on R-loops present only in CIN 612 cells as determined by RNase H1 overexpression. Pathway analysis was performed using the Hallmark database in Shiny GO 0.80 of R-loop regulated genes that contain R-loops unique to the CIN 612 cells. (B) Innate immune response genes whose expression is regulated by R-loops that are present only in CIN 612 cells. Fold changes are normalized to the mRNA counts in normal keratinocytes (n = 2). Linkages to H3K36me3 and γH2AX histones are shown on the right.

Histone modifications are differentially deposited on host chromatin within HPV positive cells

The linkage of enhanced levels of R-loops with coordinated expression of genes in multiple specific pathways suggested that additional factors act to facilitate this specificity. One way R-loops could coordinate the expression of genes in distinct pathways might be through association with different sets of modified histones that configure chromatin around these structures [44]. In addition to histones linked to chromatin states, R-loops are also associated with the modified histone, γH2AX, which is coupled with DNA break formation and may be linked to gene expression [12,45,46]. Therefore, we investigated whether there are associations between specific sets of histones and R-loop-directed gene expression that vary between normal and HPV-positive cells.

For this analysis, chromatin immunoprecipitation was performed on CIN 612 cells and normal keratinocytes for three histone marks: H3K36me3, H3K9me3, and γH2AX. H3K36me3 is typically associated with transcription, while H3K9me3 marks areas of heterochromatin [47–50]. We performed peak calling algorithms on each of the modified histones pulldown experiments using the MACS peak calling algorithm. We controlled for off-target pulldowns by background subtracting respective input control samples isolated from each of the cell lines (HFK and CIN 612 cells). The called peaks were then assigned a relative genomic location using HOMER. Peak calling analysis for these histone marks focused on the regions 2kb upstream, 2kb downstream, or in the gene body in both CIN 612 and normal keratinocytes. This analysis identified an overlap of these histones with unique and common sets of genes associated with R-loops. Overall, H3K36me3 marks were approximately 4-fold more prevalent in CIN 612 cells than in normal cells (S4 Fig). In contrast, the opposite was found for H3K9me3 marks which were reduced in CIN 612 cells relative to normal keratinocytes (S4 Fig).

It was next important to determine whether the presence of H3K36me3 and H3K9me3 formation correlated with the enhanced formation of R-loops and transcription in HPV positive cells by comparing ChIP-seq analysis for these histones to DRIP-seq and RNA-seq data, respectively. RNA-seq analysis of genes containing H3K36me3 marks identified an over 2-fold enrichment in mRNA levels of these genes over those that did not contain H3K36me3 in both normal and CIN 612 cells (Fig 5A, left). In contrast, the H3K9me3 mark on genes that were linked to R-loops in normal keratinocytes correlated with a modest 0.33-fold increase in mRNA levells, with no significant difference seen in the CIN 612 cells (Fig 5A, right). Comparing the ChIP-seq analyses with DRIP-seq, over 50% of the R-loop containing genes in CIN 612 cells were also positive for H3K36me3 (Fig 5B, bottom left), while less than 8% were associated in normal keratinocytes (Fig 5B, top left). In contrast, less than 2% of genes with R-loops in CIN 612 cells were H3K9me3 positive as compared to over 10% in normal keratinocytes (Fig 5B, right). Gene ontology analysis of the H3K36me3 and R-loop positive genes in CIN 612 cells identified pathways associated with cell cycle progression and innate immune surveillance (S5 Fig, top left). Genes in the immune response pathways identified above exhibited a strong linkage between the presence of both R-loops and H3K36me3 (Fig 4B). In contrast, in normal keratinocytes, genes associated with H3K36me3 and R-loops were found in distinctly different pathways that regulate membrane potential, protein localization, and neurogenesis (S5 Fig, top right).

10.1371/journal.ppat.1012454.g005 Fig 5 Modified histones are differentially associated with actively transcribed genes and R-loops in HPV positive cells.

(A) Total mRNA levels of genes that are H3K36me3 positive or negative (left) as well as H3K9me3 positive or negative (right) in HFKs (top panels) and CIN 612 cells (bottom two panels). The red line represents the mean. The error bars are SEM (ns, not significant; p<0.0001, ****). (B) Venn diagrams of the genomic regions containing H3K36me3 (left) or H3K9me3 (right) and R-loop peaks (MACS) overlapping between the HPV negative (HFKs) and positive (CIN612) cells (n = 2, a representative image is shown). Genes containing H3K36me3 or H3K9me3 peaks (red), R-loop peaks (green), and both (brown). P-values represent that the overlap between genic marks is above or below that expected from random distribution. RF represents representation factors. RF values greater than 1 suggest more overlap than expected from random distribution, while RF values less than 1 imply less overlap.

HPV positive cells contain high levels of modified H2AX (γH2AX) and DNA breaks, which results in the constitutive activation of DNA damage repair pathways [28,51]. Consistent with these findings, peak calling analysis identified ~4-fold more γH2AX marks (21,941 vs. 4,870) in the CIN 612 cells than in the normal cells (S6 Fig). In addition, the profile of averaged γH2AX reads across genic regions of precancerous CIN612 cells differed substantially from that seen in the normal keratinocytes. While normal keratinocytes exhibited no significant increases in γH2AX reads across genic regions, γH2AX reads in HPV positive cells increased, peaking at the TES (S6 Fig). Pathway enrichment analysis of genes containing γH2AX marks in the precancerous CIN 612 cells identified those responsible for cell cycle control, regulation of RNA biosynthetic processes, and transcription. This enrichment was only seen in HPV positive cells and not in normal keratinocytes. Interestingly, γH2AX was associated with genes with ~2-fold higher transcript levels than those without γH2AX (Fig 6A). Approximately 37% of all genic R-loops in CIN 612 cells were found to be associated with γH2AX in contrast to normal keratinocytes, where negligible levels were detected (Fig 6B). It was next important to determine whether H3K36me3 or H3K9me3 were preferentially associated with γH2AX and R-loops. Interestingly, about 25% of all R-loop-containing genes in CIN 612 cells were found to be positive for the combination of R-loops, H3K36me3, and γH2AX (Fig 6C). This level of correlation was not seen in normal cells nor with H3K9me3 and γH2AX. The R-loop associated genes which were H3K36me3 and γH2AX positive, are involved in pathways essential for viral replication, including DNA break repair and cell cycle control (Fig 6C, right).

10.1371/journal.ppat.1012454.g006 Fig 6 γH2AX is associated with R-loop formation and H3K36me3 deposition in CIN 612 cells.

(A) mRNA levels of genes γH2AX that are negative or positive in HFKs (left) and CIN 612 (right) cells. The red line represents the mean. The error bars are SEM (p<0.05, *; p<0.0001, ****). (B) Venn diagrams of the genes containing γH2AX and R-loop peaks (MACS) overlapping between the HFKs (left) and CIN 612 (right) cells (n = 2). Genes containing γH2AX peaks (red), R-loop peaks (green), and both (brown). P-values represent that the overlap between genic marks is above or below that expected from random distribution. RF represents representation factors. RF values greater than 1 suggest more overlap than expected from random distribution, while RF values less than 1 imply less overlap. (C) Venn diagrams of the genes containing γH2AX, H3K36me3, and R-loops in CIN 612 cells (left) and the GO biological processes pathway analysis of these genes (right). (D) A table showing DNA repair and metabolism genes that are differentially expressed in CIN 612 cells with corresponding marks of R-loops, γH2AX, and H3K36me3, which are present on these genes only in the precancerous cells.

The association of H3K36me3, γH2AX, and enhanced R-loops was particularly significant for genes in the DNA repair pathway. HPV proteins activate the ATM and ATR DNA repair pathways, which is critical for differentiation-dependent amplification. Our studies show that genes like ATM, ATRX, RAD51C, along with members of the Fanconi Anemia pathway (FANC-B, C, E, I, L, and M), and SETD2, the methyltransferase regulating H3K36me3, were all associated with the combination of H3K36me3, γH2AX and enhanced R-loops (Fig 6D). While a significant linkage was found between innate immune regulatory genes and the presence of both R-loops and H3K36me3, only a minimal association was found for the combination of H3K36me3, γH2AX, and enhanced R-loops. This suggests there may be a preferred linkage of H3K36me3, γH2AX, and enhanced R-loops with genes in the DNA damage repair pathway. These results indicate that the linkage between all three factors, γH2AX, H3K36me3, and R-loops, is critical for HPV pathogenesis and cancer progression.

Discussion

The levels of R-loops are increased in many cancers, and how the distribution, as well as the function of these structures, change due to the presence of high-risk HPV genomes was examined by comparing cells derived from an HPV 31 positive precancerous lesion of the cervix (CIN I) to normal keratinocytes. The levels of R-loops were found to be enriched by ~5–10 fold on individual cellular genes in CIN 612 cells in comparison to normal keratinocytes. The largest enrichment of R-loops identified in HPV positive cells was, however, associated with repetitive ALU elements, which exhibited over a 500-fold increase compared to that seen in normal keratinocytes. While the levels of R-loops are significantly increased in HPV positive cells, the overall pattern of where R-loops form on cellular genes is very similar to that detected in normal keratinocytes, with peak levels located within 2kb upstream of start sites, within gene bodies, as well as 2kb downstream of termination sequences. Approximately one-third of the R-loops identified in CIN 612 cells are located at sites similar to those found in normal keratinocytes, while about two-thirds of the R-loops are associated with unique genes only in the HPV positive cells and not in normal keratinocytes. Interestingly, the expression of genes with R-loops associated only in CIN 612 cells is divided equally between those with increased or decreased transcript levels. While no global increase in expression is associated with enhanced R-loop levels, genes in specific pathways were found to be coordinately regulated. This includes pathways associated with DNA repair, DNA replication and cell cycle, whose expression is coordinately increased. Equally interesting is the identification of genes involved in innate immune surveillance and keratinocyte differentiation, which are suppressed. All these changes may contribute to progression from normal to precancerous states as well as for the pathogenesis of high-risk HPVs, which are the etiological agents responsible for cervical intraepithelial neoplasia. This indicates that the directed formation of R-loops on specific groups of genes may provide an important function in the HPV life cycle.

The repression of genes in the innate immune surveillance pathway in the CIN 612 cells is particularly sensitive to enhanced levels of R-loops. In wild type CIN 612 cells, the expression of many innate immune regulatory genes is reduced by 2 to 5-fold from that seen in normal keratinocytes (Fig 4B). The stable overexpression of RNase H1 in CIN 612 cells resulted in increased expression of innate immune genes, including STAT1 (8-fold), NLRP3 (14-fold), AIM2 (34-fold), IL6 (190-fold), and IFNβ (over 50-fold), indicating their repression may be functionally linked to R-loop levels. Importantly, R-loops are only found to be associated with these cellular genes in the CIN 612 cells. Some of these increases exceed the amounts seen in normal keratinocytes, suggesting that multiple upstream regulators of these factors also depend on R-loop formation. Furthermore, the expression of several genes, such as RIG-I and TRIM 25, are increased by RNase H1 overexpression despite not being linked to the presence of R-loops. These are interferon stimulated genes, and the increases in expression are likely the result of the enhanced levels of IFNβ that are induced when R-loop levels are reduced [52–57]. Along with increased expression of innate immune genes, RNase H1 overexpression also reduces viral gene expression and episome levels by ~70% and 50%, respectively. Whether these reductions are due to increased expression of innate immune regulators or a direct effect due to loss of R-loops on viral episomes is unclear. In contrast to the repression of the innate immune regulatory pathway, genes in the DNA damage repair pathway are activated by enhanced levels of R-loops in CIN 612 cells. This includes genes such as ATM, Top2A, Lig1, and RAD51C. While the differences seen with the DNA damage repair genes are not as dramatic as seen with the innate immune regulated genes, the activation of the DNA damage repair pathway in HPV positive cells has been shown to be critical for viral pathogenesis and cancer progression.

Our observation that R-loops are found in association with specific sets of genes that are linked with both increased and decreased expression indicates that their formation is not merely an accidental byproduct of increased transcription but is instead the result of a directed process. One way that expression could be linked with enhanced levels of R-loops is through altered chromatin states associated with specific sets of modified histones. Previous studies have suggested an association of H3K36me3 and H3K9me3 with certain classes of R-loops, but only a limited correlation with altered expression has been described [2,58]. Our studies demonstrated that over half of R-loop associated genes in CIN I derived cells are associated with H3K36me3 marks, while only 8% are positive in normal keratinocytes. H3K36me3 has been linked with increased transcription; however, in our study, equal numbers of dually H3K36me3 and R-loop positive genes exhibit increased expression as decreased expression compared to normal keratinocytes [59–62]. This indicates that this histone mark is more likely associated with an accessible chromatin configuration rather than increased transcription alone. Both innate immune regulatory genes, as well as those in DNA damage repair, are linked with high levels of H3K36me3, and this is not seen in normal keratinocytes, demonstrating that this effect is specific to HPV positive cells. SETD2 is the methyltransferase that regulates the deposition of methyl groups on lysine 36 of histone 3 (H3K36me3), and its levels are increased in CIN 612 cells as well as other HPV positive cells [63–65]. Knockdown of SETD2 in HPV positive cells has been shown to lead to significant reductions in viral episomes, identifying it as an important regulator of viral persistence. While H3K36me3 has been identified as a mark of transcription elongation, recent studies have also linked it with DNA repair suggesting a potential link with genomic instability and DNA breaks [66–68]. A previous study linked cells with high R-loop levels to concomitant decreases in H3K9me3 levels, consistent with our studies, as CIN 612 cells contained far fewer of these marks than normal keratinocytes [69]. In contrast, no strong linkage was found between H3K9me3 and R-loop regulated gene expression in CIN 612 cells. Only 2% of R-loop associated genes were also positive for H3K9me3 as compared to 10% in normal keratinocytes.

The failure to resolve R-loops leads to the formation of DNA breaks and genomic instability [70]. HPV positive precancers, as well as other cancers, exhibit high levels of DNA break formation as indicated by enhanced amounts of γH2AX, which is often used as a surrogate marker [71]. In CIN 612 cells, high levels of γH2AX are associated with increased levels of R-loops at genes whose expression is altered. Over one third of the genes associated with R-loops in the HPV positive cells were also positive for γH2AX. In addition, 67% of the genes positive for both γH2AX and R-loops were also linked to H3K36me3. No such associations are seen in normal keratinocytes. Approximately 700 of the genes that are differentially expressed in the CIN 612 cells are linked to the combined presence of γH2AX, H3K36me3, and R-loops. Genes whose expression is positively regulated by R-loop formation and associated with both γH2AX and H3K36me3 include ATM, ATRX, ATR, Top2A, and RPA3. At the same time, genes negatively regulated by R-loops that are also H3K36me3 and γH2AX positive include JAK2 and TRIM 14. The association of γH2AX and R-loops with DNA damage repair genes may be important but the mechanism responsible is unclear. Recent studies have suggested that γH2AX might not only interact with sites of endogenous DNA breaks but also associate with DNA intermediates that form upon chromatin opening during transcription initiation [72,73]. The increased expression of genes linked with the combination of γH2AX, R-loops, and H3K36me3 in HPV positive cells compared to normal keratinocytes supports this model.

These studies identify a potential linkage between R-loops, specific histone marks, and altered transcription. However, additional factors must act to determine how genes in specific pathways are targeted. One such possibility may be the association with other non-β DNA structures like G-quadruplexes and GC skew. The relationship between G-quadruplex formation and stability of R-loops has been noted in multiple reports, and may contribute to effects in HPV positive cells [11,74,75]. Similarly, a GC skew has been reported in a number of R-loops, and a preliminary screening indicates that some but not all R-loops associated with innate immune genes have this skew, identifying an important area for future studies. In addition to structural motifs in DNAs, we have shown that inhibition of p53 leads to increased levels of R-loops in HPV positive cells, cells and has been reported in other tumor cell lines that have mutated p53 [20,69]. This indicates that factors downstream of p53 play important roles in regulating R-loop formation and that this occurs at specific sites on cellular genes. Transient inhibition of p53 in normal keratinocytes alone is, however, not sufficient to induce increased R-loop formation but our studies have shown the requirement of HPV E7 co-expression, which implicates inhibition of Rb proteins as a possible contributing factor. Additional factors that could be downstream mediators of the p53 effects on R-loop formation include members of the p21-DREAM complex, long non-coding RNAs, and APOBEC 3B proteins. In embryonic mouse stem cells, a subset of polycomb group genes was shown to be linked with R-loop formation, and overexpression of RNase H1 increased their expression, indicating a repressive effect of R-loops [76]. At the same time, RNase H1 overexpression led to decreased expression of other polycomb genes and this differential regulation is similar to our results. This R-loop dependent activity requires the cooperative action of cellular proteins, and we believe that additional factors, including modified chromatin as well as transcription factors, can provide comparable functions in HPV positive cells. It is also possible that viral proteins can contribute to regulating the expression of R-loop associated genes. Overexpression of RNAse H1 in HPV-positive cells decreased the expression of viral early genes [20], and this reduction in viral proteins could potentially impact the expression of cellular genes that are linked to the presence of R-loops at these loci.

Overall, these observations demonstrate that R-loop levels are significantly elevated within HPV positive cells compared to normal keratinocytes. While no global effect on gene expression is seen due to increased levels of R-loops, genes in pathways that are important for viral replication and cellular transformation are coordinately activated or repressed by these structures, possibly in cooperation with the recruitment of specific types of modified histones. Our studies indicate that in HPV-positive cells, R-loops contribute to regulating cellular and viral gene expression during HPV pathogenesis, including those involved in the innate immune response and DNA damage repair.

Materials and methods

Reagents

Antibodies used in these experiments were as follows: S9.6 (Millipore), Anti-Histone H3 (tri methyl K36) antibody—ChIP Grade (Abcam), Anti-Histone H3 (tri methyl K9) antibody—ChIP Grade (Abcam), Phospho-Histone H2A.X (Ser139) (D7T2V) Mouse mAb (Cell Signaling), and Mouse IgG (Diagenode). Methlyene Blue Hydrate (Sigma) was used for staining nucleic acids in dot blot assays. RNase H (ThermoFisher) was used to remove R-loops from nucleic acid extracts to determine specificity of the S9.6 antibody. Mung Bean Nuclease was purchased from New England Biologicals and was used for enzymatic digestion of samples during chromatin immunoprecipitation- and DNA:RNA immunoprecipitation-sequencing.

Cell culture and reagents

Isolation of HFKs

Neonatal human epidermis was supplied by the Skin Disease and Research Core at Northwestern University. These de-identified tissues were suspended in Hanks’ balanced salt solution (HBSS), and isolations were performed within 3 to 4 days of circumcision. The foreskins were washed in phosphate-buffered saline (PBS) before being processed. Excess blood vessels, tissue, and fat were cleaned away before being incubated overnight at 4 C in 2.4 U/ml Dispase. The following day, the epidermis was removed and incubated with 4 ml of 0.25% trypsin for 15 min. The epidermis was then scraped vigorously for 2 to 3 min before quenching the trypsin with bovine serum. The resulting suspension was then pipetted through a 40mm pore cell sieve. The cells were then spun down and resuspended in E-medium supplemented with 5ng/ml of mouse epidermal growth factor (EGF). NIH 3T3-J2 fibroblasts, growth-arrested through treatment with mitomycin-c, were seeded with the newly collected human foreskin keratinocytes (HFKs), and media was changed as required until the proliferation of the keratinocytes was achieved.

Cell culture

HFKs, HFK-31, and CIN 612 cells were all cultured in E-medium supplemented with 5ng/ml of mouse EGF. Each of these cell lines were co-cultured with NIH 3T3-J2 fibroblasts, which were growth arrested using 0.4mg/ml of mitomycin-c for at least 2 hr. To remove J2 fibroblasts prior to downstream analyses, cells were washed with Versene (0.05mM EDTA PBS) for 5 min before 2 sequential PBS washes. J2 feeders were cultured in DMEM containing 1% penicillin-streptomycin and 10% bovine serum. Cells stably overexpressing RNase H1 were generated previously [20].

Generation of cells stably maintaining HPV 31 episomes

The pBR-322min-HPV31 plasmid was digested such that the pBR-322 backbone was removed, leaving the HPV 31 genome which was recircularized. 1 μg of recircularized HPV 31 DNA was contransfected with a selection plasmid expressing a neomycin resistance cassette (PSV2neo) into around one million freshly isolated HFKs at ~60% confluence. The following day, cells were selected using 200 mg/ml G418. J2 feeders were changed on alternating days as the G418 selection. Stable maintenance of HPV 31 episomal DNA was assessed by Southern blot before expanding and performing downstream analyses on these cells.

S9.6 dot blot assay

DNA was purified from cell lysates using PhenolChloroform extractions. Samples were either left untreated or treated with 1U of RNase H for at least 1.5 hr at 37° C. DNA was then spotted onto a positively charged membrane (Zeta-probe). Membranes were then stained with Methylene blue for ~15 min before being washed with di-H2O 3 times for 5 min. Images of the Methylene blue staining were acquired to normalize to total nucleic acid content using an Odyssey Fc LiCor (LiCor BioSciences). Methylene blue staining was removed through washing with 100% ethanol for 5 min before washing with di-H2O 3 times for 10 min. Membranes were then blocked with 5% Bovine Serum Albumin (BSA) in TBST (Tris-buffered saline Tween 20) before being probed with the S9.6 anti-RNA:DNA hybrid antibody (Millipore) overnight at 4 °C. The following day, membranes were washed with TBST, probed with secondary antibody for 1 h at RT, and developed using enhanced chemiluminescence (ECL) (Fisher, 4500085). Images were taken using an Odyssey Fc LiCor.

DNA:RNA immunoprecipitation (DRIP)–qPCR

1 × 107 cells were harvested and collected in Southern lysis buffer before being treated with RNase A (5 ng/mL) and Proteinase K (7.5 ng/mL) at 37 °C overnight. DNA was purified from these samples using phenol-chloroform extractions, and 25 to 50 mg of DNA was used for each sample. DNA was sheared using a Bioruptor (Diagenode) on high power, 30 s on/90 s off cycles for 20 min or digested using 1U of mung bean nuclease for 1 h at 37 °C. Input DNA was removed before loading the samples into preblocked magnetic beads in IP buffer containing 2 mg of the RNA:DNA hybrid antibody. Immunoprecipitations were allowed to incubate overnight at 4 °C while rotating. The next day, samples were washed 8 times with RIPA buffer for 5 min while rotating. One wash in TE buffer was performed before samples were eluted for 10 min at 65 °C in 10% sodium dodecyl sulfate (SDS), 10 mM Tris pH 7.4, 50 mM ethylenediaminetetraacetic acid (EDTA). DNA was purified from these elutions using a PCR purification kit (Qiagen) and stored at −20 °C. Primer sets used to analyze S9.6 immunoprecipitated sequences are listed in the Key Resources table (Table 1).

10.1371/journal.ppat.1012454.t001 Table 1 Primer sets used in this study.

Primers	
DNA:RNA immunoprecipitation—qPCR primers			
Cellular regions			
MYADM	5’ CGT AGG TGC CCT AGT TGG GAG 3’	5’ TCC ATT CTC ATT CCC AAA CC 3’	
RPL13a	5’ AAT GTG GCA TTT CCT TCT CG 3’	5’ CCA ATT CGG CCA AGA CTC TA 3’	
EGR1	5’ GAA CGT TCA GCC TCG TTC TC 3’	5’ GGA AGG TGG AAG GAA ACA CA 3’	
SLC35B2	5’ AAG TCT TGC CCT AGC TGT GCT 3’	5’ GCC TAC ACC GCT TGT GCT TTT 3’	
SNRPN	5’ GCC AAA TGA GTG AGG ATG GT 3’	5’ TCC TCT CTG CCT GAC TCC AT 3’	
LGALS2	5’ TGA CCT CAC CTT GAC CTC TGA 3’	5’ AGC TGA ACC TGC ATT TCA ACC 3’	
ALU elements	5’ ACG AGG TCA GGA GAT CGA GA 3’	5’ CTC AGC CTC CCA AGT AGC TG 3’	
HPV 31 genomic regions			
Early PolyA	5’ GGT ATT GGT ATT GGT ATT GG 3’	5’ ACC CAT ACT ACC ATA CCT TA 3’	
Late PolyA	5’ GCG TGT GTA CTT GTA 3’	5’ GCA ACC GAA AAC GGT TAG G 3’	
Upstream regulatory regions (URR)	5’ GAT GCA GTA GTT CTG CGG TTT 3’	5’ TAT GTT GGC AAG GTG TGT TAG G 3’	

DRIP-sequencing

The same protocol was used to prepare samples for DRIP-sequencing as listed above for DRIP-qPCR. Samples were stored at -80° C until being shipped to Admera Biosciences (NJ), who performed the sequencing experiments. Briefly, the library was prepared using a KAPA HyperPrep Kit (Kapa Biosystems) following the manufacturer’s recommendation. Input DNA was end-repaired and 3’-dA tailed. Adapter was then ligated to the DNA, and the ligated product was PCR amplified and cleaned up using the SPRIselect Reagent (Beckman Coulter). Quality control was then performed for the final library, followed by sequencing.

DRIP-sequencing data analysis

Admera Biosciences (NJ) performed most of the bioinformatic analyses from our DRIP-sequencing experiments. Their bioinformatics methods are as follows: An in-house bioinformatics pipeline was used to analyze DRIP-Seq data. First, FastQC (v0.11.8) was used to check the quality of raw and trimmed reads. Trimmomatic (v0.38) was used to cut adapters and trim low-quality bases with a default setting. BWA (v0.7.10-r789) was used to map the trimmed reads to the reference genome* using the Burrows-Wheeler Alignment algorithm (BWA-MEM). Mapped reads that have low-quality MAPQ score (MAPQ < 10), not-properly-paired, or duplicated (assessed with Picard tools (v 2.20.4)) were removed. BAM was used to generate BW format (normalized by RPKM) for visualization. MACS (v2.2.4) was chosen to call peaks. If there was no replicate, the R package MAnorm (v2.2.6) was used for sample comparison. On the other hand, if there were replicates, their called peaks were merged and the DiffBind package (v2.14.0) was then used for differential analysis. Peak annotation and combined density profiles were performed by the ChIPseeker package (v1.22.1) and deepTools, respectively.

We performed the profile analysis of multiple DRIP-seq replicates from HFKs and CIN 612 cells (Fig 1D) and R-loop read distribution (Fig 2C) using the open-source Galaxy servers https://usegalaxy.org/. BAM Compare (Galaxy Version 3.5.4+galaxy0) was used to normalize either log2 IP to input ratios, or input subtracted IP reads for both biological replicates of DRIP-seq in HFK and CIN 612 cells. ComputeMatrix was used to prepare files for visualization via plotProfile and plotHeatmap (Fig 1D). CHIPseeker (Galaxy Version 1.28.3+galaxy0) was used on the BED files generated by Admera Biosciences to determine the genomic distribution of R-loop peaks. To confirm agreement between our biological replicates (S7 Fig), correlations between replicates were assessed using multiBAMSummary (Galaxy Version 3.5.4+galaxy0), and then plotting principal component analyses using plotPCA (Galaxy Version 3.5.4+galaxy0) and plotting Pearson coefficients as a heatmap using plotCorrelation (Galaxy Version 3.5.4+galaxy0). plotFingerprint (Galaxy Version 3.5.4+galaxy0) was used to determine narrow versus broad distributions of S9.6 reads across the genome.

RNA-sequencing

HFKs and CIN 612 cells were grown to confluency on 10cm dishes before removing J2 fibroblasts. Cells were scraped and centrifuged before being stored at -80° C before shipping to Admera Biosciences (NJ).

RNA-sequencing data analysis

FastQC (version v0.11.8) was applied to check the quality of raw reads. Trimmomatic (version v0.38) was applied to cut adaptors and trim low-quality bases with default setting. STAR Aligner version 2.7.1a was used to align the reads. Picard tools (version 2.20.4) was applied to mark duplicates of mapping. The StringTie version 2.0.4 was used to assemble the RNA-Seq alignments into potential transcripts. The featureCounts (version 1.6.0)/HTSeq was used to count mapped reads for genomic features such as genes, exons, promoter, gene bodies, genomic bins and chromosomal locations. The De-Seq2 (version 1.14.1) was used to do the differential analysis. Pathway analyses were performed using Shiny GO http://bioinformatics.sdstate.edu/go/.

Chromatin immunoprecipitation (ChIP)-sequencing

Formaldehyde was added to 1 × 107 cells to a final concentration of 1% for 10 min at room temperature. Excess formaldehyde was quenched upon adding 0.125M glycine before then washing samples with PBS. Cells were then incubated in collection buffer (0.1M TrisHCl pH 9.4 and 10mM DTT containing Roche Protease Inhibitor Cocktail) for 10min on ice. Cells were then collected and spun down before being sequentially washed and incubated with NCP1 (10mM EDTA, 0.5mM EGTA, 10mM HEPES pH 6.5, 0.25% Triton X100) and NCP2 (1mM EDTA, 0.5mM EGTA, 10mM HEPES, and 200mM NaCl) before being lysed in 0.5% Empigen BB, 1% SDS, 10mM EDTA, 50mM Tris HCl pH 8.0 containing Roche Protease Inhibitor Cocktail for 30 min on ice. Samples were then sonicated using a Bioruptor (Diagenode) on high power, 30 s on/90 s off cycles for 20 min. After sonication, samples were prepared exactly as described above in the DRIP-qPCR protocol. Samples were stored at -80 C before being sent off for sequencing either by Admera Biosciences (NJ) or the NUseq facility at Northwestern University.

ChIP-sequencing data analysis

Samples were either analyzed as described above in the DRIP-sequencing analysis section or the NU seq core delivered BAM files. Agreement between biological replicates was assessed using multiBAMSummary (Galaxy Version 3.5.4+galaxy0), and then plotting principal component analyses using plotPCA (Galaxy Version 3.5.4+galaxy0) and plotting Pearson coefficients as a heatmap using plotCorrelation (Galaxy Version 3.5.4+galaxy0) (S8 Fig). From the BAM files, BAM Compare (Galaxy Version 3.5.4+galaxy0) was used to normalize log2 IP to input ratios for both biological replicates of H3K36me3, H3K9me3, and γH2AX ChIPs from HFK and CIN 612 cells. ComputeMatrix was used to prepare files for visualization via plotProfile and plotHeatmap. CHIPseeker (Galaxy Version 1.28.3+galaxy0) was used on the BED files generated by Admera Biosciences or NUseq to determine the genomic distribution of each modified histone. MACS (v2.2.4) was used to call peaks. HOMER (Galaxy Version 4.11+galaxy0) was used to annotate where peaks occurred relative to their genomic location (intron, exon, etc.) and the corresponding gene name. Gene lists from HOMER were compared between DRIP, H3K36me3 ChIP, H3K9me3 ChIP, and γH2AX ChIP to obtain the overlap depicted in the Venn Diagrams. RNA counts from RNA sequencing experiments were compared to genes containing the corresponding mark (R-loops, H3K36me3, γH2AX or H3K9me3) to determine the association between mRNA levels seen in Figs 3D, 5A and 6A.

Quantification and statistical analysis

GraphPad prism was used for all statistical analyses, and all data are represented as mean +/- standard error (SEM). Two-way ANOVA and two-tailed T-tests were used to calculate p-values. Calculation of the representation factor and the associated probability of Venn diagram overlaps in Figs 5B and 6B were performed using http://nemates.org/MA/progs/overlap_stats.html from the Lund Lab. A genome size of 63,755 (CHESS database, http://ccb.jhu.edu/chess) was used to determine representation factors for the Venn Diagrams in Figs 5B and 6B. The maximum value represented as statistically significant was p = 0.05. Additional details on quantifications like replicates are specifically stated in the figure legends and methods.

Software and algorithms

GraphPad Prism was used to generate all graphs and statistical analyses of said graphs. Adobe Photoshop and Illustrator were used for the organization and preparation of digital figures. Integrated Genome Browser (BioViz) generated depth graphs of S9.6 coverage in HFK and CIN 612 cells (Fig 1E). Galaxy community servers were used to perform many of the sequencing analyses [77].

Supporting information

S1 Fig Input normalized S9.6 reads of two regions associated with R-loops in HFK and CIN 612 cells.

(A) Depth graphs of S9.6 reads normalized to the corresponding cell line’s input control reads. Reads were binned into 300bp regions during quantification using deeptools2 (BAMcompare). Two regions are pictured: Lig4 and CALML5. The red represents CIN 612 cells, while the black represents HFKs (n = 2, mean is shown).

(EPS)

S2 Fig R-loop association with RNA at different genomic locations does not demonstrate any preferential increases or decreases in either HFKs or CIN 612 cells.

Fold enrichment of S9.6 reads over input was taken from two independent DRIP-sequencing experiments. Genes that contained R-loops in the coding sequence or the corresponding 2kb flanking regions were then analyzed for their mRNA levels from two independent RNA sequencing experiments. These genes were plotted as average mRNA read counts (y-axis) versus fold enrichment of R-loops over input (x-axis). The line of best fit was calculated using Least Squares Regression (GraphPad Prism), where dashed lines represent the 95% confidence intervals (Q = 1% for detection of outliers, red). The equation for the line of best fit is shown for each graph in the top right corner.

(EPS)

S3 Fig Pathway analysis of differentially expressed genes in CIN 612 cells overexpressing (o/e) RNase H1 using the Hallmark MSigDB database.

(A) Cumulative numbers of differentially expressed genes in the CIN 612 cells overexpressing RNase H1 versus parental CIN 612 cells (red = downregulated, green = upregulated). The corresponding pathway analysis used Shiny GO 0.80 (http://bioinformatics.sdstate.edu/go/) and the Hallmark MSigDB database (downregulated = left, upregulated = right). Pathways of particular interest were those upregulated upon overexpression of RNase H1, including those responsible for an interferon alpha/gamma response, IL6 JAK STAT3 signaling, and IL2 STAT5 signaling (right). (B) Flow chart of how R-loop dependent gene expression was determined in CIN 612 cells. Negatively regulated genes were identified as those with reduced expression in parental CIN 612 cells compared to normal keratinocytes, which then increased in expression upon loss of R-loops through RNase H1 overexpression (left). 542 genes were identified, most of which were involved in immune surveillance and signaling. Positively regulated genes were identified as those with increased expression in parental CIN 612 cells compared to normal keratinocytes, which then decreased in expression upon loss of R-loops upon RNase H1 overexpression (right). 722 genes were identified, many of which were involved in DNA metabolism. Of the 1,264 genes identified as being R-loop regulated in CIN 612 cells, 833 of them contained R-loops only in the CIN 612 cells. These genes were deemed as being functionally regulated by R-loop formation for the analyses performed in Fig 4A.

(EPS)

S4 Fig H3K36me3 and H3K9me3 are differentially present on host chromatin within CIN 612 cells compared to HFKs.

(A) Venn diagram of the genomic regions containing H3K36me3 peaks (left) and H3K9me3 peaks (right) (MACS) overlapping between the HPV negative (HFKs) and positive (CIN612) cells (n = 2, a representative image is shown). (B) Fingerprint plot of H3K36me3 and H3K9me3 distribution compared to input control in HFKs and CIN 612 cells. (C) Distribution of H3K36me3 and H3K9me3 reads relative to genomic locations in HFKs and CIN 612 cells. CHIPSEEKER was used to analyze the location of each modified histones’ reads within each sample (left and right). Deeptools2 (BAMcompare) was used to input normalize H3K36me3 (middle top) and H3K9me3 (middle bottom). Regions analyzed were set as 500bp, flanking the coding sequence, and the average genic profile was visualized (ComputeMatrix). (D) HOMER was used to identify the location of where H3K36me3 (left) or H3K9me3 (right) peaks occurred within HFK and CIN 612 cells. Intergenic histone marks were filtered out, leaving only histone marks that fell within introns, exons, TES, TSS, 3’UTR, and 5’UTR. Common genes found in both HFKs and CIN 612 cells were also filtered out. Pathway analysis was then performed on the genes to which these histone marks were assigned in the CIN 612 cells or the HFKs using Shiny GO 0.80. The GO biological process database was used for all analyses.

(EPS)

S5 Fig Pathway analyses of genes enriched with modified histones (γH2AX, H3K36me3, or H3K9me3) and R-loops in CIN 612 cells.

Genes containing a modified histone mark and R-loops unique to CIN 612 cells (left) or HFKs (right) were analyzed using Shiny GO 0.80 and the GO biological process database. Due to the lack of H3K9me3 or γH2AX and R-loop containing genes in HFKs, no analysis of those genes is depicted.

(EPS)

S6 Fig γH2AX is significantly enriched on host chromatin and on genes responsible for important processes during HPV infection in CIN 612 cells.

(A) Distribution of γH2AX reads relative to genomic locations in HFKs and CIN 612 cells. CHIPSEEKER was used to analyze the location of R-loop reads within each sample (left and right). (B) Fingerprint plot of γH2AX distribution compared to input control in HFKs and CIN 612 cells. (C) Deeptools2 (BAMcompare) was used to input normalize γH2AX reads. Regions analyzed were set as 500bp, flanking the coding sequence, and the average genic profile was visualized (ComputeMatrix). (D) Venn diagram of the genomic regions containing γH2AX peaks (right) (MACS) overlapping between HFK and CIN 612 cells (n = 2, a representative image is shown). (E) HOMER was used to identify the location of where γH2AX peaks occurred within HFK (bottom) and CIN 612 cells (top). Intergenic histone marks were filtered out, leaving only histone marks that fell within introns, exons, TES, TSS, 3’UTR, and 5’UTR. Common genes found in both HFKs and CIN 612 cells were also filtered out. Pathway analysis was then performed on the genes to which these histone marks were assigned in the CIN 612 cells or the HFKs using Shiny GO 0.80. The GO biological process database was used for all analyses.

(EPS)

S7 Fig Validation of DRIP-sequencing replicates.

(A) XY correlation plots of S9.6 reads in biological replicates from HFKs (top left) and CIN 612 cells (top right). Pearson’s coefficient is labeled on each respective scatter plot. XY correlation plot comparing HFK S9.6 reads to CIN 612 S9.6 reads with Spearman’s coefficient labeled (bottom middle). These data support that there is a strong agreement between the S9.6 pulldown replicates in HFK and CIN 612 cells and that there are higher read counts in similar genomic regions in the CIN 612 cells. (B) Heatmap of Pearson coefficient values between input controls and S9.6 pulldown assays in HFK and CIN 612 cells. A strong correlation is seen among the S9.6 pulldown assays, suggesting that a majority of the reads are located in similar genomic regions. (C) Fingerprint plot analysis of input control samples and S9.6 pulldown replicates in HFK and CIN 612 cells. Input DNA reads are broadly distributed across the genome, while S9.6 reads are enriched on a much smaller proportion of DNA. (D) Principal component analysis of input control samples and S9.6 pulldown replicates in HFK and CIN 612 cells. A high degree of clustering is seen between the S9.6 replicates from the HFK and CIN 612 cells.

(EPS)

S8 Fig Validation of Modified Histones ChIP-sequencing replicates.

(A) Heatmap of Pearson coefficient values between input controls and H3K36me3 pulldown assays in HFK and CIN 612 cells (left). Principal component analysis of input control samples and H3K36me3 pulldown replicates in HFK and CIN 612 cells (right). (B) Heatmap of Pearson coefficient values between input controls and H3K9me3 pulldown assays in HFK and CIN 612 cells (left). Principal component analysis of input control samples and H3K9me3 pulldown replicates in HFK and CIN 612 cells (right). (C) Heatmap of Pearson coefficient values between input controls and γH2AX pulldown assays in HFK and CIN 612 cells (left). Principal component analysis of input control samples and γH2AX pulldown replicates in HFK and CIN 612 cells (right).

(EPS)

The Skin Biology and Diseases Research Core at Northwestern University provided foreskin keratinocytes.

10.1371/journal.ppat.1012454.r001
Decision Letter 0
Kalejta Robert F. Section Editor
Lambert Paul F. Academic Editor
© 2024 Kalejta, Lambert
2024
Kalejta, Lambert
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version0
21 Jun 2024

Dear Dr. Laimins,

Thank you very much for submitting your manuscript "HPV induced R-loop formation represses innate immune gene expression while activating DNA damage repair pathways" for consideration at PLOS Pathogens. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Paul F. Lambert

Academic Editor

PLOS Pathogens

Robert Kalejta

Section Editor

PLOS Pathogens

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064

***********************

Reviewer's Responses to Questions

Part I - Summary

Please use this section to discuss strengths/weaknesses of study, novelty/significance, general execution and scholarship.

Reviewer #1: This report by Templeton and Laimins is the follow-up to their recent PNAS paper in which they initially describe a role for R-loops in HPV infection. Overall, this study fills in some of the remaining gaps of their PNAS study. Using an assortment of NEXTgen sequencing, the authors confirm that R-loops are significantly upregulated in HPV(+) cells and are associated with different genes compared to control HFKs. There appears to be an overlap between differentially expressed genes and R-loop formation. To show that R-loops regulate the differences in expression, the authors overexpress RNase H1 in the CIN612 cells. The authors see the expected decrease in R-loops and also observe a difference in overall gene expression compared to the vector expressing cells. However, in the previous report, the authors show that overexpression of RNAse H1 impacts the viral genome. So it is not clear to me how the authors can distinguish between effects due to R-loop resolution and reduced viral gene expression. I would like to see this addressed by the authors. For example, how does RNAase H1 alter viral gene expression?

In the last part of the manuscript, the authors correlate gene expression, R-loops, and chromatin marks. The authors use venn diagrams to calculate overlap. Previous studies formally compared the R-loop and chromatin mark ChIP-Seq peaks to demonstrate overlap (e.g., 10.1016/j.molcel.2016.05.032). The authors should consider using that approach as a more sensitive and relevant option.

Reviewer #2: In the manuscript entitled “ HPV induced R-loop formation represses innate immune gene expression while activating DNA damage repair pathways”, the authors investigated the levels of co-transcriptional RNA:DNA hybrid (R-loop) structures in normal non-cancerous keratinocytes and pre-cancerous cells containing human papilloma virus (HPV) genomes. The authors discover higher R-loop levels in HPV-positive cells both at coding genes and repetitive sequences. Interestingly, genes associated with higher R-loops showed increased H3K36me3 levels and both positive and negative impact on gene expression, with upregulated genes having a function in DNA damage repair and metabolism, and downregulated genes having a function in innate immuneresponses. Thus, R-loops may provide a critical link in HPV positive cells to allow HPV pathogenesis in the host cells.

Overall, I find this an interesting study that shows clearly differences between HPV infected and non-infected cells and the suggested role of R-loops to differentially up or downregulate certain gene categories and thereby allows the virus to exploit the host genome is an interesting concept that should be reported. However, I have a few major comments that should be addressed prior to publication:

**********

Part II – Major Issues: Key Experiments Required for Acceptance

Please use this section to detail the key new experiments or modifications of existing experiments that should be absolutely required to validate study conclusions.

Generally, there should be no more than 3 such required experiments or major modifications for a "Major Revision" recommendation. If more than 3 experiments are necessary to validate the study conclusions, then you are encouraged to recommend "Reject".

Reviewer #1: 1) RNase H1 overexpression impacts viral maintenance and transcription. This should be considered in the interpretation of the RNAsH1 overexpression experiments.

2) Consider using peak-calling algorithms to calculate overlap between R-loops and chromatin marks.

Reviewer #2: 1) The DRIP-Seq data presented in Figure 1D show 3 major peaks of R-loop accumulation: i) ~1kb upstream of TSS, ii) over the TES and iii) ~1kb downstream of TES. This seems an unusual distribution of co-transcriptional structures that should primarily track with transcriptional activity. Especially peaks about 1kb upstream of TSS is rather unusual and typically not seen in other DRIP-Seq datasets. I think it would be important to further analyze this data and compare it directly with other DRIP-Seq datasets (e.g how the distribution of peaks with promoters, gene bodies, exons, introns, intergenic regions, etc. (Figure 2C) compares with other datasets.

2) The authors conclude from the DRIP-Seq data set that “only a minority of genes exhibited distinct patterns of R-loop formation in CIN 612 compared to normal keratinocytes”. Even if only a minority, this subset of genes would be very interesting to look at in more detail but the authors didn’t further explore this category and only focused on the genes that showed enhanced R-loop formation but at the same position.

3) The DRIP-Seq is a key dataset of the paper and the authors find quantitative but arguably small differences between the CIN612 and HFK cell lines. As I could not find this information, did the authors use a spike-in for normalization between samples? How do the two biological replicates correlate with each other? These would be important quality controls to assure that the differences are biologically significant and relevant.

4) The authors show in Figure 2C that ~ 60% of the DRIP-Seq reads are distributed within genic regions in both cell lines which agrees well with current literature. However, this is not in agreement with the heatmap provided in Figure 1D where it seems that the majority of peaks are upstream and downstream of genic regions. Why?

5) It is very interesting that the genes that increase R-loop levels are associated with both up and downregulated gene expression (Figure 3), which begs the question whether R-loops are cause or consequence for this transcriptional deregulation. For example, the genes could be upregulated and R-loops form as a consequence of this higher transcriptional activity. On the other hand, more stable R-loop formation could impede transcription elongation and therefore be the cause of downregulation of genes. Can the authors further analyze the subset of genes that are either up-or downregulated and see whether these genes have for example specific sequence features such as higher GC content, GC-skew, or potential for forming other secondary structures such as G4s?

6) A surprising result is presented in Figure 5A (right panels) where the authors show that genes associated with H3K9me3 (a marker for heterochromatin) in HFKs is associated with increased mRNA levels. This is not further commented. Can the authors check again their method to call H3K9me3+ versus H3K9me3- regions and make sure the thresholds they use to binarize the data are appropriate?

**********

Part III – Minor Issues: Editorial and Data Presentation Modifications

Please use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity.

Reviewer #1: 1) Please provide the details on how the venn diagram based overlaps are calculated. Specifically, what is the 'genome' size that is used as background.

2) Venn diagrams rely on arbitrary cut-offs as to what genes to include. See point 2 above about peak calling.

3) The authors demonstrate correlations between R-loops, gene expression, and chromatin marks. They should be careful not to overstate the importance of this. Without more data, it is hard to draw strong conclusions.

Reviewer #2: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example see here on PLOS Biology: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

10.1371/journal.ppat.1012454.r002
Author response to Decision Letter 0
Submission Version1
8 Jul 2024

Attachment Submitted filename: Letter to Editors Revision 7-2.docx

10.1371/journal.ppat.1012454.r003
Decision Letter 1
Kalejta Robert F. Section Editor
Lambert Paul F. Academic Editor
© 2024 Kalejta, Lambert
2024
Kalejta, Lambert
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version1
28 Jul 2024

Dear Dr. Laimins,

We are pleased to inform you that your manuscript 'HPV induced R-loop formation represses innate immune gene expression while activating DNA damage repair pathways' has been provisionally accepted for publication in PLOS Pathogens.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Pathogens.

Best regards,

Paul F. Lambert

Academic Editor

PLOS Pathogens

Robert Kalejta

Section Editor

PLOS Pathogens

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064

***********************************************************

Reviewer Comments (if any, and for reference):

Reviewer's Responses to Questions

Part I - Summary

Please use this section to discuss strengths/weaknesses of study, novelty/significance, general execution and scholarship.

Reviewer #1: Authors have addressed the concerns adequately

Reviewer #2: The authors have done a fair amount of additional analysis and provided sufficient explanations to justify their results and conclusions. I agree that this study can be published in its current form and don't have additional questions or issues.

**********

Part II – Major Issues: Key Experiments Required for Acceptance

Please use this section to detail the key new experiments or modifications of existing experiments that should be absolutely required to validate study conclusions.

Generally, there should be no more than 3 such required experiments or major modifications for a "Major Revision" recommendation. If more than 3 experiments are necessary to validate the study conclusions, then you are encouraged to recommend "Reject".

Reviewer #1: N/A

Reviewer #2: See above

**********

Part III – Minor Issues: Editorial and Data Presentation Modifications

Please use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity.

Reviewer #1: N/A

Reviewer #2: See above

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

10.1371/journal.ppat.1012454.r004
Acceptance letter
Kalejta Robert F. Section Editor
Lambert Paul F. Academic Editor
© 2024 Kalejta, Lambert
2024
Kalejta, Lambert
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
14 Aug 2024

Dear Dr. Laimins,

We are delighted to inform you that your manuscript, "HPV induced R-loop formation represses innate immune gene expression while activating DNA damage repair pathways," has been formally accepted for publication in PLOS Pathogens.

We have now passed your article onto the PLOS Production Department who will complete the rest of the pre-publication process. All authors will receive a confirmation email upon publication.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any scientific or type-setting errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript. Note: Proofs for Front Matter articles (Pearls, Reviews, Opinions, etc...) are generated on a different schedule and may not be made available as quickly.

Soon after your final files are uploaded, the early version of your manuscript, if you opted to have an early version of your article, will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Pathogens.

Best regards,

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064
==== Refs
References

1 Aguilera A , García-Muse T . R Loops: From Transcription Byproducts to Threats to Genome Stability. Molecular cell. 2012;46 (2 ):115–24. doi: 10.1016/j.molcel.2012.04.009 22541554
2 Sanz Lionel A , Hartono Stella R , Lim Yoong W , Steyaert S , Rajpurkar A , Ginno Paul A , et al . Prevalent, Dynamic, and Conserved R-Loop Structures Associate with Specific Epigenomic Signatures in Mammals. Molecular cell. 2016;63 (1 ):167–78. doi: 10.1016/j.molcel.2016.05.032 27373332
3 Yu K , Chedin F , Hsieh CL , Wilson TE , Lieber MR . R-loops at immunoglobulin class switch regions in the chromosomes of stimulated B cells. Nat Immunol. 2003;4 (5 ):442–51. doi: 10.1038/ni919 12679812
4 Petermann E , Lan L , Zou L . Sources, resolution and physiological relevance of R-loops and RNA—DNA hybrids. Nature Reviews Molecular Cell Biology. 2022;23 (8 ):521–40. doi: 10.1038/s41580-022-00474-x 35459910
5 Crossley MP , Bocek M , Cimprich KA . R-Loops as Cellular Regulators and Genomic Threats. Molecular cell. 2019;73 (3 ):398–411. doi: 10.1016/j.molcel.2019.01.024 30735654
6 Allison DF , Wang GG . R-loops: formation, function, and relevance to cell stress. Cell stress. 2019;3 (2 ):38–46. doi: 10.15698/cst2019.02.175 31225499
7 Chedin F , Benham CJ . Emerging roles for R-loop structures in the management of topological stress. The Journal of biological chemistry. 2020;295 (14 ):4684–95. doi: 10.1074/jbc.REV119.006364 32107311
8 Choo JAMY , Schlösser D , Manzini V , Magerhans A , Dobbelstein M . The integrated stress response induces R-loops and hinders replication fork progression. Cell Death & Disease. 2020;11 (7 ):538.32678076
9 Edwards DS , Maganti R , Tanksley JP , Luo J , Park JJH , Balkanska-Sinclair E , et al . BRD4 Prevents R-Loop Formation and Transcription-Replication Conflicts by Ensuring Efficient Transcription Elongation. Cell Rep. 2020;32 (12 ):108166. doi: 10.1016/j.celrep.2020.108166 32966794
10 Lam FC , Kong YW , Huang Q , Vu Han T-L , Maffa AD , Kasper EM , et al . BRD4 prevents the accumulation of R-loops and protects against transcription—replication collision events and DNA damage. Nature communications. 2020;11 (1 ):4083. doi: 10.1038/s41467-020-17503-y 32796829
11 Lee C-Y , McNerney C , Ma K , Zhao W , Wang A , Myong S . R-loop induced G-quadruplex in non-template promotes transcription by successive R-loop formation. Nature communications. 2020;11 (1 ):3392. doi: 10.1038/s41467-020-17176-7 32636376
12 Promonet A , Padioleau I , Liu Y , Sanz L , Biernacka A , Schmitz A-L , et al . Topoisomerase 1 prevents replication stress at R-loop-enriched transcription termination sites. Nature communications. 2020;11 (1 ):3940. doi: 10.1038/s41467-020-17858-2 32769985
13 Yoon J , Hwang Y , Yun H , Chung JM , Kim S , Kim G , et al . LC3B drives transcription-associated homologous recombination via direct interaction with R-loops. Nucleic acids research. 2024. doi: 10.1093/nar/gkae156 38412240
14 Chakraborty P , Huang JTJ , Hiom K . DHX9 helicase promotes R-loop formation in cells with impaired RNA splicing. Nature communications. 2018;9 (1 ):4346. doi: 10.1038/s41467-018-06677-1 30341290
15 Jaiswal AS , Dutta A , Srinivasan G , Yuan Y , Zhou D , Shaheen M , et al . TATDN2 resolution of R-loops is required for survival of BRCA1-mutant cancer cells. Nucleic acids research. 2023;51 (22 ):12224–41. doi: 10.1093/nar/gkad952 37953292
16 Prendergast L , McClurg UL , Hristova R , Berlinguer-Palmini R , Greener S , Veitch K , et al . Resolution of R-loops by INO80 promotes DNA replication and maintains cancer cell proliferation and viability. Nature communications. 2020;11 (1 ):4534. doi: 10.1038/s41467-020-18306-x 32913330
17 Stork CT , Bocek M , Crossley MP , Sollier J , Sanz LA , Chédin F , et al . Co-transcriptional R-loops are the main cause of estrogen-induced DNA damage. eLife. 2016;5 :e17548. doi: 10.7554/eLife.17548 27552054
18 Kotsantis P , Silva LM , Irmscher S , Jones RM , Folkes L , Gromak N , et al . Increased global transcription activity as a mechanism of replication stress in cancer. Nature communications. 2016;7 (1 ):13087. doi: 10.1038/ncomms13087 27725641
19 Tan SLW , Chadha S , Liu Y , Gabasova E , Perera D , Ahmed K , et al . A Class of Environmental and Endogenous Toxins Induces BRCA2 Haploinsufficiency and Genome Instability. Cell. 2017;169 (6 ):1105–18.e15. doi: 10.1016/j.cell.2017.05.010 28575672
20 Templeton CW , Laimins LA . p53-dependent R-loop formation and HPV pathogenesis. Proceedings of the National Academy of Sciences of the United States of America. 2023;120 (35 ):e2305907120. doi: 10.1073/pnas.2305907120 37611058
21 de Martel C , Plummer M , Vignat J , Franceschi S . Worldwide burden of cancer attributable to HPV by site, country and HPV type. Int J Cancer. 2017;141 (4 ):664–70. doi: 10.1002/ijc.30716 28369882
22 Kahn JA , Brown DR , Ding L , Widdice LE , Shew ML , Glynn S , et al . Vaccine-type human papillomavirus and evidence of herd protection after vaccine introduction. Pediatrics. 2012;130 (2 ):e249–56. doi: 10.1542/peds.2011-3587 22778297
23 Vu M , Yu J , Awolude OA , Chuang L . Cervical cancer worldwide. Current problems in cancer. 2018;42 (5 ):457–65. doi: 10.1016/j.currproblcancer.2018.06.003 30064936
24 Castle PE , Murokora D , Perez C , Alvarez M , Quek SC , Campbell C . Treatment of cervical intraepithelial lesions. International Journal of Gynecology & Obstetrics. 2017;138 :20–5. doi: 10.1002/ijgo.12191 28691333
25 Bruno MT , Cassaro N , Bica F , Boemi S . Progression of CIN1/LSIL HPV Persistent of the Cervix: Actual Progression or CIN3 Coexistence. Infect Dis Obstet Gynecol. 2021;2021 :6627531. doi: 10.1155/2021/6627531 33776406
26 Nedjai B , Reuter C , Ahmad A , Banwait R , Warman R , Carton J , et al . Molecular progression to cervical precancer, epigenetic switch or sequential model? Int J Cancer. 2018;143 (7 ):1720–30. doi: 10.1002/ijc.31549 29679470
27 Moody C . Mechanisms by which HPV Induces a Replication Competent Environment in Differentiating Keratinocytes. Viruses. 2017;9 (9 ). doi: 10.3390/v9090261 28925973
28 Moody CA , Laimins LA . Human Papillomaviruses Activate the ATM DNA Damage Pathway for Viral Genome Amplification upon Differentiation. PLOS Pathogens. 2009;5 (10 ):e1000605. doi: 10.1371/journal.ppat.1000605 19798429
29 Pyeon D , Pearce SM , Lank SM , Ahlquist P , Lambert PF . Establishment of Human Papillomavirus Infection Requires Cell Cycle Progression. PLOS Pathogens. 2009;5 (2 ):e1000318. doi: 10.1371/journal.ppat.1000318 19247434
30 De Geest K , Turyk ME , Hosken MI , Hudson JB , Laimins LA , Wilbanks GD . Growth and differentiation of human papillomavirus type 31b positive human cervical cell lines. Gynecol Oncol. 1993;49 (3 ):303–10. doi: 10.1006/gyno.1993.1131 8390960
31 Fehrmann F , Laimins LA . Human papillomavirus type 31 life cycle: methods for study using tissue culture models. Methods in molecular biology (Clifton, NJ). 2005;292 :317–30. doi: 10.1385/1-59259-848-x:317 15507718
32 Mehta K , Laimins L , Imperiale MJ , Munger K , Androphy E . Human Papillomaviruses Preferentially Recruit DNA Repair Factors to Viral Genomes for Rapid Repair and Amplification. mBio. 2018;9 (1 ):e00064–18. doi: 10.1128/mBio.00064-18 29440569
33 Crane H , Carr I , Hunter KD , El-Khamisy SF . Senataxin modulates resistance to cisplatin through an R-loop mediated mechanism in HPV-associated Head and Neck Squamous Cell Carcinoma. bioRxiv. 2024:2024.02.22.581374.
34 Ozbun MA , Patterson NA . Using organotypic (raft) epithelial tissue cultures for the biosynthesis and isolation of infectious human papillomaviruses. Curr Protoc Microbiol. 2014;34 :14b.3.1–8. doi: 10.1002/9780471729259.mc14b03s34 25082004
35 Loomis EW , Sanz LA , Chédin F , Hagerman PJ . Transcription-associated R-loop formation across the human FMR1 CGG-repeat region. PLoS Genet. 2014;10 (4 ):e1004294. doi: 10.1371/journal.pgen.1004294 24743386
36 Jurga M , Abugable AA , Goldman ASH , El-Khamisy SF . USP11 controls R-loops by regulating senataxin proteostasis. Nature communications. 2021;12 (1 ):5156. doi: 10.1038/s41467-021-25459-w 34526504
37 Li L , Matsui M , Corey DR . Activating frataxin expression by repeat-targeted nucleic acids. Nature communications. 2016;7 (1 ):10606. doi: 10.1038/ncomms10606 26842135
38 Ginno PA , Lott PL , Christensen HC , Korf I , Chédin F . R-loop formation is a distinctive characteristic of unmethylated human CpG island promoters. Molecular cell. 2012;45 (6 ):814–25. doi: 10.1016/j.molcel.2012.01.017 22387027
39 Ge SX , Jung D , Yao R . ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2019;36 (8 ):2628–9.
40 Bedard MC , Chihanga T , Carlile A , Jackson R , Brusadelli MG , Lee D , et al . Single cell transcriptomic analysis of HPV16-infected epithelium identifies a keratinocyte subpopulation implicated in cancer. Nature communications. 2023;14 (1 ):1975. doi: 10.1038/s41467-023-37377-0 37031202
41 Bienkowska-Haba M , Luszczek W , Zwolinska K , Scott RS , Sapp M . Genome-Wide Transcriptome Analysis of Human Papillomavirus 16-Infected Primary Keratinocytes Reveals Subtle Perturbations Mostly due to E7 Protein Expression. J Virol. 2020;94 (3 ). doi: 10.1128/JVI.01360-19 31748387
42 Bubeck D , Reijns MAM , Graham SC , Astell KR , Jones EY , Jackson AP . PCNA directs type 2 RNase H activity on DNA replication and repair substrates. Nucleic acids research. 2011;39 (9 ):3652–66. doi: 10.1093/nar/gkq980 21245041
43 Cerritelli SM , Sakhuja K , Crouch RJ . RNase H1, the Gold Standard for R-Loop Detection. Methods in molecular biology (Clifton, NJ). 2022;2528 :91–114. doi: 10.1007/978-1-0716-2477-7_7 35704187
44 Chédin F . Nascent Connections: R-Loops and Chromatin Patterning. Trends Genet. 2016;32 (12 ):828–38. doi: 10.1016/j.tig.2016.10.002 27793359
45 Jayakumar S , Patel M , Boulet F , Aziz H , Brooke GN , Tummala H , et al . PSIP1/LEDGF reduces R-loops at transcription sites to maintain genome integrity. Nature communications. 2024;15 (1 ):361. doi: 10.1038/s41467-023-44544-w 38191578
46 Scalera C , Ticli G , Dutto I , Cazzalini O , Stivala LA , Prosperi E . Transcriptional Stress Induces Chromatin Relocation of the Nucleotide Excision Repair Factor XPG. International journal of molecular sciences. 2021;22 (12 ). doi: 10.3390/ijms22126589 34205418
47 Becker JS , Nicetto D , Zaret KS . H3K9me3-dependent heterochromatin: barrier to cell fate changes. Trends in Genetics. 2016;32 (1 ):29–41. doi: 10.1016/j.tig.2015.11.001 26675384
48 Bulut-Karslioglu A , Inti A , Ramirez F , Barenboim M , Onishi-Seebacher M , Arand J , et al . Suv39h-dependent H3K9me3 marks intact retrotransposons and silences LINE elements in mouse embryonic stem cells. Molecular cell. 2014;55 (2 ):277–90. doi: 10.1016/j.molcel.2014.05.029 24981170
49 Becker JS , McCarthy RL , Sidoli S , Donahue G , Kaeding KE , He Z , et al . Genomic and proteomic resolution of heterochromatin and its restriction of alternate fate genes. Molecular cell. 2017;68 (6 ):1023–37. e15. doi: 10.1016/j.molcel.2017.11.030 29272703
50 Nicetto D , Zaret KS . Role of H3K9me3 heterochromatin in cell identity establishment and maintenance. Curr Opin Genet Dev. 2019;55 :1–10. doi: 10.1016/j.gde.2019.04.013 31103921
51 Gillespie KA , Mehta KP , Laimins LA , Moody CA . Human papillomaviruses recruit cellular DNA repair and homologous recombination factors to viral replication centers. J Virol. 2012;86 (17 ):9520–6. doi: 10.1128/JVI.00247-12 22740399
52 Cui X-F , Imaizumi T , Yoshida H , Borden EC , Satoh K . Retinoic acid-inducible gene-I is induced by interferon-γ and regulates the expression of interferon-γ stimulated gene 15 in MCF-7 cells. Biochemistry and cell biology. 2004;82 (3 ):401–5.15181474
53 Imaizumi T , Hatakeyama M , Yamashita K , Yoshida H , Ishikawa A , Taima K , et al . Interferon-γ induces retinoic acid—inducible gene-I in endothelial cells. Taylor & Francis; 2004. p. 169–73.
54 Imaizumi T , Yagihashi N , Hatakeyama M , Yamashita K , Ishikawa A , Taima K , et al . Expression of retinoic acid-inducible gene-I in vascular smooth muscle cells stimulated with interferon-gamma. Life Sci. 2004;75 (10 ):1171–80. doi: 10.1016/j.lfs.2004.01.030 15219805
55 Imaizumi T , Yagihashi N , Hatakeyama M , Yamashita K , Ishikawa A , Taima K , et al . Upregulation of retinoic acid-inducible gene-I in T24 urinary bladder carcinoma cells stimulated with interferon-gamma. Tohoku J Exp Med. 2004;203 (4 ):313–8. doi: 10.1620/tjem.203.313 15297736
56 Martín-Vicente M , Medrano LM , Resino S , García-Sastre A , Martínez I . TRIM25 in the Regulation of the Antiviral Innate Immunity. Front Immunol. 2017;8 :1187. doi: 10.3389/fimmu.2017.01187 29018447
57 Yang C , Shu J , Miao Y , Liu X , Zheng T , Hou R , et al . TRIM25 negatively regulates IKKε-mediated interferon signaling in black carp. Fish & Shellfish Immunology. 2023;142 :109095.37730077
58 García-Pichardo D , Cañas JC , García-Rubio ML , Gómez-González B , Rondón AG , Aguilera A . Histone Mutants Separate R Loop Formation from Genome Instability Induction. Molecular cell. 2017;66 (5 ):597–609.e5. doi: 10.1016/j.molcel.2017.05.014 28575656
59 Huang C , Zhu B . Roles of H3K36-specific histone methyltransferases in transcription: antagonizing silencing and safeguarding transcription fidelity. Biophys Rep. 2018;4 (4 ):170–7. doi: 10.1007/s41048-018-0063-1 30310854
60 Krogan NJ , Kim M , Tong A , Golshani A , Cagney G , Canadien V , et al . Methylation of Histone H3 by Set2 in Saccharomyces cerevisiae Is Linked to Transcriptional Elongation by RNA Polymerase II. Molecular and Cellular Biology. 2003;23 (12 ):4207–18. doi: 10.1128/MCB.23.12.4207-4218.2003 12773564
61 Li J , Moazed D , Gygi SP . Association of the Histone Methyltransferase Set2 with RNA Polymerase II Plays a Role in Transcription Elongation*. Journal of Biological Chemistry. 2002;277 (51 ):49383–8. doi: 10.1074/jbc.M209294200 12381723
62 Neri F , Rapelli S , Krepelova A , Incarnato D , Parlato C , Basile G , et al . Intragenic DNA methylation prevents spurious transcription initiation. Nature. 2017;543 (7643 ):72–7. doi: 10.1038/nature21373 28225755
63 Gautam D , Johnson BA , Mac M , Moody CA . SETD2-dependent H3K36me3 plays a critical role in epigenetic regulation of the HPV31 life cycle. PLoS Pathog. 2018;14 (10 ):e1007367. doi: 10.1371/journal.ppat.1007367 30312361
64 Mac M , DeVico BM , Raspanti SM , Moody CA . The SETD2 Methyltransferase Supports Productive HPV31 Replication through the LEDGF/CtIP/Rad51 Pathway. J Virol. 2023;97 (5 ):e0020123. doi: 10.1128/jvi.00201-23 37154769
65 Lam UTF , Tan BKY , Poh JJX , Chen ES . Structural and functional specificity of H3K36 methylation. Epigenetics & Chromatin. 2022;15 (1 ):17. doi: 10.1186/s13072-022-00446-7 35581654
66 Sun Z , Zhang Y , Jia J , Fang Y , Tang Y , Wu H , et al . H3K36me3, message from chromatin to DNA damage repair. Cell & Bioscience. 2020;10 (1 ):9. doi: 10.1186/s13578-020-0374-z 32021684
67 Li F , Mao G , Tong D , Huang J , Gu L , Yang W , et al . The histone mark H3K36me3 regulates human DNA mismatch repair through its interaction with MutSα. Cell. 2013;153 (3 ):590–600.23622243
68 Sharda A , Humphrey TC . The role of histone H3K36me3 writers, readers and erasers in maintaining genome stability. DNA Repair. 2022;119 :103407. doi: 10.1016/j.dnarep.2022.103407 36155242
69 Panatta E , Butera A , Mammarella E , Pitolli C , Mauriello A , Leist M , et al . Metabolic regulation by p53 prevents R-loop-associated genomic instability. Cell Reports. 2022;41 (5 ):111568. doi: 10.1016/j.celrep.2022.111568 36323249
70 Marnef A , Legube G . R-loops as Janus-faced modulators of DNA repair. Nature Cell Biology. 2021;23 (4 ):305–13. doi: 10.1038/s41556-021-00663-4 33837288
71 Alhmoud JF , Woolley JF , Al Moustafa AE , Malki MI . DNA Damage/Repair Management in Cancers. Cancers (Basel). 2020;12 (4 ). doi: 10.3390/cancers12041050 32340362
72 Singh I , Ozturk N , Cordero J , Mehta A , Hasan D , Cosentino C , et al . High mobility group protein-mediated transcription requires DNA damage marker γ-H2AX. Cell Research. 2015;25 (7 ):837–50.26045162
73 Dobersch S , Rubio K , Singh I , Günther S , Graumann J , Cordero J , et al . Positioning of nucleosomes containing γ-H2AX precedes active DNA demethylation and transcription initiation. Nature communications. 2021;12 (1 ):1072.
74 Ribeiro de Almeida C , Dhir S , Dhir A , Moghaddam AE , Sattentau Q , Meinhart A , et al . RNA Helicase DDX1 Converts RNA G-Quadruplex Structures into R-Loops to Promote IgH Class Switch Recombination. Molecular cell. 2018;70 (4 ):650–62.e8. doi: 10.1016/j.molcel.2018.04.001 29731414
75 Wulfridge P , Yan Q , Rell N , Doherty J , Jacobson S , Offley S , et al . G-quadruplexes associated with R-loops promote CTCF binding. Molecular cell. 2023;83 (17 ):3064–79.e5. doi: 10.1016/j.molcel.2023.07.009 37552993
76 Skourti-Stathaki K , Torlai Triglia E , Warburton M , Voigt P , Bird A , Pombo A . R-Loops Enhance Polycomb Repression at a Subset of Developmental Regulator Genes. Molecular cell. 2019;73 (5 ):930–45.e4. doi: 10.1016/j.molcel.2018.12.016 30709709
77 Community TG . The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2022 update. Nucleic acids research. 2022;50 (W1 ):W345–W51. doi: 10.1093/nar/gkac247 35446428
