1. Developmental Biology

Genome-wide discovery of cis-regulatory elements in a large genome

  1. Gillian Forbes
  2. Emilia Skafida
  3. Irene Karapidaki
  4. Savannah Moinet
  5. Mowgli Dandamudi
  6. Çağrı Çevrim
  7. Farzaneh Momtazi
  8. Chryssa Anastasiadou
  9. Sabrina Lo Brutto
  10. Michalis Averof  Is a corresponding author
  11. Mathilde Paris  Is a corresponding author
  1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, France
  2. Centre National de la Recherche Scientifique (CNRS), France
  3. Iranian National Institute for Oceanography and Atmospheric Science (INIOAS), Marine Bioscience Department, Iran
  4. Potsdam Institute for Climate Impact Research (PIK), Germany
  5. Fisheries Research Institute, Hellenic Agricultural Organization, Greece
  6. Department of Earth and Marine Sciences (DiSTeM), University of Palermo, Italy
  7. National Biodiversity Future Center (NBFC), Italy
https://doi.org/10.7554/eLife.111378.3
Share this article
Cite this article
  1. Gillian Forbes
  2. Emilia Skafida
  3. Irene Karapidaki
  4. Savannah Moinet
  5. Mowgli Dandamudi
  6. Çağrı Çevrim
  7. Farzaneh Momtazi
  8. Chryssa Anastasiadou
  9. Sabrina Lo Brutto
  10. Michalis Averof
  11. Mathilde Paris
(2026)
Genome-wide discovery of cis-regulatory elements in a large genome
eLife 15:RP111378.
https://doi.org/10.7554/eLife.111378.3
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

This important study combines chromatin accessibility and genomic DNA sequence conservation data from low-coverage genome sequencing of related species (without assembly), for the in silico identification of cis-regulatory elements in large genomes. The approach and results are compelling and well supported by the experimental validations. The work will be of interest to researchers working in the field of gene regulation and evolution, particularly because the methodology proposed can be applied to a large variety of experimental organisms.

https://doi.org/10.7554/eLife.111378.3.sa0

Abstract

Identifying cis-regulatory elements typically relies on trial and error to test the activity of DNA fragments using reporter constructs. This approach is particularly challenging in large eukaryotic genomes, where cis-regulatory elements can be dispersed over long distances, separated by large stretches of non-functional DNA. Here, we generate two types of resources that can be used to narrow the search for such elements in the 3.6 Gbp genome of Parhyale hawaiensis. First, we use bulk ATAC-seq to uncover genome-wide patterns of chromatin accessibility in Parhyale embryonic and adult tissues, and single-nucleus ATAC-seq to identify regions of open chromatin in diverse cell types. Second, by sequencing the genomes of three congeneric species – P. darvishi, P. aquilina, and P. plumicornis – we identify islands of sequence conservation across the genome, likely corresponding to functionally constrained DNA. We find that low-coverage (10–15×) short-read genome sequencing, without genome assembly, is sufficient to provide reliable maps of sequence conservation. This approach cuts the cost and labour required to generate these maps, making the identification of cis-regulatory elements more widely accessible. We demonstrate the utility of these resources by identifying cis-regulatory elements that drive robust expression of fluorescent reporters ubiquitously and in specific cell types.

Introduction

Identifying cis-regulatory elements (CREs) that regulate transcription is key for understanding how gene expression is regulated in space and time, and for targeting the expression of genetic markers and tools to specific populations of cells. The classic approach for identifying CREs is based on trial and error: DNA fragments flanking the gene of interest are selected and their activity is tested using reporter constructs (e.g. Goto et al., 1989; Arnold et al., 2013; Pfeiffer et al., 2008). In organisms with compact genomes, CREs can often be found in proximity to promoters. In large eukaryotic genomes, however, CREs can be located tens or even hundreds of kilobases away from coding sequences (e.g. Sanyal et al., 2012; Levo et al., 2022). In the latter case, the task of identifying active cis elements within a sea of non-functional DNA is laborious and inefficient. Large research communities working in established model organisms (e.g. flies and mice) or humans have historically invested substantial resources in identifying diverse CREs to support tool development (e.g. Arnold et al., 2013; Pfeiffer et al., 2008), but smaller communities working on non-conventional model organisms usually lack these resources (e.g. Sun et al., 2022; Lai et al., 2018).

A number of methods can be used to narrow down the candidate fragments of DNA to be tested. One approach is based on chromatin profiling. Transcriptional enhancers tend to be located in regions of ‘open’ chromatin, allowing access to transcription factors, they are often enriched in specific histone marks (e.g. H3K4me1 and H3K27ac) and maintain contacts with their target promoters. These features can be revealed by experimental techniques such as ATAC-seq, ChIP-seq, and Hi-C, respectively (Buenrostro et al., 2013; Lieberman-Aiden et al., 2009; Barski et al., 2007). Given that each genomic locus is found in just two copies per diploid nucleus, a major constraint in using these methods is the number of cells required to obtain a robust signal. ChIP-seq and Hi-C typically require tens of thousands to millions of nuclei per experiment (Kidder et al., 2011; Belton et al., 2012). ATAC-seq requires lower amounts of starting material – hundreds to tens of thousands of nuclei per experiment (Buenrostro et al., 2013) – and has recently even been used in the context of single-cell profiling (Buenrostro et al., 2015). Chromatin accessibility profiles can therefore be obtained in the majority of experimental systems via ATAC-seq.

A second approach for identifying candidate CREs is to look for islands of DNA that show a high degree of sequence conservation during evolution, an approach known as phylogenetic footprinting (e.g. Duret and Bucher, 1997; Blanchette and Tompa, 2002; Pennacchio et al., 2006; Woolfe et al., 2005). The slower rate of evolution, resulting from functional constraints, helps to distinguish CREs from surrounding non-coding sequences that have no function. Comparisons between species that are separated by tens (and sometimes hundreds) of millions of years of evolution are effective in differentiating putative CREs from their surrounding non-functional sequences. Sequence conservation is often lost entirely over longer evolutionary timescales. To perform these comparisons, species at the appropriate evolutionary distances must therefore be collected and sequenced. A common limitation of this approach is the labour and cost involved in genome sequencing and assembly, particularly in species with very large genomes. We propose a workaround to this problem in this paper.

A third approach for identifying putative CREs has emerged recently with the development of deep learning models capable of predicting enhancer sequences (de Almeida et al., 2022; Barbadilla-Martínez et al., 2025). To generate reliable predictions, these models must be trained on large experimental datasets, which have so far been available only for Drosophila and mammals. It is not yet clear whether these models can perform well in other, distant species.

Here, we present resources for the discovery of CREs in the crustacean Parhyale hawaiensis, an emerging experimental model in comparative developmental biology, regeneration, chronobiology and ecotoxicology (reviewed in Paris et al., 2022; Averof, 2022). Parhyale presents the typical challenges for CRE discovery in organisms with large genomes, having a genome of approximately 3.6 billion base pairs, similar in size to the human genome (Kao et al., 2016). Previous efforts to identify CREs relied on testing the activity of DNA fragments lying within a few kb from promoters (Pavlopoulos and Averof, 2005; Pavlopoulos et al., 2009; Sun et al., 2022; Ramos et al., 2019). These efforts succeeded in identifying CREs from hsp70 and opsin genes, but failed to identify additional CREs, particularly ones associated with development or cell differentiation.

We combined two orthogonal approaches – chromatin profiling and sequence conservation – to identify putative CREs in the Parhyale genome. For this purpose, we produced the following resources: (1) bulk ATAC-seq data from embryonic and adult legs of Parhyale, (2) single-nucleus ATAC-seq data from adult legs of Parhyale, which yield the chromatin profiles of >15 cell clusters, including epidermis, muscles, neurons and blood, and (3) genome-wide maps of sequence conservation, highlighting islands of conservation in relation to other Parhyale species.

We show that it is not necessary to sequence genomes at high coverage, nor to generate genome assemblies, in order to identify islands of sequence conservation. Instead, we present a method that relies on mapping short sequence reads from new species (in this case, from P. darvishi, P. aquilina, and P. plumicornis) with low stringency onto the genome assembly of the species of interest (P. hawaiensis; Kao et al., 2016). This approach cuts the cost and labour required to generate such comparative data very significantly.

To demonstrate the utility of these resources, we have used them to discover both ubiquitous and cell type-specific CREs in the genome of Parhyale. First, we used the single-nucleus ATAC-seq data to identify regions of open chromatin that are shared by all cells, or regions that are specific to particular cell types such as neurons or muscles. Among these regions, we selected ones that show sequence conservation among related Parhyale species. Putative CREs were then tested using transgenic reporters, identifying two CREs that drive ubiquitous expression (out of two tested), two that drive expression specifically in neurons (out of seven tested), and two that drive expression in muscles (out of two tested). Besides the value of these resources for the Parhyale research community, the labour- and cost-efficient approach we present here can be used to generate resources for CRE discovery in other species with large genomes.

Results and discussion

Genome-wide profiles of chromatin accessibility in embryonic and adult tissues

To explore chromatin accessibility profiles across the Parhyale genome, we performed bulk ATAC-seq on nuclei isolated from three types of tissue: (1) single Parhyale embryos collected at stages S20, S23, and S24, covering embryonic stages when many organ systems are being formed (datasets E20, E23, and E24, respectively), (2) developing legs collected at embryonic stages S25 to 26 (datasets EL1 to EL3), and (3) fully differentiated adult T4 and T5 thoracic legs (datasets L1 to L3). These datasets identify tens of thousands of peaks of open chromatin distributed across 2.8 Gbp of the assembled Parhyale genome (results summarised in Table 1). These ATAC-seq data show a relatively low abundance of accessible chromatin in exons, introns, and intergenic regions, and a strong enrichment of accessible chromatin surrounding transcription start sites (TSS; Figure 1B and D). The quality of ATAC-seq data is often measured by the fraction of reads that are found in peaks of accessible chromatin (FRiP); based on this criterion, our highest quality datasets are EL1 to EL3, L1, and E24 (see Table 1).

Figure 1
Download asset Open asset
Comparison of ATAC-seq data from embryonic and adult tissues of Parhyale hawaiensis.

(A) Genome browser plot of ATAC-seq data, focused on the 3.5 Mb spanning Hox gene cluster of Parhyale hawaiensis. In whole embryos (E20-24), ATAC-seq peaks are distributed across the entire Hox cluster. In embryonic thoracic legs (EL1-3), ATAC-seq peaks in the region of anterior and posterior Hox genes are partly suppressed. In adult T4 and T5 legs (L1-3), peaks of open chromatin are mostly visible in the region of Ubx and suppressed in other parts of the Hox cluster. (B) Histogram showing the proportion of transcription start sites (TSS), exonic, intronic, and intergenic sequences in the Parhyale genome that overlap with an ATAC-seq peak in each dataset. Only TSS and exons of intron-containing genes were used in this analysis, to exclude poorly annotated transcripts. (C) Principal component analysis of the ATAC-seq datasets, based on the number of reads assigned to each peak. PC1 and PC2 (accounting together for 44% of the variation) show the datasets clustering according to tissue of origin (whole embryos, embryonic legs, adult legs, marked in different colours). (D) Mapping of ATAC-seq reads to the sequences surrounding transcription start sites across the genome; 1 kb upstream to 1 kb downstream of each TSS is depicted in each line (half of the 27,955 TSS are shown, keeping the same order per column; see ‘Materials and methods’). Colours represent the number of mapped reads (see ‘Materials and methods’). We observe a clear enrichment of open chromatin surrounding the TSS.

Table 1
Overview of bulk and single-nuclei ATAC-seq datasets from Parhyale hawaiensis.

Further information on mapping statistics is given in Tables 5 and 6.

DatasetTotal mapped reads after filteringFraction of reads in peaks (FRiP)Number of peaks
Whole embryoE2016,431,0700.1260,435
E238,726,1740.1556,606
E247,879,1660.29119,994
Embryonic legsEL119,186,7580.3787,308
EL218,000,3840.3677,357
EL323,424,3880.42119,487
Adult legs (bulk)L120,233,1840.3276,268
L216,508,4560.1949,502
L310,206,7440.2248,324
Adult legs (single-nuclei)SN1319,174,0930.41404,519
SN2217,929,6070.43331,355
Naked DNAnDNA102,914,3080.15161,533

Principal component analysis, based on the number of reads assigned to each peak, shows these datasets clustering according to the origin of the samples (whole embryos, embryonic legs, adult legs), with the tightest clustering seen in the embryonic leg datasets (Figure 1C).

Some of the variation between samples is likely to reflect differences in chromatin accessibility between tissues and developmental stages. To test this, we focused on the genomic region that encompasses the Parhyale Hox genes. Based on the mechanisms of Hox gene regulation known to operate in Drosophila (Bowman et al., 2014; Peifer et al., 1987) and the expression of Hox genes in Parhyale (Serano et al., 2016), we expected open chromatin peaks to be distributed across the entire Hox cluster in embryos (E20, E23, and E24), but only the region surrounding the gene Ubx in adult T4 and T5 thoracic legs (L1-3). The patterns we observe partly match that expectation (Figure 1A).

Genome-wide profiles of chromatin accessibility from diverse cell types

Many research projects on development and regeneration focus on gene regulation and expression profiles that are specific to particular cell types. We have previously shown that Parhyale legs contain more than 15 transcriptionally distinct cell types, including epidermal and neuronal cells, muscle, blood, and several yet unidentified cell types (Almazán et al., 2022). To probe the cellular specificity of chromatin landscapes, we performed single-nucleus ATAC-seq (snATAC-seq) on adult Parhyale T4 and T5 legs (datasets SN1 and SN2, comprising a total of 15,969 nuclei).

By pooling the ATAC-seq reads from all the cells captured in these datasets, we generated ‘pseudo-bulk’ SN1 and SN2 datasets which are comparable with the bulk ATAC-seq datasets (Table 1, Figure 1A and B). Principal component analysis shows that these datasets cluster with the bulk ATAC-seq datasets generated from adult legs (Figure 1C).

The snATAC-seq data identify clusters of cells with distinct chromatin landscapes (Figure 2A and B). Based on the ATAC-seq peaks surrounding cell type-specific marker genes identified by snRNA-seq (Almazán et al., 2022), we are able to identify some of these cell clusters as epidermal, muscle, neuronal, and blood cells (Figure 2C and D; see ‘Materials and methods’). The cells of each cluster share ATAC-seq peaks that are specific to the cluster (e.g. Figure 2E). Thus, our snATAC-seq data identify genomic sequences that are differentially accessed in different cell types.

Figure 2 with 1 supplement see all
Download asset Open asset
Single nuclei ATAC-seq and differential chromatin accessibility across cell types.

(A–C) UMAP of integrated snATAC-seq experiments SN1 and SN2, including 15,969 cells from adult uninjured Parhyale legs (7951 cells from SN1 and 8018 cells from SN2). UMAP was colour-coded based on experiment (A), cell clusters identified from the ATAC-seq signal (B), or cell types identified from previously published snRNA-seq data (Almazán et al., 2022; see ‘Materials and methods’, Figure 2—figure supplement 1) (C). (D) UMAP of snRNA-seq data published in Almazán et al., 2022, using the same colour code. (E) Dot plots of the ATAC-seq signal at the putative CREs that we tested in vivo, per cell cluster.

Parhyale species available for comparative genomics

Our initial efforts to identify regions of sequence conservation in the Parhyale genome focused on comparisons with Hyalella azteca, a hyalid amphipod whose genome was sequenced in 2018 (Poynton et al., 2018). Comparing non-coding sequences between Parhyale and Hyalella, however, failed to reveal islands of sequence conservation, even in loci known to harbour functional CREs (see Appendix 1—figure 1). Phylogeographic studies estimate that the lineages of Parhyale and Hyalella likely diverged 110–140 million years ago (Cannizzaro and Berg, 2022), likely too long for detecting sequence conservation in CREs.

To identify regions of sequence conservation over shorter evolutionary distances, we collected three species that belong to the same genus as P. hawaiensis – P. aquilina, P. darvishi, and P. plumicornis – and we sequenced their genomes to an estimated genome coverage of 16×, 11×, and 10×, respectively (Table 2, Figure 3—figure supplement 1; these values may represent slight overestimates, see ‘Materials and methods’). Based on our analysis, the genome sizes of P. aquilina, P. darvishi, and P. plumicornis are estimated to be 1.1 Gbp, 2.8 Gbp, and 3.0 Gbp, respectively.

Table 2
Overview of short-read genome sequencing on three Parhyale species.
SpeciesEstimated genome size (Gbp)Sequenced readsGenome coverage% reads mapping to P. hawaiensisNon-exonic regions conserved in P. hawaiensis
P. aquilina1.1281,429,56316×6.4605,996
P. darvishi2.8318,968,91611×4.1525,822
P. plumicornis3.0424,305,64510×0.64156,063

The reads from each species were then mapped with low stringency to a genome assembly of P. hawaiensis (Phaw_5.0, GenBank GCA_001587735.2; updated in Almazán et al., 2022). Using this approach, 6.4% of reads from P. aquilina, 4.1% of reads from P. darvishi, and 0.46% of reads from P. plumicornis could be mapped unambiguously to the P. hawaiensis genome (Table 2), suggesting that sequence comparisons within the genus Parhyale could be informative for CRE discovery in these animals. These levels of cross-species mapping likely reflect the relative phylogenetic proximity and genome sizes of these species. A molecular phylogeny based on the sequences of conserved single-copy genes (BUSCO gene set, 29,097 aligned nucleotides across all species; see ‘Materials and methods’) reveals the phylogenetic relationships of these species (Figure 3A).

Figure 3 with 1 supplement see all
Download asset Open asset
Cross-species sequence conservation and its relation with chromatin.

(A) Molecular phylogeny depicting the evolutionary relationships of the Parhyale species included in this study (P. hawaiensis, P. aquilina, P. darvishi, and P. plumicornis), with Hyalella azteca as the outgroup. The phylogeny was based on 29,097 aligned nucleotides from conserved single-copy genes (BUSCO gene set; see ‘Materials and methods’). Divergence estimates (in million years) were calibrated against the estimated divergence between Parhyale and Hyalella (Cannizzaro and Berg, 2022). Scale bar represents number of substitutions per site. (B) Overlap between regions of sequence conservation in non-exonic sequences (TSS, introns and intergenic regions) identified by mapping of reads from P. aquilina, P. darvishi, and P. plumicornis onto the genome of P. hawaiensis (see ‘Materials and methods’). Values in parentheses indicate the number of overlapping regions that would be expected if the conserved regions were distributed randomly. (C) Quantification of the fraction of ATAC-seq peaks in non-exonic sequences that overlap a conserved region identified in P. aquilina, P. darvishi, and P. plumicornis (see ‘Materials and methods’). (D) Quantification of the fraction of ATAC-seq peaks in non-exonic sequences that overlap a conserved region in either one, two, or three of the other Parhyale species (see ‘Materials and methods’). Colour code corresponds to panel (B). Black dots in panels (C) and (D) represent the overlap that would be expected if the conserved regions were distributed randomly.

An estimate of divergence times based on the relative branch lengths in this tree, calibrated against H. azteca (Cannizzaro and Berg, 2022), suggests that P. hawaiensis has likely diverged for 20–25 million years from P. aquilina and P. darvishi, and 75–95 million years from P. plumicornis.

Genome-wide patterns of sequence conservation

The cross-species read mapping identifies hundreds of thousands of islands of sequence conservation in the genome of P. hawaiensis (Table 2, Figure 3B). Of these, 8% are found in exons, likely reflecting conservation in protein coding sequences, while 37% are found in introns, 54% in intergenic regions, and 1% overlap with promoters (TSS), marking regions that evolve at a lower rate than surrounding non-coding sequences. The peaks in non-exonic sequences include 54,950 that show conservation across all four species. If these peaks of sequence conservation were distributed randomly in the genome, only 949 would be expected to overlap between all four species (Figure 3B). This suggests that a large proportion of these sequences is conserved due to constraints that are maintained over a scale of 100 million years of divergent evolution.

The islands of sequence conservation found in non-coding sequences (promoter/TSS regions, introns, and intergenic regions) coincide with peaks of open chromatin identified by ATAC-seq more frequently than would be expected by chance (Figure 3C and D), consistent with the hypothesis that many of these regions represent cis-regulatory elements. The binding of regulatory proteins and RNA to these regions likely explains their chromatin accessibility and constrains their rate of sequence evolution.

Identification of cis-regulatory elements driving ubiquitous, neuron- and muscle-specific expression

In the past, most attempts to identify ubiquitous or cell type-specific CREs in Parhyale using DNA fragments encompassing putative promoters were unsuccessful. For example, in a recent effort to identify neuron-specific CREs, we generated reporter constructs carrying putative promoter regions from genes expressed in neurons (bruchpilot, VAChT, synaptotagmin 4) or involved in neuronal differentiation (achaete/scute, senseless, cut, prospero, neuralised). The 10 fragments we tested extended 1.5–5.5 kb upstream of the start codon, including a 3'UTR, putative TSS, and 1–5 kb of upstream sequences (see ‘Materials and methods’). None of these reporters was found to be expressed in embryos. In retrospect, it appears that some of these fragments lacked the ATAC-seq peak that is usually associated with the TSS, suggesting that the promoter region had been incorrectly identified in the absence of ATAC-seq data.

To test the utility of the ATAC-seq and sequence conservation data we generated, we used these resources to identify CREs with ubiquitous and cell type-specific activities. To identify ubiquitous promoters, we searched for genes that are robustly expressed in all cell clusters in adult Parhyale legs (Almazán et al., 2022); we then identified those with promoters that exhibit a clear ATAC-seq peak in all cell types and overlap with sequences that are conserved in at least one other Parhyale species. Based on this approach, we identified two promoter regions and an associated putative CRE belonging to the Parhyale orthologues of the genes headcase (hdc, referred to as ubi1) and muscleblind (mbl, referred to as ubi2) (Figures 2E and 4A and B, and Table 3). To test the activity of these sequences, we generated reporter constructs carrying these fragments upstream of the coding sequence of mNeonGreen and cloned these into the Minos transgenesis vector. We microinjected these constructs in early Parhyale embryos and we screened the resulting transgenic mosaic embryos by fluorescence microscopy (see ‘Materials and methods’). We found that both ubi1 and ubi2 give strong widespread expression in at least 30% of the embryos and hatchlings we screened (Figure 4A', A'', B' and B'').

Figure 4
Download asset Open asset
Identification of CREs driving ubiquitous expression.

(A, B) Genome tracks showing the ATAC-seq and sequence conservation profiles that we used to identify the ubi1 and ubi2 promoters and associated CREs. (A', B') Fluorescence observed with the ubi1.P (A') and ubi2.CRE+P (B') reporters in live late stage embryos. Image (A') was captured on a fluorescence stereoscope, image (B') is a max projection of an image stack captured by confocal microscopy. Side views of embryos, with the head located at the top right, dorsal side towards the top-left (see Browne et al., 2005). Arrowheads mark the eye; asterisks in (B') mark the base of the antennae. (A'', B'') Fluorescence observed with the ubi1.P (A'') and ubi2.CRE+P (B'') reporters in live hatchlings. These images show ventro-lateral views of unilateral genetic mosaics, with expression only on one half of the animal captured on a fluorescence stereoscope. The head is located at the top, dorsal side on the left. Arrowheads mark the eye. Scale bars, 50 µm. The raw image data are available in Suppl. Data 6, available at https://doi.org/10.5281/zenodo.19020963.

Table 3
Putative CREs with ubiquitous, neuron- and muscle-specific activity.

The transgenic reporters that we tested carry the promoter region (P) of the gene of interest and additional putative CREs upstream of the EGFP or mNeonGreen coding sequences. The fragments tested are shown in Figures 4A and B and 5A–C. The number of embryos screened, the number of embryos showing unilateral or bilateral expression of the Opsin1 transgenesis marker, and the number of embryos showing (mosaic) reporter expression are indicated. Genetic mosaicism results in different numbers of positives for the Opsin1 marker (eyes) and reporter expression (other tissues). CNS, central nervous system. The P and CRE sequences are given in Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963.

ActivityGene name, aliasPutative CRELength (bp)Embryos screenedOpsin1- positiveReporter expression
Ubiquitoushdc (MSTRG.41953)
ubi1		P102712338All tissues
(38 embryos)
mbl (MSTRG.26075)
ubi2		CRE +P501+890260-All tissues
(118 embryos)
Neuronsfutsch (MSTRG.441)
neuro1		CRE +P615+49614710No expression
jeb (MSTRG.26302)
neuro2		CRE +P976+23314825No expression
ChT (mikado.phaw_ 50.283866G186)
neuro3		P114023241No expression
αTub (MSTRG.41851)
neuro5		P954542166CNS, peripheral neurons
(108 embryos)
Cdk5α (MSTRG.7309)
neuro6		P89620890CNS, peripheral neurons
(55 embryos)
CRE +P375+89639267CNS, peripheral neurons
(81 embryos)
mav (MSTRG.827)
neuro7		CRE +P1871+230308No expression
MSTRG.27182
neuro8		CRE +P358+75122088No expression
MusclesMhc (MSTRG.35247)
Mhc		P100372-Muscles
(9 embryos)
CRE1+P735+1,003485Muscles
(6 embryos)
CRE6+P1194+1,0039212Muscles
(22 embryos)

To identify neuron-specific CREs, we used the snATAC-seq data to identify 37 regions of the genome that are preferentially accessible in neurons (relative to other cell types) and located near genes that are preferentially expressed in neurons based on snRNA-seq data (Almazán et al., 2022). Among these regions, we selected ones that show sequence conservation in at least one other Parhyale species. Based on this analysis, we selected candidate promoters and associated putative CREs from seven genes (Table 3, Figure 5A and B, Figure 5—figure supplement 1). To test their activities, we generated reporter constructs carrying these promoters and CREs upstream of the coding sequence of EGFP or mNeonGreen, cloned these into the Minos transgenesis vector (see ‘Materials and methods’) and microinjected each construct in early Parhyale embryos. The resulting transgenic mosaics were screened at late embryonic stages. Two of the seven reporters we tested, named neuro5 and neuro6 (carrying CREs associated with the genes αTub and Cdk5α, respectively), showed strong fluorescence in the brain, ventral nerve cord, and peripheral neurons in the body, legs, and antennae (Figure 5A', A'', B' and B''). We observed these patterns of fluorescence in approximately 20% of the surviving injected embryos. The neuro6 promoter alone gave expression in the nervous system, suggesting that at least part of the neuro6 reporter activity resides in the promoter. Transgenic animals carrying both neuro5 and neuro6 reporters show that these two CREs drive expression in distinct subsets of neurons (Figure 5—figure supplement 2). The embryos injected with the other five constructs showed no significant fluorescence above background levels.

Figure 5 with 2 supplements see all
Download asset Open asset
Identification of CREs driving expression in neurons and muscles.

(A–C) Genome tracks showing the ATAC-seq and sequence conservation profiles that we used to identify the neuro5, neuro6, and Mhc neuron- and muscle-specific promoters and associated CREs. (A'–C') Fluorescence observed with the neuro5-Src64B-mNeonGreen (A'), neuro6-Src64B-mNeonGreen (B'), and Mhc.CRE6+P-mNeonGreen (C') reporters in late embryos. Max projections of image stacks captured by confocal microscopy. Side views of embryos, with the head located at the top right, dorsal side towards the top-left. Arrowheads mark the eye; asterisks in (C') mark the base of thoracic legs. (A''–C'') Fluorescence observed with the of the neuro5.P (A''), neuro6.P+CRE (B''), and Mhc.CRE6+P (C'') reporters in the thoracic legs of juveniles (proximal part of the leg including coxal plate in A'', distal tip of leg in B'', proximal part of several legs in C''). All panels show max projections of image stacks captured by confocal microscopy. Asterisks in (C'') mark the base of thoracic legs. Scale bars, 50 µm. The raw image data are available in Suppl. Data 6, available at https://doi.org/10.5281/zenodo.19020963.

To identify muscle-specific CREs, we focused on the genomic region surrounding the Myosin heavy chain gene (Mhc, MSTRG.35247), which is strongly and specifically expressed in nuclei of the muscle cluster of Parhyale (Almazán et al., 2022). We identified regions in which chromatin is specifically accessible in the muscle cluster and show sequence conservation across Parhyale species (Figure 5C) and tested the activity of two of these regions using the reporter approach described earlier. In both reporters, we observed specific expression in the muscles in 12–23% of the surviving injected embryos (Figure 5C' and C''). Reporters with the Mhc promoter alone gave similar expression, suggesting that at least part of this activity resides in the promoter.

Evidence from diverse sources suggests there are different classes of cis-regulatory elements associated with different types of genes, namely with ubiquitously expressed 'housekeeping' genes, genes expressed in terminally differentiated cells and developmentally regulated genes. Each class tends to be associated with different sequence motifs and exhibits differences in nucleosome positioning, precision of transcription initiation and the use of distantly located enhancers (reviewed by Haberle and Lenhard, 2016). To test whether our approach for identifying CREs can be effective in identifying the enhancers of developmental genes, we used the ATAC-seq data from embryonic legs (EL1-3) combined with cross-species sequence conservation to identify putative enhancers in three leg patterning genes, Dll-e, dac1, and dac2 (Bruce and Patel, 2020; Kao et al., 2016; Pavlopoulos et al., 2009). We tested the activity of 13 distant CREs using fluorescent reporter constructs, carrying each CRE with the promoter region of the corresponding gene (Table 4 and Figure 6). Among these constructs, we could only detect weak ubiquitous expression driven by promoter-proximal sequences of Dll-e (Table 4). We did not find constructs that recapitulate the endogenous pattern of Dll-e, dac1, and dac2 expression in the legs of Parhyale. This result may reflect either our inability to identify the enhancers regulating these developmental genes or expression levels that are below the detection threshold of our fluorescent reporters.

Figure 6
Download asset Open asset
Putative CREs of Parhyale developmental genes tested using transgenic reporters.

Genome browser plots for Dll-e, dac1, and dac2, highlighting the promoter regions (P) and putative CREs that were tested (in grey). ATAC-seq, nDNA, and sequence conservation (P. darvishi, P. aquilina, and P. plumicornis) tracks were normalised by fragments per million and autoscaled to the highest peak in each region. For dac2, the y-axis maximum for whole embryo, embryo, and adult leg tracks was halved to account for the large TSS peak, which otherwise masks the signal at putative CREs. Regions annotated as Ns are indicated by the N stretches track.

Table 4
Putative CREs of Parhyale developmental genes tested using transgenic reporters.

The transgenic reporters tested carried the promoter region (P) of the gene of interest and additional putative CREs upstream of the mNeonGreen coding sequence. The fragments tested are shown in Figure 6. In 6–10% of the embryos screened for Dll-e activity, we observed expression patterns that were not reproducible; we interpret these as enhancer traps associated with specific sites of transgene insertion, rather than the activity of sequences carried in the reporter. The P and CRE sequences are given in Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963.

GenePutative CRELength (bp)Embryos screenedOpsin1- positiveReporter expression
Dll-e
(MSTRG.27394)		P100424345Weak ubiquitous after stage S20
Pshort74138874Weak ubiquitous after stage S20
CRE3+Pshort1209+741453107Weak ubiquitous after stage S20
CRE6+P1517+100448063Weak ubiquitous after stage S20
dac1
(MSTRG.33907)		P944607166No expression
CRE1+P2032+94416737No expression
CRE2+P819+94429247No expression
dac2
(MSTRG.33908)		P1070386112No expression
Pshort444698185No expression
CRE1+Pshort1033+444456162No expression
CRE2+P1048+1070319102No expression
CRE3+P987+10708028No expression
CRE4+P1099+10709926No expression
CRE5+P585+1070346116No expression
CRE8+P957+107010330No expression
CRE9+P916+10704714No expression
CRE10+P2358+107014528No expression
CRE11+P2973+10709217No expression

Conclusions

Previous efforts to identify cis-regulatory elements in Parhyale relied on reporter constructs carrying a few kb of sequence upstream of selected target genes, an approach that has worked well in animals and plants with relatively small genomes. As presented earlier, however, this approach was often unsuccessful in Parhyale: from tens of reporters tested, only four robust native drivers had so far been identified (Pavlopoulos and Averof, 2005; Pavlopoulos et al., 2009; Sun et al., 2022; Ramos et al., 2019). The present work adds five new drivers to that collection, including ones with ubiquitous, neuron- and muscle-specific activities. This result comes from combining information on genome-wide chromatin accessibility and evolutionary conservation profiles.

We cannot at this point distinguish the relative contributions of chromatin profiling and sequence conservation to CRE discovery, because these sources of information were not tested separately. At minimum, we can state that (1) ATAC-seq profiles serve to identify robustly the promoters of candidate genes that are active in a particular cellular context (cell type and stage) and (2) coupling this information with sequence conservation narrows down candidate promoters and distant CREs by a factor of 4–10, since only a fraction of ATAC-seq peaks show sequence conservation (Figure 3D). This represents a great improvement in our ability to select candidates to test by transgenesis, the most labour-intensive step in the process.

Our approach relies on ATAC-seq and low-coverage short-read genome sequencing, which are today widely accessible and affordable. Both methods can be applied to small samples: a few thousand cells can be sufficient for ATAC-seq (as in our single-embryo samples) and a few micrograms of DNA are sufficient for low-coverage short read sequencing of most genomes. These methods can be combined with existing genome annotations and information on gene expression or function to select relevant candidate genes. There are therefore few barriers to adopting this approach.

For ATAC-seq, good quality nuclear preparations are necessary, which could limit the adoption of this method in cases where cells are embedded within a tough extracellular matrix or a highly resistant cell wall. To detect sequence conservation by low-coverage genome sequencing, the phylogenetic distance between species is a key factor. The approach is most powerful when evolution has had sufficient time to erase sequence similarities in non-functional (neutrally evolving) sequences, but not to eliminate conservation within functionally constrained sequences. This degree of sequence divergence is usually found in species separated by tens of millions of years of evolution. In Parhyale, we have found that sequence comparisons between P. hawaiensis, P. aquilina, P. darvishi, and to some extent P. plumicornis reveal islands of sequence conservation in non-exonic sequences that cannot be found in comparisons with more distant relatives such as Hyalella (e.g. Appendix 1—figures 1 and 2).

Our strategy of mapping regions of sequence conservation by direct mapping of short sequence reads across species is much more accessible than conventional strategies that rely on genome assembly. The latter require much higher sequence coverage (>50×) from multiple libraries, long-read sequencing or other scaffolding methods, and complex bioinformatic pipelines to assemble large genomes. Moreover, these approaches are often compromised by high levels of polymorphism found in natural populations. We estimate that our approach is 5- to 10-fold cheaper than assembly based methods, even excluding labour costs.

In spite of our success in identifying CREs that are active ubiquitously or in differentiated cell types, we have not yet demonstrated that this approach can efficiently identify the CREs of developmentally regulated genes. The reasons for this failure could be trivial, such as an unfortunate selection of genes (or CREs) expressed at low levels. But there could also be substantive reasons, such as the very long range dispersal of the CREs of developmental genes, extending beyond the tens of kilobases that we scanned in this study, the presence of multiple shadow enhancers, each contributing a low level of activity in an additive manner, or the rapid evolution of these enhancers leading to low sequence conservation. Further studies could address this question by bringing in additional information, for example, by looking at the long-range contacts of promoters with other parts of the genome using Hi-C (Lieberman-Aiden et al., 2009).

Despite these caveats, the resources we present here allow us to focus our search for CREs in the vast genome of Parhyale to a small subset of DNA fragments that are both accessible and functionally constrained, thus drastically reducing the search space. This approach has led us to discover robust ubiquitous, neuron- and muscle-specific regulatory elements, which will be valuable for building new genetic and imaging tools in this species. The same resources can now be used to identify CREs that are active in a range of other cell types in Parhyale. And the same approach should be applicable, with modest investment, to most organisms.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Biological sample (Parhyale hawaiensis)Parhyale hawaiensisKao et al., 2016NCBI species ID 317513Chicago-F line, laboratory culture
Biological sample (Parhyale aquilina)Parhyale aquilinaThis paperNCBI species ID 1799758Collected in Nea Peramos, Greece
Biological sample (Parhyale darvishi)Parhyale darvishiThis paperNCBI species ID 3695942Collected in Chabahar Bay, Iran
Biological sample (Parhyale plumicornis)Parhyale plumicornisThis paperNCBI species ID 1799757Collected in Milos, Greece

Parhyale hawaiensis culture and handling

Request a detailed protocol

We used the Chicago-F inbred line of P. hawaiensis (Dana, 1853) in order to minimise polymorphisms and maximise mapping to the genome that has been sequenced from that line (Kao et al., 2016). The animals were cultured and embryos were collected in artificial sea water (specific gravity 1.022) using standard methods (Browne et al., 2005; Paris et al., 2022).

Bulk ATAC-seq on whole embryos, embryonic, and adult legs

Request a detailed protocol

The ATAC-seq protocols were based on Buenrostro et al., 2013 (for embryos) and (Corces et al., 2017) (for embryonic and adult legs), with elements from the 10x Genomics Multiome protocol.

For whole embryo datasets E20, E23, and E24, single embryos at stages S20, S23, and S24 (Browne et al., 2005), respectively, were collected and placed in the lid of an eppendorf tube, in 50 µl ice-cold lysis buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% Igepal CA-630). The egg membranes were removed, the cells were dissociated by gentle pipetting, and collected by centrifuging at 800 rcf for 10 minutes, at 4°C. The supernatant was discarded, the pellet was placed on ice and resuspended in 50 µl tagmentation mix, containing 25 μl 2x TD Buffer (Illumina #FC-121-1030) and 1.5 μl Tn5 Transposase (Illumina #FC-121-1030). Tagmentation was performed at 37°C for 30 minutes, shaking at 300 RPM.

For embryonic leg datasets EL1, EL2, and EL3, developing embryonic T4 to T8 legs were collected from pools of stage S25-26 embryos, when white eyes are just starting to be visible (Browne et al., 2005). The embryos were dissected in artificial sea water on a cold block; egg membranes were removed and 30–40 legs were collected in ice-cold buffer consisting of 10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, and 1% bovine serum albumin (BSA). The legs were then placed in 500 µl ice-cold lysis buffer consisting of 10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% Igepal CA-630 (Merck #I3021), 1 mM DTT (Roche #10197777001), 0.01% Tween-20 (Sigma #P1379), and 0.001% digitonin (Promega #G944A), and the leg tissue was mechanically dissociated by five passes through a 25Gx5/8” needle using a 2.5 ml syringe. The lysate was incubated for 5 minutes on ice and then cuticle debris was removed by filtering with a 10 μm pluriStrainer cell strainer (pluriSelect #43-10010). The nuclei were collected by centrifuging at 4°C for 5 minutes at 500 rcf, washed in 400 µl ice-cold wash buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 1% BSA, 0.1% Tween-20, 1 mM DTT), and collected again by centrifuging at 4°C for 5 minutes at 500 rcf. 18,000 to 25,000 nuclei were collected per sample, measured by staining a fraction of the suspension with trypan blue and counting the stained nuclei on a hemocytometer. Tagmentation was performed in 50 µl tagmentation mix, containing 25 µl 2X tagmentation buffer (Diagenode #C01019043) and 2.5 µl Tn5 loaded transposase (Diagenode #C01070012), at 37°C for 30 minutes, shaking at 300 RPM.

For adult leg datasets L1, L2, and L3, uninjured T4 and T5 legs were collected by anaesthetising adult animals in 0.02% clove oil in artificial seawater and cutting the legs at the proximal end of the merus (sample L1 and L2) or the basis (sample L3); 28 legs were collected for L1 and L3, 12 legs were collected for L2. Nuclei were isolated on a cold block, by opening each leg longitudinally and scraping the tissue from the cuticle using a surgical knife (Fine Science Tools #10316-14), followed by pipetting 25 times through a P200 micropipette tip in 250 µl of a lysis buffer consisting of 10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 1% BSA, 0.1% Igepal CA-630, 0.1% Tween-20, and 0.01% digitonin (for samples L1 and L3), or by five passes through a 25Gx5/8” needle using a 2.5 ml syringe (for sample L2). An equal volume of wash buffer (10 mM Tris-HCl pH 7.4; 10 mM NaCl; 3 mM MgCl2; 1% BSA; 0.1% Tween-20) was then added and cuticle pieces were removed by passing the lysate through a 10 μm pluriStrainer cell strainer (pluriSelect #43-10010). Nuclei were collected by centrifuging at 4°C for 5 minutes at 500 rcf (samples L1 and L2) or 10 minutes at 600 rcf centrifugation (sample L3). 50–90 thousand nuclei were recovered per sample, measured on a hemocytometer as described earlier. Tagmentation was performed in 50 µl tagmentation mix, containing tagmentation buffer (Illumina #15027866 for sample L1 and L3, or Diagenode #C01019043 for sample L2) and 2.5 µl Tn5 loaded transposase (Illumina #15027865 for samples L1 and L3, or Diagenode #C01070012 for sample L2) at 37°C for 30 minutes, shaking at 300 RPM.

In parallel, we generated ATAC-seq data on genomic DNA from P. hawaiensis that had been purified to remove histones and other proteins (‘naked DNA’), to reveal intrinsic biases associated with tagmentation of DNA sequence rather than chromatin. The genomic DNA was purified as described for P. aquilina and plumicornis (below). Tagmentation was performed using the Nextera XT DNA Library kit (Illumina #FC-131-1024) according to the manufacturer’s recommendations and then sequenced on a NextSeq500 (PE 2x76).

For all samples, the tagmented DNA was purified using the QIAGEN MinElute kit (QIAGEN #28204) and amplified by a total of 9–17 PCR cycles using Nextera XT Index kit v2 primers (five initial cycles followed by qPCR reaction to determine the total number of PCR cycles needed for each library). After PCR amplification, the libraries were purified by AMPure XP beads (Beckman Coulter) and were resuspended in 20 μl nuclease-free water. Sequencing was carried out on an Illumina NextSeq 500 sequencer; we carried out single-end 76 bp sequencing for the first sample we generated (E24), and then switched to paired-end 76 bp sequencing for the other samples, because this leads to more specific read mapping. The raw data (sequencing reads) of all the bulk ATAC-seq experiments are available at NCBI GEO accession GSE325175.

Reads were mapped to the P. hawaiensis genome (Phaw_5.0, GenBank GCA_001587735.2; updated in Almazán et al., 2022) using bowtie2 version 2.5.4 (Langmead and Salzberg, 2012), using very-sensitive-local parameters (Table 5). Mapped reads were filtered to remove mitochondrial reads and low-quality mapped reads (MAPQ ≥20), using SAMtools version 1.21 (Danecek et al., 2021) and BEDtools version 2.31.1 (Quinlan and Hall, 2010), and duplicates using Picard MarkDuplicates version 2.27.5 (https://broadinstitute.github.io/picard/). ATAC-seq peaks were then identified using MACS2 version 2.2.7.1 (Zhang et al., 2008).

Table 5
Read mapping of ATAC-seq data to the Parhyale hawaiensis genome.

Statistics of ATAC-seq mapping using Bowtie2 to the P. hawaiensis genome, and subsequent filtering to remove mitochondrial reads, unpaired reads (except E24, which was sequenced single end), low mapping quality reads, and PCR duplicated reads. Paired-end reads were counted as two reads.

Sample typeDatasetSequencing output - number of readsNumber of reads mapped
(%)		Number of reads after removing mitochondrial reads (mitochondrial reads removed)Number of reads after selection of properly paired reads
(non-paired reads removed)		Number of reads with mapping quality MAPQ => 20 (multimapping reads removed)Number of unique reads
(duplicate reads removed)
Whole embryoE2049,478,72642,725,420 (86.4%)42,006,520 (718,900 removed, 1.7%)40,592,846
(1,413,674 removed, 3.4%)		20,700,434 (19,892,412 removed, 49%)16,431,070 (4,269,364 removed, 20.6 %)
E2320,106,90618,672,300 (92.9%)17,169,710 (1,502,590 removed, 8%)16,543,142 (626,568 removed, 3.6%)9,007,870 (7,535,272 removed, 45.5%)8,726,174 (281,696 removed, 3.1 %)
E2425,718,86218,755,694 (72.9%)18,344,294 (411,400 removed, 2.2 %)NA (single-end sequencing)10,330,352 (8,013,942 removed, 43.7%)7,879,166 (2,451,186 removed, 23.7%)
Embryo legsEL145,551,82239,820,151 (87.4%)39,764,350 (55,801 removed, 0.1%)38,266,738 (1,497,612 removed, 3.8%)22,644,786 (15,621,952 removed, 40.8%)19,186,758 (3,458,028 removed, 15.3%)
EL244,018,06633,413,764 (75.9%)33,360,865 (52,899 removed, 0.2%)32,183,656 (1,177,209 removed, 3.5%)19,641,888 (12,541,768 removed, 39%)18,000,384 (1,641,504 removed, 8.4%)
EL359,349,07249,263,308 (83%)49,198,262 (65,046 removed, 0.1%)47,578,690 (1,619,572 removed, 3.3%)28,970,962 (18,607,728 removed, 39.1%)23,424,388 (5,546,574 removed, 19.1%)
Adult legs
(bulk)		L152,289,20836,305,142 (69.4%)35,870,798 (434,344 removed, 1.2%)34,455,566 (1,415,232 removed, 3.9%)21,191,302 (13,264,264 removed, 38.5%)20,233,184 (958,118 removed, 4.5%)
L244,194,07437,508,542 (84.9%)37,231,495 (277,047 removed, 0.7%)35,668,376 (1,563,119 removed, 4.2%)19,080,056 (16,588,320 removed, 46.5%)16,508,456 (2,571,600 removed, 13.5%)
L339,923,47419,802,704 (49.6%)19,414,574 (388,130 removed, 2%)18,575,758 (838,816 removed, 4.3%)10,637,496 (7,938,262 removed, 42.7%)10,206,744 (430,752 removed, 4%)
Adult legs (single- nuclei)SN11,351,387,9101,006,599,327 (74.5%)1,003,753,519 (2,845,808 removed, 0.3%)955,229,436 (48,524,083 removed, 4.8%)705,294,901 (249,934,535 removed, 26.2%)319,174,093 (386,120,808 removed, 54.7%)
SN2809,094,758649,723,680 (80.3%)648,556,600 (1,167,080 removed, 0.2%)615,591,128 (32,965,472 removed, 5.1%)454,077,969 (161,513,159 removed, 26.2%)217,929,607 (236,148,362 removed, 52%)
Naked DNAnDNA246,377,216236,589,795 (96%)236,456,206 (133,589 removed, 0.1%)215,062,750 (21,393,456 removed, 9%)118,945,820 (96,116,930 removed, 44.7%)102,914,308 (16,031,512 removed, 13.5%)

Single-nucleus ATAC-seq on adult legs

Request a detailed protocol

Single-nucleus ATAC-seq was performed using the 10x Genomics single-cell capturing system. Libraries were prepared using the Chromium Next GEM Single Cell ATAC Reagent Kits v1.1 according to the manufacturer’s instructions (User Guide Rev. C). Nuclear suspensions from adult Parhyale legs were prepared as described above for sample L3. After lysis and filtering, nuclei were collected by centrifuging at 4°C for 10 minutes at 600 rcf. Then, the nuclear pellets were resuspended in 7 μl of 1X Nuclei Isolation Buffer provided by the kit. Nuclei suspensions containing 15,000 nuclei were used for the tagmentation reactions, which were then loaded on the 10x Genomics Chromium Controller. Libraries were sequenced on Illumina NextSeq500 (72 bp for R1 and R2, according to 10x Genomics User Guide recommendations). SN1 data were generated with two sequencing rounds (one NextSeq Mid Output and one NextSeq High Output), resulting in a total of 675 million reads. For SN2, data were generated with one NextSeq High Output sequencing round, obtaining a total of 404 million reads. Sequencing data were processed using the Cell Ranger ATAC software (v1.2.0, 10x Genomics, Pleasanton, CA) to perform alignment, barcode filtering, and peak calling following the standard pipeline. The raw data (sequencing reads) of the snATAC-seq experiments are available at NCBI GEO accession GSE325176. Reads were mapped to the same genome assembly as the bulk ATAC-seq data (Table 5).

Analysis of snATAC-seq data

Request a detailed protocol

Analyses were performed using the Signac v1.16.0 and Seurat v5.4.0 R packages (Stuart et al., 2021; Hao et al., 2024). Normalisation and linear dimension reduction were performed using reciprocal latent semantic indexing (LSI) (Cusanovich et al., 2015), with the Seurat functions RunTFIDF, FindTopFeatures (with the min.cutoff parameter set to ‘q0’), and RunSVD, using default parameters unless otherwise specified. Datasets were integrated using reciprocal LSI projection (Stuart et al., 2019). Anchor sets were identified using the first 30 reduction axes, except the first one that strongly correlates with sequencing coverage. UMAP and cell clustering were done on the integrated dataset using (LSI) reduction (again using the 2:30 first axes of reduction). The two snATAC-seq datasets (SN1 and SN2) were processed in the same fashion.

Differentially accessible (DA) peaks were identified using the Seurat function FindAllMarkers with the parameter min.pct=0.1. We noticed that all DA regions from cluster 19 were also open in cluster 7 and most DA regions from clusters 4, 5, and 7 were also open in cluster 19. We therefore suspect that cluster 19 may contain barcodes with doublets from clusters 4, 5, and/or 7; since it contains only 95 barcodes, we decided to remove it from the dataset in subsequent analyses.

Cell type labels on the snATAC-seq clusters were obtained by comparing snATAC-seq signal in gene bodies and promoter regions with gene expression in snRNA-seq (Almazán et al., 2022), using the Seurat function TransferData. Cell labels were transferred with good confidence (prediction score above 80%) for several cell types including neurons, muscles, or blood (Figure 2—figure supplement 1).

ATAC-seq data analysis and visualisation

Request a detailed protocol

Genome browser plots for bulk ATAC-seq, snATAC-seq and naked DNA were made with pyGenomeTracks (version 3.9; Lopez-Delisle et al., 2021) and normalised by counts per million. Cell type-specific snATAC-seq tracks were made by combining the data of the two snATAC-seq experiments. The genome browser tracks of bulk and single nucleus ATAC-seq experiments are provided in Suppl. Data 1 and 2, available at https://doi.org/10.5281/zenodo.19020963.

To assess the distribution of ATAC-seq peaks across genomic features (Figure 1B, Table 6), genomic regions were classified as transcription start sites (TSS), exons, introns and intergenic regions. TSS was defined as the first base of annotated transcripts that contain at least one intron. Intronless transcripts were omitted to exclude poorly annotated transcripts from the analysis. Exonic regions were extracted from the same set of transcripts, excluding the TSS coordinate. Intronic regions were defined as the remaining regions of these transcripts that did not overlap with either the exons or TSS. Lastly, intergenic regions were defined as all remaining genomic regions, not assigned to TSS, exon, or intron. Unknown/unresolved nucleotides (Ns) in the genome and intron-less transcript regions were excluded from all feature sets.

Table 6
Number of nucleotides in ATAC-seq peaks in TSS, exon, intron, and intergenic sequences.

We determined the nucleotides of ATAC-seq peaks belonging to each of these genomic features: transcription start sites (TSS), exons, introns, and intergenic regions. Counts on TSS and exons were also determined for intron-containing genes only. Undetermined nucleotides (Ns) in the genome assembly were excluded from the counts.

Sample typeDatasetAll TSS
(56,307 nucl)		TSS in intron- containing genes
(27,955 nucl)		All exons
(129,980,737 nucl)		Exons in intron- containing genes
(94,685,505 nucl)		Introns
(673,377,575 nucl)		Intergenic
(1,390,661,853 nucl)
Whole embryoE20282021423,434,3181,551,8954,184,80610,599,571
E23315826692,368,1081,256,9333,680,69910,594,791
E24584250014,480,0143,142,23710,417,66722,338,694
Embryo legsEL1516342604,608,4932,682,9779,419,74215,523,581
EL2516144213,808,5322,242,1788,236,83913,758,216
EL3583248655,816,6243,698,36912,442,01120,580,140
Adult legsL1619353604,212,1822,470,1459,391,66615,678,750
L2406633013,265,8581,521,8174,817,6879,726,381
L3434338212,558,4981,378,5144,231,1078,936,821
Naked DNAnDNA5551249612,311,3726,304,43033,198,61770,905,248

ATAC-seq peak coordinates were determined from the MACS2 narrowPeak files for each dataset. A consolidated list of ATAC-seq peaks was generated by merging the peaks from all bulk ATAC-seq dataset using the Bedtools merge function version 2.31.1 (Quinlan and Hall, 2010). Fragment (paired-read) counts per peak were quantified using featureCounts version 2.0.6 (Liao et al., 2014). Quantification was performed separately on the single-end sample (E24) and the paired-end sequenced samples (all other samples), and the resulting count matrices were combined. This consolidated peak set was used for assessing the overlap between ATAC-seq and cross-species islands of sequence conservation (Figure 3C–D).

The merged count matrix was normalised using variance stabilising transformation (VST) from DESeq2 version 1.48.2 (Love et al., 2014). Principal component analysis (Figure 1C) was performed on the VST-transformed counts using prcomp (stats version 4.5.2) with scaling enabled, and visualised using ggplot2 version 4.0.1 (Wickham, 2016).

To generate the TSS enrichment heatmaps (Figure 1D), we normalised the ATAC-seq signal across datasets, the MACS2 bedgraph files representing the fragment coverage across the genome were scaled by fragments per million. For each dataset, the total number of fragments after mapping and filtering was divided by one million and the bedgraph signal values were divided by this scaling factor. Scaled bedgraph files were converted to bigwig format using bedGraphToBigWig version 377 (Kent et al., 2010). A signal matrix was generated using computeMatrix from deeptools version 3.5.6 (Ramírez et al., 2016) in reference-point mode, centred on TSS, with a bin size of 10 bp and ±1 kb flanking regions. Heatmaps were visualised using deeptools plotHeatmap, with the colour scale capped at the 98th percentile across all datasets. TSSs were sorted by mean signal intensity for each row in dataset EL3, which was the dataset with the highest FRiP score.

Field collection and identification of Parhyale species

Request a detailed protocol

P. plumicornis (Heller, 1866) were collected on the island of Milos, Greece, in rocky shores near Klima (coordinates 36.736594°N 24.417607°E, and N 36.732361°N 24.422056°E). Live adults were examined under a Nikon field microscope for preliminary identification based on the presence of the characteristic dense arrays of plumose setae on the ventral side of the second antennae, notably on articles 4 and 5 of the peduncle and on the first articles of the flagellum. Identification was then corroborated by additional characters scored on adults cultured in the laboratory, including a smooth body, length of second antenna about 2/3 of body length, inflated peduncle of the second antennae, subchelate first and second gnathopods, robust and spinose pereopods, first uropod carrying a robust seta on the distal-lateral side of the peduncle, well-defined inner ramus of the third uropod, which is not fused to the peduncle, and deeply cleft telson.

P. aquilina (Costa, 1857) were collected in shallow water, in a well-protected bay near Nea Peramos (Kavala), in Northern Greece (coordinates 40.825880°N 24.309377°E). Preliminary identification of adult males was based on the characteristic stout and inflated shape of the dactylus of the first pair of gnathopods, observed on a Nikon SMZ800N stereoscope. Females were selected based on the observation that they formed mating couples with those males. They could be identified more definitively as P. aquilina after genome sequencing (see below), based on sequence similarity with a previously sequenced gene fragment from this species (Iaciofano and Lo Brutto, 2017).

P. darvishi (Momtazi and Maghsoudlou, 2016) were collected by hand between cobbles in the intertidal zone of Chabahar Bay (25.352533°N 60.598494°E), Gulf of Oman, placed in 99% ethanol and stored at –20°C. Preliminary identification of the adult males was conducted on a Nikon stereomicroscope based on the following diagnostic features: first antenna lacking dense setae, a stout second gnathopod, and an anterodistal tuft of setae on the maxillipedal palp. The identification was confirmed microscopically by comparing diagnostic characters with the original species description (Momtazi and Maghsoudlou, 2016), including the absence of setae on the posterior margin of the propodus of pereopods 6 and 7, the presence of five apical robust setae on the outer ramus of uropod 3, and two robust setae on the palm of the first male gnathopod, with the proximal seta being distinctly incurved. Female specimens were selected based on their association with males in the field and confirmed in the laboratory by examining their diagnostic features under a microscope: a transverse palm defined by two stout setae, a medial process on the posterior margin of the palm of the second gnathopod, and oostegites bearing curl-tipped setae.

Low-coverage genome sequencing and cross-species read mapping

Request a detailed protocol

Genomic DNA was extracted from single individuals that were flash frozen in liquid nitrogen, using the NucleoSpin Tissue kit (Macherey Nagel #740952) according to the manufacturer’s recommendations. Whole genome libraries were built for each species using the Nextera XT DNA Library kit (Illumina #FC-131-1024) according to the manufacturer’s recommendations and then sequenced on a NextSeq500 (PE 2x76). The raw data (sequencing reads) obtained from P. aquilina, P. darvishi, and P. plumicornis are available at NCBI BioSample accessions SAMN56533921, SAMN56533922 and SAMN56533923, respectively.

We noticed that one sequencing round on P. plumicornis included reads corresponding to a plasmid that was used in our laboratory, which includes three fragments from the genome of P. hawaiensis (the PhOpsin1 promoter and two fragments from the Dll-e locus), representing a contamination of our library. We removed the contaminating reads by mapping the sequenced reads onto the plasmid sequence, using bowtie with the parameters "--end-to-end --very-sensitive". We discarded the reads with matches to the plasmid sequence, except ones that had two or more mismatches to the plasmid (found in ≥2 reads). Overall a total of 37,385 reads were discarded from the sequencing output file, leaving 424,305,645 reads.

To estimate genome coverage, the reads sequenced from P. aquilina, P. darvishi, and P. plumicornis were used to generate a genome assembly for each species using ABySS-pe v4.2.1 (Jackman et al., 2017), with parameters kc = 2 k=25 n=2. The resulting assemblies were highly fragmented (N50<1 kb), as expected considering the low sequencing coverage. Sequences corresponding to the crustacean BUSCO gene set (version odb12) were identified in these assemblies using BUSCO v6.0.0 (parameter --metaeuk) (Tegenfeldt et al., 2025). To estimate sequencing coverage for each species, we mapped the reads back to the genome assembly using bowtie2 v2.3.5.1 (Langmead and Salzberg, 2012) with stringent parameters (--end-to-end --very-sensitive --score-min L,−0.2,–0.2). We discarded reads that harboured over two mismatches or a mapping quality below 30, and duplicate reads. We then computed the per base coverage on the BUSCO genes (distributions shown in Figure 3—figure supplement 1). Estimates of genome size were obtained by dividing the number of sequenced bases by the mean coverage. The estimated coverage was 16×, 11×, and 10×, corresponding to estimated genome sizes of 1.1 Gbp, 2.8 Gbp, and 3.0 Gbp for P. aquilina, P. darvishi, and P. plumicornis, respectively. As a control, we followed the same procedure for P. hawaiensis, for which genome size is known (Parchem et al., 2010; Kao et al., 2016), and found a genome size of 3.0 Gbp instead of 3.6 Gbp (with a genome coverage of 5.8×). We cannot determine whether our predicted genome size is more accurate than the published estimate or whether our method slightly overestimates coverage (and thus underestimates genome size). In any case, these results suggest that the genome sizes of P. darvishi and P. plumicornis are similar to that of P. hawaiensis, but the genome of P. aquilina is about three times smaller.

Reads from all species were mapped to a genome assembly of P. hawaiensis (Phaw_5.0, GenBank GCA_001587735.2; updated in Almazán et al., 2022), using bowtie2 v2.3.5.1 with the parameters "--very-sensitive-local --local –3 10", and low-quality mapping were removed using the function “view” from samtools v1.21 with the parameter '-q 10". Conservation tracks were produced using the function pileup from macs2 with the parameter "--extsize 75".

Genome browser plots for sequence conservation were generated with pyGenomeTracks (version 3.9; Lopez-Delisle et al., 2021) and normalised by counts per million, as previously described for ATAC-seq. The genome browser tracks showing the cross-species sequence conservation are provided in Suppl. Data 3, available at https://doi.org/10.5281/zenodo.19020963.

Islands of sequence conservation were determined for each species as non-exonic genomic regions (TSS, introns and intergenic regions) with at least one mapped read; their lengths varied from 23 to 5771 nucleotides, with a median of 59–69 nucleotides. We considered that conserved sequence islands overlap between species if their coordinates on the P. hawaiensis genome are shared over at least 50% of the length of the shortest island (Figure 2B). We followed the same logic for assessing the overlap of ATAC-seq peaks with conserved regions (Figure 2C and D). For the latter analysis we used the consolidated list of ATAC-seq peaks (see above).

Phylogeny and estimate of divergence times

Request a detailed protocol

The molecular phylogeny of the sequenced Parhyale species (Figure 3A) was generated based on the sequences of conserved single-copy BUSCO genes. We ran BUSCO with the same parameters as described earlier on three additional species with available genome assemblies: P. hawaiensis (Phaw_5.0, GenBank GCA_001587735.2; updated in Almazán et al., 2022), H. azteca (NCBI RefSeq assembly GCF_000764305.2, Poynton et al., 2018), and Daphnia pulex (NCBI RefSeq assembly GCF_021134715.1). We selected the 23 BUSCO genes that were found in all six species, aligned each gene set separately using macse v2.07 that preserves coding frames (Ranwez et al., 2018), and concatenated the gene alignments using AMAS (Borowiec, 2016). We selected the 29,097 sites that were represented in all four Parhyale species and in H. azteca. The concatenated alignment is available in Suppl. Data 4, available at https://doi.org/10.5281/zenodo.19020963. The phylogeny was generated by the Maximum Likelihood approach using iqtree v3.0.1 (Wong et al., 2025) with the parameters -m MFP -bb 1000 -nt AUTO. Divergence times were estimated using iqtree v2.4.0 on the output tree and the same input alignment, calibrated against an estimated divergence of 110–140 million years between P. hawaiensis and H. azteca (Cannizzaro and Berg, 2022).

Identification of CREs using transgenic reporters

Request a detailed protocol

Our initial efforts to identify neuron-specific CREs (prior to availability of ATAC-seq and sequence conservation data) focused on candidate genes expressed in neurons or involved in neurogenesis. The reporter constructs included 1.5–5.5 kb of sequence immediately upstream of the start codon of each gene: reporters with 1.5 and 3.6 kb upstream sequence for bruchpilot (mikado.phaw_50.283862bG132); 2.3 and 4 kb for VAChT (MSTRG.41555); 2 kb for synaptotagmin 4 (mikado.phaw_50.283817fG572); 4.2 kb for achaete/scute (MSTRG.45479); 4.8 kb for senseless (MSTRG.40789); 3.9 kb for cut (MSTRG.41143); 5.5 kb for prospero (MSTRG.45422); 5 kb for neuralised (MSTRG.5439). The sequences of these fragments are available in Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963. Based on our transcriptome assemblies (Kao et al., 2016; Sinigaglia et al., 2022), these fragments include the 5'UTR, the 5' end of the annotated transcript (putative TSS) and 1–5 kb of sequences further upstream. These putative CREs were amplified by PCR on Parhyale genomic DNA and cloned immediately upstream of the coding sequence of EGFP and the SV40 early transcriptional terminator, and placed in a Minos transposon vector (Pavlopoulos and Averof, 2005) that also carries PhOpsin1-EGFP as a transgenesis marker (Ramos et al., 2019). These reporter constructs were co-injected with Minos transposase mRNA into early Parhyale embryos, as described previously (Pavlopoulos and Averof, 2005; Kontarakis and Pavlopoulos, 2014); typically, the injection mix contained 250 ng/µl donor plasmid and 200 ng/µl Minos transposase mRNA. The injected embryos were raised in filtered artificial sea water containing 100 units/mL penicillin, 100 μg/ml streptomycin and 0.25 μg/ml amphotericin B (Gibco #15240096). The surviving embryos were screened for fluorescence on a Leica MZ16F stereoscope during mid and late embryonic stages. At least 50 embryos expressing the PhOpsin1-EGFP marker were scored per construct, but no additional expression was detected in the nervous system.

Candidate ubiquitous CREs were identified by combining information from snRNA-seq, snATAC-seq, and sequence conservation. First, using the integrated snATAC-seq dataset, we selected the peaks that are accessible in at least 20% of the cells within every cell cluster and do not show differential accessibility across clusters (Seurat function FindAllMarkers p_val_adj >0.5). Second, among those peaks, we selected the ones that are nearest a gene (using the Signac function ClosestFeature) that is expressed in at least 20% of the cells within each cell cluster, and is not differentially expressed across clusters (Seurat function FindAllMarkers p_val_adj >0.5, or pct.1>0.5 and pct.2>0.5 and abs(pct.1-pct.2)<0.1), based on the snRNA-seq data generated by Almazán et al., 2022. Third, among the resulting peaks, we selected those that included regions of cross-species sequence conservation. Among the resulting peaks, we focused on the orthologues of Drosophila genes headcase/hdc (MSTRG.41953) and muscleblind/mbl (MSTRG.26075) (Figure 4A and B). To the mbl promoter region, we added an ATAC-seq peak with conserved sequences that is located within an intron >50 kb downstream of the promoter (Figure 4B). These candidate CREs were placed upstream of the coding sequence of the mNeonGreen fluorescent protein fused with histone H2B to drive nuclear localisation (Wolff et al., 2018), followed by the P10 transcriptional terminator (Pfeiffer et al., 2012), designated as H2B-mNeonGreen-P10, by gene synthesis (Twist Bioscience, USA), in the context of a Minos vector carrying the PhOpsin1-EGFP transgenesis marker (described above). The mbl promoter fragment includes a splice donor sequence and the 5' end of the gene’s first intron, which carries ATG triplets that could be misread as start codons. To ensure that this region is spliced out, we added the 3' end of an intron and splice acceptor from the Parhyale Distal locus (Kontarakis et al., 2011), prior to the mNeonGreen coding sequence. The annotated sequences of both CREs are available in Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963. These reporter constructs were injected into Parhyale embryos as described above, and screened in mid and late embryonic stages and in juveniles using a Zeiss Axiozoom V16 stereoscope. The image shown in Figure 4B' was captured on a Zeiss LSM 800 confocal microscope equipped with a Plan-Apochromat 20×/0.8 M27 objective (Zeiss 420650-9901-000).

Candidate neuron-specific CREs were identified by combining information from RNA-seq, ATAC-seq, and sequence conservation. First, we selected the peaks in the integrated snATAC-seq dataset that are enriched in the cells predicted to be neurons compared to the cells predicted to be non-neuronal (with a confidence score ≥0.8, based on the label-transfer approach described earlier), using the Seurat function FindAllMarkers (test.use = 'LR',min.pct=0.2) with the following filtering parameters: corresponding ATAC-seq reads should be found in at least 20% of the predicted neuronal cells, and in at least 10% more cells than in the other cell types combined, the log2 fold change enrichment should be over 3 and the adjusted p-value below 10^(–100). We only kept the peaks that were nearest to a gene labeled as neuronal-specific, using our snRNA-seq data (Almazán et al., 2022) and the ClosestFeature function from Signac and FindMarkers from Seurat. A list of 37 peaks near 33 genes was manually inspected on the IGV genome browser using our bulk and snATAC-seq and sequence conservation data and reduced to a set of 7 genes, including 7 promoter regions and 5 additional associated CREs. Promoter regions included the sequence from the start of the ATAC-seq peak up to the nucleotide before the start codon; associated CREs included the sequence of the entire ATAC-seq peak. Reporter constructs were generated by placing the candidate CRE and promoter sequences upstream of the coding sequences of beta galactosidase and EGFP, separated by the T2A ribosome skipping peptide, and the SV40 early transcriptional terminator. These reporters were placed in a Minos vector carrying the PhOpsin1-EGFP transgenesis marker (as described above). The reporter constructs were generated by gene synthesis (Twist Bioscience) and conventional subcloning. The annotated sequences are available in the Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963. The activity of these reporters was tested by microinjecting Parhyale embryos as described above. Approximately 400–700 embryos were injected per construct. The surviving embryos were screened for fluorescence on Leica MZ16F and Zeiss AxioZoom V16 fluorescence stereoscopes during stages S20 to S28. For each construct, we found that 20–30% of injected embryos surviving to late embryonic stages expressed the PhOpsin1-EGFP transgenesis marker, either unilaterally or bilaterally (10–80 positive embryos per construct). In addition to this marker, approximately 20% of the embryos injected with the neuro5 and neuro6 constructs showed strong fluorescence in the brain and ventral nerve cord.

A second set of constructs were made in which the neuro5 and neuro6 regulatory elements are used to express membrane-localised fluorescent proteins (using the Src64B and HRas tags; Karapidaki et al., 2026), to label neuronal projections. This includes two constructs carrying neuro5-Src64B-mNeonGreen and neuro6-Src64B-mNeonGreen (followed by the P10 transcriptional terminator; Pfeiffer et al., 2012), and a construct carrying both neuro6-mScarlet3-HRas (followed by the SV40 early transcriptional terminator) and neuro5-Src64B-mNeonGreen (followed by the P10 transcriptional terminator). The first two were cloned in a Minos vector carrying the PhOpsin1-EGFP transgenesis marker (see above); the latter was cloned in a Minos vector carrying no additional marker. These constructs were microinjected and screened as described above. Mosaic G0 embryos were imaged on a Zeiss LSM 800 confocal microscope equipped with a Plan-Apochromat 20×/0.8 M27 objective (Zeiss 420650-9901-000).

To identify muscle-specific transgenic line, we focused on the Myosin heavy chain (Mhc) gene of Parhyale (MSTRG.35247), which is specifically and strongly expressed in all muscle cells in the snRNA-seq data of Almazán et al., 2022. We identified the promoter of Mhc and two putative CREs that are accessible specifically in muscle nuclei and contain regions that are conserved across Parhyale species. The sequences of these putative CREs are available in Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963. These sequences were placed upstream of the coding sequence of mNeonGreen followed by the P10 transcriptional terminator, by gene synthesis (Twist Bioscience), in the context of a Minos vector carrying the PhOpsin1-EGFP transgenesis marker (as described above). The activity of these reporters was tested by microinjection into Parhyale embryos, as described above. Late embryos (from S26) were screened for both Mhc and PhOpsin1 activity using a Leica MZ16F fluorescence stereoscope.

To identify CREs regulating the activity of leg patterning genes, we explored ATAC-seq peaks with sequence conservation around the genes dll-e, dac1 and dac2, corresponding to genes MSTRG.27394, MSTRG.33907, and MSTRG.33908 in the genome annotation of Sinigaglia et al., 2022. Putative CREs and promoter regions, identified by the TSS peak, were placed upstream of H2B-mNeonGreen-P10, within the Minos transposon vector containing the PhOpsin1-EGFP transgenesis marker, as described earlier. The CREs were tested by microinjecting Parhyale embryos, as described above. Embryos were screened for CRE activity every day from the start of leg development (S19) to hatching, and PhOpsin1 activity was scored at the final developmental stages. CREs were considered inactive, or below our detection threshold, when at least 10 PhOpsin1-positive embryos were identified while embryos did not consistently show the expected distal or mid leg expression pattern. Long and short versions of the promoter region were tested for Dlle and dac2. Not all obvious ATAC-seq peaks around Dlle were tested since both long and short versions of the promoter alone gave weak ubiquitous expression. The sequences of these putative CREs are available in Suppl. Data 5, available at https://doi.org/10.5281/zenodo.19020963.

Appendix 1

Appendix 1—figure 1
Download asset Open asset
Sequence comparison of known Parhyale CREs with homologous regions from Hyalella.

Sequence comparisons of known Parhyale CREs with homologous regions from Hyalella azteca. Dot plot alignments between previously characterised Parhyale regulatory fragments from the hsc70 (PhMS, Pavlopoulos and Averof, 2005), opsin 1 and opsin 2 genes (Ramos et al., 2019) and homologous regions of the Hyalella azteca genome (accessions NW_025942174, NW_025945614 and NW_025930472, respectively; Poynton et al., 2018). Two putative opsin 1 genes were found in Hyallella. The dot plots were made using the EMBOSS Dotmatcher tool (Rice et al., 2000) with low stringency settings (window size 10, threshold 30). Regions of sequence conservation appear as diagonal lines. The only region of significant sequence similarity, in opsin 2, coincides with part of the coding sequence.

Appendix 1—figure 2
Download asset Open asset
Chromatin accessibility and sequence conservation in previously identified Parhyale CREs.

(A–C) Chromatin accessibility and sequence conservation profiles of previously identified Parhyale CREs. Genome browser plots showing the patterns of chromatin accessibility and cross-species sequence conservation in the three genomic loci of Parhyale hawaiensis that harbour the previously known cis-regulatory elements PhMS, PhOpsin1, and PhOpsin2 (highlighted in grey). We observe ATAC-seq peaks and overlapping islands of sequence conservation among Parhyale species in all three fragments. Note that the same elements show no significant sequence conservation when compared with homologous regions from the more distant species Hyalella azteca.

Data availability

The sequencing data of the bulk and single-nuclei ATACseq experiments were deposited in GEO (accessions GSE325175 and GSE325176). The NGS sequencing data from the genomes of P. aquilina, P. darvishi and P. plumicornis were deposited in NCBI (BioSample accessions SAMN56533921, SAMN56533922 and SAMN56533923). The supplementary data files are available at https://doi.org/10.5281/zenodo.19020963.

The following data sets were generated
    1. Forbes G
    2. Skafida E
    3. Cevrim C
    4. Averof M
    5. Paris M
    (2026) NCBI Gene Expression Omnibus
    ID GSE325175. Genome-wide discovery of cis-regulatory elements in a large genome [bulk ATAC-seq].
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE325175
    1. Forbes G
    2. Skafida E
    3. Averof M
    4. Paris M
    (2026) NCBI Gene Expression Omnibus
    ID GSE325176. Genome-wide discovery of cis-regulatory elements in a large genome [scATAC-seq].
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE325176
    1. Averof M
    2. Anastasiadou C
    3. Paris M
    (2026) NCBI BioSample
    ID SAMN56533921. Short genomic reads from a Parhyale aquilina individual.
    https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN56533921
    1. Momtazi F
    2. Averof M
    3. Paris M
    (2026) NCBI BioSample
    ID SAMN56533922. Short genomic reads from a Parhyale darvishi individual.
    https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN56533922
    1. Averof M
    2. Loa B
    3. Paris M
    (2026) NCBI BioSample
    ID SAMN56533923. Short genomic reads from a Parhyale plumicornis individual.
    https://www.ncbi.nlm.nih.gov/biosample/?term=SAMN56533923
    1. Forbes G
    2. SkafidaIrene E
    3. Savannah K
    4. Mowgli M
    5. Çağrı D
    6. Farzaneh C
    7. Chryssa M
    8. Sabrina A
    9. Lo B
    10. Averof M
    11. Paris M
    (2026) Zenodo
    Supplementary data for 'Genome-wide discovery of cis-regulatory elements in a large genome'.
    https://doi.org/10.5281/zenodo.19020963

References

  60. Book
    1. Wickham H
    (2016) Ggplot2
    Springer International Publishing.
    https://doi.org/10.1007/978-3-319-24277-4

Article and author information

Author details

  1. Gillian Forbes

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Generated and analysed the bulk ATAC-seq data, Optimised the transgenic reporter approach, Tested the CREs of developmental genes, Writing – review and editing
    Contributed equally with
    Emilia Skafida and Irene Karapidaki
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0007-7737-5926

  2. Emilia Skafida

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Generated bulk ATAC-seq data, Generated and analysed the single-nucleus ATAC-seq data
    Contributed equally with
    Gillian Forbes and Irene Karapidaki
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1671-5487

  3. Irene Karapidaki

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Optimised the transgenic reporter approach, Identified the neuron-specific CREs
    Contributed equally with
    Gillian Forbes and Emilia Skafida
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0008-8913-3219

  4. Savannah Moinet

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Identified the muscle-specific CREs
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0000-8883-3522

  5. Mowgli Dandamudi

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Identified the ubiquitous CREs
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3649-3882

  6. Çağrı Çevrim

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Generated bulk ATAC-seq data
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4720-7944

  7. Farzaneh Momtazi

    1. Iranian National Institute for Oceanography and Atmospheric Science (INIOAS), Marine Bioscience Department, Tehran, Iran
    2. Potsdam Institute for Climate Impact Research (PIK), Potsdam, Germany
    Contribution
    Collected and identified P. darvishi
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3074-0865

  8. Chryssa Anastasiadou

    Fisheries Research Institute, Hellenic Agricultural Organization, Nea Peramos, Greece
    Contribution
    Collected P. aquilina
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1780-887X

  9. Sabrina Lo Brutto

    1. Department of Earth and Marine Sciences (DiSTeM), University of Palermo, Palermo, Italy
    2. National Biodiversity Future Center (NBFC), Palermo, Italy
    Contribution
    Provided advice for field collections and identified P. plumicornis
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9964-904X

  10. Michalis Averof

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Conceived and supervised the project, Collected P. aquilina and P. plumicornis, Writing – original draft
    For correspondence
    michalis.averof@ens-lyon.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6803-7251

  11. Mathilde Paris

    1. Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, Lyon, France
    2. Centre National de la Recherche Scientifique (CNRS), Paris, France
    Contribution
    Conceived and supervised the project, Sequenced the genomes, Determined the genome-wide patterns of sequence conservation, Performed the phylogenetic analysis, Writing – review and editing
    For correspondence
    mathilde.paris@ens-lyon.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6166-6672

Funding

Fondation pour la Recherche Médicale (EQU202303016278)

  • Gillian Forbes
  • Irene Karapidaki
  • Mowgli Dandamudi
  • Michalis Averof

European Research Council

https://doi.org/10.3030/694918
  • Michalis Averof

Boehringer Ingelheim Fonds

  • Çağrı Çevrim

European Commission

https://doi.org/10.3030/766053
  • Emilia Skafida
  • Mathilde Paris

European Commission (CN_00000033, CUP B73C22000790001)

  • Sabrina Lo Brutto

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Oliver Coleman of the Natural History Museum in Berlin for advice on collecting Parhyale, Pavlos Vidoris for support during field work, Chiara Sinigaglia for advice on leg staging and dissection, Alba Almazan for advice on nuclear isolation, and Benjamin Gillet and Sandrine Hughes of the IGFL sequencing platform for support with Illumina sequencing. This research was supported by grants from the Fondation pour la Recherche Médicale (grant EQU202303016278) and the European Research Council (grant ERC-2015-AdG 694918), a doctoral fellowship from Boehringer Ingelheim Fonds, and a fellowship from the Marie Curie ITN programme EvoCell (H2020-MSCA-ITN-2017 #766053), under the European Union Horizon 2020 programme. SLB was supported by the National Recovery and Resilience Plan of the European Union (project CN_00000033, CUP B73C22000790001).

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.111378. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2026, Forbes, Skafida, Karapidaki et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 643
    views
  • 29
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Gillian Forbes
  2. Emilia Skafida
  3. Irene Karapidaki
  4. Savannah Moinet
  5. Mowgli Dandamudi
  6. Çağrı Çevrim
  7. Farzaneh Momtazi
  8. Chryssa Anastasiadou
  9. Sabrina Lo Brutto
  10. Michalis Averof
  11. Mathilde Paris
(2026)
Genome-wide discovery of cis-regulatory elements in a large genome
eLife 15:RP111378.
https://doi.org/10.7554/eLife.111378.3

Share this article

https://doi.org/10.7554/eLife.111378

Categories and tags