extfieldsofstudy: Medicine Biology

text:
DDX17 is required for efficient DSB repair at DNA:RNA hybrid deficient loci

Abstract Proteins with RNA-binding activity are increasingly being implicated in DNA damage responses (DDR). Additionally, DNA:RNA-hybrids are rapidly generated around DNA double-strand breaks (DSBs), and are essential for effective repair. Here, using a meta-analysis of proteomic data, we identify novel DNA repair proteins and characterise a novel role for DDX17 in DNA repair. We found DDX17 to be required for both cell survival and DNA repair in response to numerous agents that induce DSBs. Analysis of DSB repair factor recruitment to damage sites suggested a role for DDX17 early in the DSB ubiquitin cascade. Genome-wide mapping of R-loops revealed that while DDX17 promotes the formation of DNA:RNA-hybrids around DSB sites, this role is specific to loci that have low levels of pre-existing hybrids. We propose that DDX17 facilitates DSB repair at loci that are inefficient at forming DNA:RNA-hybrids by catalysing the formation of DSB-induced hybrids, thereby allowing propagation of the damage response.

INTRODUCTION
Accurate repair of DNA damage prevents the mutations that are the driving force behind carcinogenesis. To maintain fidelity our cells have evolved a network of repair pathways that detect and resolve the various forms of DNA damage including: nucleotide adducts, inter-strand crosslinks, collapsed replication forks and strand breaks (1). Double-strand breaks (DSBs) are considered the most toxic form of DNA damage due to their relatively high propensity for causing mutations (2) and, mutations in multiple DSB repair genes are associated with an increased cancer risk (3,4).
Initial signalling in response to DSBs is carried out by the three apical DDR kinases: ataxia telangiectasia mutated (ATM), ataxia telangiectasia and rad3-related (ATR) and DNA-dependent protein kinase (DNA-PK) (5). An early and essentially universal marker of DSBs is the phosphorylation of Serine 129 of histone H2AX to form γ-H2AX which directly recruits phosphorylated MDC1 to DNA breaks. The E3-ubiquitin ligase RNF8 then binds to MDC1 and mono-ubiquitylates other factors, including L3MBTL2 (6) and possibly the H1 linker histone, thereby creating a marker that recruits the E3 ligase RNF168 (7,8). RNF168 is currently thought to poly-ubiquitylate H1, contributing to decompaction of chromatin, and to mono-ubiquitylate H2A and H2AX at two positions; K13 and K15 (8,9). These H2A ubiquitylations facilitate the recruitment of 53BP1 and BRCA1 that determine pathway choice, either non-homologous end-joining (NHEJ) or homologous recombination (HR), respectively.
HR is a high-fidelity replication-based process, whereas NHEJ involves rapid processing and ligation of the broken DNA ends that can be error prone. There is growing evidence that the DDR signalling cascade is more complex as a number of additional factors have been implicated, including additional ubiquitin ligases, such as HUWE1 (10), as well as RNA binding proteins, such as DROSHA and WRAP53β (11,12). Chromatin is not only a barrier to DNA repair but also provides the substrate for post-translational modifications and chromatin remodelling enzymes required for ordered recruitment of both repair factors and the chromatin remodelling enzymes essential for repair of and recovery from DNA damage.
There is currently growing evidence for a role for RNA and RNA binding proteins in the DDR and specifically in DSB repair (13)(14)(15). A wide variety of RNA binding proteins have already been characterised as DNA repair genes including the splicing factor THRAP3 (16) and the DEAD-box helicase Senataxin (17) . In addition, both on-going transcription and RNA itself have been implicated in the damage response, with transcription inhibitors (18)(19)(20) and RNase treatments (18,21) both being shown to significantly impair the damage response. While DNA:RNA-hybrids have been detected around active DSBs and shown to be required for efficient repair (12,17,19,22), more recently, RNA polymerase III has been shown to be recruited to DSBs where it catalyses transcription templated from the 3' ssDNA overhang produced upon resection (15). In addition, many RNA binding proteins have been identified that promote their formation and resolution around break sites (12,17,23,24). Defining the roles of RNA and RNA binding proteins in the DDR is an emerging area requiring characterisation.
Here, we use a meta-analysis of proteomic screens to identify novel DNA repair proteins and highlight DDX17 as a potential repair factor. Our investigations provide evidence that DDX17 is important for maintaining genome stability via DSB repair, specifically via facilitating RNF168 recruitment and subsequent histone ubiquitylation at DSBs. Consistent with these observations, cells deficient for DDX17 function display impaired DSB repair in both NHEJ and HR pathways resulting in DNA damage accumulating to toxic levels. Additionally, we characterised a requirement for DDX17 in promoting DSB-induced DNA:RNA-hybrid formation. Interestingly, this role is preferentially required around DNA breaks in loci that are naturally deficient for DNA:RNA-hybrids. Our data suggest that DDX17 specifically facilitates DSB repair by promoting DSB-induced hybrid formation at regions of the genome which are inefficient at forming hybrids.

Transfection and Drug treatments
For DIvA, A549 and BJ-5ta cells, Dharmafect 1 was used at a final dilution of 1/1000 and siRNA targeting either DDX17 (Thermo, s2062 for A549 cells, s20623 for all other experiments), RAD54 () or an untargeted control (Horizon Discovery, D-001810-01-05) were used at a final dilution of 20nM.
Cells were cultured for 24h prior to transfection and treated 48h after transfection. siRNA and Dharmafect 1 (Horizon Discovery, T-2001-03) were first diluted separately in optimem to a volume of 5% the final desired medium volume, and incubated for 5min at RT. siRNA and Dharmafect 1 were then combined in 1:1 ratio and incubated for 20 min at RT. The remaining 90% of the final volume of medium was then added and the medium on the cells was immediately replaced with this. For siRNA transfections in hTERT-RPE1, U2OS and HR-U2OS, cells were transfected with Oligofectamine Reagent (LifeTechnologies) according to manufacturer's instructions. Briefly, 1.2-1.5x10 5 cells were plated on a 35mm cell culture dish. After 24h cells were transfected with 40pmol of negative control siRNA (Dharmacon) or siRNA targeting the gene of interest. Cells were treated and harvested 48h post transfection. ATMi (KU55933, SelleckChem) and DNA-PKi (NU7026) were used at 10uM unless otherwise stated. Olaparib (AZD2281, SelleckChem) and ICRF-193 (I4659, Sigma) were used at indicated doses. Irradiation was performed using a Mainance Millenium Sample Irradiator containing a Cs-137 sealed source.

Western blotting
Cells were harvested by lysing in 1.2x sample loading buffer and sonicated with a Diagenode Biorupter for 5 minutes on high to shear the DNA. Samples were run on either 6% or 10% polyacrylamide gels and transferred in tris-glycine transfer buffer containing 20% methanol and 0.01% SDS onto nitrocellulose membranes. Membranes were blocked using 5% BSA in TBST and primary antibody probing was done overnight at 4°C with 5% BSA in TBST. Antibodies and dilutions are listed in Supplementary Table 1. Membranes were then washed for 10 min in TBST 3 times at RT and incubated with secondary antibodies (Li-COR Biosciences) at RT for 1 h at a dilution of 1/10,000.
Membranes were then washed for 10 min in TBST 3 times at RT and imaged with a Li-COR Odyssey.

Clonogenic survival assay
For U2OS or hTERT-RPE1 cell lines, cells were trypsinized 48h after siRNA transfection and counted. 500 cells were plated onto a 60mm dish. For IR sensitivity assays, cells were allowed to adhere for 1h prior to irradiation at indicated doses (Mainance Millenium Sample Irradiator containing a Cs-137 sealed source), while for Olaparib and ICRF-193 treatments, cells were directly plated in media containing the indicated dose for the duration of the experiment. Cells were grown at 37°C for 10-14 days until colonies were an average of 1-2mm in diameter. Colonies were stained with DMMB (0.25% Dimethylmethylene blue in 50% Methanol) and counted.

Neutral Comet Assay
Cells were treated with siRNA 48h prior to IR (3Gy, Mainance Millenium Sample Irradiator containing a Cs-137 sealed source) and harvested at indicated timepoints. Neutral comet assays were carried out according to manufacturer's guidelines (Trevigen). Briefly, cells were harvested and combined with LMAgarose (Trevigen) at a final concentration of 1x10 5 cells/ml and loaded onto polylysine slides.
The agarose plugs were allowed to solidify at 4°C for 1h before being immersed in lysis buffer (Trevigen) overnight at 4°C. Slides were then equilibrated in cold Electrophoresis buffer (100mM Tris pH9.0, 300mM NaAz) for 30min prior to electrophoresis for 1h at 24V. The DNA was precipitated for 30min at RT in DNA precipitation buffer (1M NH4Ac in EtOH) and washed with 70%EtOH for a further 30min. Slides were allowed to dry overnight at 37°C prior to staining with SYBR green (Roche, S7563). Images were acquired on a DeltaVision integrated microscope system using the Applied Precision SoftWoRx acquisition software mounted on an IX71 Olympus microscope with a 10x air objective (Imsol). All images were taken as single slices using a CoolSNAP HQ2 ICX-285 CCD camera. Comet analysis was performed using the CometScore software from Tritek Corporation.

Metaphase spreads
Cells were treated with 2Gy IR (Xstrahl RS320) and allowed to recover for 48h. Cells were then captured in metaphase by treating with 200ng/mL colcemid (Sigma, D7385) for 2h and were then harvested by trypsinisation as well as their growth medium. Cells were pelleted and resuspended in 10mL 75mM potassium chloride and allowed to swell by incubating at 37°C for 30 min. 5mL of icecold fixative (75% methanol, 25% acetic acid) was then slowly added. Cells were then twice pelleted, resuspended in 10mL fixative and incubated on ice for 2 min. Cells were pelleted and resuspended in 4mL of fixative, then dropped onto glass slides from a height of 15-20cm using a p200 pipette. Slides were then steamed over a water bath for 10 s and allowed to dry overnight. The slides were then stained by immersing them in water containing 0.1µg/mL DAPI and then washed by re-immersing in water and again allowed to dry overnight before coverslips were mounted with Vectashield anti-fade mounting medium (Vector Laboratories, H-1000). Spreads were imaged using a Carl Zeiss LSM 710 confocal microscope and counted manually in ImageJ.

Immunofluorescence
Cells were cultured on glass coverslips and either treated with IR, for U2OS, A549 and RPE-1 cells, or with 300nM hydroxytamoxifen (Sigma, H6278) for 4 hours, for DIvA cells. Cells were washed once in PBS and then pre-extracted at RT for 3 min in CSK buffer (100mM NaCl, 300mM sucrose, 3mM MgCl2, 10mM PIPES pH 7.0, 50mM NaF, 5mM Sodium orthovanadate, 10mM β-glycerol phosphate and 0.7% Triton X-100). Cells were then washed once in CSK buffer, once in PBS and then fixed using 4% paraformaldehyde in PBS for 20 min at RT. Cells were then washed once in PBS, once in TBST and then blocked with TBST containing 10% goat serum (Merck, G9023). Cells were then washed twice in TBST for 5 min at RT and incubated overnight at 4°C with primary antibodies diluted in TBST containing 1% goat serum. All primary antibodies and dilutions are listed in Supplementary   Table 1. Cells were then washed 4 times with TBST for 5 min at RT before incubating with Alexa-Fluor conjugated secondary antibodies at RT for 1h. Cells were then washed 4 times in TBST for 5 min, dipped in water to remove residual buffer and then mounted on glass slides using Vectashield anti-fade hard-set mounting medium containing DAPI (Vector Laboratories, H-1500). Slides were imaged using a Carl Zeiss LSM 710 confocal microscope and images were analysed using the Fiji distribution of ImageJ(33) via the FindFoci plugin (34). FindFoci settings were calibrated on 4 randomly selected images and these settings were then used to count the foci in all imaged nuclei.
For U2OS and hTERT RPE1 WT cells, cells were transfected 48h prior to treatment. Cells were then fixed with 4%PFA (EMS) for 10min at RT before permeabilisation with 0.25% Triton-X100 in PBS for 2min at RT. After blocking with 1%BSA in PBS for 1h, cells were incubated for 1h with primary antibody at 37°C and subsequent secondary antibody for 1h at 37°C. Slides were mounted using Vectashield mounting media with DAPI (VectorLaboratories, H1200). Images were acquired on a DeltaVision integrated microscope system using the Applied Precision SoftWoRx acquisition software mounted on an IX71 Olympus microscope with a UPLFLN 40x objective (numerical aperture [NA] 1.3) (Imsol). All images were taken as Z-slices (0.5μm thickness) using a CoolSNAP HQ2 ICX-285 CCD camera. Parameters for image acquisition were kept constant throughout each experiment. Images were deconvolved using SoftWoRx conservative deconvolution. Quantification was carried out using FIJI software (33). In brief, after deconvolution, the images were projected using 'sum slices' (FIJI software) to prevent data loss, and pixel clusters above a defined threshold were counted as individual foci. For representative images, the slices were projected using 'max projection'.

DRIP-Seq
5ng of DRIP DNA was library prepped using the Ultra II DNA library prep kit for Illumina (NEB, E7645L) according to the manufacturer protocol using 6 PCR cycles. Libraries were then subjected to paired-end 75 cycle sequencing on an Illumina NextSeq 500 an the raw data is available at Array Express, accession number E-MTAB-10814. Analysis was done using custom Bash and R scripts, all scripts are available at https://github.com/Bushell-lab/drip_seq . Briefly, fastq files were aligned using Bowtie2 (35) and then sorted and indexed using Samtools (36). Paired alignments were converted to bedpe files using Bedtools and then cut to standard bed files and converted back to bam files. Read coverage was calculated using Samtools Depth and then normalised to total library size. These coverage files were then imported into R for further processing, statistical analysis and plotting.

Proteomic meta-analysis
Published data from proteomic studies was accessed via publication supplementary materials and filtered for significant hits mostly using the parameters from the original publication. These datasets were compiled into an SQL database using custom python scripts allowing for analysis of hit frequency and dataset overlaps. For protein interaction networks, interaction strengths were generated via the STRING database(37) and were then modelled in Cytoscape (38) with groups being created via statistical clustering using the ClusterONE plugin (39).
Protein-protein interaction networking of the 119 proteins identified by over one third of the studies in our meta-analysis revealed 2 major protein clusters; DNA repair and transcription/splicing as well as several other proteins that are partly associated with ribosomal RNA biogenesis and chromatin organisation (Fig. 1C). Gene ontology analysis of the 119 enriched genes confirmed this result, as DNA repair and splicing GO terms were among the most significant along with terms associated with DNA replication and chromosome organisation (Supplementary Figure S1C). Ranking all genes by the frequency in which they were identified in these studies again highlighted commonly identified DNA repair factors, such as MDC1 and MSH6. However, some genes which have not been reported to be DNA repair factors were also significantly enriched, such as DDX17 (Fig. 1B). Interestingly, several of these enriched genes have also been shown to interact with DNA:RNA-hybrids (50), and DDX17 was found to be the most significantly enriched hybrid interactor without a previously described role in DNA repair (Fig. 1B) (51,52).
A further analysis of the phosphoproteomics studies in the meta-analysis found that only 33 genes are shared between all three studies (Fig. 1D). The low overlap between these studies is likely due to methodological differences, although those identified in all three studies might be expected to be those most likely to be core DNA repair factors and, as expected, a majority of these 33 proteins are indeed canonical DNA repair factors. However, eight were associated with transcription and splicing processes, DDX17 among them (Fig. 1E). Interestingly, all three studies identified the same DDX17 phospho-peptide which harboured 2 "SQ" motifs that are the target of the DDR kinases ATM and ATR (Supplementary Figure S1D). Collectively, this meta-analysis highlighted not only a strong enrichment of transcription and splicing factors in DNA repair processes but, specifically, highlighted DDX17 as a potentially novel DNA repair factor.

DDX17 is required for cell survival in response to DNA damage via double-strand break repair
To investigate this potential role of DDX17 in the DDR, we initially conducted clonogenic survival assays with RPE-1 cells in response to an increasing dose of irradiation (IR) upon siRNA mediated knockdown of DDX17. In RPE-1 cells, DDX17 knockdown resulted in a significant reduction in cell survival in response to as little as 0.5 Gy IR relative to control cells ( Fig. 2A). This reduction in cell survival was comparable to that observed for treatment with the ATM inhibitor KU55933, suggesting an important role for DDX17 in the DDR. IR sensitivity upon DDX17 depletion was also observed in an independent cell line (Supplementary Figure S2A). These U20S cells were also treated with Olaparib and the topoisomerase II inhibitor ICRF-193, agents specific for HR or NHEJ, respectively ( Fig. 2B-C). In all of these experiments, DDX17 knockdown was found to significantly reduce cell survival in response to all the DSB-inducing agents used demonstrating a broad role for DDX17 in the DDR. To directly assess the persistence and resolution of DSBs, we next used neutral comet assays.
Comet tail moments, the measure of DNA breaks, significantly increased with both control and DDX17 knockdown cells 15 minutes post 3Gy IR. However, 24 hours after irradiation cells in which DDX17 had been depleted had significantly higher comet tail moments than in control cells (Fig.2D-E, Supplementary Figure S2B). This is indicative of a compromised DSB repair and demonstrates that DDX17 depletion significantly reduces DSB repair capacity. To further investigate the direct effect of DDX17 depletion upon maintenance of genome stability, we employed metaphase spreads. DDX17 depleted cells treated with 2 Gy IR and allowed to recover for 2 days (Supplementary Figure S2C) displayed substantial loss of chromosomes (Fig. 2F-G), consistent with sustained, unresolved DNA damage (53,54). Combined with the comet assay results, this provides evidence for the direct involvement of DDX17 in DSB repair and the maintenance of genome stability.
Given these pronounced DNA repair phenotypes upon DDX17 knockdown, we quantified the focal recruitment of phosphorylated H2AX (γ-H2AX), 53BP1 and BRCA1 using immunofluorescence. For these experiments, we used the "Damage-Induced via AsiSI (DIvA)" cell system pioneered by the Legube lab that uses an AsiSI restriction enzyme fused to an oestrogen receptor to induce doublestrand breaks at known genomic loci in response to 4-hydroxytamoxifen (OHT) treatment. DDX17 depletion resulted in a small but significant increase in yH2AX foci per cell (Fig. 2H), again suggesting that DDX17 knockdown results in an accumulation of unresolved DNA damage. However, immunofluorescence of both 53BP1 and BRCA1 showed substantially reduced focal recruitment of both after DDX17 knockdown (Fig. 2I-J), which cannot be explained by reduced levels of either 53BP1 or BRCA1 upon DDX17 depletion (Supplemental Fig. 2D). We further confirmed these results with IR-induced foci (IRIF) quantifying both γ-H2AX and 53BP1 in U2OS cells. DDX17 depletion once more results in a subtle increase in γ-H2AX IRIF but a significant decrease in 53BP1 foci (Supplementary Figure S2E-G). These results indicate that DDX17 contributes to DDR signalling at an early stage in the cascade. To further quantify the role of DDX17 in DSB repair, we conducted repair assays, using the well-established GFP reporter cell lines DR-GFP and EJ5, which quantify the efficiency of HR and NHEJ, respectively. DDX17 depletion significantly reduced the efficiency of both HR and NHEJ repair efficiency, with the NHEJ reporter efficiency being reduced to a level comparable to 53BP1 knockdown (Supplementary Figure S2H-I). Thus, our data have identified a crucial role for DDX17 in promoting the recruitment of downstream DSB repair factors to facilitate the repair process and ultimately maintain genome stability.

Immunofluorescence identifies the ubiquitin cascade as the point of DDX17 action in DSB repair
The recruitment of 53BP1 or BRCA1 to DSBs is a relatively early event in the DDR that follows a cascade of histone modifications and protein recruitment. Therefore, we analysed the recruitment of several components of the cascade to determine DDX17's point of action.
As our previous immunofluorescence experiments found no decrease in γ-H2AX focus formation upon DDX17 depletion in DIvA cells with induced DSBs (Fig. 2H), we next probed for recruitment of the ubiquitin ligase RNF8 which occurs through direct interaction with MDC1. DDX17 depletion did not significantly perturb the formation of RNF8 damage-induced foci (Fig. 3A), suggesting a function downstream of RNF8. Recruitment of the subsequent ubiquitin ligase, RNF168, is dependent upon prior RNF8 recruitment and, interestingly, we observed a significant decrease in RNF168 damageinduced foci upon DDX17 depletion (Fig. 3B). In addition, the abundance of conjugated ubiquitin, detected with a specific antibody, was also similarly reduced upon DDX17 depletion (Fig. 3C). Thus, DDX17 functions downstream of RNF8 to facilitate RNF168 focal recruitment and polyubiquitylation at sites of DNA damage. We used quantitative western blotting of mono-ubiquitylated γ-H2AX, which is dependent upon RNF168 and required for subsequent 53BP1 and BRCA1 recruitment, to support these immunofluorescence data and similarly observed reduced damage-dependent monoubiquitylated γ-H2AX upon DDX17 knockdown (Fig. 3D-E). As there is no change in either RNF8 or RNF168 protein levels (Fig. 3D), the observed defects in recruitment and ubiquitylation result from regulation of protein recruitment to chromatin proximal to DSBs sites, rather than any changes in gene expression. Consistent with RNF168 activity being required for their recruitment, both 53BP1 and BRCA1 display defective focal recruitment to sites of DNA damage (Fig. 3F). We independently validated key results using a different DNA damaging agent and siRNA in A549 cells (Supplemental Fig. 3A-D). Thus, DDX17 is required for the recruitment and activity of RNF168, but not RNF8, to chromatin in the vicinity of DSBs, identifying DDX17 as a crucial factor in the regulation of DSB signalling functioning between these two E3 ubiquitin ligases.

DDX17 promotes double-strand break induced DNA:RNA-hybrids at hybrid deficient genomic loci
Loss of RNF168 recruitment and reduced ubiquitylation at sites of DNA damage is consistent with the phenotype previously observed for the microRNA biogenesis factor DROSHA(12), a known interactor of DDX17 (32,55). This study found DROSHA to be required for the formation of DNA:RNA-hybrids at DSBs, and that these hybrids were critical for effective DSB repair. Given that DDX17 was also previously reported to be a DNA:RNA-hybrid interactor (12) As previously reported (12,17,24,56), DRIP-qPCR found that DNA:RNA-hybrid levels increased around the DSB sites in response to induction, but not around an uncut control locus, and that all observed signal was sensitive to RNase-H pre-treatment, indicating that the enrichment is specifically from DNA:RNA-hybrids (Supplementary Figure S4B). Upon DDX17 knockdown, we found that DSBinduced DNA:RNA-hybrids were dramatically reduced at all 4 loci we quantified by qPCR, suggesting that DDX17 is required for the formation of these structures in response to break formation (Fig. 4A,   Supplementary Figure S4B). To examine this phenotype at higher resolution, we also conducted highthroughput sequencing of our DRIP samples (i.e. DRIP-Seq), allowing us to analyse the formation of these hybrids across all the sites cut by the AsiSI enzyme in DIvA cells. An overall analysis of the 99 frequently cut AsiSI loci confirmed that DDX17 depletion causes a significant decrease in DSBinduced DNA:RNA-hybrids (Fig. 4B, Supplementary Figure S4C). Comparing AsiSI sites previously reported to be prone to either NHEJ or HR repair(57) found a significant decrease in DSB-induced DNA:RNA-hybrids at both groups upon DDX17 knockdown (Supplementary Figure S4D). In addition, comparing DSB sites with high transcriptional activity to those with low transcriptional activity found that DDX17 again significantly reduced DSB-induced DNA:RNA-hybrids in both groups (Supplementary Figure S4E).
A site-by-site analysis showed that DDX17 knockdown only significantly affects a subset of AsiSI DSB sites; approximately half of these sites show little to no change (Fig. 4C). This allowed us to categorise the sites into either DDX17-dependent, where DDX17 knockdown reduces DSB-induced hybrid formation, or DDX17-independent, where DDX17 knockdown has no significant effect on DSBinduced hybrid formation (Fig. 4C). This classification separates the effects of DDX17 depletion, as the DDX17-dependent sites show no DSB-induced hybrids after DDX17 knockdown whereas the DDX17-independent sites still have a strong induction of DNA:RNA-hybrids levels compared to CTRLsi (Supplementary Figure S4F-H). Classification of these groups allows us to elucidate the factors that contribute to DDX17-dependence of DSB-induced hybrids by comparing the genomic features of these two groups. Interestingly, there were eight sites that showed increased DSB-induced hybrids upon DDX17 knockdown (Fig. 4C), though this increase is quite subtle and there are too few of these sites to accurately investigate this mechanism. We have previously reported that DSBinduced DNA:RNA hybrid formation positively correlates with the pre-existing transcriptional activity of the damaged locus (22). We therefore initially compared the relative transcriptional activity of the DDX17-dependent and independent groups and determined that, although the difference is subtle, DDX17-dependent loci tend to have lower transcriptional activity than the independent sites (Fig. 4D).
Since there is a subtle relationship between DDX17-dependency and transcriptional activity, we next chose to investigate features related to transcription. We examined the levels of pre-existing DNA:RNA-hybrid between the groups which showed that DDX17-dependent sites have substantially lower levels of pre-existing DNA:RNA-hybrids than DDX17-independent sites (Fig. 4E) and this differential was proportionally greater than their decreased relative transcriptional activity.
While it is known that transcriptional activity and DNA:RNA-hybrids are strongly linked, they are not directly proportional due to complex regulatory features governing DNA:RNA-hybrids. Recent studies have shown that DNA and RNA-binding proteins, chromatin modifications, chromatin remodelling and processes such as splicing and replication have significant influence on DNA:RNA-hybrid levels (50,(58)(59)(60)(61)(62)(63). In the presence of DDX17 we found no correlation between pre-existing DNA:RNA-hybrids and DSB-induced hybrids (Fig. 4F). However, upon DDX17 depletion we see a dramatic shift resulting in a positive correlation between these features (Fig. 4G). Therefore, in the presence of DDX17, DSBinduced hybrids do not rely on pre-existing hybrid levels, but with loss of DDX17 the generation of DSB-induced hybrids becomes dependent on the innate ability of the locus to form these hybrid structures. Interestingly, a re-analysis of our previous DRIP-seq data conducted with DROSHA knockdown did not show a significant positive correlation between DSB-induced hybrids and preexisting hybrids (Supplementary Figure S4I Figure S4J). The RBXML1 locus shows DSBinduced DNA:RNA-hybrids in both control and DDX17 knockdowns, whereas the CAPZB locus has little to no induction upon DDX17 knockdown (Fig. 4H-I, Supplementary Figure S4K).
We have determined that DDX17 facilitates DSB-induced DNA:RNA-hybrid formation in the vicinity of DSBs and that this effect is most notable at genomic loci with innately low levels of DNA:RNA-hybrids.
In fact, DDX17 depletion leads to DSB-induced hybrid levels being directly proportional to the loci's pre-existing hybrid levels, suggesting that without DDX17 this DSB-dependent induction of hybrids relies on the loci's natural ability to form hybrids. As a result, given that they are inefficient at forming hybrids normally, we suggest that DDX17-dependent loci require enzymatic assistance to form DSBinduced DNA:RNA-hybrids, whereas loci of high pre-existing hybrid levels are naturally more efficient at forming DNA:RNA hybrids and are therefore not as dependent upon DDX17.

DISCUSSION
In this study we used a meta-analysis of 11 published large-scale proteomic datasets investigating DNA repair mechanisms to identify potentially novel DNA repair factors. A number of these publications noted that RNA related gene groups were significantly enriched in their results (40,42,44,64,65), and some specifically noted DDX17 as a top hit (44,66). Our analysis yielded a significant number of genes associated with gene expression and splicing, alongside DNA repair genes, as those most commonly identified. In particular, DDX17 was identified with high confidence, as a potential DNA repair gene. Therefore, we chose to investigate the possible role of DDX17 in the DNA damage response.
Clonogenic survival assays revealed that depletion of DDX17 sensitises cells to IR, Olaparib and ICRF-193 indicating a broad role for DDX17 in the response to agents that induce DSBs. Indeed, neutral comet assays found DDX17 depletion to significantly decrease repair of IR-induced DSBs.
Consistent with defective DSB repair, metaphase spreads revealed that DDX17 is necessary to prevent chromosomal loss post IR. Thus, DDX17 is important for the maintenance of genome stability.
Localisation of certain DSB repair factors or specific posttranslational modification to focal repair structures demonstrated that while H2AX phosphorylation was not affected by DDX17 depletion, recruitment of both 53BP1 and BRCA1 was significantly reduced. Defective recruitment of these regulators of DSB repair pathway choice suggests defective NHEJ and HR pathways of repair, which was confirmed in NHEJ-and HR-specific reporter assays. Consistent with a function upstream of 53BP1 and BRCA1, recruitment of RNF168 was significantly DDX17-dependent while RNF8 was DDX17-independent, thus placing the function of DDX17 in signalling DNA damage between these two E3 ubiquitin ligases. As 53BP1 and BRCA1 compete for recruitment to the ubiquitylated H2A and H2AX generated by RNF168, a role for DDX17 in regulating RNF168 can also be inferred from DDX17-dependent formation of conjugated ubiquitin and mono-ubiquitylated γ-H2AX at sites of DNA damage. Competition between 53BP1 and BRCA1 recruitment to RNF168-dependent ubiquitylation of histones is a key regulatory point in the choice between HR and NHEJ and subject to complex regulation (10)(11)(12)(67)(68)(69). It is likely that this stage of the DDR is subject to even more complex regulation than initially believed due to many fine-tuning controls over downstream pathway choice.
DNA:RNA-hybrids are now known to rapidly form at DSBs and are required for efficient repair (12,17,19). Furthermore, several RNA binding proteins having been reported to regulate the formation of DNA:RNA-hybrids around DSBs (12,17,23,24,70). As DDX17 is known to interact with DNA:RNA-hybrids(50), we examined a potential role for DDX17 in regulating DSB-induced DNA:RNA hybrids by quantifying and mapping these hybrids across the genome. We found that DDX17 promotes the formation of DSB-induced DNA:RNA-hybrids and this role is particularly important at loci that normally have low levels of pre-existing DNA:RNA-hybrids. Upon DDX17 depletion, DSB sites in loci with low pre-existing DNA:RNA-hybrids were unable to form the DSB-induced hybrids necessary to propagate an efficient damage signal (12,20,71). However, the role of DDX17 in stimulating these hybrids at loci normally expressing high levels of DNA:RNA hybrids was less important as such loci could form sufficient DSB-induced hybrids even upon depletion of DDX17. We hypothesise that loci with naturally low DNA:RNA-hybrid levels require enzymatic assistance from DDX17 to efficiently form DSB-induced hybrids as they lack features, such as preferential DNA sequence and chromatin structure (72,73), required for high levels of hybrid formation. Since it has been shown previously that these DSBinduced hybrids are necessary for efficient DSB signalling and repair (12,20,71), we propose that DDX17's role in the DDR is to facilitate the formation of DSB-induced hybrids, a role that is particularly important at loci with normally low levels of DNA:RNA hybrids.
Several canonical DNA repair genes have also been shown to function via RNA and DNA:RNA-hybrid related processes. For example, the HR factor RAD52 was found to facilitate RNA-dependent DNA repair via strand-invasion of RNA into double-stranded DNA to form a DNA:RNA-hybrid (20,74,75). A previous study also reported that the NHEJ complex can associate with DNA:RNA-hybrid at DSBs and that this was necessary for error-free repair (19). Future studies will be required to determine whether DDX17 has functions that interplay with these known DNA repair genes and DNA:RNAhybrids. However, as DDR ubiquitylation of histones and signalling is DDX17-dependent, it is likely that DDX17-dependent regulation of DNA:RNA hybrids contributes to the remodelling of chromatin flanking DSBs, which is known to be essential for DSB signalling (76)(77)(78). This is consistent with the well described ability of DNA:RNA-hybrids to regulate gene expression and chromatin architecture (79,80). A role for DNA:RNA hybrids in chromatin remodelling at DSBs would explain how the loss of hybrid promoting proteins such as DDX17 and DROSHA impact upon the early DDR cascade.
Elucidating the role of DSB-induced DNA:RNA hybrids in signalling and repairing DNA damage will require identification and full characterisation of all the proteins regulating these structures. The identification of DDX17 as a new contributor to efficient DNA double strand break repair via both major repair pathways, NHEJ and HR, advances our knowledge of the emerging roles for RNA in DNA-repair.   Same as (H) but for BRCA1 immunofluorescence. All statistics were done using unpaired, directional Wilcoxon tests, *p<0.05, **p<0.01, ***p<0.001.


The paternal genetic legacy of Hungarian-speaking Rétköz (Hungary) and Váh valley (Slovakia) populations

One hundred and six Rétköz and 48 Váh valley samples were collected from the contact zones of Hungarian-Slovakian territories and were genotyped for Y-chromosomal haplotypes and haplogroups. The results were compared with contemporary and archaic data from published sources. The genetic composition of the Rétköz population from Hungary and the Váh valley population from Slovakia indicates different histories. In the Rétköz population, the paternal lineages that were also found in the Hungarian Conquerors, such as R1a-Z93, N-M46, Q-M242, and R1b-L23, were better preserved. These haplogroups occurred in 10% of the population. The population of the Váh valley, however, is characterized by the complete absence of these haplogroups. Our study did not detect a genetic link between the Váh valley population and the Hungarian Conquerors; the genetic composition of the Váh valley population is similar to that of the surrounding Indo-European populations. The Hungarian Rétköz males shared common haplotypes with ancient Xiongnu, ancient Avar, Caucasian Avar, Abkhazian, Balkarian, and Circassian males within haplogroups R1a-Z93, N1c-M46, and R1b-L23, indicating a common genetic footprint. Another difference between the two studied Hungarian populations can be concluded from the Fst-based MDS plot. The Váh valley, in the western part of the Hungarian-Slovakian contact zone, is genetically closer to the Western Europeans. In contrast, Rétköz is in the eastern part of that zone and therefore closer to the Eastern Europeans.

Introduction
The Carpathian Basin was historically the destination for several nomadic tribes that migrated westwards from Inner and Central Asia towards Europe. The ancient Hungarians (Steppe Magyars) entered the Carpathian Basin from the east in the ninth century CE and settled there (Róna-Tas, 1999). The genetic legacy of the populations in the Carpathian Basin can be studied only if the genetic structure of them is compared to that of ancient populations. Fortunately, in recent years, the number of aDNA genetic results based on both uniparental and whole genome data has increased in Hungary (Neparáczki et al., 2017a;Neparáczki et al., 2017b;Neparáczki et al., 2018;Olasz et al., 2018;Neparáczki et al., 2019;Csáky et al., 2020;Fóthi et al., 2020;Maár et al., 2021;Nagy et al., 2021). Several Y-STR and Y-SNP results from contemporary Hungarian-speaking populations have been published (Völgyi et al., 2009;Bíró et al., 2015;Pamjav et al., 2017). However, fewer genetic results exist from populations living in geographically isolated areas in the Carpathian Basin, such as the Bodrogköz .
Based on published uniparental genetic data, the contemporary Hungarians possess some genetic similarities to ancient conquerors in the Carpathian Basin, as well as to Central/Inner Asian populations and populations in the Ural Mountains (Bíró et al., 2015;Fehér et al., 2015;Neparáczki et al., 2017a;Pamjav et al., 2017;Huang et al., 2018;Dudás et al., 2019;Csáky et al., 2020). The core of the Hungarian Conquerors may have originated from Inner Asia/Southern Siberia, and then likely admixed with the peoples they encountered during their migration westwards. Examples of this may be found in the Ural and the Caucasus Mountains. The Y-chromosomal relationship for Inner Asian origin is the presence of haplogroups C-M86, R1a-Z93, and N-M46 in ancient or present-day Hungarian populations Fóthi et al., 2020). The Y-chromosomal relationship for Ural or Caucasian origin is the presence of Y-haplogroups N-L1034, R1b-L23, Q-M242, and G2a-L156 in Hungarian conquerors or modern Hungarian populations (Biró et al., Pamjav et al., 2017;Fóthi et al., 2020). The mtDNA relationship for Inner Asia/Southern Siberia origin is the presence of mtDNA haplogroups A, B, H6, and T1a* in the Hungarian Conquerors or present-day Hungarians (Neparáczki et al., 2017a). The mtDNA haplogroup X2f represents a genetic link to Caucasian populations (Neparáczki et al., 2017a). Other mtDNA studies of early Hungarians revealed a diverse composition of haplogroups with significant Asian affinity (Bogács-Szabó et al., 2005;Branstatter et al., 2007;Tömöry et al., 2007) and concluded that the influence of Asian migration into Eastern Europe is detectable in the Székely (or Sekler) population (Romania) and the relatively high proportion of Asian mtDNA haplogroups A, B, C, G, and Y, totaling 7.9% altogether, may be an indicate an influx of Asian nomads into the Székely population in medieval Transylvania (Branstatter et al., 2007).
Thus, the present-day Hungarian speakers have a very complex genetic history. Moreover, the population consists of several ethnic groups in the Carpathian Basin, such as the Székely, Csángó, Palóc, Jász, and Kun. However, Y-chromosomal studies about them and other regional population groups are still limited.
In our previous study, we studied paternal genetic composition of the Bodrogköz population in Eastern Hungary and found genetic similarities to that of ancient Hungarians . To continue the study, we analyzed two Hungarian speaking populations, the Rétköz population in Eastern Hungary, as well as the Váh valley population in Slovakia. To explore their population histories, all Rétköz and Váh valley samples were surveyed for 23 Y-STRs and over 40 Y-SNPs. The resulting data were compared to our published Y-chromosomal data from Eurasian populations and contemporary Hungarians, as well as to Eurasian populations from other published studies.
The Rétköz is a region in Eastern Hungary, whereas the Váh valley is in Western Slovakia. The populations of both places live within the Hungarian-Slovakian contact zone and are Hungarian-speakers. The Rétköz is near the Bodrogköz , but genetic studies have not been conducted in this region or in the Váh valley. The common features of these three regions (Bodrogköz, Rétköz and Váh valley) are that ancient Hungarian Conquerors lived there and that many of their cemeteries have been excavated.
In this study, we present new Y-STR and Y-SNP data of two Hungarian-speaking populations in the contact zone of Hungarian-Slovakian territories and compare them to contemporary Eurasian and available aDNA data to gain further knowledge about the genetic history of these populations.

Materials
We collected samples from 106 unrelated males from Rétköz, Hungary and 48 males from the Váh valley, Slovakia ( Figure 1). Informed consent was obtained from all live participants included in the study.

Testing of Y-STR and Y-SNP markers
Genomic DNA was extracted from buccal swabs using the Investigator kit and EZ1 robotic system (Qiagen, Germany), as described in the manufacturer's instructions. The samples were quantified using the Quantifiler Human kit and the ABI 7500 Realtime PCR System (Thermo Fisher Scientific, Waltham, MA, USA).
DNAfrom the Rétköz and Váh valley populations was surveyed for genetic variation using the Promega PowerPlex Y23 kit. Allele sizing and calling were determined with the ABI3500 Genetic Analyzer and GeneMapper ID-X v.1.4 software. To test for Y-SNP markers (Y chromosomal  Single Nucleotide Polymorphism), we performed amplifications of 1-2 ng genomic DNA with Custom TaqMan probes and analyzed the relative fluorescence of the PCR productsin an ABI 7500 Real-time PCR instrument using SDS.1.2.3 software. The SNP markers tested were CT-M168, R1b-Z2103, D-M174, N-M231, N1-LLy22g, N-L708, N-M46, N-L1034, N-VL29, N-Z1936, N-F4205, R1a-Y2633, and N-Y24365. The haplogroups are described in accordance with the generally accepted nomenclature, as is common practice (Karafet et al., 2008 and the ISOGG).
The Y-STR (Y chromosomal Short Tandem Repeat) haplotypes in this study were sent to the YHRD (accession numbers: YA004754 for the Rétköz and YA004755 for the Váh valley populations).

Phylogenetic analysis
To examine the STR variation within the haplogroups, relational networks were constructed using the Network 10.2.0.0 program (Bandelt et al., 1999). Repeats of the DYS389I locus were subtracted from the DYS389II locus, and the DYS385 locus was excluded because the Network program cannot handle the duplicated locus. To put the results into a more extensive geographical context, we included haplotypes of 10 overlapping evolutionarily stabile loci from other Eurasian populations. The rho statistic in the network program was used to estimate the time to the most recent common ancestor (TMRCA) of haplotypes among the compared haplogroups (Bandelt et al., 1999). Evolutionary time estimates were calculated according to Zhivotovsky et al. and STR mutation rate was assumed to be 6.9 × 10 −4 /locus/25 years (Zhivotovsky et al., 2004) as is common practice. STR-based TMRCA estimates were for information only and are not discussed in this paper due to their unreliability and insignificance to the primary purpose of the study. Frontiers in Genetics frontiersin.org 03

Genetic structure
Based on the Y-STR haplotypes, pairwise Rst (stepwise mutation) genetic distances were computed with YHRD. org's online AMOVA, and the MDS plot constructed (release 66) (Willuweit and Roewer 2015). Pairwise Fst genetic distances were calculated based on haplogroup frequencies using Arlequin 3.5 software (Excoffier and Lischer, 2010). MDS plot from the Fst values was constructed with the cmdscale function of R (R Core Team, 2021). Rst is an analogy of Fst based on allele size differences; it is defined as the correlation of allele sizes between tested markers within populations. Fst (fixation index) is more efficient when there are high levels of gene flow, whereas Rst reflects population differentiation better under low gene flow (Balloux and Goudet, 2002). That's why we used both Rst and Fst-based metric MDS analyses for population comparisons.
Haplotype and haplogroup frequencies and their diversity values were calculated using the formula from Nei 1973.
The haplotype and haplogroup diversities of the Rétköz group were 1.00 and 0.901, respectively, whereas these values

Phylogenetic analysis
Based on the Y-STR haplotypes, networks were constructed for haplogroups that may be linked to those of the Hungarian Conquerors, because these haplogroups can be found in ancient Hungarians and originated from Inner/ Central Asia. We also included 10 Y-STR haplotypes from other Eurasian populations and aDNA results from published sources to widen the geographical context. We have included six networks (R1a-Z93, N-Tat, Q-M242, R1b-P25/M343, E1b1-M78, and G2a-L156) that are potentially helpful in uncovering the genetic legacy of the populations being studied. Other networks, which were constructed for haplogroups of European origin (R1a-M458, R1a-Z280, R1b-P312, R1b-U106 and I2a-P37) and could not be related to Hungarian Conquerors, were not included in the study. Supplementary Table S2 summarizes data for all the published aDNA samples used in the study including haplotypes, haplogroups, ages and geographic origins. Figure 2 depicts an MJ network of 121 R1a-Z93 haplotypes from the 15 populations tested by us for this study or previously published (Bíró et al., 2015;Underhill et al., 2015;Dudás et al., 2019) and aDNA samples (Olasz et al., 2018;Fóthi et al., 2020;Keyser et al., 2021). Three modern Hungarian and one aDNA samples (II54 from the Hungarian Royal Basilica of Székesfehérvár), including 1 Rétköz sample (Rétköz 40), formed a common branch with three ancient Xiongnu samples (TUK45, TUK04, TUK25) on the right side of Figure 2. On this branch, one modern Hungarian sample shared two haplotypes (1 TUK25 and one Bashkirian Mari) in haplotype cluster 3. Two Uzbek samples (1 Khwarizm and one Fergana sample) can be derived from cluster 3 (see Figure 2). The five Bashkirian Mari and two Uzbek (Khwarizm) samples are clustered one molecular step from cluster three at the end of this branch. Cluster five includes one Bashkirian Mari and two Uzbek males from Tashkent in Uzbekistan and located at one molecular step from cluster 1. Cluster one includes three haplotypes: one Hungarian aDNA (Nagykörös Gr2), one Xiongnu aDNA (TUK09A) and one modern Armenian sample, separated by one molecular step (DYS389I) from cluster 5. The founding haplotype, which may have arisen in the common ancestor of Frontiers in Genetics frontiersin.org 05 populations or males and is shared by them. Based on this, cluster 1 may be the founding haplotype, as it contains two ancient haplotypes from males that lived 1,000 to 2,500 years ago. All the samples and branches are derived from this haplotype. Cluster two includes one Rétköz Hungarian (Rétköz = R48), one Mongolian, one Altaian, and one Andronovo aDNA (S10) samples. Cluster four can be separated one molecular step from cluster five and includes one Hungarian and two Uzbek males from Tashkent in Uzbekistan. The Hungarian King Bela III is found three molecular steps away from cluster 1 (see B3 in Figure 2). The remaining three Hungarian haplotypes are outliers in the network and are not shared by any of the samples.

Median joining network of 121 R1a-Z93 haplotypes
The other studied samples in the network either form independent clusters, such as Altaians, Khakassians, Khanties, and Uzbek Madjars, or are scattered within the network. Other aDNA samples (S26, ARZT1, ARZT28, and MN376 from Keyser et al., 2021) in the network form a different branch, with almost all Khakassian samples located on the left side of Figure 2. The age of accumulated STR variation (TMRCA = Time to Most Recent Common Ancestor) within the R1a-Z93* lineage for 121 samples is estimated as 13 ± 3 kya (95% CI = 10.0-16.0 kya), considering that cluster one is the founder haplotype, which is higher than that of the SNP-based calculation (4.6 kya, 95% CI: 4.2-5.0 kya) (Adamov et al., 2015; www.yfull.com).

Median-joining network of 179 N-M46 haplotypes
A median-joining network of 179 N-M46 (previously N1c, Karafet et al., 2008) haplotypes was generated using populations we previously studied (Bíró et al., 2015;Fehér et al., 2015), as well as those researched by others (Pimenoff et al., 2008;Ilumäe et al., Figure 3 (haplotype cluster 1). Cluster one is the founding haplotype. Cluster 2, another large cluster, includes 43 haplotypes (32 Buryat; two Mongolian; one Hungarian; two Northern Mansi; and six ancient Avars from Csáky et al., 2020) from five populations, as shown in Figure 3. Cluster two is located one molecular step from cluster 1. The only difference is in the DYS391 locus, with allele 11 in cluster one and allele 10 in cluster 2. As seen in Figure 3, 73% of Buryats belonged to cluster 2, indicating that the Buryats we studied belong to a young and isolated population (Dudás et al., 2019). Two Hungarian males, including Rétköz 19 (R19) sample, derives from cluster two via one Buryat haplotype (see three in Figure 3). One Hungarian aDNA haplotype (Ö52/50, Fóthi et al., 2020) was positioned six mutational steps away from cluster one and formed a haplotype branch with two present-day Hungarian males (see aH in Figure 3). Three modern Hungarian males shared 1-1 haplotypes with Finnish, Northern Mansi, or Bashkirian males (see black arrows in Figure 3).
The other population samples included in the network formed independent clusters, such as Bashkirian Mari, Finns, Khanties, Bashkirians, and Khanties from another research (Ilumäe et al. (2016). They also shared haplotype clusters or were scattered in the network.

Median-joining network of 153 R1b-M343 haplotypes
A median-joining network of 153 R1b-M343/P25 (Eastern R1b-L23) haplotypes was generated using samples from FTDNA, populations we previously studied (Balanovsky et al., 2011;Bíró et al., 2015) (Figure 4). The founder R1b*-M343 (L23) haplotype was shared by nine samples, including the Rétköz 41 (R41) sample, from two populations (7 Hungarians and two Avars from the Caucasus), as shown in Figure 4. The majority Hungarian haplotypes, including four Rétköz (R28, R52, R39 and R37) and one Váh valley samples (V24), appear to be descended from the founding haplotype (F = founder), because these haplotypes differ one molecular step from the founding haplotype. The pattern of these haplotype clusters is starlike, representing a set of closely related haplotypes of Hungarian males. Three Rétköz R1b samples (R54, R59 and R75) clustered with Ossetian samples from the Caucasus (indicated by the red circles in a dashed circle under the founding haplotype cluster in Figure 4.). In addition to the founding haplotype, some Hungarian haplotypes were shared with Abkhazian (turquoise and red circles) or with Hungarians and Avars from the Caucasus (red and black circles), which are marked with an arrow. All other Hungarian haplotypes were scattered within the network. It is interesting that either Hungarian haplotypes shared a common haplotype with Avar males or clustered with Avar haplotypes in the network (red and black circles).
The other population samples in the network formed independent clusters, such as Lezgian, Ossetian, and Circassian from the Caucasus or were scattered in the network. Primarily Balkarian (Caucasus) and Uzbek males form an independent branch (R1b-M73) on the right side of the network (Figure 4). The age of accumulated STR variation within the R1b-M343/P25 lineage for 153 samples is estimated to be 22.8 ± 3.8 kya (95% CI = 19.0-26.6 kya), considering the founder haplotype is the ancestral one (F in Figure 4), which is in agreement with the time calculated on sequence data (yfull.com).

Median-joining network of 167 Q-M242 haplotypes
The median-joining network of 167 Q-M242 haplotypes from 12 populations does not show a star-like pattern, but rather a more diverse one. However, it was simple to analyze the genetic relationship of the Hungarian samples ( Figure 5). The largest haplotype cluster, to which a relatively large number of individuals belong, is shared by four populations (see M on Figure 5). This cluster includes 22 Tuvinian, seven Todjin, seven Altaian, and one Mongolian haplotypes. Another large haplotype cluster (cluster 1) includes 14 males from three populations (2 Sojot, one Tojin, and 11 Tuvinian males), whereas almost all other haplotypes are scattered within the network.
The Rétköz 80 (R80) sample shares the same haplotype as an Uzbek male, and the Rétköz 101 (R101) sample clusters with a Mongolian male. The third Rétköz male (R58) is located six steps away from an Uzbek sample and is an outlier. It is noteworthy that some Hungarian samples, primarily Hungarian-speaking Csángó and Székely males in Transylvania, Romania form a separate branch with the ancient Xiongnu males (TUK01, TUK20, and TUK26) (see the black and red circles within the dashed circle in the lower left part of Figure 5), indicating a close genetic affinity with each other. Two ancient Xiongnu (TUK43 and TUK44) samples and one ancient Avar sample from published sources (Csáky et al., 2020;Keyser et al., 2021) are more closely related to a European (Lithuanian) and an Altai sample (see the top of Figure 5.), which indicate different histories. We constructed another network, including Jewish and Balkarian samples, with 192 haplotypes. Two Hungarian males (Hu423 and E699) formed a cluster with the Jewish samples, suggesting a genetic link (data not shown). These Hungarian samples form a cluster with European samples on the right side of Figure 5 (see red H and pink circles). Three Hungarian samples (Csangó 68, Székely 180, and Hu1769) form a branch with Balkarian males (data not shown). These three Hungarian samples form a branch with other Hungarian and Xiongnu samples in Figure 5. Based on this network, we have distinguished the Hungarian samples that have a genetic relationship with the Jewish and Balkarian samples.
The age of accumulated STR variation within the Q-M242 lineage for 167 samples is estimated to be 17.7 ± 3.6 kya (95% CI = 14.1-21.3 kya), considering that the modal haplotype ( Figure 5) is the founder, which agrees with the previous calculated 15-25 ky (Karafet et al. (2008).

Median-joining network of 189 G2a-L156 haplotypes
The MJ network of 189 G2a-L156 haplotypes is depicted in Figure 6. The samples belong to four subgroups (P303, L497, M406 and P16). The model haplotype cluster (cluster 1 = M in Figure 6) is shared by three populations, including two Hungarian, one German and one Balkarian males. The haplotype cluster is likely to be the founder haplotype, as all branches are derived from this and is the median center of the branches. The biggest haplotype cluster is shared by four populations, including six Hungarians (Rétköz 11 = R11), one German, one Balkarian, and one Circassian chromosome fromthe Caucasus (cluster two in Figure 6). Cluster three is shared by three populations, including one Hungarian, one Avar and two Circassian males from the Caucasus. One Hungarian male matches a German haplotype (red-purple Frontiers in Genetics frontiersin.org 09 circle), and another Hungarian male matches with a Balkarian male (red-dark blue circle in Figure 6), indicating a common origin with Hungarians. We marked these haplotypes with an arrow in Figure 6, as well. All other Hungarian males derive from clusters 1-3, and a shared haplotype was not observed. It appears that the Hungarian haplotypes, including Rétköz (R11, R23, R33, R68 and R70) and Váh valley (V34, V43) samples, are generally clustered together on a common branch with Balkarian/Karachay, Avar, or German haplotypes, suggesting that these chromosomes have common origin. Ossetian, Lezgian, and Abkhazian males from the Caucasus form independent clusters, except for the Circassians, which are scattered within the network.
The age of accumulated STR variation within the G2a1-L156 lineage for 189 samples is estimated to be 19.2 ± 3.7 kya (95% CI = 15.5-22.9 kya), considering the modal haplotype ( Figure 6) is the founder one. This is in line with the 15.0 ky calculated by Rootsi et al. (2012). Figure 7 depicts the median-joining network of 183 E1b1-M78 haplotypes. The network shows a star-like pattern. The biggest cluster (cluster 1 = M) is the modal haplotype shared by six populations, including 43 males. The model haplotype includes seven Hungarian (R32 and R85), 12 Albanian, five Ukrainian, 10 Serbian, seven Bulgarian, and two Slavic haplotypes. There are several haplotype clusters showing different admixture of the seven populations researched in the study, such as cluster 2 (5 Hungarian, one Albanian, one Ukrainian, one Serbian, and three Bulgarian males), cluster 3 (1 Hungarian, two Albanian, one Ukrainian, four Serbian, and three Bulgarian males), cluster 4 (3 Hungarian, one Bulgarian, one Serbian, and one Slavic males), cluster 5 (4 Albanian, one Rétköz Hungarian, and one Bulgarian males), cluster 6 (2 Váh valley Hungarian, two Albanian, and one Serbian males), cluster 7 (1 Hungarian, one Rétköz and one Albanian males), and cluster 8 (1 Hungarian, one Váh valley, and two Albanian males). The remaining Hungarian and other samples are scattered within the network or have shared haplotypes, like Ukrainian-Hungarian or Austrian-Bulgarian (see arrows in Figure 7).

Median-joining network of 183 E1b1-M78 haplotypes
The age of accumulated STR variation within the E1b1-M78 lineage for 183 samples is estimated to be 6.7 ± 1.3 kya (95% CI = 5.4-8.0 kya), considering the modal haplotype (see M in Figure 7) is the founder, which is in line with the arrival of Neolithic farmers in Europe (Battaglia et al., 2009).

Genetic structure
We constructed a non-metric multidimensional scaling (MDS) based on Y-chromosomal haplotypes (2,405 haplotypes) that consisted of 23 STR loci available from 14 populations (www. yhrd.org). The Rst-genetic distances and Rst p-values of the The cluster one populations showed close genetic affinities (<0.05) with most of the European populations that were compared, such as Croatian, Czech, and Polish. Interestingly, cluster one population is closely related to the Xinjiang Uighur population in China. Finnish (0.2844), Bashkirian Mari (0.1722), and Estonian (0.1181) populations, however, are genetically farthest from cluster 1.
The haplogroup frequency data used for population comparisons, as well as their corresponding references are included in Supplementary Table S4. Pairwise Fst-distances and p-values for 61 populations, including Rétköz, Váh valley, and other Eurasian populations from published sources were calculated as shown in Supplementary Table S5 and presented in a metric MDS plot (Figure 9). Between the two studied populations, the pairwise Fst-distance was insignificant (p > 0.05) (Supplementary Table S5). Furthermore, Rétköz and Váh valley had insignificant Fst-values, with 5 and 11 populations, respectively. Among them, three overlapped: two Hungarian (Bodrogköz and Szeged) and a Slovenian population. Czech and Western Slovak have an insignificant Fst-distance from Rétköz, but not from Váh valley. Váh valley showed an insignificant Fst-distance from eight other populations: Bulgarian, Gagauz, Greek, Moldavian, Serbian, and three Hungarian (Csángó, Székely, and a representative Hungarian population).
The location of the studied populations on the MDS plot is consistent with the geographical distances between them (Figure 9). Populations from the same geographic region were clustered together. Hungarian populations, however, overlapped with Balkanian populations. The Rétköz and Váh valley samples were shown to be relatively far from the Hungarian Conqueror samples.

Phylogenetic analysis
The primary objectives of this study were to create a male phylogeny of the Hungarian-speakers in Rétköz and the Váh valley and to compare them with other Eurasian populations, with aDNA samples from their probable ancestral geographical origin, and with previously studied populations.
To examine the genetic variation within the Hungarian groups, we used evolutionarily stable binary markers (SNPs) to define the haplogroup of each Y-chromosome, and then examined the STR-defined variation within each haplogroup. The established geographical specificity of Y haplogroups means that the haplogroups observed in the Hungarian-speakers may be attributed to their geographical origins or may be related to the Hungarian Conquerors. These conclusions, however, remain uncertain.
The overall pattern of haplogroup distributions in the studied populations was similar, but haplogroups R1a-Z93, N1c-L1034, N1c-VL29, N1c-Z1936, and Q-M242 appeared only in the Rétköz population (Table 1). Haplogroups G2a-L156 and R1b-M343/P25 (L23) were observed more frequently in the Rétköz population. Both lineages range predominantly across the Caucasus and Asia. As such we focused on the genetic history of these haplogroups.

Median-joining network of 121 R1a-Z93 haplotypes
Underhill et al. state that the paragroup R1a-Z93 was most common in the South Siberian Altai region but that it also occurred in Kyrgyzstan and, in all Iranian populations. R1a-Z2125 occurred at highest frequencies in Kyrgyzstan and among Afghan Pashtuns. Additionally, the R1a-Z93 haplogroup was also common in Afghan (Tajik) and Caucasian ethnic groups (Tajik). As such, we included the populations (Z93 haplotypes) from the regions analyzed by other researchers (Underhill et al., 2015), ancient Hungarian, Xiongnu, and Avar samples (Olasz et al., 2018;Csáky et al., 2020;Fóthi et al., 2020;Keyser et al., 2021) in the network analysis.  It is surprising that the present-day Hungarian, ancient Hungarian and Xiongnu samples are located adjacent to each other on this branch, indicating that 1-1 mutation step separates them. Keyser et al., 2021 demonstrated that ancient Xiongnu (TUK45, TUK04, TUK25 and TUK09A) samples belonged to haplogroups R1a-Z93 (Z2125), which are also included in our network. These samples clustered with our Hungarian aDNA (II54) and three present-day Hungarian samples, including Rétköz 40 (R40), in the branch mentioned above (see Figure 2). Furthermore, a Hungarian aDNA haplotype (Nagykörös Gr2) shared the same haplotype as the Xiongnu TUK09A sample (cluster one in Figure 2.). Likewise, a modern Hungarian sample shared the same haplotype as the Xiongnu TUK25 sample (cluster three in Figure 2). aDNA studies showed that the Hungarian King Béla III, and another sample (II54) from Royal Basillica belonged to haplogroup R1a-Z93 (Olasz et al., 2018), and two R1a-Z93 samples were found among Hungarian Conqueror population as well (Fóthi et al., 2020). Additionally, some Xiongnu R1a-Z93 haplotypes matched those of the Hungarian Conquerors (Keyser et al., 2021). These observations confirm a genetic relationship between Hungarian Conqueror and the ancient Xiongnu people, suggesting descent from a shared common ancestor. These observations also reflect that the descendants of the Xiongnu people and the ancient Hungarians are found among modern Hungarian populations, indicating the presence of a common genetic trace (i.e., genetic continuity).

Median-joining network of 153 N-M46 haplotypes
The Y-chromosomal haplogroup N-M231 is widespread from Scandinavia to the Kamchatka in North Eurasia and is the most frequent single haplogroup in Siberia (Karafet et al., 2008;Rootsi et al., 2012). Based on the geographic distribution of the most basal N*-M231, its lineage most likely originated in Southeast Asia, whereas its most widespread subgroup is N-M46 (Rootsi et al., 2012). Other authors consider Southern Siberia as its geographical origin (Derenko et al., 2007).
The Hungarian Y-chromosomal gene pool has only a small percentage of N-M46 (1%; Tat = M46) and has a distribution typical of East-Central Europe (Völgyi et al., 2009). However, its incidence is higher among the Hungarian-speaking Bodrogköz (6.2%) and Székely (6.3%) (Transylvania, Romania) populations (Biró et al., Pamjav et al., 2017). Based on Hungarian aDNA studies, N-M46 is detected at a higher frequency in Hungarian Conquerors (Neparaczki et al., Fóthi et al., 2020). These results support genetic continuity between the ancient Hungarian and the present-day Hungarian populations.
The Rétköz N-M46 sample 1 (R1) from this study and four other Hungarian samples belong to haplotype cluster one and share 31 haplotypes with 8 Eurasian populations (Hungarian, Mongolian, Southern Mansi, Avar, Northern Mansi, Finnish, Bashkirian, and Khanty), indicating they may all share a common genetic history (Figure 3). The second Rétköz (R19) N-M46 sample shared a haplotype with the general Hungarian population (see three in Figure 3.). The third Rétköz (R4) haplotype does not share a haplotype and instead clusters with Finnish haplotypes. This is consistent with the N-Z1936 subgroup, as this SNP is common in the Finnish population. Other modern Hungarian N-M46 samples, including Hungarian aDNA (aH: Őrkút Gr50), are located on the top or the right side of the network, along with Bashkirian, Finnish, and samples originating in the Ural Region (Southern and Northern Mansi, Khanty). Three Hungarian males share a haplotype with the Finnish, Northern Mansi, and Bashkirian samples, suggesting a shared ancestor (see the black arrows). One aDNA Xiongnu (X: TUK30A) sample from Mongolia clusters with Mongolian samples on the left side of the network, indicating a genetic relationship between these groups.

Median-joining network of 153 R1b-M343 haplotypes
The most frequent Western European lineage, haplogroup R1b-M269, was originally thought to have originated in the Palaeolithic. Recent analysis, however, suggests a Neolithic origin (Batini et al., 2017). Most R1b-M412 chromosomes belong to Western European males, but another subgroup, R1b-L23, is commonly referred to as "Eastern European R1b". Its frequency among Turkish, Caucasian, and some SE European and Circum-Uralic populations is about 10% (Myres et al., 2011).
Our network analysis showed that the R1b*-M343 (xM412) haplotypes are divided into two subclades: R1b-L23 (xM412) and R1b-M73. The Rétköz and other Hungarian haplotypes cluster with the L23 samples included in the study and share haplotypes with Hungarian, Abkhazian, and Avar males from the Caucasus, suggesting a genetic link between them (see Figure 4). These Avars are unlikely the same as the Pannonian Avars in the Carpathian Basin because present-day examples live in the Northeastern Caucasus, primarily in Dagestan, Kalmykia, and Chechnya in Russia. The earliest mention of the Avars in European history is by Priscus, who reported in 463 C E that a joint delegation of Saragurs, Urogs, and Unogurs requested an alliance with Byzantium. The delegation claimed that in 461 CE, their peoples were displaced by the Sabirs due to pressure from the Avars (Priscus, 463 CE). The German researcher Karl Heinrich Menges stated that these Avars have nothing to do linguistically with the Proto-Mongolian Avars of the Great Migrations, but as some of the latter may have taken refuge in this region, their name has become the common name of the Avars of present-day Dagestan, who have no name of their own (Menges, 2011).
The haplotype analysis revealed a genetic relationship between the contemporary Hungarian-speakers and Caucasian populations (Avars, Abkhazians, and Ossetians) included in this study. It should be noted, however, that the peoples of East-Central Europe received Central-Inner Asian and Caucasian Frontiers in Genetics frontiersin.org 13 genes before and after the ancient Hungarians settled in the Carpathian Basin (pre-Hungarians: Sarmatian-Alans, Huns, Avars, and Onogur-Bulgars; post-Hungarians: Pechenegs, Jassic peoples, and Cumans) (Lipták 1979). Based on aDNA studies, haplogroup R1b-L23 was found in Hungarian Conquerors (Neparáczki et al., 2019;Fóthi et al., 2020), supporting the genetic footprint observed in our results, but it was not possible to determine when and where this genetic trace was introduced into the gene pool of the ancient Hungarians.
Based on the network analysis of haplogroup Q-M242, is that Hungarian-speaking populations living in more isolated areas (no admixture), such as Csángó and Székely males in Transylvania, Romania, clustered on the same branch as the ancient Xiongnu people lived in Inner Asia about 2000 years ago, suggesting these males are genetically related (see Figure 5). An Uzbek sample from Khwarezm (Biró et al., 2015) is also on this branch, in addition to this study's R80 sample, which shares the same haplotype as an Uzbek sample from China (Huang et al., 2018), implying a common genetic trace. Neparáczki et al., 2019 reported that, in Hungary, haplogroups Q-M25 andQ-F1096 (xM25) were found in the Hunnic period (fifth century CE) and the Hungarian Conquest period (9th-11th centuries CE), respectively. They state that although Q-M25 is rare in Europe its highest frequency is among Hungarian speaking Székely population in Transylvania, Romania. Furthermore, the authors noted that ancient samples with haplogroup Q1a2-M25 are known from the Bronze Age Okunevo and Karasuk cultures, as well as from Middle Age Tian Shan Huns and Hunnic-Sarmatians, suggesting that this lineage may be of Hunnic origin in Europe. This is confirmed by the Hun/1 sample, derived from Transylvania (Neparáczki et al., 2019). These observations also confirm the genetic relationship between contemporary Hungarians (Székely and Csángó) and the ancient Xiongnu people of Inner Asia that we found in our study.
The remaining two Rétköz Q samples (R101 and R58) cluster with Central (Uzbek, Tajik, and Turkmen) and Inner Asian (Mongolian) samples in a subbranch of the main branch, which indicates that these samples are distant from the Hungarian-Xiongnu subbranch but that they nevertheless originate from a common branch ( Figure 5).
Another interesting observation is that three Hungarian males on the Hungarian-Xiongnu branch ( Figure 5) form a cluster with Balkarian males in another network (see Results, data not shown). The Balkars identify as a Turkic people and speak the same language as the Karachays from Karachay-Cherkessia. Balkars and Karachay are sometimes referred to as a single ethnicity (Ehaotufca 2010). Our comparative phylogenetic study of folk music and genetics, indicate that the Volga-Sicilian-Turkish-Karachay-Hungarian-Finnish-Dakota and the Chinese-Mongol-Volga-Sicilian -Turkish-Karachay-Hungarian cultures have the largest sets of common melody types, suggesting the existence of a common "parent language" from which their music evolved. The results may also show that the populations that incorporate this musical style into their culture have common genetic roots, in particular that the development of this musical culture in their early history may be attributed to their common genetic ancestors (Pamjav et al., 2012). In this study, the similar melodies in Hungarian and Karachay folk music may suggest their common genetic traces. Unfortunately, however, comprehensive genetic studies on Balkars/Karachays are not available.

Median -joining network of 189 G2a-L156 haplotypes
Haplogroup G isassociated with the spread of agriculture, particularly in Europe. Haplogroup G was first discovered in Europe and Georgia (Semino et al., 2000) and was later detected in Caucasian and Hungarian populations (Nasidze et al., 2004;Völgyi et al., 2009). The frequency of haplogroup G2a-P15 (L156) is about 4% in the general Hungarian population (Völgyi et al., 2009), 6% in the Hungarian-speaking Csángó population (Transylvania, Romania), and 4% in the Hungarian-speaking Székely population (Transylvania, Romania) (Bíró et al., 2015). As the frequency of the G haplogroup is low among Hungarian-speakers, only the L156 SNP was tested from the downstream SNPs, which is phylogenetically equivalent to SNP P287 marker (www. phylotree.org/Y/tree/G.htm).
Rootsi et al. analyzed 113 Hungarian males, of which 2 belonged to haplogroup G-M201 (1.8%) and to subgroups G-L497 (0.9%) and G-M406 (0.9%). These authors showed that SNP P303 defines the most frequent and widespread subhaplogroup G, whereas P303-related L497 lineages occur in Europe, where they likely Frontiers in Genetics frontiersin.org 14 originated. The highest frequency of G2a-P303 is detected in populations from the Caucasus, specifically among South Caucasian Abkhazians (24%), Northwest Caucasian Adyghe (39.7%), and Cherkessians (36.5%) (Rootsi et al., 2012). Another frequent subclade is M406, which is the sister clade of P303. The G2a-M406 has a peak frequency in the Mediterranean and Central Anatolian (6-7%) populations, as well as in Greek (4%) and Italian (3%) populations. It is not detected in many other regions with high P303 frequency (Rootsi et al., 2012). The G2a-P16 (P18) lineage is specific to the Caucasus, accounts for a third of the Caucasian male gene pool and has a high frequency in the Southern and Northwestern Caucasus, with the highest frequency among North Ossetians (63.6%). Outside the Caucasus (Anatolia, Armenia, Russia, and Spain), the P16 lineage is either present at less than 1% or is absent (Rootsi et al., 2012).
Based on our network analysis, four subclades can be distinguished as L497, P303, M406, and P18 (under P16), as the samples with known haplogroups in the network are clustered on the same branches. Hungarians are included in L497, P303, and M406 subclades, except for P18 (1 Csángo sample), indicating that the gene flow from the P18 subclade of the Caucasus has been negligible. We included nine samples from the Rétköz and Váh valley populations, each of which falls into these three subgroups. The results indicate a close genetic relationship between the Hungarian-Balkar/Karachay and Hungarian-Avar males, because they share common haplotypes (cluster 1-3) or cluster in similar haplotypes (P303, L97, and M406). Thus, these observations likely suggest a common genetic origin, rather than coicidence.
The results of the aDNA studies in the Carpathian Basin show the haplogroup G2a-L156. One ancient Avar and four Hungarian Conqueror samples belong to the subgroup G2a-L293, two Hungarian Conqueror samples belong to haplogroup G2a-U1, and one ancient Hungarian sample belongs to haplogroup G2a-L30 (Neparáczki et al., 2019;Fóthi et al., 2020). SNPs L293 and L30 are downstream of L156 but upstream of M406. U1 and L497 are sister clades of P303 under M406 SNP (http://www. phylotree.org/Y/tree/G.htm). This indicates that SNPs L293 and L30 are older than SNPs M406 and P303, whereas SNPs U1 and L497 are the youngest.
As such, the modern and ancient DNA results support the case for a common genetic origin of the Hungarian and Caucasion populations.
As shown in Figure 7, the E1b1-M78 chromosomes of Hungarian-speaking populations originated in the Balkans and were introduced into Hungarian gene pool. Their location on the outer edge of the network, indicating many mutational steps, suggests that the M78 chromosomes in the Váh valley population appear to have separated from the common European ancestor earlier than those of the Rétköz population.
In aDNA studies, E1b1-M78 was observed among the Middle and Late period Avar (650-710 CE), as well as among the Hungarian Conquerors (895-mid Xth century) (Neparáczki et al., 2019). However, the absence of STR precludes its inclusion in the network analysis.
Based on the network topologies, the most diverse haplogroups are Q-M242, G2a-L156 and R1b-L23, which can be seen from the pattern of the networks, that is, the number of unique haplotypes of a male within the network is significant as well as the age of the accumulated STR variation (TMRCA) also supports it. The divergence times of haplogroups N-M46 and R1a-Z93 are almost the same (TMRCA) and within the network there are many clusters where several males share the same haplotype. Haplogroup E-M78 is the youngest, as the number of unique haplotypes is less than in the other networks.

Genetic structure
Based on Fst analysis, the Váh valley population shows a stronger genetic similarity to Balkan populations than does the Rétköz population. This is primarily due to the relatively high E and I2a haplogroup frequencies. Both haplogroups are common in Balkan populations but are less frequently found in Eastern Europeans (Regueiro et al., 2012;Doğan et al., 2016). The Rstbased MDS plot supports this observation, as the Váh valley popular was closer to a Balkan population (Croatia), whereas the nearest population to Rétköz is a Central European population (Czech). A further difference between the two studied Hungarian populations Frontiers in Genetics frontiersin.org 15 can be concluded from the Fst-based MDS plot. It shows that the Váh valley, which is at the western part of the Hungarian-Slovakian contact zone, is genetically closer to the Western Europeans and that Rétköz is at the eastern part of that zone and is therefore genetically closer to the Eastern Europeans. This study did not detect a genetic linkage between the Váh valley and the Hungarian Conquest period populations, and the results also reveal the Váh valley's relative isolation from neighboring populations.
In summary, the results obtained by us, which show that the genetic relationships between modern Hungarians, Hungarian Conquerors, Asian Huns (Xiongnu) and ancient Avars are continuous, are fully supported by the results of the whole genome sequencing data performed by Hungarian researchers (Maróthi et al., 2022).

Conclusion
The genetic composition of the Rétköz (Hungary) and Váh valley (Slovakia) populations indicate different histories. In the Rétköz population, the paternal lineages that were also found in the Hungarian Conquerors, such as haplogroups R1a-Z93, N-M46, Q-M242, R1b-L23, and G2a-L156, were better preserved. The genetic composition of the Váh valley population is similar to that of the surrounding Indo-European populations.
Further studies are needed to clarify if the common genetic footprints were acquired directly or indirectly. Comprehensive studies from European, Central Asian, and Caucasian populations should be conducted for the haplogroups using more downstream SNPs and NGS sequencing to learn more about the origins, expansion, and ethno-linguistic affiliations of the populations.

Data availability statement
The data presented in the study can be found in the YHRD respository, accession numbers: YA004754 and YA0047555.

Ethics statement
The studies involving human participants were reviewed and approved by the according to the recommenndations of the Data Protection Code the Humanities Research Center (MTA BTK-KP/450-17/2018). The patients/participants provided their written informed consent to participate in this study.

Author contributions
HP and EF designed the study and conducted the experiment. HP wrote the manuscript. EF and AF collected the samples. AT, DD, EF, VK, and HP analyzed the data. AT, DD and AF visualized data. All authors reviewed the manuscript. HP and EF contributed equally to this work.

Funding
This work was supported by the Thematic Excellence Program (TKP20-NKA-020) of the Ministry of Innovation and Technology in Hungary.