Introduction

Apart from linear chromosome, extrachromosomal DNA in circular form also exists in the nuclei of eukaryotes (Noer et al., 2022). Circular DNAs could be roughly classified into two groups based on their cell origins (Chiu et al., 2020). The one exclusively present in cancerous cells is usually referred to as ecDNA (extrachromosomal DNA), which is of megabase long in average and plays key roles in tumorigenesis (T. Wang et al., 2021). The ecDNA biogenesis is linked to ‘episome model’ (Carroll et al., 1988), chromothripsis (Shoshani et al., 2021), breakage-fusion-bridge (Coquelle et al., 2002), or translocation-deletion-amplification model (Van Roy et al., 2006). The other class found in somatic and germ cells is usually called eccDNA (extrachromosomal circular DNA). The size of eccDNA ranges from dozens of bases to hundreds of kilobases (Paulsen et al., 2018). Contrary to ecDNAs, the biogenesis mechanisms and biological functions of eccDNAs are relatively less experimentally characterized, and current studies show inconclusive or even contradictory results.

The genomic origins of eccDNAs have been extensively investigated in different cells and conditions with the application of Circle-seq and its refined derivatives (Mehta et al., 2020; Moller et al., 2015; Mouakkad-Montoya et al., 2021; Y. Wang et al., 2021), where eccDNAs are detected via rolling circle amplification and deep sequencing. While CpG islands (Shibata et al., 2012), gene-rich regions (Moller et al., 2018), and repeat elements, e.g., LTR (long terminal repeat) (Moller et al., 2015), LINE-1 (Dillon et al., 2015), segmental duplication (Mouakkad-Montoya et al., 2021), or satellite DNA (Mouakkad-Montoya et al., 2021), are hotspots for eccDNA formation, others found that eccDNAs are nearly random with regard to genomic distribution (Moller et al., 2020), or even made opposite observations (Henriksen et al., 2022). Epigenomically, the overall higher GC content, the periodicities of dinucleotide and eccDNA size convergently point to that nucleosome wrapping of DNA might contribute to the formation of small-sized eccDNAs (Shibata et al., 2012; Y. Wang et al., 2021), an intriguing starting point for mechanistic understanding of eccDNA origination. However, direct evidence for coincident positioning of eccDNAs and nucleosomes is still lacking, not to mention specific epigenetic marks on nucleosomes that are tightly associated with eccDNA formation.

EccDNAs are increased upon DNA damages (Moller et al., 2015; Paulsen et al., 2021), suggesting them as byproducts of successive DNA repairs. Among diverse repair pathways, it was reported that eccDNA levels particularly depend on resection after double-strand DNA break (DSB) and repair by microhomology mediated end joining (MMEJ) (Paulsen et al., 2021; Y. Wang et al., 2021). In further support of the involvement of MMEJ, microhomology is found around eccDNA break-points (Lukaszewicz et al., 2021; Moller et al., 2015). However, in all studies using short-read sequencing technologies, eccDNA break-points are mis-annotated if microhomologous sequences are present around due to their interference to precise break-point detection (see our main text) (Dillon et al., 2015; Henriksen et al., 2022; Kumar et al., 2017; Lv et al., 2022; Mann et al., 2022; Moller et al., 2018; Moller et al., 2020; Paulsen et al., 2019; Prada-Luengo et al., 2019; Shibata et al., 2012; Sin et al., 2020; Y. Wang et al., 2021; P. Zhang et al., 2021), or the eccDNA identification doesn’t depend on precise break-point detections at all (Moller et al., 2015; Mouakkad-Montoya et al., 2021). The contribution of microhomology to eccDNA generations thus needs to be revisited with precise mapping of break-points.

Alternative or additional mechanisms might be involved in germline eccDNA formation. During meiosis, two spatially-closed break sites catalyzed by SPO11 at recombination hotspots may release eccDNAs accompanied by de novo deletions on linear chromosomes (Lukaszewicz et al., 2021). Consistently, germline microdeletions display similar sequence features with eccDNAs (Shibata et al., 2012). However, a recent study reported that the creation of germline eccDNAs negatively correlate with meiotic recombination rates (Henriksen et al., 2022). Therefore, it remains to be determined whehter meiosis might significantly contribute to eccDNA biogenesis.

We envision that eccDNA landscape during spermatogenesis is ideal for clarifying the abovementioned issues and so better understand the biogenesis mechanisms and biological implications of eccDNAs. Only a small fraction of histones will survive from the histone-to-protamine transition in mature sperms, allowing us to more specifically correlate eccDNA origination with histones. Studying eccDNAs in germline cells rather than somatic cells could help reveal to what extent meiosis might contribute to eccDNA generation and de novo structural variations that can be passed to offspring. Therefore, in this study, we profiled eccDNAs via Circle-seq in human sperms and different developmental stages of mouse germ cells with an improved analysis pipeline to identify eccDNAs at nucleotide resolution. We conclude that germline eccDNAs are formed largely by microhomology mediated ligation of nucleosome-protected fragments, and barely contribute to de novo genomic deletions at meiotic recombination hotspots.

Results

Widespread eccDNA formation in human and mouse germline cells

Because it was reported that high content of eccDNAs existed in sperms (Henriksen et al., 2022), we examined the genome-wide eccDNA landscape in two human sperm samples by Circle-seq (Moller et al., 2018; Moller et al., 2015; Moller et al., 2020) (see Materials and methods), and indeed found that there were widespread eccDNAs across the human genome (Figure 1A, Figure 1—figure supplement 1A). This motivated us to further investigate the biogenesis mechanism, particularly whether it might be linked to specific spermatogenesis processes. Given that it is ethically prohibited and technically challenging to collect pure spermatogenic cell types from human individuals, we therefore turned to use mouse model to study the eccDNA formation during spermatogenesis.

Overview of eccDNA formation during mouse spermatogenesis. (A) Schematic representation of Circle-seq in human sperm cells and mouse SPA, RST, EST and sperm cells validated with immunochemistry. SYCP3: a component of the synaptonemal complex; γH2AX: a marker for double-strand breaks; SP65: a marker for acrosome organelle; TUBULIN: structural component of manchette in EST and flagellum axoneme in sperm cells. (B) Number of eccDNAs detected in different cell types. *two-sided t-test p-value < 0.05; ** two-sided t-test p-value <0.01. (C) Size distribution of eccDNAs during mouse spermatogenesis. Dotted lines indicate multiplies of 180 bp. (D) A representative genomic locus showing the gene annotation, Circle-seq signals, detected eccDNAs, and SINE and DNA repeat elements. Highlighted in red rectangle is a large-sized eccDNA. (E) Enrichment of eccDNAs at given genomic regions relative to randomly-selected control regions. (F) Enrichment of eccDNAs at given repeat elements relative to randomly-selected control regions.

A series of cell divisions and morphological changes are involved in spermatogenesis, where spermatogonial stem cells develop into spermatocytes (SPA) via mitosis, and SPA then undergo meiosis to produce haploid round spermatids (RST), which will take a dramatic morphological change and chromatin compaction to produce elongated spermatids (EST) and finally matured sperms (Hess et al., 2008; Roosen-Runge, 1962). We isolated SPA, RST and EST using flow cytometry and collected sperms from mouse cauda epididymis (Hayama et al., 2016) for subsequent Circle-seq (Moller et al., 2018; Moller et al., 2015; Moller et al., 2020) (see Materials and methods). All four cell types were validated with known markers and cell morphology (Figure 1A). To account for sample variations, up to five biological replicates and ~150 million reads for each cell type were sequenced for eccDNA detection. From ~1,500 to ~180,000 high-confidence eccDNAs were identified, suggesting widespread circular DNA formation during mouse spermatogenesis (Figure 1B; see Materials and methods). The reproducible rate of eccDNAs with 50% reciprocal overlap between biological replicates was only ~2.4% in average, a level comparable to previous studies (Henriksen et al., 2022; Moller et al., 2018) (Figure 1—figure supplement 1B). As noted earlier (Moller et al., 2018), the detected eccDNAs seemed not saturated (Figure 1—figure supplement 1C; see Discussion and Materials and methods), which might underlie the observed low reproducibility. Nevertheless, principal component analysis suggested that the within-group similarity was marginally higher than the between-group similarity (Figure 1—figure supplement 1D), allowing investigation of stage-specific eccDNA features during mouse spermatogenesis.

The detected germline eccDNAs verified known genomic features of eccDNAs. First, the natural size distribution of eccDNA is usually distorted in Circle-seq as smaller eccDNAs tend to be overrepresented in rolling circle amplification (Mohsen et al., 2016; Moller et al., 2018; Moller et al., 2015). As expected, the detected eccDNA population was dominated by small-sized eccDNAs, most of which were ~180bp or ~360bp long (Figure 1C). However, eccDNA size could occasionally reach to tens of kilobases (Figure 1D). Second, eccDNAs from different cell types were all enriched at gene-rich regions, especially 5’UTR (Figure 1E, Figure 1—figure supplement 1E), corroborating the reported association between eccDNA frequency and gene density in somatic cells (Dillon et al., 2015; Shibata et al., 2012). EccDNAs were also highly associated with SINE but not LINE elements (Figure 1F), and quantitative analysis revealed that eccDNA biogenesis was positively correlated with SINE density (Figure 1—figure supplement 1F), but negatively correlated with LINE density (Figure 1—figure supplement 1G). Given that SINE and LINE elements function to orchestrate chromosomes into gene-rich A compartment and gene-poor B compartment respectively (Lu et al., 2021), the positive correlation between eccDNAs and SINE elements might further support that eccDNAs are overall highly associated with the gene-rich regions. Interestingly, we also noticed a strong association between eccDNAs and DNA transposons (Figure 1F), suggesting that DNA transposons might get circularized rather than or in addition to re-integrated into the genome, an interesting possibility awaiting further investigations. Altogether, the genome-wide eccDNA landscape during mouse spermatogenesis allow us to further study the biogenesis mechanism and function of eccDNAs.

High eccDNA load and periodic eccDNA size distribution in mouse sperm cells

Notably, sperm cells had 97,372 eccDNAs detected in average, a number significantly higher than those in SPA (15,246), RST (18,426) and EST (3,591) cells (Figure 1B). SPA cells did not show higher eccDNA numbers (Figure 1B), suggesting that meiosis does not seem to contribute significantly to eccDNA biogenesis. Since the same amount of eccDNAs (10ng) were used for library construction and all samples were sequenced in comparable and sufficiently-deep depth, it suggests that eccDNA species in sperm cells has higher complexity. However, the higher starting cell number for sperm cells might account for the larger diversity of sperm eccDNA species (see Discussion and Materials and methods), otherwise, it would be interesting to explore any specific features of sperm cells underlying the higher load of eccDNAs.

In contrast to SPA, EST and RST eccDNAs showing the unimodal distribution that was centered at ~180bp, sperm-derived eccDNAs showed a multimodal distribution with a pronounced periodicity of ~180bp (Figure 1C), which was readily seen in individual samples (Figure 1—figure supplement 2). Given that each nucleosome consists of 147bp DNA wrapping itself around a histone core, the ~180bp-long fragments likely corresponded to histone core region plus ~20-30bp linker regions, as observed in apoptotic cells (Matassov et al., 2004). Although the identified eccDNAs in all spermatogenesis stages were likely related to nucleosomes, the different modes of size distribution might be due to distinct nucleosome compositions and structures between sperm and other spermatogenic cells.

Mouse sperm eccDNAs come from DNA fragments protected by histones

Only a small fraction of histones will be retained in mouse sperm cells after histone-to-protamine transition (Torres-Flores et al., 2020), permitting us to more specifically correlate eccDNAs with histones. We were therefore motivated to see whether the detected eccDNAs were derived from the retained histones in mature sperm cells. We noted that sperm eccDNAs had higher GC content than surrounding regions as well as control regions randomly selected across the genome (Figure 2A), resembling the sequence feature of nucleosome-protected DNA fragments. Consistently, sequence-based prediction revealed significantly higher nucleosome occupancy probability for ~180 bp (from 175bp to 185bp) and ~360bp (from 355bp to 365bp) sperm eccDNA regions (Figure 2B; see Materials and methods). A small dip was observed at the center of ~360bp eccDNA regions, which likely corresponded to the linker region between two nucleosomes (Figure 2B, right).

Association between sperm eccDNAs and nucleosome positioning. (A) GC contents of sperm eccDNAs, regions upstream and downstream of eccDNAs, and randomly-selected length-matched control regions. ***two-sided Wilcoxon test p-value <0.001. (B) Predicted probability of nucleosome occupancy for eccDNA and randomly-selected length-matched control regions (highlighted by red-shaded area), and surrounding regions. Boxplots showing the probability distribution of individual eccDNAs and control regions. ***two-sided Wilcoxon test p-value <0.001. (C) Enrichment of eccDNAs at different histones and histone modifications. (D and E) ChIP-seq signal distribution at [-1.8kb, +1.8kb] of ~180bp (D) and ~360 bp (E) eccDNA regions. ChIP-seq signals quantified as reads density are color-coded below heatmaps.

It is a common practice to re-use publicly-available genomics data generated in the same cell types for integrative analysis. Taking advantage of public ChIP-seq data for histones and their modifications in mouse sperm cells (Jung et al., 2019; Jung et al., 2017; Singh et al., 2021), we found that eccDNAs were significantly enriched with certain histone variants and modifications (Figure 2C), and 7.46% of sperm eccDNAs in total were intersected with at least one ChlP-seq peaks. Considering that histones occasionally retained in sperms might not generate strong ChIP-seq signals exceeding the peak calling cutoff, a meta-gene analysis of ChIP-seq signals at and around sperm eccDNA regions will likely provide more insights. Interestingly, enrichment of H3 histone and H2A.Z, TH2A and TH2B histone variants but depletion of H3.3 variant was observed at ~180bp sperm eccDNA regions (Figure 2D). These eccDNAs also showed strong associations with H3K27ac, H3K4me1, H3K9ac and H3K27me3 modifications, however, no enrichment was seen for H3K4me3, and H3K36me3 and H3K9me3 signals were comparable with or even lower than randomly-selected regions as control (Figure 2D).

We next examined ~360bp sperm eccDNAs which supposedly correspond to two nucleosomes, and made similar observations. Centers of ~360bp eccDNAs were well positioned between two adjacent nucleosomes consisting H3 histone and H2A.Z histone variants, and H3K27ac, H3K4me1, H3K9ac and H3K27me3 histone modifications (Figure 2E). Similar to ~180bp eccDNAs, ~360bp eccDNAs did not show association with H3.3 or H3K4me3, or stronger association than randomly-selected regions with H3K36me3 and H3K9me3 either (Figure 2E). Although H3.3 variant coincides with active transcription, it is also well known for its localization at heterochromatin region and its roles in promoting heterochromatin formation by inhibiting H3K9/K36 histone demethylase (Udugama et al., 2022). Together, euchromatin is generally more preferred than heterochromatin for eccDNA biogenesis, which is consistent with the enrichment of sperm eccDNAs at gene-rich regions (Figure 1D).

Large-sized eccDNAs are preferentially generated from heterochromatin regions

Intriguingly, periodic distribution of nucleosomes was observed for ~360bp eccDNAs but not for ~180bp ones, indicating that eccDNAs from di-nucleosomes but not mono-nucleosomes preferentially originate from well-positioned nucleosome arrays (Figure 2E). We were further prompted to ask whether eccDNAs of different sizes are originated from different genomic regions. Indeed, small-sized eccDNAs (<3kb) were more enriched at H3K27ac-marked euchromatin regions, while large-sized ones (≥3kb) at H3K9me3-marked heterochromatin regions (Figure 3A). Accordingly, small-sized eccDNAs were generally more associated with genic regions, while large-sized ones with non-genic regions (Figure 3B). Since sperm eccDNAs in this study were dominantly small-sized ones (Figure 1C), strong enrichment of eccDNAs at genic regions was observed (Figure 1E). However, strong depletion at genic regions was reported for human sperm eccDNAs in a recent study (Henriksen et al., 2022). Close inspection suggests that the discrepancy is partially reconciled in the light of two eccDNA groups of different sizes. Henriksen et al studied eccDNAs with the size largely ranging from ~3kb to 50kb (Henriksen et al., 2022), rather than small-sized ones reported by us and many others (Dillon et al., 2015; Moller et al., 2018; Moller et al., 2020; Paulsen et al., 2019; Shibata et al., 2012; Y. Wang et al., 2021). In contrast to small-sized eccDNAs, the large-sized sperm eccDNAs detected in this study also displayed weak but apparent negative correlation with gene density and Alu elements (Figure 3C and D).

Large-sized eccDNAs are preferentially generated from heterochromatin regions. (A) Distribution at H3K27ac- and H3K9me3-marked regions for eccDNAs of different sizes. (B) Distribution at different genomic regions for eccDNAs of different sizes. (C) Number of small (<3kb) vs large (≥3kb) eccDNAs per Mb as a function of gene number per Mb. Pearson correlation coefficients and two-sided t-test p-values are indicated. (D) Number of small (<3kb) vs large (≥3kb) eccDNAs per Mb as a function of Alu number per Mb. Pearson correlation coefficients and two-sided t-test p-values are indicated.

Germline eccDNAs as apoptotic products are not asscociated with meiotic recombination hotspots

The observed association between eccDNA and oligonucleosomal DNA fragmentation (Figure 2) is a typical feature of cell death. The spontaneous death of germ cells has been observed during the normal spermatogenesis (Liu et al., 2017; Shaha et al., 2010; Weinbauer et al., 2001; Young et al., 2001), however, it is still debatable whether spermatids and sperm can undergo apoptosis (Lachaud et al., 2004). Thus, sperm-derived eccDNAs might be associated with apoptosis (if exists) or un-programed cell death of germ cells during the normal spermatogenesis (see also Discussion). In support of this hypothesis, all features associated with mouse germline eccDNAs identified in this study (Figure 1C, E and F) closely matched with those of eccDNAs whose generation is dependent on apoptotic DNA fragmentation in somatic cells (Figure 4—figure supplement 1) (Y. Wang et al., 2021).

During meiosis, two spatially-closed cleavage sites catalyzed by SPO11 at recombination hotspots could release eccDNAs and generate de novo genomic deletions (Lukaszewicz et al., 2021), which, if transmitted to offspring, might contribute to structural variations within population. Since most sperm eccDNAs likely result from oligonucleosomal DNA fragments in sperm cells undergoing cell death (Figures 2 and 3) and SPA cells undergoing meiosis does give rise to more eccDNAs than other cells (Figure 1B), meiotic recombination is unlikely the major mechanism for germline eccDNA generation. To test this hypothesis, we first investigated to what extent eccDNA break-points well correspond to recombination hotspots defined as SPO11 or PRDM9 binding sites (Alleva et al., 2021; Lange et al., 2016). We noted that there was only a small number of eccDNAs with both break-points located in one recombination hotspot or two different hotspots (Figure 4A). These eccDNAs only constituted <0.15% (or <350) of mouse germline eccDNAs, suggesting a very low level of coincidence between eccDNA generation and meiotic recombination (Figure 4A). Consistently, only dozens of, or a few hundred eccDNAs in mouse germline cells coincided with known genomic deletions within mouse population (Figure 4B). Altogher, germline eccDNAs are likely apoptotic products that are not asscociated with meiotic recombination hotspots and heritable genomic deletions.

Microhomology directed ligation accounts for emergence of most eccDNAs. (A) Numbers of eccDNAs or randomly-selected control regions overlapped with recombination hotspots in mouse. (B) Shown are numbers of mouse sperm eccDNAs or randomly-selected control regions having 95% reciprocal overlap with different types of structural variations. (C) Illustrated are how an eccDNA with homologous sequences (CGA) at two ends is identified from short-read sequencing data by our methods vs other methods. (D) Percentages of homologous sequences of different lengths (coded by different color saturation levels) are shown for eccDNAs and randomly-selected control regions. (E) GC content of homologous sequences and randomly-selected control regions. (F) Percentages of homologous sequences of different lengths (coded by different color saturation levels) are shown for small-sized (<3kb) and large-sized (≥3kb) eccDNAs. (G) Length of microhomologous sequences as a function of the eccDNA size. Data points are shown as median plus lower (25%) and upper (75%) quartiles. The shaded area is 95% confidence interval of linear regression line. Pearson correlation coefficient and two-sided t-test p-value are indicated. (H) Sequencing motif analysis for ±10bp leftmost left ends and ±10bp leftmost right ends of eccDNAs with no perfectly matched homologous sequences observed. (I) Model for MMEJ-directed eccDNA biogenesis.

Microhomology directed ligation is the major biogenesis mechanism of germline eccDNAs

We therefore further explored how nucleosome-protected DNA fragments get circularized into eccDNAs. As suggested by previous studies, MMEJ is implicated in eccDNA biogenesis (Lukaszewicz et al., 2021; Moller et al., 2015; Paulsen et al., 2021; Y. Wang et al., 2021). The precise distribution of microhomologous sequences relative to eccDNA break-points will help better understand how and to what extent MMEJ might contribute to eccDNA biogenesis. However, we noted that the presence of microhomologous sequences will hinder precision eccDNA break-point identification (Figure 4C), which is not well dealt with by all existing methods (Figure 4—figure supplement 2A) (Kumar et al., 2017; Mann et al., 2022; Prada-Luengo et al., 2019; P. Zhang et al., 2021). Short sequencing reads spanning eccDNA break-points will be mapped to the genome as split reads, with its first part mapped to the right end of eccDNA, and the second part to the left end. If the sequence in front of the left eccDNA end is homologous to the right eccDNA end, or if the sequence following the right eccDNA end is homologous to the left eccDNA end, the homologous regions will be included in both parts of split reads to reach to a maximal length of matches, and many existing methods will mistake the eccDNA plus two homologous regions as the whole eccDNA region (Figure 4C). Being aware of it, we developed a base-resolution method for eccDNA identification on the basis of previous efforts (Figure 4—figure supplement 2B) (Kumar et al., 2017; Moller et al., 2018). When homologous sequences are present, we record the coordinates of the leftmost form of eccDNA and an offset corresponding to the length of homologous sequences to represent all possible eccDNA variants (Figure 4C). Our method has comparable sensitivity and specificity with existing methods (Figure 4—figure supplement 2C), and is easy-to-use and freely-available.

In contrast to simulated controls (15%), more than one third of eccDNAs had ≥1bp homologous sequences, most of which were shorter than 5bp (Figure 4D), suggesting the involvement of MMEJ in eccDNA biogenesis. The GC content of homologous sequences was higher than that of simulated control regions, permitting stronger base-pairing for efficient MMEJ (Figure 4E). Considering that two free-ends of long DNA fragments might be not as spatially close as those of short DNA fragments, formation of longer eccDNA should more rely on longer homologous sequences for stable base pairing. Indeed, large-sized eccDNAs in SPA, RST and EST cells did show higher percentage of ≥2bp homology than small-sized eccDNAs, and large-sized eccDNAs in sperm cells showed higher percentage of >5bp homology (Figure 4F). A significant positive correlation between lengths of homologous sequences and eccDNA sizes was observed (Figure 4G). We further reasoned that for the remaining two thirds of eccDNAs that were lack of perfectly-matched homologous sequences, imperfect homologous sequences might be present. Accordingly, we noted the same sequence motifs between eccDNA starts and sequences following eccDNA ends, and between eccDNA ends and sequences in front of eccDNA starts (Figure 4H). Similar observations have been made also by others before (Sin et al., 2020), however they failed to precisely locate the homologous sequences relative to eccDNA break-points. Altogether, we propose that sticky ends of DNA fragments with homologous sequences might base pair with each other, and then be ligated by DNA ligase, e.g., DNA ligase III (Y. Wang et al., 2021), to form eccDNAs (Figure 4I). In sum, among all proposed mechanisms (Chiu et al., 2020; Dillon et al., 2015; Moller et al., 2015; Paulsen et al., 2021; Sin et al., 2020), MMEJ-mediated ligation accounts for emergence of most eccDNAs at least in germline cells.

The biogenesis mechanism of germline eccDNAs is conserved between human and mouse

Upon revealing the major biogenesis mechanism of mouse germline eccDNAs, we further examined whether the mechanism is conserved in human by analyzing the features of human sperm eccDNAs detected in this study (see Materials and methods). Similarly, human sperm eccDNAs were originated from oligonucleosomal fragmentation as well, as suggested by the pronounced size periodicity of ~180bp (Figure 5A), higher GC content and nucleosome occupancy probability (Figure 5B and C). We further performed base-level comparison between eccDNAs and meiotic recombination hotspots. Consistent with what we observed in mouse, only about 600 eccDNAs (<0.7% of human sperm eccDNAs) were located in one recombination hotspot or two different hotspots (Figure 5D), and only dozens of, or a few hundred eccDNAs in human sperm cells coincided with known genomic deletions (Figure 5E). Overall, our analysis disfavors the meiotic recombination to eccDNA biologenesis. Instead, we observed higher frequency of micro-homologous sequences for eccDNAs than simulated control regions (Figure 5F). Moreover, large-sized eccDNAs had longer homologous sequences possibly for stable base pairing, as indicated by the strong positive correlation between lengths of homologous sequences and eccDNA sizes (Figure 5G). Altogether, these results suggest microhomology directed ligation of nucleosome protected DNA fragments is a conserved pathway for germline eccDNA generation in both human and mouse.

The biogenesis mechanism of germline eccDNAs is conserved between human and mouse. (A) Size distribution of sperm eccDNAs in two biological replicates. (B) GC contents of sperm eccDNAs, regions upstream and downstream of eccDNAs, and randomly-selected length-matched control regions. ***two-sided Wilcoxon test p-value <0.001. (C) Predicted probability of nucleosome occupancy for eccDNA and randomly-selected length-matched control regions (highlighted by red-shaded area), and surrounding regions. Boxplots showing the probability distribution of individual eccDNAs and control regions. ***two-sided Wilcoxon test p-value <0.001. (D) Numbers of eccDNAs or randomly-selected control regions overlapped with recombination hotspots in human. EccDNAs located completely within a hotspot (Intra-), or with both ends overlapped with two different hotspots (Inter-) are shown separately. (E) Shown are numbers of human sperm eccDNAs or randomly-selected control regions having 95% reciprocal overlap with different types of structural variations. (F) Percentages of homologous sequences of different lengths (coded by different color saturation levels) are shown for eccDNAs and randomly-selected control regions. (G) Length of microhomologous sequences as a function of the eccDNA size. Data points are shown as median plus lower (25%) and upper (75%) quartiles. The shaded area is 95% confidence interval of linear regression line. Pearson correlation coefficient and two-sided t-test p-value are indicated.

Discussion

The biogenesis mechanisms and biological implications of eccDNAs are big puzzles. Increasingly evidence has emerged to relate eccDNA formation with nucleosome-protected DNA fragments, however, to our knowledge, no study provides direct evidences. The comparison between high number of sperm eccDNAs and retained histones in sperm cells after histone-to-protamine transition (Torres-Flores et al., 2020) provide a good way to clarify the abovementioned issues. Here, we profiled eccDNA landscape during mouse spermatogenesis, and explicitly linked eccDNAs with nucleosome-protected DNA fragments (Figure 2). According to these results, we speculate that any oligonucleosomal fragments might have a chance of eccDNA formation as long as microhomologous sequences are present at the ends. In other words, any single cells might generate distinct sets of eccDNAs and thus the theoretical number of eccDNA is extremely large. This might explain why sperm cells, which have a higher starting cell number, have more eccDNAs detected, even though the same amount of eccDNAs (10 ng) for each cell type was used for library preparation and high-coverage sequencing (~150 million reads).

Multiple biogenesis mechanisms of eccDNA have been proposed (Chiu et al., 2020; Dillon et al., 2015; Moller et al., 2015; Paulsen et al., 2021; Sin et al., 2020), it is however unclear which one plays the dominant role. With our nucleotide-resolution eccDNA break-point detection method, we clearly demonstrated that microhomologous sequences are present at boundaries of most eccDNAs (Figure 4), strongly suggesting the MMEJ-medicated eccDNA formation as the major mechanism for eccDNA biogenesis. The biogenesis mechanism of germline eccDNAs is conserved between human and mouse. However, it remains unclear whether eccDNA generation from healthy cells not undergoing apoptosis is also medicated by MMEJ.

Spontaneous germ cell death has been observed during the normal spermatogenesis (Liu et al., 2017; Shaha et al., 2010; Weinbauer et al., 2001; Young et al., 2001). Given the periodic size distribution (Figure 1C) and similar genomic features with eccDNAs from apoptotic DNA fragmentation in somatic cells (Figure 1 and Figure 4—figure supplement 1), we reason that mouse germline eccDNAs may be also the germ cell death products. As the final place for sperm maturation and storage, epididymis might contain unhealthy sperm cells undergoing cell death and so DNA fragments, which might additionally contribute to the high eccDNA loads of sperm cells than others. It is possible that failure of histone-to-protamine exchange might account for some sperm cell death. It has been shown that H3.3, generally linked to H3K4me3 marks, plays key roles in modulating TP1 removal and PRM1 incorporation for nucleosome eviction and replacement by protamine (Wang et al., 2019), and loss of H3.3 will increase the cell death rate (Yuen et al., 2014). H3K36me3 and H3K9me3, whose demethylase is inhibited by H3.3 (Udugama et al., 2022), are also involved in histone replacement (Wang et al., 2019), failure of which might underlie the eccDNA generation as well. In line with it, sperm eccDNA regions correspond to depletion of H3.3 variant and H3K4me3, H3K36me3 and H3K9me3 histone modifications (Figure 2D and E). This suggests the number of sperm eccDNAs might be highly associated with cell death, and sperm eccDNA may serve as a clinical biomarker for quality assessment of healthy sperms.

It is still debatable whether meiotic recombination is an important source of eccDNA generation. Both positive and negative correlations between eccDNAs and meiotic recombination rates at chromosomal level have been reported previously (Henriksen et al., 2022; Lukaszewicz et al., 2021). Our results suggest that this discrepancy can be largely reconciled in the light of two eccDNA groups of different sizes. We noted that small-sized and large-sized eccDNAs are preferentially derived from euchromatin and heterochromatins respectively (Figure 3). Since meiotic recombination hotspots are enriched at euchromatin regions (Lange et al., 2016), small-sized eccDNAs are therefore positively correlated with recombination rates but the opposite is true for large-sized eccDNAs (Figure 3—figure supplement 1). Therefore, the observed correlations between eccDNA density and recombination meiotic rate at chromosomal level are simply indirect reflections of their preferences to different chromatin regions.

We further investigated the coincidence between eccDNAs and meiotic recombination hotspots and structural variations at base level, and found that sperm cell eccDNAs does seem to be a major source or by-products of de novo deletions, likely because that extensive eccDNA formation might only occur at dead sperm cells that having no chance to be transmitted to the next generation. However, we cannot fully exclude the possibility that some eccDNAs emerging in normal sperm cells and do not affect their viability will be accompanied by de novo deletions or re-integrated into genome to create structural variations. Furthermore, they might stimulate innate immune responses (Y. Wang et al., 2021), server as extracellular vesicles encoding small RNAs (Paulsen et al., 2019), and be taken as biomarkers (Lv et al., 2022), which represent exiting research directions.

Methods

Isolation of mouse germ cells

Testes from adult C57BL/6 mice were decapsulated, and the seminiferous tubules were torn into small pieces. After incubating in 8 ml PBS containing 1 mg/ml collagenase (Sigma, C5138, St. Louis, MO, USA) and 1 mg/ml hyaluronidase (Sigma, H3506, St. Louis, MO) at 37°C for 6 min, the dispersed seminiferous tubules and cells were incubated at 37°C for 5 min with gentle shaking. Cells were collected by centrifugation at 200 × g for 5 min at 4°C, and washed once with PBS, resuspended in 15 ml PBS containing 0.25% Trypsin (Gibco, 25200072) and 1 mg/ml DNase I (Sigma, AMPD1), and incubated at 37°C for 5 min with gentle shaking. Thereafter, cells were collected and washed once with PBS. After filtration through a 70 μm Nylon Cell Strainer (BD Falcon, 352350), the cells were suspended in 30 mL of PBS and stained at 37°C with 5 mg/mL Hoechst 33342 (Thermo Scientific, 62249) for 30 min. Hoechst fluorescence were detected with a 450-nm band-pass filter for blue fluorescent (Hoechst Blue) or a 675-nm band-pass filter for red fluorescent (Hoechst Red). Mouse spermatocytes, round spermatids and elongating spermatids were sorted by BD FACSAria™ Fusion (BD Biosciences). Mouse mature spermatozoa were collected form mouse cauda epididymis. The cauda epididymis was quickly cut into pieces and incubated in 1 ml pre-warmed human tubal fluid (HTF) (Millipore, MR-070-D) for 15 min at 37°C, thus allowing the mature spermatozoa to release from the tissue. After filtration through a 40 μm Nylon Cell Strainer (BD Falcon, 352340), mature spermatozoa were collected and washed three times with PBS. About 106 SPA, RST and EST cells and 107 sperm cells were obtained and used for the eccDNA extraction.

Immunofluorescence

Mouse germ cells were spread on glass slides that were air dried. The slides were fixed in 4% PFA at room temperature for 5 min, washed with PBS three times. After blocking with 5% bovine serum albumin, the primary antibodies (rabbit anti-SYCP3 pAb (Abcam, ab15093, 1:400); mouse anti-γH2AX pAb (Millipore, 05-636, 1:400); Mouse anti-sp56 mAb (QED Bioscience, 55101, 1:200); Rabbit Tubulin pAb (ABclonal, AC007, 1:100)) were added to the slides, followed by 14-16 h incubation at 4°C. After washing the slides with PBS, the secondary antibodies were added, followed by 1.5 h incubation. The cell nuclei were stained with DAPI for 5 min. Images were observed using a fluorescence microscope SP8 microscope (Leica).

Human adult sperm sample preparation

The sperm donation candidates in the present study were healthy young Chinese men. Each candidate completed a medical examination and extensive medical/social questionnaire to exclude any potential individuals with genetic or major medical problems (such as cardiovascular diseases and sexually transmitted diseases) listed in the Basic Standard and Technical Norms of Human Sperm Bank published by Chinese Ministry of Health. Smokers, drug abusers, and heavy drinkers were also excluded. The rest of the candidates signed a voluntary sperm donation informed consent and agreed to live in Beijing for at least 6 months. The sperm bank also recorded the candidates’ age, date of birth, and date of semen collection. The ethical approval in this study were provided by the Reproductive Study Ethics Committee of Peking University Third Hospital (2017SZ-035). Semen samples were selected through 40% density gradient of PureSperm (Nidacon International, Molndal, Sweden) by centrifugation (500g, 30 min) at room temperature and washed with phosphate-buffered saline (PBS) for three times; the obtained spermatozoa were used for the eccDNA extraction.

Circle-Seq

Purification of mouse and human eccDNAs was performed as previously described, with minor modifications (Moller et al., 2018). In brief, samples were resuspended in 500 μl Lysis solution (10mM Tris pH7.4, 100mM NaCl, 1% SDS, 1% Sarkosyl, 150mM DTT) with 10μL Proteinase K (Thermo Scientific, EO0491), and incubated overnight at 55 °C. After cell lysis, phenol: chloroform: isoamyl alcohol was added and mixed, centrifuged at 13,000 g for 15min at 4°C. The supernatant was moved to a new tube, incubated with 500 μL isopropanol at room temperature for 10 min, and centrifuged at 13,000 g for 15min at 4 °C. The resulted DNA pellet was washed with 1 ml 70% ethanol, treated with alkaline to separate chromosomal DNA from eccDNAs by rapid DNA denaturing–renaturing, followed by column chromatography on an ion exchange membrane column (TIANprep Mini Plasmid Kit, DP103). Column-bound DNA was eluted in TE buffer (10 mM Tris-Cl, pH 8.0; 1 mM EDTA, pH 8.0). The remaining linear DNA was treated by exonuclease (Plasmid-Safe ATP-dependent DNase, Lucigen, E3101K) at 37°C for 1 week, during which additional ATP and DNase was added every 24 h (30 units per day) according to the manufacturer’s protocol (Plasmid-Safe ATP-dependent DNase, Lucigen, E3101K). The eccDNA samples were cleaned by phenol: chloroform: isoamyl alcohol once, followed by ethanol precipitation. EeccDNAs were then amplified by phi29 polymerase (NEB, M0269L) at 30°C for 16 h. Phi29-amplified DNA was cleaned by phenol: chloroform: isoamyl alcohol once, followed by ethanol precipitation. The DNA samples were sonicated to a set size of 250 bp with an M220 Focused-ultra sonicator (Covaris, Woburn). Sequencing libraries were generated with NEBNext® UltraTM II DNA Library Prep Kit for Illumina (NEB, E7645S), according to the manufacturer’s instructions. 10 ng of eccDNA samples were used for library preparation. NEBNext Multiplex Oligos for Illumina (Set 1, NEB #E7600) were used for PCR amplification of adaptor-ligated DNA. Libraries were purified with SPRIselect® Reagent Kit (Beckman Coulter, Inc #B23317). Paired-end 150 bp sequencing was performed on Illumina NovaSeq 6000 System.

EccDNA detection

We developed a nucleotide-resolution eccDNA detection pipeline on the basis of previous efforts (Figure 4—figure supplement 1B) (Kumar et al., 2017; Moller et al., 2018). Briefly, SeqPrep (https://github.com/jstjohn/SeqPrep) was used to trim adapter sequences and merge the overlapping paired-end reads into singleton long reads, followed by reads mapping to GRCm38 reference genome using BWA MEM version 0.7.17-r1188 (Li et al., 2009). Samblaster version 0.1.26 (Faust et al., 2014) or an in-house Perl script was used to remove PCR duplicates and separate alignments into split reads, discordant and concordant reads. Candidate eccDNAs are firstly identified based on split reads (high-confidence ones). When homologous sequences are present at eccDNA ends, we recorded the coordinates of the leftmost form of eccDNA and an offset corresponding to the length of homologous sequences to represent all possible eccDNA variants. Potential split reads that failed to be mapped as split reads in the first place (low-confidence ones) as well as discordant reads were identified and counted using in-house Perl scripts. The average coverages (in terms of RPK) for candidate eccDNAs and surrounding regions were then calculated based on all different types of reads. Any eccDNA supported by at least two high-confidence split reads or discordant reads, with its 95% region covered by at least one read, and with its average coverage twice of that of its surrounding region, is considered as a high-confidence eccDNA.

EccDNA characterization

Gene structure annotations in mouse and human are based on Ensembl GRCm38 Release 102 and Ensembl GRCh38 Release 104, respectively. Repeat elements are annotated by RepeatMasker database. Random region generation, sequence composition calculation, and overlapping region determination were all done with Bedtools version 2.30.0 (Quinlan et al., 2010). The nucleosome occupancy probability of the eccDNAs and the surrounding regions was predicted using predNuPoP function of NuPoP version 2.2.0 R package (Xi et al., 2010). Motif analysis and visualization were done by Two Sample Logo web server (Vacic et al., 2006). Public ChIP-seq datasets (Jung et al., 2019; Jung et al., 2017; Singh et al., 2021) were aligned to GRCm38 reference genome using Bowtie, and sorted and indexed using SAMtools version 1.7 (Danecek et al., 2021). Picard MarkDuplicates version 2.18.14 was used to remove PCR duplicates. Peak calling was done with Macs2 version 2.2.7.1 (Y. Zhang et al., 2008). Coverage file was generated from BAM file and visualized using deeptools version 3.5.1 (Ramirez et al., 2016).

Data and software access

The raw sequencing data reported in this study have been deposited in the Genome Sequence Archive (Chen et al., 2021) in National Genomics Data Center (Members et al., 2022) (GSA accession numbers: CRA008015 and HRA002957) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa. The eccDNA mapping and detection workflows are available at https://github.com/NjuChenlab/eccDNA_detector_tools.

Competing interest statement

The authors declare no competing financial interests.

Acknowledgements

The authors are grateful to members of the J.-Y.C., W.L. and C.L. labs for cooperation, reagent sharing, and insightful discussions during the course of this investigation. This study is supported by the MOST [2021YFF1201500 to J.-Y.C] and NSFC [32170653 to J.-Y.C, 81925015 and 32230029 to W.L., and 32270898 to C.L.].

Author contributions

J.-Y.C., W.L, C.L. and H.J. conceived and designed the experiments. C.L. and Z.Z. performed the majority of the experiments with the help from S.X., Y.C., Y.C. and J.W., J.H. was responsible for the processing and analysis of most genomic data. J.-Y.C. interpreted the data and wrote the paper.

Detection and characterization of eccDNAs. (A) Genome-wide distribution of eccDNAs detected in two human sperm samples. (B) Histograms showing distributions of reproducible rates between any two replicates of the same biological sample. Left: data generated in this study; Right: public datasets [Moller et al, Nat Commun (2018) 9:1069; Henriksen et al, Mol Cell (2022) 82: 209-217.e7]. (C) Number of detected eccDNAs as a function of non-redundant uniquely-mapped reads (useful reads) randomly sampled from indicated Circle-seq sequencing libraries. (D) Principal component analysis of biological replicates using Jaccard indexes with all other replicates as features. 75% confidence ellipses are indicated. (E-G) eccDNA number per 2 mega bases as a function of gene density (E), density of SINE (F) and LINE (G) elements.

Length distribution of eccDNAs in different cells. Dotted lines indicate multiplies of 180 bp.

Correlation between the density of small-sized (A) vs large-sized (B) eccDNAs and the meiotic recombination rate [Jensen-Seaman et al., Genome Res (2004) 14, 528-538].

High similarity between sperm eccDNAs detected in this study and those from apoptotic DNA fragmentation reported previously. (A) Enrichment of eccDNAs reported in Nature (2021) 599: 308-314 at given genomic regions relative to randomly-selected control regions. (B) Enrichment of eccDNAs reported in Nature (2021) 599: 308-314 at given repeat elements relative to randomly-selected control regions. (C) EccDNA numbers per 2Mb region in this study as a function of eccDNA numbers per 2 Mb reported in Nature (2021) 599: 308-314.

Evaluation of our nucleotide-resolution eccDNA detection method. (A) A representative genomic locus showing Circle-seq signals and eccDNAs detected by individual methods. Enlarged are the first eccDNA as well as its sequence information. Homologous sequences (TAC) at the eccDNA end are highlighted in red. (B) Schematic representation of our eccDNA detection pipeline. (C) Fraction of eccDNAs detected by individual methods (x-axis) that are co-detected by other methods (y-axis).