Abstract

Background

   The Pacific white shrimp (Litopenaeus vannamei) is the most widely
   farmed shrimp species globally, yet the epigenetic regulation
   underlying its embryonic development remains largely unexplored.
   Histone modifications are known to orchestrate gene expression during
   early development in model organisms, but their role in crustaceans is
   poorly understood.

Results

   In this study, we present the first comprehensive histone modification
   landscape during L. vannamei embryogenesis using CUT&Tag (Cleavage
   Under Targets and Tagmentation). We profiled high-resolution landscapes
   of four histone marks (H3K4me1, H3K4me3, H3K27ac, H3K27me3) across
   seven developmental stages from blastula to nauplius, revealing dynamic
   chromatin state transitions associated with developmental progression.
   Integration with transcriptomic data uncovered a strong temporal
   correlation between chromatin states and gene expression, particularly
   during zygotic genome activation (ZGA). Furthermore, our analysis
   uncovered key developmental genes associated with critical biological
   processes such as molting, body segmentation, and neurogenesis,
   providing novel insights into the epigenetic regulation of these
   events. Functional annotation of cis-regulatory elements based on
   histone marks identified candidate enhancers and regulatory loci linked
   to these key genes.

Conclusions

   Our study provides the first epigenomic framework of shrimp
   embryogenesis, uncovering chromatin-based regulatory mechanisms during
   early development. The identification of stage-specific enhancers and
   active chromatin regions offers valuable resources for functional
   genomics in crustaceans and sheds light on conserved and divergent
   aspects of ZGA regulation beyond model systems.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s13072-025-00615-4.

   Keywords: Litopenaeus vannamei, CUT&Tag, Histone modifications,
   Epigenetic regulation, Zygotic genome activation

Background

   The Pacific white shrimp (Litopenaeus vannamei) is one of the most
   economically significant farmed shrimp species globally, owing to its
   rapid growth, high yield, and superior meat quality. According to the
   2022 statistics from the Food and Agriculture Organization (FAO), the
   aquaculture production of L. vannamei accounted for 51.7% of the
   world’s total crustacean production, reaching 5.8 million tons, ranking
   it first worldwide [[36]1]. As a member of the Decapoda within
   Crustacea, L. vannamei undergoes a unique early developmental pattern,
   which includes several distinct stages: embryo, nauplius, zoea, mysis,
   and postlarvae. During the embryonic phase, development progresses from
   a zygote to a multicellular structure, and through key stages such as
   blastula, gastrula, limb bud embryo, and larva in membrane. Following
   hatching, the shrimp progresses through successive nauplius, zoea,
   mysis, and post larval stages, marked by distinct morphological changes
   [[37]2, [38]3]. These early developmental stages are marked by
   significant morphological and physiological changes, which play a
   crucial role in the overall life cycle. The process of metamorphosis
   represents a critical transition in L. vannamei’s development and
   serves as an example of life cycle modularity [[39]3, [40]4]. However,
   these pronounced physiological changes during early development can
   lead to high variability and difficulty in larval rearing, posing a
   significant challenge in shrimp aquaculture [[41]3, [42]5].

   The early embryonic development of animals is primarily regulated by
   maternal gene products. As development progresses, control gradually
   shifts to the zygotic genome during the maternal-to-zygotic transition
   (MZT), a process characterized by precise and coordinated changes in
   gene expression. Studies have shown that the timing and intensity of
   ZGA vary across species, reflecting the diversity of regulatory
   mechanisms involved [[43]6]. Despite its critical role in development,
   ZGA in L. vannamei remains unexplored, leaving a significant gap in
   understanding how gene expression is regulated during this pivotal
   phase. Notably, L. vannamei exhibits unique developmental features,
   such as extensive physiological changes in tissue formation, dietary
   transitions, muscle development, and exoskeleton remodeling. These
   processes, while essential for its life cycle, are still not fully
   understood. Current studies on L. vannamei have primarily focused on
   its morphological, biochemical, and transcriptomic characteristics
   [[44]3], with limited attention to the epigenetic mechanisms,
   particularly histone modifications, that underlie gene expression
   regulation during early embryogenesis.

   Epigenetic regulation, including DNA methylation, histone
   modifications, and chromatin remodeling, modulates gene expression
   without altering the DNA sequence [[45]7]. Histones, as fundamental
   components of chromatin, undergo various post-translational
   modifications (PTM), such as methylation, acetylation, and
   phosphorylation, which impact chromatin dynamics and facilitate the
   binding of transcriptional regulators [[46]8, [47]9].
   Post-translational modifications of histones have been widely
   demonstrated to play critical roles in regulating gene expression
   during animal embryonic development [[48]8, [49]10]. Moreover, specific
   types of histone modifications can serve as predictors of gene
   expression activity [[50]11, [51]12]. Among these, four
   well-characterized histone modifications include H3K4me1
   (monomethylation of histone 3 lysine 4), H3K4me3 (trimethylation of
   histone 3 lysine 4), H3K27me3 (trimethylation of histone 3 lysine 27),
   and H3K27ac (acetylation of histone 3 lysine 27). H3K4me1 is associated
   with active enhancers, H3K4me3 marks promoters of actively transcribed
   genes, H3K27me3 is linked to transcriptional repression, and H3K27ac
   indicates regions of open chromatin. These histone modifications play
   pivotal roles in regulating gene expression and chromatin structure
   during developmental processes [[52]11, [53]12]. However, their roles
   in the early embryonic development of L. vannamei remain largely
   unexplored.

   Cleavage Under Targets & Tagmentation (CUT&Tag) is an advanced
   genome-wide chromatin profiling technique that has emerged as a
   powerful tool for studying histone modifications [[54]13]. Compared to
   traditional chromatin immunoprecipitation followed by sequencing
   (ChIP-seq), CUT&Tag offers several advantages, including a higher
   signal-to-noise ratio, greater sensitivity, reduced experimental time,
   and lower cell requirements [[55]13]. This method has proven effective
   for analyzing chromatin modifications and transcription factor binding
   sites [[56]14]. This innovative approach provides a powerful tool to
   investigate the epigenetic landscape regulating early embryonic
   development in L. vannamei.

   In this study, we employed CUT&Tag and RNA-seq technologies to generate
   high-resolution maps of the histone modifications H3K4me1, H3K4me3,
   H3K27me3, and H3K27ac during critical stages of early embryonic
   development in L. vannamei. Our primary goal is to investigate the role
   of these histone modifications in the MZT and to capture the dynamic
   epigenetic changes that regulate gene expression during this period.
   This study advances our understanding of epigenetic regulation during
   L. vannamei embryogenesis, offering insights to improve shrimp
   aquaculture. Furthermore, by annotating the genome-wide cis-regulatory
   elements of L. vannamei using multiple histone modifications, this
   research lays the groundwork for further functional genomics studies in
   this economically important species.

Methods

Experimental animal

   The samples of L. vannamei used in this study were collected on-site at
   Hairen Aquatic Seed Industry Technology Co., Ltd (Hebei, China).
   High-quality male and female broodstock were selected by experienced
   technicians for artificial mating. The reproductive behavior of the
   shrimp was closely monitored, and once spawning was observed,
   fertilized eggs were immediately collected. The eggs were washed with
   clean seawater and incubated under strictly controlled conditions
   (temperature: 30 ± 0.5 °C; salinity: 33‰). Embryonic development was
   monitored in real time under a microscope. At each target developmental
   stage, embryos were promptly collected using a nylon mesh to ensure
   stage-specific sampling.

   Samples for CUT&Tag were processed immediately after collection. For
   RNA-seq, embryos at the same developmental stages were washed three
   times with sterile seawater, flash-frozen in liquid nitrogen, and
   transported on dry ice to our laboratory, where they were stored at
   − 80 °C until further processing. This study included seven
   developmental stages: blastula (blas), gastrula (gast), limb bud embryo
   (lbe), larva in membrane (lim), nauplius I (N1), nauplius III (N3) and
   nauplius VI (N6) stages.

CUT&Tag experiment and data analysis

   The CUT&Tag libraries were prepared using the Illumina Hyperactive
   Universal CUT&Tag Assay Kit (Vazyme, China, #TD903), with nuclei
   extracted from the embryos following a previously described protocol
   [[57]14]. Approximately 200 embryos were collected at each
   developmental stage for each biological replicate using a nylon mesh.
   The freshly collected embryos were then washed 2–3 times with Wash
   Buffer to remove residual seawater and impurities. To isolate nuclei,
   300 µL of pre-chilled NE Buffer was added to the embryos, which were
   then homogenized on ice using sterile surgical scissors. An additional
   500 µL of pre-chilled NE Buffer was added, and the homogenate was
   pipetted thoroughly to generate a single-cell suspension. The
   suspension was incubated on ice for 5 min and filtered through a 40 μm
   strainer to remove debris. The filtrate was centrifuged at 1300 × g for
   4 min, and the resulting pellet was washed 2–3 times with 500 µL of
   Wash Buffer. Finally, the nuclei were resuspended in 100 µL of Wash
   Buffer, mixed gently, and prepared for downstream applications. All
   CUT&Tag procedures were performed in accordance with the manufacturer’s
   instructions and the published protocol [[58]14]. Two biological
   replicates were included for each developmental stage, and
   approximately 50,000 nuclei were used per sample. Antibodies used in
   this study are listed in Table [59]S1. The CUT&Tag libraries were
   sequenced on the Illumina NovaSeq 6000 platform (Novogene, Tianjin,
   China).

   Raw sequencing data were subjected to quality control using the NGS QC
   Toolkit v2.3 to generate clean reads [[60]15]. The clean reads were
   aligned to the L. vannamei reference genome (NCBI accession:
   ASM378908v1) using the BWA-MEM v0.7.17 algorithm [[61]16], and mapping
   statistics are summarized in Table [62]S2. Using awk, reads with
   abnormal pairing relationships (e.g., no distance or distances greater
   than 1000 bp) were removed, retaining only fragments with lengths
   between 100 and 1000 bp. These filtered fragments were subsequently
   used for peak calling with MACS2 v2.2.7.1 [[63]17], applying the
   following parameters: macs2 callpeak -t A.bam -c IgG.bam -f BAMPE -g
   1.66e9 --keep-dup auto -q 0.05.

   Annotation and visualization of the peaks were conducted using R
   packages GenomicFeatures [[64]18], ChIPseeker [[65]19], and ggplot2
   [[66]20]. Additional analyses included fragment size distribution,
   reproducibility assessment of replicate samples, signal heatmaps of
   H3K4me3 at transcription start sites (TSS), peak signal enrichment
   heatmaps, and peak distribution across genomic features. The analyses
   followed R scripts from a previously published pipeline [[67]21].

Chromatin state analysis

   Chromatin state analysis was performed using the ChromHMM v1.26 model,
   which identifies distinct chromatin states and annotates their
   genome-wide distribution based on combinatorial patterns of histone
   modifications at each locus [[68]22]. The model was trained using
   dynamic profiles of four key histone modifications (H3K4me1, H3K4me3,
   H3K27me3, H3K27ac) obtained from the above mentioned seven embryonic
   developmental stages. To determine the optimal number of chromatin
   states, we performed parameter tuning and model evaluation. The results
   indicated that models with more than 10 states exhibited substantial
   functional redundancy among states, while models with fewer than 10
   states failed to effectively distinguish certain biologically
   meaningful features. Therefore, the final model employed the LearnModel
   function to define 10 chromatin states, balancing the need to reduce
   redundancy while preserving the resolution required to capture critical
   biological signatures. The genome-wide distribution patterns of
   chromatin states across distinct embryonic developmental stages are
   presented in Supplementary File [69]1. The chromatin state maps were
   visualized using Integrative Genomics Viewer (IGV) v2.15.2 to
   facilitate interpretation [[70]23].

RNA-seq experiment and data analysis

   RNA-seq samples were collected in parallel with those used for the
   CUT&Tag experiment, from the same batch of embryos, with three
   biological replicates for each developmental stage. Total RNA was
   extracted using the RNAprep Pure Tissue Kit DP431 (TIANGEN, Beijing,
   China), followed by DNase I (Takara Bio, Inc.) treatment to remove DNA
   contamination. RNA quality and purity were assessed using a Nanovue
   Plus spectrophotometer (GE Healthcare, New Jersey, USA) and agarose gel
   electrophoresis. Transcriptome libraries were prepared using the
   NEBNext^® Ultra™ RNA library prep kit for Illumina^® (NEB, Ipswich, MA,
   USA) and sequenced on the Illumina NovaSeq PE150 platform (Illumina,
   San Diego, CA, USA), generating approximately 6 GB of raw data per
   sample.

   The raw sequencing data underwent quality control using the NGS QC
   Toolkit v2.3 [[71]15] to produce clean reads (quality control results
   are summarized in Table [72]S3). The clean reads were aligned to the L.
   vannamei reference genome using STAR v2.7.11b [[73]24]. Gene
   quantification was performed using featureCounts [[74]25], and
   expression levels were normalized to transcripts per million mapped
   reads (TPM).

   For data exploration, principal component analysis (PCA) was conducted
   after z-score normalization. Differential gene expression analysis was
   carried out using DESeq2 [[75]26] with thresholds of|fold change| >2
   and p < 0.05. KEGG Pathway enrichment analysis of differentially
   expressed genes (DEGs) was conducted using the R package
   clusterProfiler. The expression levels of key genes during embryonic
   development were quantified in TPM and visualized using GraphPad Prism
   v8.0.2.

Temporal analysis of gene expression dynamics

   Gene expression dynamics across developmental stages were analyzed
   using Short Time-series Expression Miner (STEM) v1.3.13. The average
   TPM values from three biological replicates per stage were
   log2-normalized, and nine model profiles were pre-defined. Significant
   trends (p < 0.05) were identified, with a minimum fold-change threshold
   of 2 for gene selection. Genes from each profile were further analyzed
   for functional enrichment using KEGG pathway analysis via the
   clusterProfiler R package.

Integrated analysis of cut & tag and RNA-seq data

   To minimize the impact of low-confidence and weak signals that might
   confound downstream functional analyses and to enhance the robustness
   of our conclusions, we employed a stringent peak filtering strategy.
   Specifically, Bedtools v2.31.1 [[76]27] was used to identify
   intersecting peaks between the two biological replicates, and only
   high-confidence peaks that were consistently detected across both
   samples were retained. These validated peaks were then used for
   subsequent analyses, including histone modification peak
   quantification, gene-level annotation, and integration with
   transcriptomic datasets.

   Histone modifications were analyzed based on their well-characterized
   regulatory roles: H3K4me3 is enriched at the promoters of actively
   transcribed genes, H3K27me3 is associated with transcriptional
   repression mediated by the Polycomb Repressive Complex 2 (PRC2),
   H3K4me1 serves as a marker of enhancers, and H3K27ac identifies
   accessible chromatin regions [[77]11, [78]12]. Based on these
   definitions, genomic regions concurrently marked by H3K4me3 and
   H3K27ac, but devoid of H3K27me3, were designated as “active open
   sites”. Genes associated with these regions and exhibiting an
   expression level of TPM > 0.5 were defined as actively expressed genes.
   KEGG pathway enrichment analysis was performed on this gene set to
   identify significantly overrepresented pathways (p < 0.05). To assess
   dynamic expression trends, we employed STEM [[79]28] on log₂-normalized
   expression data. Trends were considered significant if p < 0.05, with a
   minimum fold change ≥ 2 used as a selection criterion. In parallel,
   genomic regions showing an overlap of at least 1 bp between H3K4me1 and
   H3K27ac peaks, while lacking H3K27me3, were defined as active enhancers
   at specific developmental stages.

   Visualization of selected gene loci was performed using IGV v2.15.2
   [[80]23], and statistical analyses of gene expression were conducted
   using GraphPad Prism v8.0.2. Gene expression heatmaps were standardized
   using z-score normalization to facilitate cross-sample comparison.

Results

High resolution mapping of histone modifications in L. vannamei embryos using
CUT&Tag

   To validate the application of CUT&Tag in L. vannamei embryos, we first
   assessed its efficiency. Alignment rates of sequencing reads to the L.
   vannamei reference genome consistently exceeded 90% across all samples,
   confirming high-quality data generation (Table [81]S2). Taking the N3
   stage as an example, the fragment size distribution displayed the
   typical nucleosomal ladder pattern, characteristic of successful
   CUT&Tag experiments in the histone modification samples. In contrast,
   control samples exhibited a random size distribution, further
   validating the specificity of the technique (Fig. [82]1A) [[83]29]. The
   reproducibility of the experiment was assessed by correlation analysis,
   which showed high consistency among biological replicates (Fig. [84]1B
   & S1A). To further validate the specificity of the CUT&Tag results and
   the accuracy of peak calling, we conducted signal enrichment analysis
   at peak regions using heatmaps. The results showed significant signal
   enrichment at peak centers for histone modification groups, whereas the
   control group (Immunoglobulin G antibody, IgG antibody) exhibited no
   enrichment due to the lack of specific binding at these regions
   (Fig. [85]1C). In addition, H3K4me3 showed strong enrichment around
   transcription start sites (TSS) (Fig. [86]1C), consistent with its
   known localization pattern [[87]30], confirming the specificity and
   reliability of the data. Similar patterns were observed for the other
   six developmental stages (Figure [88]S1A). Genome-wide analysis of
   histone mark distribution at the seven developmental stages revealed
   that H3K4me1, H3K4me3, H3K27ac, and H3K27me3 were predominantly located
   in promoter regions and distal intergenic regions, rather than within
   gene bodies, consistent with their regulatory roles in gene expression
   [[89]11, [90]12]. Specifically, approximately 20%-30% of H3K4me1,
   50%-70% of H3K4me3, 10% of H3K27me3, and 30%-40% of H3K27ac were
   located in promoter regions. In distal intergenic regions, around
   50%-60% of H3K4me1, 30% of H3K4me3, 60%-70% of H3K27me3, and 40%-50% of
   H3K27ac were distributed (Fig. [91]1D and E, S1B; Table [92]S4). These
   results demonstrate that CUT&Tag provides robust detection of histone
   modifications in the embryonic L. vannamei genome, offering reliable
   insights into histone mark distribution and their potential regulatory
   roles during early development.

Fig. 1.

   [93]Fig. 1
   [94]Open in a new tab

   Quality assessment of CUT&Tag at the nauplius III stage. (A) Fragment
   size distribution. Plots showing the proportion (left) and raw number
   (right) of DNA fragments by size. Histone modification samples show
   strong enrichment of mononucleosome-sized fragments, whereas control
   samples exhibit a random size distribution. (B) Replicate
   reproducibility. Correlation analysis of genome-wide read distribution
   using 500 bp bins demonstrates high reproducibility between biological
   replicates of histone modification samples. (C) Signal enrichment at
   specific loci. The left two panels show the enrichment of CUT&Tag
   signal for H3K4me3 and IgG control at transcription start sites (TSS)
   of annotated genes. The right two panels depict signal enrichment
   within ± 1 kb of the peak summit center for H3K27me3 and IgG. Heatmaps
   (bottom) and line plots (top) confirm significant enrichment for
   histone marks compared to controls. (D) Peak distribution across
   genomic features. (E) Enrichment of peaks relative to the TSS

Chromatin state annotation based on histone modification marks during early
embryonic development in L. vannamei

   The distribution patterns of histone modifications and their influence
   on chromatin accessibility enable genome-wide chromatin state
   annotation [[95]22]. Such annotations could serve as a valuable
   resource for identifying regulatory features associated with key
   developmental processes in the L. vannamei genome. To provide the L.
   vannamei research community with a comprehensive, dynamic, and
   accessible genomic annotation, we utilized the dynamic profiles of four
   informative histone modifications obtained across multiple early
   developmental stages in this study. Using ChromHMM, we performed a
   detailed analysis of chromatin states in the L. vannamei genome during
   early embryonic development. As this is the first chromatin state
   analysis conducted for this species, our study provides a novel
   perspective. Ten distinct chromatin states were identified and named
   following ChromHMM conventions [[96]22]. The states are as follows:
   Bivalent enhancer 1 (EnhBiv1), Weak transcription (TxWk), Active TSS
   (TssA), Flanking TSS downstream (TssFlnkD), Weak flanking bivalent
   TSS/Enhancer (BivFlnkWk), Bivalent enhancer 2 (EnhBiv2), Quiescent/low
   (Quies), Flanking bivalent TSS/Enhancer (BivFlnk), Weak repressed
   polycomb (ReprPCWk), and Repressed polycomb (ReprPC). These chromatin
   states encompass a range of functional patterns, including active
   transcription, enhancers, bivalent regions, repressive marks, and
   quiescent regions (Figs. [97]2A, B).

Fig. 2.

   [98]Fig. 2
   [99]Open in a new tab

   ChromHMM-based chromatin state analysis in L. vannamei. (A) Heatmap
   (left) showing emission parameters for the Roadmap Epigenomics 10-state
   expanded model, based on H3K4me1, H3K4me3, H3K27me3, and H3K27ac
   signals. Darker colors indicate higher signal enrichment, while lighter
   colors represent weaker signals. The adjacent heatmap (right) displays
   fold enrichment for each chromatin state in 200-bp bins within 2 kb of
   transcription start sites (TSSs). State descriptions and abbreviations
   are provided to the right of the heatmaps. (B) Example browser view of
   ChromHMM-defined chromatin states for the nauplius III embryo genome.
   The top track represents the annotated reference genome, followed by an
   integrated chromatin state track and individual tracks for each
   chromatin state. (C) ChromHMM genome annotations for seven embryonic
   developmental stages, illustrated using the [100]XM_027355294.1 (e75)
   gene region as an example. Each row corresponds to a specific stage,
   with consistent color coding for the chromatin states in A, B and C.
   (D) Bar chart displaying the temporal expression levels of the e75 gene
   during L. vannamei embryogenesis, quantified in TPM. (E) IGV tracks
   showing e75 expression and histone modification patterns from the
   gastrula to nauplius VI stages. The signal dynamics of active chromatin
   marks (H3K4me3 and H3K27ac) correlate well with the expression level of
   e75

   To assess the validity of chromatin state annotations, we manually
   examined the relationship between gene expression levels and chromatin
   state categories. Overall, gene expression levels corresponded well
   with the annotated chromatin states. As a representative example, we
   investigated the e75 gene ([101]XM_027355294.1), known to be associated
   with molting. From the gastrula to N6 stages, the promoter region of
   e75 was annotated as an Active TSS (TssA) chromatin state,
   characterized by strong enrichment of H3K4me1, H3K4me3, and H3K27ac.
   Visualization using IGV and gene expression analysis revealed that e75
   expression begins during the gastrula stage, consistent with its
   chromatin state and histone modification patterns (Fig. [102]2C-E).

   These findings support the reliability and biological relevance of the
   chromatin state predictions. The chromatin state analysis results,
   including the stage-specific tracks of each individual chromatin state
   as well as composite tracks integrating all identified states, as
   illustrated in Fig. [103]2B and C, are detailed in Supplementary file
   [104]1. This is the first chromatin state annotation study for L.
   vannamei and it offers a foundational resource for future studies
   investigating the molecular mechanisms of L. vannamei embryonic
   development.

Dynamic changes in histone modifications and transcriptome during ZGA

   The precise and orderly activation of the zygotic genome is crucial for
   normal embryonic development, and histone modifications have been shown
   to play a key role in this process [[105]6]. To uncover the dynamic
   changes in histone modification and transcriptome during ZGA in L.
   vannamei, we analyzed the four histone modifications (H3K4me1, H3K4me3,
   H3K27me3, and H3K27ac) alongside transcriptomic profiles across seven
   key developmental stages. Principal component analysis (PCA)
   demonstrated strong similarity between the transcriptomes of the
   replicate samples across developmental stages, confirming the
   reliability of the transcriptomic data (Figure [106]S2).

   The number of significantly enriched histone modification peaks
   (Fig. [107]3A) and the heatmap of signals surrounding the peaks
   (Fig. [108]3C) revealed that H3K4me1 and H3K4me3 were already highly
   enriched starting from the blastula stage, while H3K27me3 and H3K27ac
   exhibited significant enrichment beginning at the gastrula stage.
   Notably, H3K4me1, H3K27me3, and H3K27ac reached their highest levels of
   enrichment at the N1 stage, whereas H3K4me3 peaked during the lim
   stage. The number of genes corresponding to these histone modification
   peaks also followed the same trend (Fig. [109]3B), indicating that the
   signal changes of these histone modifications were relatively evenly
   distributed across these different developmental stages rather than
   being concentrated in some specific regions.

Fig. 3.

   [110]Fig. 3
   [111]Open in a new tab

   Dynamics changes of histone modifications and transcriptome during ZGA.
   (A) The number of histone mark-associated peaks across seven
   developmental stages. (B) The number of histone mark-associated genes.
   (C) Heatmaps of signal enrichment. Enrichment patterns of H3K4me1,
   H3K4me3, H3K27me3, and H3K27ac, with signals centered on peak summits
   and extended to ± 1,000 bp regions across the seven developmental
   stages. (D) The proportion of expressed genes (with an average TPM
   value > 0.5) relative to the total gene count in the L. vannamei genome
   at each developmental stage. (E) A schematic model depicting the
   temporal dynamics of histone marks (red curves) and transcription
   levels (blue curves) during early embryonic development. Developmental
   stage-specific embryonic schematics in this figure were adapted from a
   previously published study [[112]3]

   During the MZT, maternal mRNA undergoes degradation before the
   activation of zygotic gene expression. Consequently, the number of
   expressed genes in the early embryonic transcriptome decreases
   initially, followed by a gradual increase upon ZGA [[113]6]. In L.
   vannamei, our study showed that expressed genes during the blastula
   stage accounted for 23.27% of the total genes, with a steady increase
   observed across developmental stages until expression plateaued during
   the lim stage (Fig. [114]3D). This suggests that ZGA in L. vannamei
   likely initiates before the blastula stage. This interpretation is
   supported by the enrichment of H3K4me1 and H3K4me3 signals, along with
   detectable levels of H3K27me3 and H3K27ac, as early as the blastula
   stage. This is consistent with observations in other species, where ZGA
   is typically accompanied by the re-establishment of histone
   modifications [[115]6].

   Overall, our data demonstrate that histone modifications and
   transcriptomic changes are tightly coordinated during early embryonic
   development in L. vannamei, displaying a strong temporal dependency and
   correlation (Fig. [116]3E).

Differential gene expression analysis across developmental stages of L.
vannamei embryos

   To investigate gene expression dynamics across developmental stages of
   L. vannamei embryos, differential gene expression (DEG) analysis was
   performed between adjacent stages, followed by KEGG pathway enrichment
   analysis. The number of upregulated and downregulated DEGs at each
   stage is shown in Fig. [117]4A.

Fig. 4.

   [118]Fig. 4
   [119]Open in a new tab

   Analysis of differentially expressed genes (DEGs) by RNA-seq. (A) Bar
   chart of DEGs between neighboring developmental stages. (B) Bubble
   chart of KEGG pathway enrichment analysis of DEGs across developmental
   stages. Bubble size represents the number of DEGs associated with each
   pathway, and bubble color intensity reflects the significance level of
   enrichment

   Notably, a comparison between the gastrula and blastula stages revealed
   a predominance of upregulated DEGs, likely reflecting the large-scale
   activation of the zygotic genome during this phase. As development
   progressed, both upregulated and downregulated genes increased in
   number, correlating with the heightened levels of repressive H3K27me3
   and active H3K27ac histone modifications observed at later stages. This
   pattern aligns with the increasing complexity of embryonic tissues and
   organs, necessitating precise gene expression regulation for proper
   cell differentiation [[120]3].

   KEGG pathway enrichment analysis identified stage-specific biological
   processes and signaling pathways. In the gastrula stage, DEGs were
   significantly enriched in pathways related to cell status and signal
   transduction, including tight junctions, apoptosis, and the cGMP-PKG
   pathway. During the lbe stage, DEGs were enriched in pathways involved
   in transcription, translation, and signal transduction, such as
   ribosome biogenesis, mRNA surveillance, and the Hippo signaling
   pathway. For the lim stage, DEGs were associated with pathways related
   to the nervous system, immune system, digestive system, endocrine
   system, and sensory systems, including Wnt signaling, Hippo signaling,
   neuroactive ligand-receptor interaction, platelet activation, fat
   digestion and absorption, thyroid hormone signaling, and
   phototransduction. At the nauplius stages (N1, N3, N6), DEGs were
   enriched in pathways related to metabolism (amino acids, lipids,
   carbohydrate, glycan biosynthesis) and hormone biosynthesis, as well as
   pathways linked to the digestive and nervous systems (Fig. [121]4B).

   These results demonstrate that DEGs at different developmental stages
   are involved in distinct biological processes and signaling pathways,
   reflecting the dynamic changes in cell status and function during
   embryonic development in L. vannamei.

Functional characterization of actively expressing genes at each
developmental stage of L. vannamei embryo

   While transcriptomics reflects cumulative transcript levels, histone
   modifications provide real-time insights into the chromatin’s
   accessibility and regulatory status. These modifications help identify
   cis-regulatory elements, such as promoters and enhancers, and clarify
   their roles in gene expression and regulation [[122]11, [123]12].
   During embryonic stages that inherit large amounts of maternal mRNA,
   histone modification data help mitigate the confounding effects of
   maternal mRNA fluctuations, enabling the identification of key
   zygotically expressed genes essential for developmental progression
   [[124]31]. To define such regulatory elements, we identified “active
   open sites” as genomic regions marked by H3K4me3 and H3K27ac, but
   lacking H3K27me3. Genes associated with these regions and exhibiting
   nonzero expression were classified as actively expressing genes, which
   are likely to play crucial regulatory roles during specific
   developmental stages.

   To investigate this, we integrated CUT&Tag and RNA-seq data to identify
   actively expressing genes at each developmental stage of L. vannamei
   embryogenesis. The results indicated that most genes associated with
   active open sites exhibited expression at certain stages
   (Fig. [125]5A), suggesting that chromatin modification states can
   reliably predict gene expression. The absence of expression in some
   cases could be due to the fact that chromatin opening precedes
   transcriptional activation [[126]32]. To validate the functional
   significance of these genes, we performed a KEGG pathway enrichment
   analysis, revealing stage-specific functional patterns. The results
   indicated that the majority of significantly enriched pathways were
   associated with signal transduction, and the key pathways at different
   developmental stages exhibited a high degree of consistency with their
   biological characteristics. During the gastrula and lbe stages, the
   significant enrichment of pathways related to cell growth and
   differentiation (e.g., the Hippo signaling pathway) closely aligned
   with the rapid proliferation of embryonic cells and the events of germ
   layer differentiation [[127]3, [128]33]. At the lim stage, the
   significant enrichment of immune system pathways (e.g., the T cell
   receptor signaling pathway) coincided with the critical period when
   larvae lose the protection of the egg membrane, reflecting the
   biological need for the rapid establishment of immune defense
   mechanisms during this stage [[129]34]. In the nauplius stages (N1-N6),
   the significant enrichment of transport and catabolism pathways (e.g.,
   fatty acid metabolism) was closely related to the energy remodeling
   required for enhanced locomotion and dietary shifts in the larvae
   (Fig. [130]5B) [[131]3, [132]33].

Fig. 5.

   [133]Fig. 5
   [134]Open in a new tab

   Functional enrichment of actively expressed genes at each developmental
   stage. (A) Venn diagram illustrating the overlap between genes
   identified by CUT&Tag and RNA-seq. The pink circle represents genes
   expressed at a given stage (TPM > 0.5), while the blue circle
   represents genes associated with active open sites (marked by H3K4me3
   at promoter regions and overlapping with H3K27ac peaks, but not
   H3K27me3). The intersection of these two sets highlights genes that are
   actively expressed at that stage. (B) KEGG pathway enrichment analysis
   of the actively expressed genes, with significantly enriched pathways
   (p < 0.05) visualized in a bubble plot

   Further, to understand the temporal expression patterns of these genes
   and their relationship to developmental stages, we conducted gene
   expression temporal analysis. This categorized the genes into nine
   distinct profiles. These profiles exhibit expression characteristics
   correlated with key developmental events and significantly enriched
   biological pathways (Fig. [135]6). Profile 1 shows the strongest
   expression during the blastula stage, followed by a gradual
   downregulation. It is significantly enriched in pathways such as the
   Hippo signaling pathway and apoptosis-related pathways, suggesting its
   role in regulating early embryonic proliferation and morphogenesis.
   During the blastula stage, cells maintain proliferative homeostasis
   through the Hippo pathway, while the apoptosis pathway eliminates
   abnormal cells [[136]35, [137]36]. The downregulation of proliferation
   genes after the gastrula stage marks a shift in developmental focus
   toward tissue differentiation [[138]3]. Profile 3 exhibits peak
   expression during the lim to N3 stages and is enriched in pathways such
   as T-cell receptor signaling, neurotrophin signaling, and the
   phosphatidylinositol signaling system. These pathways may drive immune
   cell differentiation, synaptic formation, and neural crest migration,
   indicating that the lim stage is a critical node for the development of
   the nervous and immune systems [[139]34, [140]37]. Profile 6 exhibits a
   fluctuating expression pattern, with its dynamic expression
   significantly associated with transcriptional regulation, translation
   mechanisms, and the TGF-β signaling pathway. This dynamic regulatory
   mechanism may precisely coordinate the spatiotemporal specificity of
   embryonic development by phase-specific activation or inhibition of key
   pathways [[141]38]. In contrast, Profile 8 maintains high expression
   levels from the gastrula stage to the lim stage and from the N3 to N6
   stage, and is significantly enriched in signaling transduction pathways
   such as mTOR and HIF-1. It likely provides sustained support for global
   differentiation and organ maturation during embryonic development by
   regulating metabolic adaptability [[142]39, [143]40]. Profile 9
   exhibits a continuous upregulation trend and is enriched in pathways
   such as oxidative phosphorylation, amino acid metabolism, and lipid
   metabolism. These pathways provide energy for the later stages of
   embryonic volume expansion and post-hatching material transformation,
   while maintaining metabolic and tissue homeostasis [[144]3, [145]33].
   The detailed KEGG enrichment results for each expression profile are
   provided in Table [146]S5.

Fig. 6.

   [147]Fig. 6
   [148]Open in a new tab

   Temporal profiles of actively expressed genes and KEGG pathway
   enrichment for each profile. The left panel shows the 9 temporal
   expression profiles identified from actively expressed genes at each
   developmental stage, with the number of genes in each profile indicated
   above the corresponding line graph, which illustrates the trends in
   gene expression over time. The right panel displays bubble charts of
   KEGG pathway enrichment analysis for each gene expression profile

   Overall, this study integrates CUT&Tag and RNA-seq data to reveal the
   dynamic changes in histone modifications and gene expression during the
   embryonic development of L. vannamei, highlighting significant
   enrichment of signaling pathways at various stages. It provides
   insights into the regulatory patterns and molecular mechanisms driving
   embryonic development in L. vannamei by linking histone modifications
   to gene expression changes.

Identification and functional characterization of key enhancer-associated
genes during L. vannamei embryonic development

   Enhancers are key cis-regulatory elements that regulate gene expression
   levels. H3K4me1 is a characteristic histone modification associated
   with enhancers [[149]12]. A genomic region marked by both H3K4me1 and
   the accessible chromatin mark H3K27ac, but lacking the repressive
   histone modification H3K27me3, is indicative of an active enhancer at a
   given developmental stage [[150]11, [151]12, [152]41]. Genes associated
   with such enhancers are likely to be in an enhanced expression state.
   Using CUT&Tag data, we identified active enhancers and their
   corresponding genes, hypothesizing that genes harboring active
   enhancers and exhibiting concurrent transcriptional activity play
   pivotal roles in embryonic development. To refine our selection of
   developmentally essential genes, we intersected enhancer-associated
   genes with actively expressed genes identified in the previous section,
   revealing key regulatory genes at each developmental stage
   (Fig. [153]7A).

Fig. 7.

   [154]Fig. 7
   [155]Open in a new tab

   Identification and KEGG functional enrichment of actively expressing
   genes with active enhancer marks. (A) Venn diagrams illustrating the
   overlap between actively expressing genes and active enhancer-marked
   genes. The blue circle represents the actively expressing genes,
   defined as genes with H3K4me3 peaks in their promoter regions,
   overlapping with H3K27ac peaks (≥ 1 bp overlap) but showing no overlap
   with H3K27me3 peaks. The pink circle represents the active
   enhancer-marked genes, defined as genes exhibiting H3K4me1 and H3K27ac
   peak overlap (≥ 1 bp) while lacking H3K27me3 overlap. The intersection
   of these two groups represents genes that are both actively expressed
   and marked by active enhancers. (B) KEGG pathway enrichment analysis of
   actively expressing genes with active enhancer marks. The significantly
   enriched pathways (p < 0.05) were visualized using a bubble plot, where
   bubble size represents the number of intersecting genes associated with
   each pathway, and color intensity reflects the statistical significance
   (P-value) of enrichment

   To elucidate the biological pathways and functional processes enriched
   in genes marked by active enhancers, we performed KEGG enrichment
   analysis. The results demonstrated a significant enrichment in
   signaling pathways, with stage-specific variations (Fig. [156]7B).
   During the gastrula and lbe stages, pathways related to cell growth and
   differentiation were significantly enriched, while at the lim stage,
   pathways associated with differentiation and the immune system were
   prominently enriched. From the N1 to N6 stages, pathways related to
   transport and catabolism were significantly enriched, which was
   consistent with the KEGG results of actively expressing genes
   (Fig. [157]7B). The characteristic biological events at each stage are
   closely associated with the pathway enrichment patterns [[158]3,
   [159]33, [160]34] Notably, during the gastrula and lbe stages,
   PRC-related pathways involved in H3K27me3 establishment were
   significantly enriched, consistent with our earlier observation of
   rapid H3K27me3 deposition during these stages (Fig. [161]3A, B).
   Additionally, gastrula-stage enrichment analysis identified the
   ATP-dependent chromatin remodeling pathway, suggesting that chromatin
   reprogramming during early embryonic development is not solely dictated
   by maternal factors but also requires zygotic gene activation. These
   findings highlight a complex interplay between histone modifications,
   chromatin remodeling, and zygotic genome activation.

   Our approach successfully identified several key regulatory genes
   associated with stage-specific developmental events in L. vannamei,
   encompassing core biological processes such as cell proliferation and
   differentiation, segmental pattern formation, tissue morphogenesis,
   molting regulation, systemic development, metabolic reprogramming,
   muscle differentiation, and cuticle formation (Fig. [162]8A-F, Table
   [163]S6). For instance, the Hox family gene [164]XM_027356713.1
   (exd-like) regulates the development of body segments and the nervous
   system, and participates in establishing the proximal, medial, and
   distal regionalization patterns of appendages [[165]37, [166]42]. Its
   promoter maintains an active open state from the lbe to N3 stages and
   is regulated by enhancers (Fig. [167]8A and D). The [168]XM_027355389.1
   (60a-like) gene controls imaginal disc development and embryonic
   morphogenesis, and its mutation leads to phenotypes in Drosophila such
   as small eyes, increased thoracic bristles, abnormal appendages, and
   female sterility [[169]43]. Its promoter remains in an active open
   state from the lim to N3 stages and is regulated by enhancers
   (Fig. [170]8B and E). The [171]XM_027365316.1 (pro-resilin-like) gene
   is involved in cuticle formation, enhancing flexibility and providing
   elastic support for joints, cuticles, and appendages [[172]44]. Its
   promoter maintains an active open state from the lim to N6 stages and
   is regulated by enhancers (Fig. [173]8C and F). Further analysis
   reveals a high correlation between the histone modification patterns of
   these key genes and their dynamic expression profiles across different
   developmental stages. These findings not only validate the
   effectiveness of our active expression gene and enhancer screening
   strategy but also demonstrate that histone modification signatures can
   serve as important molecular markers for identifying key developmental
   regulatory genes.

Fig. 8.

   [174]Fig. 8
   [175]Open in a new tab

   Representative actively expressing genes with active enhancer marks at
   certain stages. (A-C) Bar chart depicting the expression levels of key
   genes during the embryonic development of L. vannamei, quantified in
   transcripts per million (TPM). (D-F) Integrated genome viewer (IGV)
   tracks displaying RNA-seq and CUT&Tag peak profiles for RNA expression,
   H3K4me1, H3K4me3, H3K27ac, and H3K27me3 signals at these gene loci,
   illustrating their expression dynamics and histone modification
   landscapes throughout embryonic development. Control tracks are
   included for comparative analysis, providing a reference for background
   signal levels

   Through an integrative analysis of histone modifications and gene
   expression dynamics, this study reveals how histone modifications
   cooperate with transcriptional regulation to drive embryonic
   development. This systematic approach not only enables the precise
   identification of developmentally essential genes but also offers novel
   insights into the molecular mechanisms underlying early embryogenesis
   in L. vannamei.

Discussion

   Embryonic development is initially governed by maternal gene products,
   and control is gradually transferred to the zygotic genome during
   zygotic genome activation (ZGA). This transition is tightly regulated
   through dynamic changes in histone modifications, which have been
   extensively studied in model organisms but remain largely unexplored in
   crustaceans [[176]6]. In this study, we provide the first
   high-resolution, genome-wide maps of four histone modifications
   (H3K4me1, H3K4me3, H3K27ac, and H3K27me3) during embryonic development
   in a crustacean species. By integrating CUT&Tag data with
   transcriptomic profiles, we systematically identified cis-regulatory
   elements and their stage-specific activity, thus offering novel
   insights into the epigenetic regulation of zygotic genome activation
   (ZGA), cell lineage specification, and morphogenesis in L. vannamei
   embryos.

Histone modifications and ZGA in L. vannamei

   Although the timing and patterns of histone modification reprogramming
   and ZGA vary among species, these processes are evolutionarily
   conserved in sexually reproducing animals [[177]31, [178]45, [179]46].
   Since germ cells are highly specialized, they must regain cellular
   totipotency after fertilization to enable proper embryonic development
   and tissue differentiation.

   Our findings suggest that ZGA in L. vannamei likely occurs before or
   around the blastula stage. This conclusion is supported by the
   significant increase in expressed genes from the blastula to gastrula
   stages (Fig. [180]3D), suggesting large-scale zygotic genome
   activation. Additionally, the presence of active histone modifications
   in the blastula stage and their substantial upregulation in the
   gastrula stage further corroborate this (Fig. [181]3B). Due to the
   reliance of L. vannamei on external fertilization mechanisms (where
   males transfer spermatophores to the female’s thelycum during mating,
   and sperm-egg binding occurs externally upon ovulation) [[182]47],
   coupled with the rapid transition from multicellularization to the
   blastula stage within minutes post-fertilization, the developmental
   process is highly sensitive to environmental fluctuations such as
   temperature and salinity variations [[183]44, [184]47, [185]48]. This
   significantly exacerbates asynchrony in embryonic development prior to
   the blastula stage. Consequently, obtaining earlier-stage samples to
   precisely determine the timing of the maternal-to-zygotic transition
   (MZT) remains challenging.

   Previous studies in model species and humans already showed tremendous
   difference in the timing of MZT and dynamics of histone modifications
   during embryogenesis across a broad range of species [[186]6, [187]49].
   Comparative analysis of L. vannamei with other marine invertebrates
   also highlights interspecies differences in MZT timing and histone
   modification dynamics. In Mulinia lateralis, MZT occurs between the
   morula and gastrula stages, with weak histone modification signals
   before the morula stage. H3K27me3 is established first, while H3K4me1
   and H3K4me3 are significantly upregulated only until the gastrula stage
   (our unpublished data). In contrast, Crassostrea gigas exhibits high
   levels of all four histone modifications as early as the zygote stage
   [[188]50]. These findings suggest that, despite evolutionary proximity,
   species-specific differences exist in histone modification
   reprogramming and MZT timing. Future studies incorporating more species
   and higher-resolution techniques such as single-cell CUT&Tag could
   provide deeper insights into chromatin accessibility and
   cell-type-specific regulatory mechanisms during early development.

Regulation of key developmental processes

   Our study identified significant enrichment of H3K4me1, H3K27me3, and
   H3K27ac at the N1 stage, while H3K4me3 levels peaked at the lim stage.
   This suggests that histone modification dynamics are closely linked to
   the developmental needs of each stage. The increased H3K4me3 during the
   lim stage may reflect extensive tissue differentiation as the embryo
   prepares to hatch, while the concurrent enrichment of activating
   (H3K4me1, H3K27ac) and repressive (H3K27me3) marks at the N1 stage
   suggests a complex interplay between activation and repression during
   metamorphosis [[189]3, [190]33].

   The time-series analysis of histone marks and the transcriptome
   revealed a high degree of temporal correlation between the two,
   indicating that histone modifications and transcriptional activities
   collaboratively drive gene regulatory processes during early embryonic
   development. This underscores the significance of epigenetic mechanisms
   in orchestrating developmental transitions.

   L. vannamei exhibits unique stages during early embryonic development,
   undergoing significant physiological changes, particularly in tissue
   formation, dietary transitions, muscle development, and exoskeleton
   remodeling. However, current research on early embryonic development of
   L. vannamei mainly focuses on morphological, biochemical, and
   transcriptomic characteristics, with limited studies on epigenetic
   aspects, particularly histone modifications [[191]3, [192]33, [193]34].
   In this study, we utilized CUT&Tag and RNA-seq analysis for the first
   time to explore potential molecular mechanisms regulating physiological
   changes during early embryonic development from the perspective of
   histone modification-mediated gene expression.

Identification of key developmental genes and their cis-regulatory elements

   Through a comprehensive analysis of histone modifications and
   transcriptomic data, we identified genes that are actively expressing
   at each stage and further screened out those with active enhancer
   marks. This approach not only circumvents the interference of maternal
   mRNA fluctuations in transcriptomic analysis but also enables precise
   identification of developmental regulatory genes and their networks
   through epigenetic features, providing new insights into the molecular
   mechanisms underlying embryonic development in invertebrates.

   By defining “active open regions” (genomic regions marked by both
   H3K4me3 and H3K27ac activation marks but lacking the repressive
   H3K27me3 modification) and integrating RNA-seq expression data
   (TPM > 0.5), we successfully identified epigenetically regulated
   actively expressed genes across various stages of embryonic development
   in L. vannamei (Fig. [194]5A). To further elucidate the regulatory
   mechanisms of key genes, we screened for genes marked by active
   enhancers by combining the co-localization signals of the enhancer
   marker H3K4me1 and the open chromatin marker H3K27ac (Fig. [195]7A).
   These genes not only exhibited significant transcriptional activity but
   also showed a strong correlation between the dynamics of histone
   modifications in their enhancer regions and the stage-specific
   functional requirements. For instance, the promoter regions of the
   segmental and nervous system development regulator gene exd-like
   [[196]37, [197]42] and the cuticle formation gene pro-resilin-like
   [[198]44] maintained an active chromatin state from the lim to the
   nauplius stages, with their enhancer modification levels significantly
   upregulated during development, showing high synchronization with gene
   expression dynamics (Fig. [199]8). This finding indicates that
   epigenetic regulation and transcriptional activation collaboratively
   and precisely drive the execution of developmental programs.

   This study represents the first successful localization of key
   developmental genes and their cis-regulatory elements in L. vannamei
   through dynamic histone modification profiling. This strategy not only
   validates the feasibility of histone modifications as predictive
   markers of gene expression but also provides valuable candidate genes
   and regulatory element resources for future functional genomics
   research [[200]11, [201]12].

   In addition, we conducted an in-depth analysis of the chromatin states
   of L. vannamei using the open-source software package ChromHMM
   [[202]22]. By annotating chromatin states based on histone
   modifications, our study provides a valuable resource for future basic
   biological research on shrimp embryonic development. However, it is
   worth noting that although our integrated CUT&Tag and RNA-seq analyses
   revealed a strong correlation between predicted regulatory elements and
   gene expression patterns, cellular heterogeneity inherent in
   whole-embryo samples may introduce confounding signals into chromatin
   state annotations. In future studies, the application of single-cell
   CUT&Tag or spatially resolved epigenomic techniques will be essential
   for generating more precise, stage- and tissue-specific chromatin state
   maps.

Conclusions

   This study provides new epigenetic resources for the economically
   significant crustacean species L. vannamei. Notably, this is the first
   successful application of CUT&Tag in L. vannamei. By analyzing histone
   marks (H3K4me1, H3K4me3, H3K27me3, and H3K27ac) and corresponding
   transcriptomes at seven critical stages of early embryonic development,
   we revealed that dynamic changes in histone modifications coincide
   closely with transcriptional activation during the ZGA of L. vannamei.
   This correlation suggests that histone modifications and
   transcriptional activity jointly drive gene regulation during early
   embryonic development, offering new insights into the ZGA process in
   shrimp. Additionally, through differential expression analysis, and
   integrated analysis of histone modification and transcriptome data, we
   explored potential molecular mechanisms underlying physiological
   changes during early embryonic development mediated by histone
   modifications. This study identified a series of potential key genes
   involved in development, providing new insights from the perspective of
   histone modification-mediated gene expression. This represents the
   first characterization of histone modifications during early
   development in shrimp. The findings of this study will serve as a
   significant reference for research in shrimp developmental biology and
   aquaculture.

Supplementary Information

   Below is the link to the electronic supplementary material.
   [203]Supplementary Material 1^ (9.4KB, xlsx)
   [204]Supplementary Material 2^ (11.2KB, xlsx)
   [205]Supplementary Material 3^ (11.1KB, xlsx)
   [206]Supplementary Material 4^ (13KB, xlsx)
   [207]Supplementary Material 5^ (29.4KB, xlsx)
   [208]Supplementary Material 6^ (14.1KB, xlsx)
   [209]Supplementary Material 7^ (2MB, docx)

Acknowledgements