A simple ATAC-seq protocol for population epigenetics [ version 2 ; peer review : 2 approved ]

We describe here a protocol for the generation of sequence-ready libraries for population epigenomics studies, and the analysis of alignment results. We show that the protocol can be used to monitor chromatin structure changes in populations when exposed to environmental cues. The protocol is a streamlined version of the Assay for transposase accessible chromatin with high-throughput sequencing (ATAC-seq) that provides a positive display of accessible, presumably euchromatic regions. The protocol is straightforward and can be used with small individuals such as daphnia and schistosome worms, and probably many other biological samples of comparable size (~10,000 cells), and it requires little molecular biology handling expertise. Open Peer Review


Introduction
Understanding the dynamic cross-talk between epigenetic mechanisms and environmental cues of animal populations is of fundamental importance for ecologists and evolutionary biologists (Shi et al., 2019). The dynamics of chromatin has long been of interest as a source of phenotypic variance within and among animal populations (Hu & Barrett, 2017;Zhang et al., 2018) and can affect their ecological performance (Augusto et al., 2019;Hawes et al., 2018). In eukaryotic cells, chromatin is a dynamic structure that provides epigenetic information to control cell and gene function (Chen & Dent, 2014). The physical organization of the chromatin landscape modulates accessibility of genomic regions and dynamically response to both external and internal stimuli. In general, accessible genomic regions are enriched in regulatory elements important for gene activity while inaccessible regions restrict binding of transcriptional regulators resulting in gene silencing (Stergachis et al., 2013). Assay for transposase accessible chromatin with high throughput sequencing (ATAC-seq) is a technique used to assess genome-wide chromatin accessibility. ATAC-seq works similarly as DNase-seq (DNase I hypersensitive sites with high-throughput sequencing) (Song & Crawford, 2010), and determines which genomic regions are accessible to Tn5 transposase (i.e. open chromatin regions), presumably the regulatory regions. Tn5 transposase inserts Illumina adapter sequences upon accessing the chromatin, which removes the need for additional steps to make the sequencing libraries later. This simple and efficient protocol reduces the starting material required, compared to DNase-seq. It also avoids many other steps such as the interaction with antibodies (e.g. ChIP-seq) or chemical treatment (e.g. FAIRE-seq, WGBS) that might introduce bias. Deep sequencing of the PCR amplified Tn5 accessible regions provides a high-resolution map of accessible chromatin regions in the genome. We reasoned that this technique can not only be used to establish functional links between chromatin structure and gene function, but also to quantify epigenetic diversity in populations. This would require generation of ATAC-seq chromatin maps in single individuals. In addition, the technique should be sufficiently robust to be used by scientists who are experts in the field of population (epi)genetics and ecology, but having potentially received little training in molecular biology.
Here, we describe a streamlined and robust method for ATACseq of individuals of the crustacea Daphnia pulex and for the trematode Schistosoma mansoni. Our procedure is based on the protocol from Buenrostro et al. (2015); Corces et al. (2016) and Nextera DNA Library Preparation Kit (2017). Besides their ecological and epidemiological importance, both abovementioned organisms show high phenotypic plasticity in response to environmental cues (e.g. the presence of predator for Daphnia) or during their development (schistosomes). There is a rich literature that has documented the amazing property of Daphnia to modify their phenotypes at the morphological, physiological, behavioral and more recently at the molecular levels in response to a large panel of environmental stressors including diet, pollution, heavy metals, and predator kairomones (reviewed in (Harris et al., 2012;Riessen, 1999)). Schistosomes also deal with a multitude of signals from the water environment as well as cues that come from their hosts, shaping morphology, metabolism, and infection success in the short-term and also their full development later in life (Augusto et al., 2017;Augusto et al., 2019;Roquis et al., 2018). We and others have characterized many aspects of epigenetic mechanisms behind the phenotypic plasticity of schistosomes and their cross-talk with environmental cues. Epigenetics of D. pulex phenotype plasticity is still poorly understood. To validate the robustness of the method, we were experimenters with different levels of expertise in molecular biology to run the experimental procedure independently using D. pulex and S. mansoni. We show here that our procedure provides robust results with individual D. pulex and with single adult worms of S. mansoni, but other organisms of similar cell number can probably also be used.

Materials and methods
Animal sampling A batch of ~300 commercial Daphnia pulex was obtained from a commercial supplier (Aqualiment: http://www.aqualiment.eu/). For schistosomes, fresh adult worms were collected from female Swiss OF1 mice (weight mean 18g) supplied by Charles River, L'Arbresle, France. Mice had been infected by peritoneal injection with 150 mixed sexes cercaria. Water and food were given ad libitum, 12h light/dark cycle, 25°C. Housing, feeding and animal care followed the national ethical standards established in the writ of 1 February 2013 (NOR: AGRG1238753A) setting the conditions for approval, planning and operation of establishments, breeders and suppliers of animals used for scientific purposes and controls. The French Ministère de l'Agriculture et de la Pêche and French Ministère de l'Éducation Nationale de la Recherche et de la Technologie provided permit A66040 to our laboratory for experiments on animals and certificate for animal experimentation (authorization 007083, decree 87-848) for the experimenters. Hepatic perfusions were performed with lethal injection of 1 mg per kg body weight of sodium pertobartial

Amendments from Version 1
In this new version, we have added information that was requested by the reviewers. There are also new authors who contributed to this new, strongly revised version. In response to the reviewer comments we added: -Environmental cues -Counting of nuclei -post-sequencing bioinformatics detection of chromatin structure differences -two population epigenetics pilot studies Daphnia and Schistosoma In particular, we added two examples of pilot population epigenetics studies to show the feasibility of the method from animal sampling to the bioinformatics analysis. However, the emphasis is on the wet-bench part.
Any further responses from the reviewers can be found at the end of the article REVISED solution (Dolehal, Vetoquinol, Lure, France) after 65 days post infection. Living adult male worm was individually transferred to a 1.5 mL Eppendorf tube and immediately processed for ATAC-seq library preparation.

Environmental cues
At their arrival, D. pulex were immediately split into two sets of equal size (~ 150 × 2) and placed in two independent experimental tanks (i.e. initial density of 75 ind/L), hereafter called the 'stress' and the 'control' tanks. Each experimental tank consisted in a 2-L plastic aquaria (L × l × h = 18 ×12 × 11 cm) supplied with clean water, inside of which a floating plastic fish breeding isolation box (L × l × h = 12.5 × 8 ×7 cm) was placed. These isolation boxes are transparent with a series of 1 mm cracks on the bottom wall to allow water connection between the tanks and inside the isolation box.
D. pulex were acclimated in their respective experimental tanks out of the isolation box for 20 days prior to starting the experiment. This lag time before the experiment also allowed the production of new D. pulex offspring born in our experimental setup. During this acclimating period, only negligible mortality was observed and newly hatched D. pulex were observed in the two experimental tanks. After this 20-day acclimating period, a predator (i.e. a guppy fish previously trained to eat D. pulex) was introduced into the isolation box of one experimental tank during 15 days (i.e. hereafter called the 'stress treatment', compared to the 'control treatment') ( Figure 1A). During the experiment the fish was fed every other day with 10 D. pulex collected alternatively from the stress and the control tank (i) to avoid subsequent biases in density between the experimental treatments and (ii) to account for a possible effect of D. pulex sampling on congeners' Figure 1. Schematic representation of experimental design. (A) Daphnia were put into a water tank and allowed to acclimate (start population). Then, two experimental tanks were set up following strictly the same design. The only difference was the presence of a predator (a guppy trained to eat daphnia) in the floating plastic fish breeding isolation box in the stress treatment. (B) S.mansoni infected snails were used to produce cercariae that were separated into two populations, one treated with Latex in well water, the other without treatment. After one hour cercaria were used to infect mice. Adult worms were recovered at 65 days post-infection by perfusion and used for ATAC-seq.
responses. D. pulex sampling for fish feeding was achieved using a sterile 3-ml plastic transfer pipet. This experimental setup allowed the D. pulex of the stress treatment to experience an indirect predation pressure (i.e. without predation risk) through a direct visual contact with the predator and an olfactory contact with environmental cues released by the predator. Over the experiment, the D. pulex and the predator were maintained at room temperature following the natural photoperiod and the former were fed ad libitum with clean phytoplancton (i.e. Chlorella sp.) reared in our lab facilities. Living D. pulex were sampled by pipetting through a 1 mL automatic pipette with enlarged openings of the pipetting tips. To avoid experimenter bias, 13 different persons sampled at least one individual each. Finally, each specimen was then individually transferred to a 1.5 mL Eppendorf tube and was immediately processed for ATAC-seq library preparation as follows. Controls were done without organisms as input.
For schistosomes, treatment of vertebrate infective larvae (cercariae) and recovery of adult worms was done as described in Augusto et al., 2017; briefly, E. milii var. hislopii latex was collected at Ilha do Governador district (22°48´09´´S/ 43°12´35´´W), Rio de Janeiro, Brazil lyophilized at -52°C on 8 ×10 -1 mBar for three 12-hour cycles in a Modulyo 4K Freeze Dryer with an acrylic chamber (Edwards High Vacuum Int., UK). The dose of the powdered lyophilized latex used to expose cercariae was 1.4 mg/L, described by Augusto et al., 2017 as LC 50 for the intermediate host Biomphalaria glabrata. Cercariae were collected from infected host B. glabrata and split into two groups as follows: one group described as 'Control treatment' ( Figure 1B) was kept in a tank of water for one hour while other group described as 'Latex treatment' was exposed to a solution of E. milii lyophilized latex in distilled water (1.4 mg/L), both for one hour. We infected 10 female mice (4 weeks-old Swiss-Webster mice, weight mean 18 g) with 150 exposed cercariae per mouse and another 10 mice were infected with 150 mock treated cercariae per mouse, all using standard percutaneous inoculation and mixed sexes. Finally, parasites couples were recovered at 65 days postinfection by perfusion. (Figure 1B). Males were manually separated from females and only male worms were used.

Counting of nuclei
Individual male S. mansoni worms and individual D. pulex were transferred into 1.5 ml Eppendorf tubes and all excess liquid was removed. Animals were resuspended in 20 -100 µL of Nuclei EZ lysis buffer (Sigma N3408-200ml) and grinded with disposable polypropylene pestles. Two µL were applied on microscope glass slides and 2 µL ProLong Diamond Antifade mountant with DAPI (Invitrogen P36966) were added to stain nuclei. Nuclei in the total volume were counted under a Zeiss fluorescence microscope.

Transposase mixture
The necessary material is listed in Table 1 and Table 2 and must be prepared in advance. In addition, nuclease-free water, high fidelity DNA polymerase for PCR and corresponding buffers, freshly prepared 80% ethanol, refrigerated centrifuge, 0.2 ml PCR tubes, 1.5 ml tubes, ThermoMixer with agitation, PCR thermal cycler, qPCR instrument, magnetic rack, 1 mL pipette, 100 µL pipette, and 10 µL pipette are needed.
An Eppendorf ThermoMixer was then set with agitation to 37°C and the following steps performed. This gives 50 µL of transposase mixture for each sample. The samples are pipetted up and down 10 times to disrupt cells. In our hands, the addition of IGEPAL directly into the tagmentation reaction, and the agitation eliminated the need for a separate cell lysis step and streamlined the protocol. Chromatin tagmentation This step uses the Nextera Tn5 transposome to 'tagment' the chromatin, which is a process that fragments the chromatin and tags the DNA with adapter sequences in a single step.
• Tagmentation reactions are incubated at 37°C for 30 min in an Eppendorf ThermoMixer with agitation at 300 rpm.
• Tagmented chromatin is immediately purified using a QIAGEN MinElute Reaction Cleanup kit or a QIAquick PCR Purification Kit, and purified DNA is eluted into 20 µL of elution buffer (10 mM Tris-HCl, pH 8).

Or
Option 2 (for NEB mix, more convenient but more expensive). Combine the following in a PCR tube for each sample: 20 µL purified transposed DNA; 2.5 µL Universal Ad1_noMX primer In both options the final volume is 50 µL. The samples are pre-amplified using a PCR machine with the program described in Table 4.
In order to reduce GC and size bias in the subsequent PCR, the PCR dynamics is monitored using qPCR to stop amplification prior to saturation. To run a qPCR side reaction, we combined the following depending on the option that had been chosen previously: Option 1: 5 µl PCR product of the initial pre-amplification reaction (keep the remaining 45 µL at 4°C); 2.5 µl 5x GoTaq G2 buffer; 0.1 µL GoTaq 2; 3.14 µl Nuclease-free MilliQ water; The samples are amplified in a qPCR machine with the program set out in Table 5.
To calculate the optimal additional number of cycles needed for the remaining 45 µL PCR, relative fluorescence is plotted against cycle number and the cycle number that corresponds to one-third of the maximum fluorescent intensity is determined ( Figure 2). In our experience, the total number of amplification cycles must not exceed 21 (Augusto et al., 2019).   The remaining 45 µL PCR reaction is run with the additional number of cycles and purified with a QIAGEN MinElute Reaction Cleanup kit or a QIAquick PCR Purification Kit, or similar, and eluted into a total of 45 µL of elution buffer (10 mM Tris-HCl, pH 8). Elution can be done in two rounds.
Fragments are separated by electrophoresis through a 1.5% agarose gel or on a Bioanalyzer chip. A ladder that corresponds to the nucleosome-free region and multiple nucleosome-size fragments should be seen (one nucleosome = about 150 bp).
A single band at around 150 bp indicates sample degradation or over-fragmentation. Ideally, five bands should be obtained, three bands are acceptable ( Figure 3A).

AMPure XP beads double-side purification
This step enriches for the nucleosome-free (~300 bp) as well as di and tri-nucleosome fragments. Removing small fragments (primer dimers) is important for optimal sequencing. First transfer 45 µL to an Eppendorf tube (or use PCR tube directly), add 22.5 µL (0.5X original volume, to remove large fragments) AMPure XP beads, pipet up and down 10 times to mix thoroughly. Incubate at room temperature for 10 minutes and place tubes in magnetic rack for 5 minutes. Transfer supernatant to new tube and add 58.5 µL (1.3X original volume, to remove small fragments) AMPure XP beads, pipet up and down 10 times to mix thoroughly. Incubate at room temperature for 10 minutes, place tubes in magnetic rack for 5 minutes and discard supernatant. Wash beads with 200 µL 80% ethanol (freshly made), pipet ethanol over beads 10 times, then discard ethanol. Ensure all ethanol is removed. Leave tube on magnetic rack with cap open for 3 to maximum 10 minutes depending on ambient humidity. The beads should be 'glowing' but not wet. Be careful not to over-dry them, which will decrease elution efficiency. Resuspend beads in 20 µL nuclease-free water, pipet up and down 10 times to mix thoroughly, place tube in magnetic rack for 1-5 minutes and transfer supernatant to new tube. This step can be replaced by Diagenode IP-Star, size selection 320 bp.
We have not systematically investigated if different purification procedures influence on the result. Purified libraries should be stored at -20°C and can be used for sequencing after up to 4 months.

Libraries check
Size profiling can be performed using an Agilent Bioanalyzer High Sensitivity DNA Assay. Expected profiles are shown in Figure 3B. Bioanalyzer profiles or KAPA library quantification kit were used to quantify libraries and proceed to sequencing. We present here data sequenced on a NextSeq550 High Output Flowcell as paired-end and 75 bp. bowtie2-align-s basic-0 -p 6 -x genome -N 1 -L 32 -i S,1,1.15 --n-ceil L,0,0.15 --dpad 15 --gbar 4 --end-to-end --score-min L,-0.6,-0.6. Uniquely aligned reads were retained by filtering the tag "XS:i:" that is absent in their alignement annotations.  (2) we performed a multivariate HMM over the combined ATAC-seq samples in each condition. For that, BAM files were processed under the differential mode, with a false discovery rate (FDR) cutoff of 0.05 and bin size of 500.

Results
The method can be used by scientists with low expert level in molecular biology The protocol described in the methods section was tested by 13 experimenters with molecular biology expert level ranging from untrained to over several 10 years of experience, or some who had retired from active wet-bench work several years ago. In only two cases ATAC-seq library production did not succeed.

ATAC-Seq can be used on individual Daphnia and individual Schistosoma adults
Our ATAC-seq procedure delivered reproducible chromatin profiles for individual D. pulex and adult S. mansoni. Projection of ATAC-seq reads on a metagene profile indicated that Tn5 accessible and thus presumably open chromatin structure occurs at the TSS and in gene bodies (Figure 4).

Exposure to predator cues leads to morphological differences in Daphnia
Our results show that on average, the (LL-SL)/SL ratio calculated for D. pulex from the stress treatment (N = 14; Mean = 0.24 ± 0.072) was significantly higher than that from the control treatment (N = 12; Mean = 0.15 ± 0.039; Mann-Whitney U Test, U = 19, Z = -3.32, P < 0.001) ( Figure 5). This result confirms the expected induction of anti-predatory morphs in the stress treatment. It is noteworthy that the quantified morphological response to predation pressure observed in the stress treatment most likely reflects a more general response of stressed D. pulex including morphological, physiological and behavioural changes (Boersma et al., 1999) Our first intention in comparing D. pulex from the two experimental  treatments was to confirm that we effectively induced a global response in stressed individuals, these responses having been otherwise much better documented previously (Riessen, 1999).
Exposure to predator cues leads to differences in chromatin structure between exposed (stressed) and unexposed (control) Daphnia pulex Using the DESeq2 procedure described above for 'start' vs. 'control' we identified 66,194 differences between 'control' and 'stressed'. This is by far too many, and indeed, shifts in MA plots (not shown) indicated that the assumption that is underlying the algorythm used in DESeq2 and the requires that most sites do not change, was violated. Metagene profiles, using the same number of aligned reads over the entire genome, lend further support to the finding that 'stressed' samples had on average fewer reads over genes than 'control' samples indicating major changes in chromatin structure ( Figure 6).
This also means that there is a large number of regions for which no reads could be recovered in the stressed samples. This is not due to a general lower accessibility of Tn5 to the cells and nuclei because of a thicker cuticle or a similar phenotypic trait because the insert size distribution of start, control and stressed populations are similar (Supplementary file 1). If DNA was more inaccessible in the stressed population we would expect longer fragments. To cope with the general decrease of ATAC-Seq reads in the stressed population, we resorted to ChromstaR, a HMM based software that was developed for ChIP-Seq analysis but that in principle can also be used for ATAC-Seq and is probably less sensitive to zero values. Under the constraints of numerous instances of an 'absence of data' (Tn5 inaccessible), ChromstaR identified 87 regions that are different between start and control, and stress. All were visually inspected using MACS2 average profiles, normalised by the same number of aligned reads over the genome. Among these 87 regions, ATAC signal was down in stressed samples compared to 'control and start' in 45 regions (52%), down in 'stress and control' compared to 'start' in 16 (18%), up in 'stress and control' in 3 (3.4%), and down in 'control' in only 1 (1.1%). Seven regions showed a heterogenous pattern on ATAC signals. In 15 regions differences were considered too weak (17%) suggesting that fine tuning of ChromstaR parameters might be necessary in the future. These results are in line with a general decrease in ATAC signal in the stressed samples, i.e. chromatin becomes less accessible and/or less heterogenous. It is interesting to note that for 20 regions adjacent ATAC signals (less than 2kb apart) were detected, lending further support to the idea that chromatin structure changes occur in a controlled fashion. Clustering of the samples clearly regroups control and stressed samples (Figure 7).
Exposure to Latex leads to differences in chromatin structure between schistosoma adults that developed from exposed (stressed) and unexposed (control) cercaria With DESeq2 we found 296 differences between schistosoma adults that developed from latex exposed cercaria and controls with adjusted p-value ≤0.05. MA plots were symmetric around 0 ( Figure 8A), and PCA plots ( Figure 8B) indicated clear segregation of both sample groups. We used in this small study only four worms to demonstrate the feasibility of the method.

Discussion
Phenotypically, plasticity plays an important role in development and evolution. The relative contribution of genetic and epigenetic components to heritable plasticity is a matter of lively scientific debate (Hu & Barrett, 2017;Roquis et al., 2018). One of the caveats of analyzing epigenetic information is that it is stored in several, very different bearers of information (e.g. DNA methylation, modification of histones, non-coding RNA and topological position in the interphase nucleus).
Nevertheless, these types of information converge towards a change in chromatin structure which can be approximated by DNA accessibility. We reasoned that a straightforward ATAC-seq method to map the chromatin accessibility  status in populations with high phenotypic plasticity would facilitate further investigations of the role of epigenetics in plasticity. This study field is also of particular importance to field ecologists. We therefore set-out to establish a robust, easy to use protocol that can be used with little molecular biology training. Our protocol was successfully used in the framework of a summer school 'Epigénétique en Ecologie et Evolution' by participants with different levels of expertise in molecular biology using D. pulex. We also used single adult S. mansoni worms as biological material in a small pilot study. We believe that our protocol is suitable for fast epigenotyping of other organisms as well. From our experience, the only parameter that might be necessary to optimize is Tn5 to chromatin ratio if over-or under-fragmentation occurs. A potential issue is contamination with microorganisms whose DNA might be present in the libraries.

Open Peer Review José María Santos-Pereira
Andalusian Center for Development Biology (CABD), Spanish National Research Council, Pablo de Olavide University, Seville, Spain In this protocol article, Augusto and colleagues describe a version of the popular ATAC-seq method for chromatin accessibility profiling in the invertebrate organisms Daphnia pulex and Schistosoma mansoni. Their main claim is that this protocol can be used by researchers with little molecular biology training. However, only the wet-lab component of the ATAC-seq experiment is covered by the manuscript, while the computational analysis of the obtained sequences is omitted. I think that this is a major limitation that has to be addressed before indexing (please, see below).
Major points: Using Tn5 transposase with the appropriate number of cells is critical for this protocol.
Authors mention that TAGmentation time should be reduced in case that libraries are over-TAGmented. However, they do not mention the cell numbers used with these organisms, nor the number of Schistosoma individuals. This is important to be clearly explained in the manuscript.

○
The authors claim that this protocol produces robust results. However, the results are not shown at all. A complete ATAC-seq protocol should cover from the sample preparation in the lab to the analysis of the generated results. I think this is a major limitation of this article and that the authors should explain, at least, the primary computational analysis of the data, the performed quality controls and examples of data visualization. An example of basic analysis of the data, such as peak calling and motif enrichment analysis in the called peaks would be much more helpful to understand how this technique may help to study epigenetic diversity at the population level. Of course, these data have to be available for the community, for example by a GEO accession code.

Minor points:
The origin of the Daphnia samples is not explained.

○
While I assume that the use of 1% Igepal in the TAGmentation reaction is what lyses the samples, this has to be stated more clearly, since usually cell lysis is performed before TAGmentation in other ATAC-seq protocols.

○
The Promega GoTaq G2 polymerase is not a high-fidelity enzyme and therefore I would not recommend its use for this purpose, since this could result in mutations in the amplified molecules and decreased alignment efficiency to the reference genome. Have the authors noted this when analyzing their data in comparison with NEB Next?

○
How was the quality of the results obtained by the 13 volunteers? If this protocol can be performed by researchers with little molecular biology training, then the authors should show that the experiments performed by their volunteers were of enough quality to be used.

○
The agarose gel showed in the supporting material has no lane labels, so it cannot be known what is shown. It would be helpful to show pictures of the agarose gels underlying Fig. 2 data in the main figure. 1-In this protocol article, Augusto and colleagues describe a version of the popular ATACseq method for chromatin accessibility profiling in the invertebrate organisms Daphnia pulex and Schistosoma mansoni. Their main claim is that this protocol can be used by researchers with little molecular biology training. However, only the wet-lab component of the ATAC-seq experiment is covered by the manuscript, while the computational analysis of the obtained sequences is omitted. I think that this is a major limitation that has to be addressed before indexing (please, see below). It was our initial intention to show the robustness of the wet-bench part. The computational analysis of the obtained sequences is now also provided.
Major points: 2-Using Tn5 transposase with the appropriate number of cells is critical for this protocol.
Authors mention that TAGmentation time should be reduced in case that libraries are over-TAGmented. However, they do not mention the cell numbers used with these organisms, nor the number of Schistosoma individuals. This is important to be clearly explained in the manuscript.

A sentence describing the number of Schistosoma individuals was added at the animal sampling and transposase mixture section. Number of nuclei were counted and added to the results section, paragraph "ATAC-Seq can be used on individual Daphnia and individual Schistosoma adults"
3-The authors claim that this protocol produces robust results. However, the results are not shown at all. A complete ATAC-seq protocol should cover from the sample preparation in the lab to the analysis of the generated results. I think this is a major limitation of this article and that the authors should explain, at least, the primary computational analysis of the data, the performed quality controls and examples of data visualization. An example of basic analysis of the data, such as peak calling and motif enrichment analysis in the called peaks would be much more helpful to understand how this technique may help to study epigenetic diversity at the population level. Of course, these data have to be available for the community, for example by a GEO accession code We have added a data analysis section and two examples of population epigenetics study to show the feasibility. The following sentence was added to the Data availability section: This article contains supporting information online on NCBI SRA (BioProject PRJNA587385).

Minor points 4-The origin of the Daphnia samples is not explained
This is now in the Animal sampling part.
5-While I assume that the use of 1% Igepal in the TAGmentation reaction is what lyses the samples, this has to be stated more clearly, since usually cell lysis is performed before TAGmentation in other ATAC-seq protocols.
We understand the reviewer's concern. When AM, RCA and CG optimized the method we realized that adding IGEPAL directly into the tagmentaion reaction did not alter the outcome and that it streamlined the protocol and made it more straightforward.
6-The Promega GoTaq G2 polymerase is not a high-fidelity enzyme and therefore I would not recommend its use for this purpose, since this could result in mutations in the amplified molecules and decreased alignment efficiency to the reference genome. Have the authors noted this when analyzing their data in comparison with NEB Next? We did not observe differences between both polymerases following our pipeline. GoTaq is cheaper and that could be an advantage when large number of libraries must be prepared.
7-How was the quality of the results obtained by the 13 volunteers? If this protocol can be performed by researchers with little molecular biology training, then the authors should show that the experiments performed by their volunteers were of enough quality to be used.
As indicated "The protocol described in the methods section was tested by 13 experimenters with molecular biology expert level ranging from untrained to over several 10 years of experience, or some who had retired from active wet-bench work several years ago. In only two cases ATAC-seq library production did not succeed. " We also added the Figure  The authors describe here a simple ATAC-seq protocol for its use in population epigenetics. ATACseq is now a common molecular biology technique used to assess chromatin accessibility in nuclei of culture cells or tissues. The authors have adapted this technique to small non-model organisms, which could be a way to determine epigenetic variation/signatures possibly associated with phenotypic variance in animal populations. Although this protocol seems really straightforward and useful for the community (even for non-molecular biologist), several points need to be clarified before its indexing.

Major points:
Would be nice to have an estimation of the size and cell number in both S.mansoni and D.pulex specimens. This could be important to further adapt this protocol with other organisms of similar size, smaller or larger.

○
There is no mention of nuclei isolation or tissue treatment before Transposase treatment. Maybe worth to mention if such steps are needed or not with these organisms, and why.

○
There is no mention at all of the sequencing part per se of this ATAC-seq protocol. The protocol closes with "library check". A paragraph should be added to present the sequencing part of this ATAC-seq protocol. In particular, mention if D.pulex and S.mansoni genomes are well annotated and/or mention if alternative procedure can be used for sequence data analysis in these cases. Page 2/ Mat&Met: I would separate "Animal sampling" and "transposase mixture" sections.

○
Page 2/Table2/Reagents produced in the laboratory: I would recommend using commercially available PBS1X, and not sure why another Tagmentation (TD) buffer is indicated here. Usually, the TD buffer is directly provided with the enzyme (as indicated in Table 1). If so, Table 2 is not needed.

○
Page 4/Figure1 could be a bit more explicit. At least indicate on the graph the one-third of the maximum fluorescent intensity to determine the additional PCR cycles.

○
Page 4/AMPure XP beads double-side purification: "This step enriches for the nucleosome-free (~300 bp)" … "regions as well as di and tri-nucleosome fragments". Maybe more correct like this. This specific point could be also mentioned in the Figure 3 legend. ○ Also mentioned: "First transfer 45 ul to an Eppendorf tube (or use PCR tube directly)". After the additional PCR cycles and purification, DNA should be in 20 ul elution buffer. Need to be corrected or clarified for consistency in the following of the protocol.

○
Page 5/Results: Most data are also presented via a zenodo interface. However, the legends and labels of these figures are not always explicit, and the correspondence with the actual figures of the protocol is not always intuitive, something that could be readily improved. , one of large fragments in the range of 10 kb. As they did for Figure 2, the authors should comment on that and propose an alternative solution when this particular situation is encountered, i.e. needs further purification or it is acceptable to process.

Is the description of the method technically sound? Yes
Are sufficient details provided to allow replication of the method development and its use by others? Partly If any results are presented, are all the source data underlying the results available to ensure full reproducibility? Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article? Partly 2-There is no mention of nuclei isolation or tissue treatment before Transposase treatment. Maybe worth to mention if such steps are needed or not with these organisms, and why We appreciate the reviewer's concern. This was one of the streamlining steps and in our hands no nuclei isolation was necessary. We tested both approaches. Maybe the small size of the organisms play a role.
3-There is no mention at all of the sequencing part per se of this ATAC-seq protocol. The protocol closes with "library check". A paragraph should be added to present the sequencing part of this ATAC-seq protocol. In particular, mention if D.pulex and S.mansoni genomes are well annotated and/or mention if alternative procedure can be used for sequence data analysis in these cases. The sequencing and analysis parts were added.
4-What would have been very interesting is to provide an illustration that this technique is indeed suitable to quantify epigenetic diversity in populations. For instance, it would have been nice to have examples of a genome browser display emanating from their multiple ATAC-seq libraries in D.pulex and/or S.mansoni. Maybe, the authors could also refer to some published or ongoing work. We have added the description of the population epigenetics part. We agree that IGV visualisation is definitely useful to get an impression of the data especially in relation to genome annotations. However, we feel that the display is not so useful for the population studies. It was not our intention to use this method to infer any functional relationships between chromatin structure and e.g. gene expression. 7-Page 2/Introduction/end of the paragraph: I would move "Controls must be done without organisms as input" in the Mat&Met/transposase mixture section. This was done.
9-Page 2/Table2/Reagents produced in the laboratory: I would recommend using commercially available PBS1X, and not sure why another Tagmentation (TD) buffer is indicated here. Usually, the TD buffer is directly provided with the enzyme (as indicated in Table 1). If so, Table 2 is not needed. We understand the reviewer concern, however we would like to keep this information on the paper. Our group is running several experiments using ATAC-seq and this is a useful information in the case of missing PBS 1X and/or TD buffer. When we use the Illumina kit in the lab, we can often dilute very much the Tn5 enzyme and what gets