Genome-wide discovery of cis-regulatory elements in a large genome
- Institut de Génomique Fonctionnelle de Lyon (IGFL), École Normale Supérieure de Lyon, CNRS and UCBL Lyon 1, France
- Centre National de la Recherche Scientifique (CNRS), France
- Iranian National Institute for Oceanography and Atmospheric Science (INIOAS), Marine Bioscience Department, Iran
- Potsdam Institute for Climate Impact Research (PIK), Germany
- Fisheries Research Institute, Hellenic Agricultural Organization, Greece
- Department of Earth and Marine Sciences (DiSTeM), University of Palermo, Italy
- National Biodiversity Future Center (NBFC), Italy
Figures
Figure 1
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.
Figure 2 with 1 supplement
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.
Figure 2—figure supplement 1
Label transfer scores from snRNA-seq to snATAC-seq.
(A) Distribution of scores for the transfer of cell labels from the snRNA-seq (Almazán et al., 2022) to snATAC-seq (SN1 and SN2) data. (B) UMAP of the snATAC-seq experiments, colour-coded according to the label transfer scores. Both panels show that cell identities of several cell types (including muscles, neurons, epidermis, and haemocytes) can be transferred from RNA-seq to ATAC-seq clusters with relatively high confidence from the RNA-seq data.
Figure 3 with 1 supplement
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.
Figure 3—figure supplement 1
Estimates of genome coverage in P. aquilina, P. darvishi, and P. plumicornis.
Distribution of per nucleotide sequence coverage, based on the sequences of single-copy BUSCO genes in (A) P. hawaiensis, (B) P. aquilina, (C) P. darvishi, and (D) P. plumicornis (see ‘Materials and methods’).
Figure 4
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.
Figure 5 with 2 supplements
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.
Figure 5—figure supplement 1
Putative CREs of Parhyale neuron-specific genes tested using transgenic reporters.
Genome browser plots for 5 loci harbouring putative neuron-specific CREs (highlighted in grey). The tracks of the other two loci that we tested, neuro5 and neuro6, are shown in Figure 5A and B.
Figure 5—figure supplement 2
Comparison of the activity of neuro5 and neuro6 reporters.
Dorso-lateral view of a transgenic embryo carrying the neuro5>Src64B-mNeonGreen (in cyan) and neuro6-mScarlet3-HRas (in red) transgenes. Dorsal side is located towards the top-left. Max projection of image stack captured by confocal microscopy. Scale bar, 50 µm.
Figure 6
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.
Appendix 1—figure 1
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
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.
Tables
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.
| Dataset | Total mapped reads after filtering | Fraction of reads in peaks (FRiP) | Number of peaks | |
|---|---|---|---|---|
| Whole embryo | E20 | 16,431,070 | 0.12 | 60,435 |
| E23 | 8,726,174 | 0.15 | 56,606 | |
| E24 | 7,879,166 | 0.29 | 119,994 | |
| Embryonic legs | EL1 | 19,186,758 | 0.37 | 87,308 |
| EL2 | 18,000,384 | 0.36 | 77,357 | |
| EL3 | 23,424,388 | 0.42 | 119,487 | |
| Adult legs (bulk) | L1 | 20,233,184 | 0.32 | 76,268 |
| L2 | 16,508,456 | 0.19 | 49,502 | |
| L3 | 10,206,744 | 0.22 | 48,324 | |
| Adult legs (single-nuclei) | SN1 | 319,174,093 | 0.41 | 404,519 |
| SN2 | 217,929,607 | 0.43 | 331,355 | |
| Naked DNA | nDNA | 102,914,308 | 0.15 | 161,533 |
Table 2
Overview of short-read genome sequencing on three Parhyale species.
| Species | Estimated genome size (Gbp) | Sequenced reads | Genome coverage | % reads mapping to P. hawaiensis | Non-exonic regions conserved in P. hawaiensis |
|---|---|---|---|---|---|
| P. aquilina | 1.1 | 281,429,563 | 16× | 6.4 | 605,996 |
| P. darvishi | 2.8 | 318,968,916 | 11× | 4.1 | 525,822 |
| P. plumicornis | 3.0 | 424,305,645 | 10× | 0.64 | 156,063 |
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.
| Activity | Gene name, alias | Putative CRE | Length (bp) | Embryos screened | Opsin1- positive | Reporter expression |
|---|---|---|---|---|---|---|
| Ubiquitous | hdc (MSTRG.41953) ubi1 | P | 1027 | 123 | 38 | All tissues (38 embryos) |
| mbl (MSTRG.26075) ubi2 | CRE +P | 501+890 | 260 | - | All tissues (118 embryos) | |
| Neurons | futsch (MSTRG.441) neuro1 | CRE +P | 615+496 | 147 | 10 | No expression |
| jeb (MSTRG.26302) neuro2 | CRE +P | 976+233 | 148 | 25 | No expression | |
| ChT (mikado.phaw_ 50.283866G186) neuro3 | P | 1140 | 232 | 41 | No expression | |
| αTub (MSTRG.41851) neuro5 | P | 954 | 542 | 166 | CNS, peripheral neurons (108 embryos) | |
| Cdk5α (MSTRG.7309) neuro6 | P | 896 | 208 | 90 | CNS, peripheral neurons (55 embryos) | |
| CRE +P | 375+896 | 392 | 67 | CNS, peripheral neurons (81 embryos) | ||
| mav (MSTRG.827) neuro7 | CRE +P | 1871+230 | 30 | 8 | No expression | |
| MSTRG.27182 neuro8 | CRE +P | 358+751 | 220 | 88 | No expression | |
| Muscles | Mhc (MSTRG.35247) Mhc | P | 1003 | 72 | - | Muscles (9 embryos) |
| CRE1+P | 735+1,003 | 48 | 5 | Muscles (6 embryos) | ||
| CRE6+P | 1194+1,003 | 92 | 12 | Muscles (22 embryos) |
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.
| Gene | Putative CRE | Length (bp) | Embryos screened | Opsin1- positive | Reporter expression |
|---|---|---|---|---|---|
| Dll-e (MSTRG.27394) | P | 1004 | 243 | 45 | Weak ubiquitous after stage S20 |
| Pshort | 741 | 388 | 74 | Weak ubiquitous after stage S20 | |
| CRE3+Pshort | 1209+741 | 453 | 107 | Weak ubiquitous after stage S20 | |
| CRE6+P | 1517+1004 | 480 | 63 | Weak ubiquitous after stage S20 | |
| dac1 (MSTRG.33907) | P | 944 | 607 | 166 | No expression |
| CRE1+P | 2032+944 | 167 | 37 | No expression | |
| CRE2+P | 819+944 | 292 | 47 | No expression | |
| dac2 (MSTRG.33908) | P | 1070 | 386 | 112 | No expression |
| Pshort | 444 | 698 | 185 | No expression | |
| CRE1+Pshort | 1033+444 | 456 | 162 | No expression | |
| CRE2+P | 1048+1070 | 319 | 102 | No expression | |
| CRE3+P | 987+1070 | 80 | 28 | No expression | |
| CRE4+P | 1099+1070 | 99 | 26 | No expression | |
| CRE5+P | 585+1070 | 346 | 116 | No expression | |
| CRE8+P | 957+1070 | 103 | 30 | No expression | |
| CRE9+P | 916+1070 | 47 | 14 | No expression | |
| CRE10+P | 2358+1070 | 145 | 28 | No expression | |
| CRE11+P | 2973+1070 | 92 | 17 | No expression |
Key resources table
| Reagent type (species) or resource | Designation | Source or reference | Identifiers | Additional information |
|---|---|---|---|---|
| Biological sample (Parhyale hawaiensis) | Parhyale hawaiensis | Kao et al., 2016 | NCBI species ID 317513 | Chicago-F line, laboratory culture |
| Biological sample (Parhyale aquilina) | Parhyale aquilina | This paper | NCBI species ID 1799758 | Collected in Nea Peramos, Greece |
| Biological sample (Parhyale darvishi) | Parhyale darvishi | This paper | NCBI species ID 3695942 | Collected in Chabahar Bay, Iran |
| Biological sample (Parhyale plumicornis) | Parhyale plumicornis | This paper | NCBI species ID 1799757 | Collected in Milos, Greece |
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 type | Dataset | Sequencing output - number of reads | Number 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 embryo | E20 | 49,478,726 | 42,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 %) |
| E23 | 20,106,906 | 18,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 %) | |
| E24 | 25,718,862 | 18,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 legs | EL1 | 45,551,822 | 39,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%) |
| EL2 | 44,018,066 | 33,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%) | |
| EL3 | 59,349,072 | 49,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) | L1 | 52,289,208 | 36,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%) |
| L2 | 44,194,074 | 37,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%) | |
| L3 | 39,923,474 | 19,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) | SN1 | 1,351,387,910 | 1,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%) |
| SN2 | 809,094,758 | 649,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 DNA | nDNA | 246,377,216 | 236,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%) |
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 type | Dataset | All 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 embryo | E20 | 2820 | 2142 | 3,434,318 | 1,551,895 | 4,184,806 | 10,599,571 |
| E23 | 3158 | 2669 | 2,368,108 | 1,256,933 | 3,680,699 | 10,594,791 | |
| E24 | 5842 | 5001 | 4,480,014 | 3,142,237 | 10,417,667 | 22,338,694 | |
| Embryo legs | EL1 | 5163 | 4260 | 4,608,493 | 2,682,977 | 9,419,742 | 15,523,581 |
| EL2 | 5161 | 4421 | 3,808,532 | 2,242,178 | 8,236,839 | 13,758,216 | |
| EL3 | 5832 | 4865 | 5,816,624 | 3,698,369 | 12,442,011 | 20,580,140 | |
| Adult legs | L1 | 6193 | 5360 | 4,212,182 | 2,470,145 | 9,391,666 | 15,678,750 |
| L2 | 4066 | 3301 | 3,265,858 | 1,521,817 | 4,817,687 | 9,726,381 | |
| L3 | 4343 | 3821 | 2,558,498 | 1,378,514 | 4,231,107 | 8,936,821 | |
| Naked DNA | nDNA | 5551 | 2496 | 12,311,372 | 6,304,430 | 33,198,617 | 70,905,248 |
Additional files
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)
Genome-wide discovery of cis-regulatory elements in a large genome
eLife 15:RP111378.
https://doi.org/10.7554/eLife.111378.3
