- Research
- Open access
- Published:
Differences in alternative splicing events in the adaptive strategies of two Daphnia galeata genotypes induced by fish kairomones
BMC Genomics volume 25, Article number: 725 (2024)
Abstract
Background
Daphnia galeata is a suitable model organism for investigating predator-induced defense. Genes and pathways exhibiting differential expression between fish kairomone-treated and untreated groups in D. galeata have been identified. However, understanding of the significance of alternative splicing, a crucial process of the regulation of gene expression in eukaryotes, to this mechanism remains limited. This study measured life-history traits and conducted short-read RNA sequencing and long-read isoform sequencing of two Korean D. galeata genotypes (KB1 and KE2) to uncover the genetic mechanism underlying their phenotypic plasticity under predation stress.
Results
KB1 exhibited strategies to enhance fertility and decrease body length when exposed to fish kairomones, while KE2 deployed an adaptive strategy to increase body length. Full-length transcriptomes from KB1 and KE2 yielded 65,736 and 57,437 transcripts, respectively, of which 32 differentially expressed transcripts (DETs) were shared under predation stress across both genotypes. Prominent DETs common to both genotypes were related to energy metabolism and the immune system. Additionally, differential alternative splicing (DAS) events were detected in both genotypes in response to fish kairomones. DAS genes shared between both genotypes may indicate their significant role in the post-transcriptional stress response to fish predation. Calpain-3, involved in digestion and nutrient absorption, was identified as a DAS gene in both genotypes when exposed to fish kairomones. In addition, the gene encoding thymosin beta, which is related to growth, was found to be a statistically significant DAS only in KB1, while that encoding ultraspiracle protein, also associated with growth, was only identified in KE2. Moreover, transcripts encoding proteins such as EGF-like domain-containing protein, vitellogenin fused with superoxide dismutase, and others were identified overlapping between DAS events and DETs and potentially elucidating their association with the observed phenotypic variation in each genotype.
Conclusions
Our findings highlight the crucial role of alternative splicing in modulating transcriptome landscape under predation stress in D. galeata, emphasizing the requirement for integrating gene expression and splicing analyses in evolutionary adaptation studies.
Background
Organisms interact with various biotic and abiotic environmental factors through phenotypic plasticity, which is evidenced by behavior or physiological changes, and adapt to environmental changes [1]. Phenotypic plasticity refers to an organism’s ability to produce various phenotypes in response to environmental conditions, ultimately aiding their survival and reproduction across various habitats [2]. These mechanisms can enhance population persistence and facilitate adaptations to novel and shifting conditions by mitigating the impact of selection on changing environments [3]. Predation is a key biotic factor in structuring whole communities, maintaining species diversity, and driving natural selection in populations [4, 5]. Aquatic predators including vertebrates and invertebrates, release kairomones into the water that trigger diverse responses after being detected by their prey, ultimately reducing prey vulnerability [6, 7]. These predator-induced responses are a typical example of phenotypic plasticity and have been extensively reported in variety of Daphnia species [8, 9].
Daphnia is an extensively employed model organism in ecology, evolution, and ecotoxicology, due to its capacity to develop distinct alternative phenotypes in response to environmental signals [10, 11]. Members of the genus Daphnia play a crucial role in linking trophic levels, bridging primary producers with consumers in freshwater ecosystems, which renders them susceptible to heightened predation pressure [10, 12]. Upon exposure to fish kairomones, variations in life-history traits and morphological, physiological, and behavioral plasticity, along with other effects, have been reported [13,14,15]. Daphnia galeata, among Daphnia species, is a relatively large zooplankton discovered across the northern hemisphere and recognized as a dominant plankton species in the Han River, Korea [16]. While this species does not show significant morphological changes or diel vertical migration in response to vertebrate predator cues, it displays noticeable phenotypic variations in life-history traits under predation risk [17]. For example, D. galeata shows an adaptive strategy for ensuring generational survival by producing more offspring when exposed to fish kairomones [15, 18]. Additionally, Daphnia is found in diverse environments across a wide ecological condition with high degree of genetic variation in its populations [19]. Genetic divergence between Daphnia genotypes led to diverse phenotypic responses to stressors [20]. Therefore, it is not possible to generalize the impact of stressors on monoclonal populations to the species level, as genetically diverse populations will exhibit variations in stress tolerance and response [21]. Our recent study revealed that fish kairomone-treated and untreated groups of D. galeata shared differentially expressed transcripts (DETs) among genotypes, although each genotype displayed different genes that exhibited altered expression under fish predation exposure [18]. This implies that stress response mechanisms might be quite uniform within a species, although genotypes with varying genetic backgrounds could employ distinct strategies for adapting to environmental changes.
Recent transcriptome studies by RNA sequencing (RNA-seq) have offered crucial insights into the molecular mechanisms underlying transcriptional responses to fish predation. The alteration of pre-mRNA splicing represents another significant source of phenotypic variation [22]. However, limited research has been conducted on the contribution of other mechanisms that regulate gene expression through RNA processing, such as alternative splicing (AS), a eukaryotic regulatory mechanism which selectively processes exons and introns in the pre-mRNA transcript. AS represents a fundamental mechanism regulating gene expression, contributing to transcriptional diversity and activity in eukaryotes [23]. Studies have revealed extensive AS events in which > 50% of Daphnia pulex and Daphnia magna genes undergo this regulatory process [24]. Notably, the Dscam gene exemplifies the impact of AS, generating at least 13,000 potential mRNA transcripts in Daphnia [25]. The main types of AS events include alternative first exon (AF), retained intron (RI), alternative 3′ splice site (A3), alternative 5′ splice site (A5), skipped exon (SE), alternative last exon (AL), and mutually exclusive exon (MX) [26].
RNA-seq, as a high-throughput next-generation sequencing technique, has become an essential tool for transcriptome-wide analysis of differential gene expression and gene regulatory networks. However, the limitations of RNA-seq, such as its short read lengths and amplification biases, restrict the acquisition of accurate full-length transcripts and information on differential mRNA splicing, especially for species that lack high-quality reference transcriptomes [27]. Recently, Isoform sequencing (ISO-seq) has been used to analyze the full-length transcriptomes of various organisms [28, 29]. Additionally, in previous studies, the integrative analysis of ISO-eq and RNA-seq data could accurately quantify transcript abundance [29]. However, ISO-seq comes with higher costs and error rates compared to RNA-seq, as its results require correction and validation using high-throughput, high-accuracy short-read data from RNA-seq [29]. Therefore, the combination of the two sequencing strategies can lead to a comprehensive characterization of transcriptomes and could accurately identify post-transcriptional modification events and DE isoforms [28].
A previous study reported the important role of AS in shaping transcriptome diversity in response to metal exposure in D. pulex [24]. Nonetheless, there is a scarcity of analyses examining the proportional impact of AS on gene expression regulation across Daphnia species, which serves as a prominent model organism in ecological genomics. Differential splicing analysis has the capacity to reveal functional diversity that cannot be detected solely through differential gene expression analysis. In elucidating the genes and molecular pathways implicated in stress responses, therefore, it serves as a valuable complement to differential gene expression analysis [24]. Furthermore, the ability of AS patterns to vary among individuals indicates that splicing is a dynamic mechanism that can influence stress responses in a genotype-specific manner [30]. Therefore, investigating the genetic basis of differential gene expression and splicing among D. galeata genotypes can help distinguish common stress response pathways from genotype-specific responses to fish predation exposure. In this study, we integrated differential expression (DE) and differential splicing analyses to identify the role of AS in mediating responses to fish kairomone-induced stress in D. galeata. This study is the first to report the application of ISO-seq coupled with RNA-seq in D. galeata, advancing our understanding of the biological significance of AS events in this species and its impact on transcriptome diversity in response to fish kairomone exposure.
Methods
Collection of samples and processing
This study used two D. galeata genotypes collected from two different environments, Paldang Lake and Baekhak Reservoir, in Gyeonggi Province, Korea, in June 2022, using a 55 μm mesh-sized Kitahara Zooplankton net (Samjee-tech, Korea). The collected samples were transferred to the laboratory and each egg-bearing female was transferred into a separate beaker for parthenogenetic reproduction. D. galeata individuals were cultured in controlled laboratory conditions (20 °C, 16 h light/8 h dark cycle, ADaM [Aachener Daphnien Medium]), and were fed with 1.0 mg C L-1 Chlorella vulgaris once daily. The two genotypes examined were determined through genotyping using 13 microsatellite markers by polymerase chain reaction (PCR) amplification in three individuals per genotype (Table S1) [31]. Thirteen fluorescein amidite (FAM) labeled primer sets were designed to target these microsatellite loci. The total 25 μL PCR products contained 2.5 μL of 10X Pfu Reaction Bufer, 0.5 μL of dNTP mix (10 mM), 0.25 μL of Add Pfu DNA Polymerase (Addbiomeditech, Daejeon, Korea), 1 μL of template DNA, 1 μL of each primer (10 pmol/μL), and 18.75 μL of distilled water. The amplification protocol consisted of an initial denaturation at 95 °C for 2 min, followed by 35 cycles of denaturation at 95 °C for 30 s, annealing at 52 °C for 1 min, extension at 72 °C for 1 min, and a final elongation step at 72 °C for 5 min. The PCR products were analyzed using 3730xl DNA Analyzer (Applied Biosystems, Foster City, USA) and genotyped using GeneMapper 4.0 (Applied Biosystems, Foster City, USA). As a result, we identified two different genotypes with sequence differences in 13 microsatellite markers (Table S2).
Experimental set-up
This study utilized two D. galeata genotypes (KB1 and KE2) that were cultured for one year prior to the experiment and had established clonal lines. The experimental design and procedures were derived from previous studies [18, 20]. In short, D. galeata was exposed to kairomone-enriched medium, obtained from aquaria housing freshwater mandarin fish (Siniperca scherzeri), to simulate fish predation. The kairomone-containing medium, filtered and mixed with ADaM medium, was optimized for the survival of D. galeata while inducing changes in life-history traits. Both genotypes were raised for two consecutive generations in kairomone-free medium (control group) and kairomone-containing medium (experimental group) before the experiment to minimize inter-individual variability. The experiment was conducted for 21 days for each experimental individual to investigate the long-term effects of fish kairomones. During this period, we documented reproductive life-history traits including “Age at First Reproduction,” “Number of Offspring First Brood,” “Total Number of Broods,” and “Total Number of Offspring.”. We also assessed growth-related traits such as “Somatic Growth Rate” and “Body Length.”. The "Somatic Growth Rate" of D. galeata was calculated by dividing the change in body length of neonates at the start (t0) to adults at the end (t21) of the experiment by the total duration of the experiment. The length of D. galeata's body was measured from the top of the head, through the center of the eye, to the base of the spine [20]. The individuals were examined with a CX22LED microscope (Olympus, Japan), and captured with an eXcope T500 microscope camera (DIXI Science, Korea). Differences in life-history traits between the presence and absence of fish kairomones in each genotype was compared by conducting one-way ANOVA statistical analyses using the aov() function in R version 4.2.2 [32, 33]. Samples from both genotypes were collected to perform ISO-seq and RNA-seq. For RNA-seq, egg-bearing females were pooled per sample (n = 20) after 21 days, resulting in three biological replicates for each experimental condition (control and fish kairomone-exposed) in the two genotypes, totaling 240 individuals (two genotypes × two conditions × 20 individuals × three biological replicates). For ISO-seq, only 20 female individuals from each genotype were pooled after exposure to fish kairomones. Afterward, the samples were stored in Qiagen RNAprotect tissue reagent (Qiagen, Valencia, CA, USA) at 4ºC overnight and then kept at − 80ºC until RNA extraction as per the manufacturer instructions.
RNA extraction, library construction, and sequencing
Since we could not obtain sufficient amounts of RNA from a single individual, we pooled the 20 experimental individuals per sample. Sampling females in the inter-molt stage ensures stable gene expression; therefore, only egg-bearing experimental females were pooled. [34]. Total RNA was isolated using the Qiagen RNeasy Plus Universal Mini Kit (Qiagen, Valencia, CA, USA) following the manufacturer’s instructions after homogenization with a disposable pestle and homogenizer. RNA quality was evaluated with an Agilent 2100 bioanalyzer (Agilent Technologies, Amstelveen, The Netherlands), and RNA quantification was conducted with an ND-2000 Spectrophotometer (Thermo Inc., DE, USA). Total RNA samples with RIN greater than 7 were used for constructing the cDNA libraries for ISO-seq and RNA-seq.
The RNA-seq library was constructed from the total RNA using the NEBNext Ultra II Directional RNA-seq Kit (New England BioLabs, Inc., UK). Briefly, The Poly(A) RNA Selection Kit was used to isolate mRNA (Lexogen, Inc., Austria). The isolated mRNAs were then utilized for cDNA synthesis and shearing according to the manufacturer’s instructions. Indexing was performed utilizing Illumina indexes 1–12, followed by PCR enrichment. Subsequently, the libraries were evaluated using the TapeStation HS D1000 Screen Tape (Agilent Technologies, Amstelveen, The Netherlands) to ascertain the mean fragment size. The library quantification kit on a StepOne Real-Time PCR System (Life Technologies, Inc., USA) was utilized for library quantification. Three libraries of the KB1 and KE2 genotypes were constructed and sequenced on the Illumina NovaSeq 6000 (Illumina, Inc., USA) platform to generate 101 bp paired-end reads. For library preparation for ISO-seq, 300 ng of RNA was used for full-length transcript sequencing with a Pacbio Sequel IIe system (Pacific Biosciences, CA, USA) according to the manufacturer’s instructions. The ISO-seq library was prepared following the ISO-seq protocol using the NEBNext Single Cell/Low Input cDNA Synthesis & Amplification Module, SMRTbell prep kit 3.0. A 10 μL library was prepared with the SMRTbell Prep Kit 3.0, followed by the annealing of SMRTbell templates using the Sequel II Binding Kit 3.1. The Sequel II sequencing kit 2.0 and SMRT cells 8 M Tray were used for sequencing. Finally, SMRT cells using 24-h movies were sequenced for each SMRT cell using the PacBio Sequel IIe (Pacific Biosciences) sequencing platform by Macrogen (Seoul, Korea). The following steps followed the PacBio Sample Net-Shared Protocol, accessible at https://www.pacb.com/.
ISO-seq long-read data analysis
A detailed overview of the experimental setup is provided in Fig. 1. Raw Pacbio sequencing data (hifi_reads.bam) was generated by performing circular consensus sequences (CCS) Analysis (–min-rq 0.9) on Sequel IIe in HiFi mode. The raw data generated from each cell was processed individually using the following procedures: 1) demultiplexing was performed using the Demultiplexing Barcodes utility on SMRTLink 11.1.0.166339 software; 2) cDNA primer and the spurious false positive signal analysis was carried out using lima (version 2.6.99) with the parameters –isoseq –peekguess; 3) PolyA tail trimming and concatemer removal were performed to generate full-length non-concatemer (FLNC) reads using isoseq3 (version 3.8.1) refined with the parameters –require-polya; 4) these FLNC reads were clustered to obtain consensus isoforms, and nonredundant transcript sequences were obtained after further clustering using isoseq3 cluster with the parameters –verbose –use-qvs. Thus, high-quality and polished full-length consensus sequences were produced. These polished consensus sequences were further subjected to correct and remove redundancy with the RNA-seq data using LoRDEC software to obtain high-quality (HQ) transcripts [35].
Collapsing redundant transcript isoforms and assessment of completeness
A two-step collapsing strategy [28] was implemented to eliminate redundant transcript isoforms. The HQ isoforms were first aligned and sorted to generate SAM files using minimap2 (version 2.18) with the parameters -ax splice, -uf, -C5, and –secondary = no, using the genome of D. galeata as reference genome (GCA_918697745.1) [36, 37]. The redundant isoforms were collapsed based on the alignment results using the cDNA cupcake software with the following parameters: –min_aln_coverage 0.95, –min_aln_identity 0.85, and –dun-merge-5-shorter [38]. Furthermore, unmapped transcripts were also collapsed using Cogent v8.0 software and cDNA cupcake software with default parameters [39]. In this process, different gene families were initially identified from unmapped transcripts. A “fake genome” was produced by concatenating all cogent unassigned contigs. Using it as a reference genome, we collapsed these unmapped transcripts using cDNA cupcake software as described in the previous steps. The fasta and gff files from two processes were merged into one file. Finally, highly similar sequences with 99% identity were removed from the merged transcript isoforms using CD-HIT-EST (version 4.8.1) and subsequently generated full-length transcriptomes for the two genotypes for further analysis [40]. The completeness of the full-length transcriptomes of KB1 and KE2 were assessed using BUSCO v5.1.3 software against the Arthropoda lineages’ gene sets (arthropoda_odb10) [41]. After assessing completeness, full-length transcripts were categorized by comparing them to the reference genome annotation using gffcompare software (version 0.12.2) [42].
Functional annotation of full-length transcripts
To predict potential coding regions and open reading frames (ORFs) of full-length transcripts, the TransDecoder pipeline with default parameters was used [43]. The resulting ORFs were searched against three databases (UniProt, TrEMBL, and Non-redundant (NR)) using the BLASTp function of diamond (version 2.1.8) with the following parameters: –evalue 1e-5, –max-target-seqs 1, and –outfmt 6 [44,45,46,47]. In addition, the resulting ORFs were identified using eggNOG-mapper (version 2.1.11) with –evalue 1e-5 parameters to obtain functional annotation information, including that from the Clusters of Orthologous Groups (COG) database, the Gene Ontology (GO) knowledgebase, the Kyoto Encyclopedia of Genes and Genomes (KEGG), and the protein families database Pfam [48]. The accession numbers of seed_ortholog obtained from the eggNOG database were converted to UniProt IDs. The top hit for each query sequence was chosen for transcriptome annotation and subsequent characterization.
RNA-seq short-read data analysis and differential transcript expression analysis
The quality control of raw RNA sequencing data was performed using FastQC (version 0.11.9) [49]. Adapter and low-quality reads were removed using Trimmomatic (version 0.36) with the default parameters, except for the TRAILING parameter set to 20 [50]. After trimming, the read quality was re-evaluated using FastQC to confirm the removal of adapters. To quantify the expression of clean reads, the align_and_estimate_abundance.pl script from the Trinity pipeline was utilized. All clean data from RNA-seq were aligned to the full-length transcriptomes quickly and accurately using RNA-seq by Expectation Maximization (RSEM) [51]. DE analysis was conducted at the isoform level using the DESeq2 package for R v4.1.3 within the Trinity pipeline [32, 52]. The results were filtered using adjusted p-value (padj) < 0.05 and a fold change of ≥ 4 (log2FC ≤ − 2 or log2FC ≥ 2) [53]. Finally, GO and KEGG pathway enrichment analyses were performed to determine the functions of DETs. Ortholog IDs of DETs were extracted from the results of eggNOG-mapper for GO and KEGG enrichment analyses, which were carried out using ShinyGO (version 0.61) [54]. Categories with corrected p-values < 0.05 were selected for further analysis. To verify the quality of the triplicated samples, Principal component analysis (PCA) was conducted using DEseq2 and ggplot2 [55]. In addition, the expression of DETs was visualized as a volcano plot using the R package ggplot2.
Differential alternative splicing (DAS) event analysis
A GTF annotation file obtained from genotypes KB1 and KE2 was used to investigate the differences in AS events in response to fish kairomones using SUPPA2 software [56]. Using the generateEvents function of SUPPA2 software with default parameters, seven AS event types were produced from the GTF annotation file, including AF (alternative first exon), RI (retained intron), A3 (alternative 3′ splice site), A5 (alternative 5′ splice site), SE (skipped exon), AL (alternative last exon), and MX (mutually exclusive exon) (Figure S1). To discover the DAS events in both genotypes, the inclusion levels (PSIs) per AS event were determined by the psiPerEvent function of SUPPA2 software using the results of transcripts per million reads (TPM) and AS events. The TPM value of each sample was generated using the align_and_estimate_abundance.pl script from the Trinity pipeline and then used to calculate additional inclusion levels (PSIs). Furthermore, significant DE analysis was performed using the diffSplice function of SUPPA2 software using the criteria |∆PSI|≥ 0.15 and p-value < 0.05. Splicing plots were obtained using ggsashimi (version 1.1.5) software [57]. Short-read RNA-seq data were utilized to validate AS events and splice sites. Finally, to determine the functions of significant DAS events, GO and KEGG pathway enrichment analyses were performed using ShinyGO (version 0.61) [54]. Only enriched categories showing corrected p-values < 0.05 were selected. Additionally, a Venn diagram was utilized to compare the commonalities and differences between DETs and DAS genes identified in both genotypes (https://bioinformatics.psb.ugent.be/webtools/Venn/).
Results
Varied life-history strategies under fish kairomones stress
To elucidate the genotype-specific responses to predator stress, six life-history traits linked to reproduction and growth were measured in two D. galeata genotypes. Reaction norms of life-history traits showed their intraspecific variation between the two D. galeata genotypes (Fig. 2). A one-way ANOVAs revealed that “fish kairomones” had a significant effect (P < 0.05) on the six life-history traits, with the exception of “Somatic Growth Rate” in genotype KB1 (Fig. 2). In contrast, fish kairomones had a significant effect (P < 0.05) only on the “Number of Offspring First Brood” and “Body Length” in genotype KE2 (Fig. 2B and F). Interestingly, contrasting patterns in terms of life-history traits associated with growth were exhibited by the two genotypes. Regarding the life-history traits related to growth, the “Body Length” of KE2 significantly increased (P < 0.01), while that of KB1 significantly decreased (P < 0.05), in the presence of fish kairomones (Fig. 2F). In summary, genotype KE2 showed an adaptive strategy by increasing its body length rather than improving fertility. However, genotype KB1 showed the adaptive strategy of improved fertility to sustain generational survival, producing many offspring, and reducing its body length to avoid being discovered and attacked by predators.
Summary of ISO-seq data
A total of 1,891,084 and 1,632,379 subreads were recovered for control and experimental groups of KB1, along with 1,477,802 and 1,611,476 subreads for control and experimental groups of KE2, respectively. These subreads yielded 1,724,362 and 1,475,596 CCS for the control and experimental groups of KB1, along with 1,346,904 and 1,464,581 CCS for the control and experimental groups of KE2, respectively. After obtaining the FLCN reads from each condition via IsoSeq3 refinement, we merged them by genotype. Finally, after clustering and polishing steps, a total of 100,250 (28) and 85,882 (24) HQ (low-quality) isoforms were obtained for KB1 and KE2 genotypes, respectively. We excluded low-quality isoforms from downstream analysis (Table S3).
After processing ISO-seq data, RNA-seq data was employed to correct and remove redundancy in the polished sequences. RNA-seq data is summarized in Table S4. Based on the results of reference genome mapping, redundant isoforms were collapsed, generating 61,967 and 54,222 unique isoforms in KB1 and KE2, respectively (Table S5). The remaining unmapped sequences were utilized to construct a "fake genome" to collapse redundant isoforms, yielding 4,979 and 5,070 unique isoforms in KB1 and KE2, respectively (Table S5). Finally, all unique isoforms were merged, and redundant isoforms were subsequently collapsed using CD-HIT. This process yielded 65,736 isoforms for KB1 and 57,437 isoforms for KE2, with average lengths of 2027.5 bp and 2053.5 bp, respectively (Table S5).
The completeness and quality of the KB1 and KE2 transcriptomes are important prerequisites for further analyses. Of the 1,013 BUSCO groups searched, 98.7% and 98.4% of complete orthologs belonging to the phylum Arthropoda were detected in KB1 and KE2 genotypes, respectively (Table S6) and includes complete single copies and duplicates (putative paralogs or complete genes with multiple copies). In comparison to the reference genome annotation, a total of 23,736 and 21,363 potentially novel isoforms (coded as j) and a total of 7,997 and 6,645 unknown isoforms (coded as u) were respectively annotated in the full-length transcriptomes of KB1 and KE2, showing that the full-length transcriptomes contain more novel isoforms (Figure S2).
Identification of coding regions and functional annotation
Functional annotation was carried out for both genotypes to elucidate the biological context of the transcripts. A total of 65,736 and 57,437 transcripts were used as query sequences for functional annotation, with TransDecoder identifying 60,606 and 53,576 transcripts having the longest ORFs in KB1 and KE2, respectively. Among the transcripts from both genotypes, a BLASTp-based homology search against the NR database revealed the highest annotation rates, with 59,768 (98.62%) and 52,626 (98.23%) transcripts annotated in the KB1 and KE2 genotypes, respectively (Table S7). The number of transcript isoforms annotated in Uniprot, TrEMBL, and eggNOG database were 43,218/38,593, 59,738/52,848, and 54,671/48,540 for KB1/KE2, respectively (Table S7). Overall, considering various databases, a total of 59,794 and 52,877 transcripts were respectively annotated in KB1 and KE2.
Differential expression analysis under fish kairomone exposure
PCA was conducted to verify the deviation of triplicated samples by visualizing variations in transcript expression (Fig. 3A and B). In genotype KB1, the first principal component (PC1) accounted for 44% of the variance, whereas PC2 explained 23% of the variance (Fig. 3A). There was no distinct clustering observed in the read counts by condition, which could be attributed to effects related to the variability among triplicate samples. In genotype KE2, the control condition showed clear clustering, while the fish-kairomone treatment condition exhibited variability among triplicate samples (Fig. 3B).
A total of 855 and 227 DETs were identified between samples exposed to fish kairomones versus untreated samples of KB1 and KE2, respectively (Fig. 3C and D). Among these, 691 DETs from KB1 and 191 from KE2 were annotated against the eggNOG database. This resulted in 426 DETs (221 up- and 205 down-regulated) in KB1, and 144 DETs (73 up- and 71 down-regulated) in KE2 (Table S8). Afterward, we analyzed the association between phenotypic variation and genotype-specific responses to fish kairomone exposure using a Venn diagram (Figure S3), which shows that the number of DETs may differ from those mentioned above due to the annotation of multiple transcripts with the same protein name. A total of 32 DETs were common to both genotypes, with 20 exhibiting consistent expression profiles (e.g., up- or downregulation) in both genotypes (Figure S3 and Table S9). Consistently up-regulated DETs encoded acetyl-CoA carboxylase, homeodomain transcription factor, and others. In contrast, consistently downregulated DETs encoded long-chain-fatty-acid–CoA ligase, a C-type lectin domain-containing protein, and others. Second, the remaining 12 DETs showed inconsistent expression (i.e., they exhibited contradictory expression patterns in both genotypes; Table S9). For example, NAD( +) kinase was upregulated in KB1 but downregulated in KE2. In contrast, muscle-specific actin and others were downregulated in KB1 but upregulated in KE2. Finally, to associate the DETs exclusively found in each genotype with observed phenotypic variation, we compiled a list of 330 (182 up- and 148 down-regulated) and 90 (57 up- and 33 down-regulated) unique DETs in KB1 and KE2, respectively (Table S9). Up-regulated DETs related to reproduction were identified only in KB1, including transcripts encoding a phosphopyruvate hydratase, H( +)-transporting two-sector ATPase, NADH-ubiquinone oxidoreductase, and others involved in energy metabolism. In addition, downregulated DETs related to growth, such as those encoding a receptor protein-tyrosine kinase, delta-like protein, and others, were identified only in KB1. Conversely, upregulated DETs related to growth, such as that encoding a CRIB domain-containing protein, were identified in KE2 and not in KB1.
Investigation of differential alternative splicing events in response to fish kairomone exposure
To explore differences in AS events between the control (without fish kairomones) and experimental (fish kairomone-exposed) groups, AS analyses were performed in both genotypes. The total number of AS events was 19,682 and 17,494 in KB1 and KE2 genotypes, respectively. The distribution of AS types was similar across the genotypes, with AF being the most abundant (33.33%–34.9%), followed by RI (27.9%–29.9%), A3 (12.3%–12.7%), A5 (11.7%–11.9%), SE (7.1%–7.6%), AL (4.0%–4.7%), and MX (1.0%) (Table S10). In addition, DAS events were examined to shed light on the relationship between phenotypic variation and genotype-specific responses to fish kairomones in both genotypes. Comparisons of the AS events between control and experimental groups identified 630 (294 up- and 336 down-regulated) and 335 (192 up- and 143 down-regulated) significantly differentially spliced genes in KB1 and KE2, respectively. Among these, 572 and 310 DAS genes were annotated against the eggNOG database in KB1 and KE2, respectively. If DAS events were identified a given gene and their event_IDs showed encoded the same protein but showed conflicting ∆PSI values (e.g., one event_ID had a ∆PSI greater than 0.15, but another had a ∆PSI less than 0.15) that gene was excluded. This resulted in 363 DAS genes (182 positively- and 181 negatively-spliced) in KB1, and 219 DAS genes (121 up- and 98 down-regulated) in KE2 (Table S11). A total of 23 common DAS genes were identified in both genotypes when exposed to fish kairomones (Table S12). Among the common DAS genes, the gene encoding calpain-3, involved in digestion and nutrient absorption, underwent AF splicing (Fig. 4A and B). Interestingly, this gene displayed an increase in exon inclusion (∆PSI: 0.46, occurring in CAKKLH010000325.1:1,894,575–1,894,925:1,895,154:1,894,575–1,897,881:1,898,336:-) in KB1 but a decrease in exon inclusion (∆PSI: − 0.22, occurring in CAKKLH010000325.1:1,894,575–1,894,929:1,895,145:1,894,575–1,897,881:1,898,603:-) in KE2 under fish kairomone exposure (Table S12). In addition, growth-related thymosin beta was found to display statistically significant DAS (∆PSI: − 0.27, occurring in CAKKLH010000342.1:111,977–112,421:111,929–112,421:-) categorized as an alternative 3′ splice site (A3) only in KB1 (Fig. 4C and Table S12). However, the gene encoding ultraspiracle, which is associated with growth, was found to display DAS categorized as an alternative 5′ splice site (A5) (ΔPSI: − 0.18, occurring in CAKKLH010000192.1:25,109–25177:25,086–25177: +) only in KE2 (Fig. 4D and Table S12). The growth-related DAS gene encoding small monomeric GTPase was only identified in KB1 (Table S12). Furthermore, DAS genes related to growth such as those encoding SSD domain-containing protein (ΔPSI: 0.39 and AF event) and CTP synthase (ΔPSI: 0.31 and RI event) were only identified in KE2, which displayed a significant increase in "Body Length" when exposed to fish kairomone (Table S12).
Comparison of the DETs and DAS genes identified in response to fish kairomone exposure
The number of DETs and DAS genes from both genotypes shown in Fig. 5 may differ from those mentioned above due to the annotation of multiple transcripts with the same protein name. Two transcripts, encoding an EGF-like domain-containing protein and vitellogenin fused with superoxide dismutase, were detected as both DETs and significant DAS events in both genotypes (Table S13). There were 39 transcripts detected as both DETs and significant DAS events in KB1 (Table S13). Among them, transcripts related to reproduction, such as those encoding peroxidase and BZIP domain-containing protein, were identified. On the other hand, five transcripts were detected as both DETs and significant DAS events in KE2 (Table S13). Growth-related transcripts such as that encoding the chitin-binding type-2 domain-containing protein, were identified in KE2. Furthermore, the transcript of the ZP domain-containing protein, which plays a pivotal role in oocyte and gamete development, was upregulated and differentially spliced (∆PSI: 0.51 and RI event) in KB1, while it was downregulated without differential splicing in KE2 (Table S13). The putative ecdysteroid 25-hydroxylas, which plays vital roles in regulating arthropod maturation, growth, and development, was upregulated and differentially spliced (∆PSI: 0.15 and RI event) in KE2, while only being differentially spliced (∆PSI: − 0.45 and A5 event) in KB1 (Table S13).
Enrichment analysis of GO and KEGG pathways of DETs and DAS genes
GO and KEGG enrichment analyses were performed using DETs detected from both genotypes to gain further insights into their underlying molecular mechanisms. We found that most of the significantly enriched GO biological processes (BPs) were related to small molecule, phosphate-containing compound metabolic processes, and macromolecule localization in both genotypes (Table S14). Most enriched GO molecular functions (MFs) were associated with the binding of various molecular compounds (e.g., small molecule, anion, carbohydrate derivative, nucleotide, ATP, etc.) in KB1. Conversely, GO MFs related to growth such as structural molecule activity, structural constituent of the cuticle, Wnt-protein binding, and others were enriched only in KE2, which displayed significantly increased "Body Length" under fish kairomone exposure (Table S14). Furthermore, downregulated DETs in KB1 were enriched in GO cellular components (CCs) such as cuticular protein and endocuticle structural glycoprotein SgAbd-2, which is likely associated with the decrease in "Body Length" under fish kairomone exposure in this genotype. DETs were used to search pathways that contribute to alterations in life-history traits. In KEGG enrichment analysis metabolic pathways, such as those regulating fatty acid biosynthesis and metabolism were identified in both genotypes (Table S15). The pathways of FoxO signaling and mucin-type O-glycan biosynthesis, which play important roles in regulating individual growth and development, were identified under fish kairomone exposure only in genotype KB1, which shows late maturation and decreased body size (Table S15). The Wnt signaling pathway, involved in neuronal development, was identified only in KE2 (Table S15). In the enrichment analysis using DAS genes, BPs related to metabolic processes were identified in both genotypes (Table S16). GO MFs including extracellular matrix structural constituent and copper ion binding were identified in KB1, while various enzymatic (oxidoreductase, carbon–oxygen lyase, etc.) and binding (L-ascorbic acid, ferrous iron, etc.) activities were enriched in KE2 (Table S16). In addition, metabolism-related pathways were only identified in KE2, with no statistically significant pathways found in KB1 (Table S17).
Discussion
Predator-induced responses in Daphnia have been extensively researched [58]. Prey species exhibit various adaptations in response to predation, including life history, behavior, and morphological and physiological trait alterations [59]. The life-history traits exhibited plasticity in response to vertebrate predation risk, with organisms showing earlier or later maturation and adjusting their body size in response to this risk [18, 20]. While many studies have explored the molecular mechanisms underlying the responses to predator stressors based on gene expression, the role of AS remains overlooked. Our previous study also focused on differences in gene expression rather than differential splicing responses to fish kairomones [18]. AS serves as a critical regulatory mechanism enhancing transcriptome and proteome diversity, potentially exerting a significant role in mediating responses to stress alongside gene expression alterations. To the best of our knowledge, this is the first study to find the association of predator-induced phenotypic variation with AS patterns and transcriptional profiles in D. galeata.
Daphnia can increase its body size in response to predator-induced stress, a phenomenon known as inducible defense. Larger body size confers advantages such as increased fecundity and improved competitive ability, which facilitate Daphnia's survival in predator-rich environments. Conversely, Daphnia species can also decrease their body size in response to predators; this adaptive response decreases their likelihood of being detected and eaten by predators. Additionally, previous studies in Daphnia have also shown that energy resources are preferentially allocated for reproduction upon exposure to fish kairomones, increasing the rate at which neonates reach adulthood and facilitating the production of more offspring [58]. Similar to a previous study [6], our results demonstrated altered life-history traits related to reproduction and growth in both genotypes under fish predator stress (Fig. 2). KB1 exhibited defense strategies by enhancing fertility to yield a large number of offspring to ensure the survival of successive generations and reducing “Body Length” to evade detection and predation (Fig. 2). A previous study also identified early maturation and decreased size in Daphnia as a response to vertebrate predators [13]. Early reproduction before predation provides an ecological advantage for maintaining population numbers [13]. In KE2, exposure to predatory fish kairomones led to no significant reproduction-related changes (Fig. 2). However, growth-related life-history traits, such as “Body Length,” increased in KE2, reflecting patterns that contrasted with those observed in KB1 (Fig. 2F). Upon detecting the presence of predatory fish, Daphnia species initiate specific physiological and behavioral mechanisms to enhance their chances of survival. One potential mechanism for phenotypic variation is energy reallocation, where organisms prioritize growth and development over reproduction, resulting in larger body sizes. Exposure to predatory fish kairomones can induce an adaptive response in Daphnia, leading to an increase in body length. Prior research has proposed that when facing threats of predation by fish, Daphnia may exhibit accelerated growth during the juvenile stages and subsequently allocate resources to high reproductive efforts and early spawning upon reaching maturity [14]. Additionally, a larger size at first reproduction could be exhibited in response to fish predation, representing an ecological cost (community trade-off) deviating from the classical prediction of reduced body size in response to fish predation risk [60, 61].
DE analysis was conducted to understand the genetic mechanisms regulating the appropriate functional response in the two D. galeata genotypes that exhibit contrasting phenotypic plasticity under predation stress. Among consistently down-regulated DETs, there were transcripts encoding long-chain-fatty-acid—CoA ligase and C-type lectin domain-containing protein. Previous studies have shown that upon exposure to invertebrate predators in D. pulex, the expression levels of proteins involved in the elongation of very-long-chain fatty acids and C-type lectins decreased [58]. In addition, the downregulation of C-type lectin-related genes was identified in D. galeata upon parasite exposure [62]. This phenomenon seems to occur when Daphnia activate their inducible defenses in response to threats from predators or parasites. The DETs that exhibited inconsistent expression patterns across genotypes include NAD( +) kinase, which was upregulated in KB1 but downregulated in KE2 (Table S5). NAD( +) kinase plays a pivotal role in various metabolic processes, encompassing the production of energy within mitochondria, which is required for cell division, growth, and development in Daphnia [63]. A previous study has shown that a higher abundance of proteins related to energy metabolism and biosynthetic processes supports the increased energy demand of the embryo for expressing Triops-induced defensive traits in D. magna [64]. In contrast, muscle-specific actin was downregulated in KB1, but upregulated in KE2, which displayed significantly increased “Body Length” under exposure to fish kairomones (Fig. 2F and Table S5). In a previous proteomic study, evidence for a higher abundance of muscle-specific actin was found in D. magna exposed to chemical cues of the predator Triops cancriformis [64]. The study also suggested that increasing muscular mass may compensate for the consequences of carapace fortification or altered hydrodynamics resulting from modified carapace morphology [64]. This could be associated with muscle-specific actin expression, although further research is required to fully understand its specific mechanisms under predation stress in Daphnia. The two D. galeata genotypes exhibited higher numbers of genotype-specific DETs than common DETs, which is in line with our previous findings [18]. The growth-related DETs, such as a receptor protein-tyrosine kinase and a delta-like protein, were down-regulated only in KB1 (Table S4), which showed significantly decreased “Body Length” in response to fish kairomones. Conserved receptor protein-tyrosine kinases play a crucial role in various aspects of cell fate determination, differentiation, patterning, proliferation, growth, and survival across metazoans [65]. The delta-like protein is involved in cell–cell notch signaling, which is critical for the development and adult homeostasis of animals [66]. Thus, KB1 seems to adopt an adaptive strategy by reducing its body size, aiming to evade detection or predation by suppressing growth through the downregulation of receptor protein-tyrosine kinases and delta-like proteins associated with growth. However, further research is required to determine the specific role of these proteins in D. galeata under predator-induced stress.
Various AS events were identified in the two genotypes, and several DAS genes appeared to be potentially associated with altered life-history traits in response to fish kairomone exposure. The gene encoding calpain-3 underwent an AF event and was identified as differentially spliced in both genotypes (Figs. 4A and B). In previous research, the calpain gene (CalpA) in Drosophila was found to be specifically expressed in the midgut, which is involved in digestion and nutrient absorption [67]. This implies that CalpA, by participating in protein degradation and regulation, can potentially enhance energy metabolism and ensure sufficient energy and resources for reproduction. Interestingly, calpain-3 displayed an increase in exon inclusion (∆PSI: 0.46) in KB1, which showed increased fertility upon exposure to fish kairomone (Table S7). In contrast, calpain-3 displayed a decrease in exon inclusion (∆PSI: − 0.22) in KE2, which showed no significant differences in life-history traits related to reproduction except for "Number of Offspring First Brood" when exposed to fish kairomone. This indicates that the differential splicing of calpain-3 may be associated with the reproductive-related life-history trait changes observed in the two genotypes in response to fish kairomones. According to a previous study, three types of AS occurred in the mouse calpain-3 gene, including an AF event [68]. These isoforms play a crucial role in muscle tissue by maintaining the structure and function of muscle cells [68]. Various calpain-3 isoforms are expressed in muscle cells, and mutations in this gene can affect their protein degradation and regulation, which are crucial factors in muscle growth and function maintenance [68]. Thus, the differential splicing of calpain-3 may also be associated with the observed life-history trait changes related to growth in the two genotypes. Additionally, growth-related genes such as thymosin beta, which displays DAS classified as an A3 event (∆PSI: − 0.27), were identified in KB1 but not in KE2 (Fig. 4C and Table S7). Thymosin beta, which is the most abundant thymosin type in the body and plays multifaceted roles, including its well-recognized participation in actin regulation by binding to actin monomers [69]. In Caenorhabditis elegans, tetraThymosinbeta plays a crucial role in specific developmental stages that require actin polymerization [70]. Actin is a crucial component in cytoskeleton dynamics, and actin-binding proteins (ABPs) are recognized as critical contributors to skeletal muscle health [71]. The dysfunction of ABP disrupts the essential signaling pathways for muscle development, resulting in muscle disorders [71]. Therefore, the negatively spliced transcript of thymosin beta may regulate the decrease in "Body Length" in genotype KB1, as a response to fish kairomones. On the other hand, ultraspiracle, an essential component of the ecdysone receptor (EcR) and heterodimer receptor complex in insects, displayed DAS classified as an A5 event (ΔPSI: − 0.18) only in KE2 (Fig. 4D and Table S7). Based on previous studies, ultraspiracle protein (USP) is known to play a significant role in insect growth [72]. It is involved in growth and development through the formation of a complex with the EcR, forming heterodimers that bind to DNA and activate specific genes [73]. This regulates the expression of genes related to growth, development, metamorphosis, and reproduction, playing a key role in the regulation of physiological processes in insects. Thus, the differential splicing of USP may also be associated with the observed life-history trait changes related to growth in the two D. galeata genotypes. However, further research is required to elucidate the specific mechanisms of the alternative transcripts of USP under predation stress in D. galeata.
Additionally, we examined the functions of genes displaying DETs and DAS upon exposure to fish kairomones, analyzing their association with phenotypic variation in the two genotypes (Table S8). Overlapping transcripts were found in the DETs and DAS genes, which could account for the association with the phenotypic variation seen in each genotype. In KB1, reproduction-related genes, such as those encoding a peroxidase and BZIP domain-containing protein were identified as common among DET and DAS genes (Table S8). Peroxidase, which plays a fundamental role in the hardening process of the chorion during oogenesis by catalyzing the formation of di- and tri-tyrosine bonds directly associated with the chorion hardening process in D. melanogaster, was identified as an upregulated DET and positive DAS gene [74]. In addition, the BZIP domain-containing protein, identified as an upregulated DET and DAS gene (with positive ∆PSI) in this study, has been found to play essential roles in both Drosophila development and reproduction [75]. In summary, based on the observed phenotypic variations, KB1 exhibited increased fecundity under predator stress, potentially mediated by these proteins. On the other hand, transcripts encoding the chitin-binding type-2 domain-containing protein related to growth were identified as DET and DAS genes in KE2 (Table S8). Chitin filaments and cuticle proteins form the essential constituents of the cuticle [76]. In a previous study, two chitin deacetylase genes linked to neckteeth induction were identified in D. pulex [77].
The aforementioned genes are consistent with the GO and KEGG pathway enrichment analyses conducted in this study, in which several processes related to development and morphogenesis were over-represented. Growth-related GO terms, analyzed using DETs from both genotypes, such as structural molecule activity, structural constituent of the cuticle, Wnt-protein binding, and others, were specifically enriched in KE2 only (Table S9). In D. pulex, the structural constituent of the cuticle was enriched in response to Chaoborus kairomones, indicating its association with cuticle remodeling [58, 78]. Moreover, genes involved in body remodeling and cuticle protein activation were found to be enriched in D. magna when exposed to kairomones from vertebrate and invertebrate predators [79]. Therefore, our results suggest that these GO terms may potentially be related to the increase in "Body Length" in genotype KE2 as a response to fish kairomones. Some reports indicate that signaling pathways regulate developmental changes upon exposure to predators, subsequently influencing reproduction and growth. In the KEGG enrichment analysis using DETs, mucin-type O-glycosylation was enriched under fish kairomone exposure only in KB1 (Table S10). Mucin-type O-glycosylation, a predominant type of O-glycosylation, is an evolutionarily conserved posttranslational modification in animals and it is often found on mucins in the mucous membranes of the digestive tract [80]. Recent studies in model organisms have demonstrated essential roles for O-glycans in embryonic development, organogenesis, and tissue homeostasis [81]. Hence, the activation of this pathway may enhance embryonic development and account for the increased fertility of KB1 in response to fish kairomones. Wnt signaling, a key pathway that affects the late developmental stage of Daphnia [82], was enriched only in KE2 (Table S10). Transcriptional profiling of D. mitsukuri exposed to fish kairomones indicated its response to predation risk through the regulation of Wnt signaling activity [82]. As observed in previous studies, a genotype exhibiting smaller body sizes in response to fish kairomones displayed downregulated DETs involved in this pathway [18]. Therefore, the upregulation of DETs within the Wnt pathway, crucial for animal development and growth, might contribute to the larger body size observed in genotype KE2 [83].
The intraspecific differences observed in response to fish kairomones may arise from differences in expressed transcripts or from splicing machinery functions induced by the kairomones [84]. AS is crucial in increasing proteome diversity, and the frequent occurrence of AS events has been suggested to be responsible for the increased phenotypic complexity of eukaryotes [25]. The abundance of a particular transcript isoform is influenced by three factors: the transcription rate of the gene, the splicing rate of the primary transcript to yield the specific transcript isoform, and a combination of both [85]. Hence, changes in the splicing ratios of exons and introns can modulate transcript abundance independently of overall gene expression. Previous studies have suggested that differences in AS patterns, rather than differences in gene expression, play a more significant role in the diversification of the pharyngeal jaws of cichlids in Lake Tanganyika [86]. Moreover, the exon–intron boundaries comprise highly conserved splice site sequences that must be accurately identified by the spliceosome for adequate splicing [22]. Splicing factors interact with both intronic and exonic sequences to dictate the frequency of exon inclusion in the final transcript [87]. If splicing is not accurately regulated, it may lead to the occurrence of premature stop codons and transcripts with modified amino acid sequences [22]. In addition, the presence of sequence variations near the splice site may affect the competitiveness of the splice acceptor sites, leading to phenotypic variation [88]. Therefore, further investigation into sequence variations at splice sites across genotypes is required to understand the genetic mechanisms of spliced genes associated with phenotypic variation. Future studies should also perform quantitative RT-PCR (qRT-PCR) on the DETs identified in this study to validate our results. Additionally, we intend to validate the AS events selected through bioinformatics analysis, using both long and short-read data in this study, by employing RT-PCR in the future. Importantly, epigenetic modifications, including cytosine methylation, may contribute to our observed results [89]. These factors play a significant role in phenotypic plasticity by regulating gene expression and shaping the organism’s response to environmental stimuli.
Conclusions
This study aimed to integrate gene expression and splicing analyses to comprehensively analyze the genetic mechanisms of phenotypic plasticity under predation stress by examining phenotypic variations induced by fish kairomones in two D. galeata genotypes from the Han River in Korea. The two genotypes showed variation in life-history traits related to reproduction and growth when exposed to fish kairomones, indicating that each genotype displays unique adaptive strategies for coping with predation-induced stress. DE analyses revealed common and unique DETs related to reproduction and growth in response to fish kairomones in both genotypes, implying that they could account for the variation in life-history traits observed in both genotypes. Moreover, DAS analysis exhibited multiple alternatively spliced transcripts associated with life-history trait variation in response to fish kairomones that were not detected in DE analysis. This suggested that splicing may have an equal or even greater impact than DE alone on the post-transcriptional stress response to fish kairomones.
Availability of data and materials
The raw transcriptome data generated and/or analyzed from this study are available on the Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/sra/PRJNA1085014).
Abbreviations
- ABP:
-
Actin-binding proteins
- ADaM:
-
Aachener Daphnien Medium
- AF:
-
Alternative first
- AL:
-
Alternative last
- AS:
-
Alternative splicing
- BP:
-
Biological processes
- BUSCO:
-
Benchmarking universal single-copy ortholog
- CC:
-
Cellular components
- CCS:
-
Circular consensus sequences
- COG:
-
Clusters of Orthologous Groups
- DAS:
-
Differential alternative splicing
- DE:
-
Differential expression
- DET:
-
Differentially expressed transcripts
- EcR:
-
Ecdysone receptor
- EGFR:
-
Epidermal growth factor receptor
- FLNC:
-
Full-length non-concatemer
- GO:
-
Gene Ontology
- HQ:
-
High-quality
- KEGG:
-
Kyoto Encyclopedia of Genes and Genomes
- MF:
-
Molecular functions
- N.S.:
-
No significance
- NR:
-
Non-redundant
- ORF:
-
Open reading frames
- PCA:
-
Principal component analysis
- RI:
-
Retained intron
- SE:
-
Skipped exon
- SMRT:
-
Single-molecule real-time
- TPM:
-
Transcripts per million
References
Stearns SC. The evolutionary significance of phenotypic plasticity. Bioscience. 1989;39:436–45.
Agrawal AA. Phenotypic plasticity in the interactions and evolution of species. Sci. 2001;294:321–6.
Forsman A. Effects of genotypic and phenotypic variation on establishment are important for conservation, invasion, and infection biology. Proc Natl Acad Sci. 2014;111(1):302–7.
Fine PV. Ecological and evolutionary drivers of geographic variation in species diversity. Annu Rev Ecol Evol. 2015;46:369–92.
Aldana M, Maturana D, Pulgar J, García-Huidobro MR. Predation and anthropogenic impact on community structure of boulder beaches. Sci Mar. 2016;80:543–51.
Stibor H, Lüning J. Predator-Induced Phenotypic Variation in the Pattern of Growth and Reproduction in Daphnia hyalina (Crustacea: Cladocera). Funct Ecol. 1994;8:97–101.
Schoeppner NM, Relyea RA. Interpreting the smells of predation: how alarm cues and kairomones induce different prey defences. Funct Ecol. 2009;23:1114–21.
Yin M, Laforsch C, Lohr JN, Wolinska J. Predator-induced defense makes Daphnia more vulnerable to parasites. Evolution. 2011;65(5):1482–8.
Herzog Q, Rabus M, Ribeiro BW, Laforsch C. Inducible Defenses with a" Twist": Daphnia barbata Abandons Bilateral Symmetry in Response to an Ancient Predator. PLoS ONE. 2016;11:e0148556.
Lampert W. Daphnia: development of a model organism in ecology and evolution. Int Ecol Inst. 2011;21:1–275.
Miner BE, de Meester L, Pfrender ME, Lampert W, Hairston NG. Linking genes to communities and ecosystems: Daphnia as an ecogenomic model. Proc Royal Soc. 2012;279:1873–82.
Shapiro J. The importance of trophic-level interactions to the abundance and species composition of algae in lakes. Hyper Eco. 1980;2:105–16.
Lampert W. Phenotypic plasticity of the size at first reproduction in Daphnia: the importance of maternal size. Ecology. 1993;74(5):1455–66.
Carter MJ, Silva-Flores P, Oyanedel JP, Ramos-Jiliberto R. Morphological and life-history shifts of the exotic cladoceran Daphnia exilis in response to predation risk and food availability. Limnologica. 2013;43(3):203–9.
Tams V, Nickel JH, Ehring A, Cordellier M. Insights into the genetic basis of predator-induced response in Daphnia galeata. Ecol Evol. 2020;10(23):13095–108.
Choi TJ, Do TD, Kim JI, Kim CB. Analysis of the complete mitogenome of Daphnia galeata from the Han River, South Korea: structure comparison and control region evolution. Funct Integr Genom. 2023;23(1):65.
Stich HB, Lampert W. Predator Evasion as an Explanation of Diurnal Vertical Migration by Zooplankton. Nature. 1981;293:396–8.
Choi TJ, Han SM, Malik A, Kim CB. Comparative transcriptome analysis of two Daphnia galeata genotypes displaying contrasting phenotypic variation induced by fish kairomones in the same environment of the Han River, Korea. BMC Genom. 2023;24(1):580.
Simon JC, Pfrender ME, Tollrian R, Tagu D, Colbourne JK. Genomics of environmentally induced phenotypes in 2 extremely plastic arthropods. J Hered. 2011;102(5):512–25.
Tams V, Lüneburg J, Seddar L, Detampel JP, Cordellier M. Intraspecific phenotypic variation in life history traits of Daphnia galeata populations in response to fish kairomones. PeerJ. 2018;6:e5746.
De Coninck DI, Janssen CR, De Schamphelaere KA. An investigation of the inter-clonal variation of the interactive effects of cadmium and Microcystis aeruginosa on the reproductive performance of Daphnia magna. Aquat Toxicol. 2013;140:425–31.
Xu S, Snyman M. The effects of mutations on gene expression and alternative splicing. Proc R Soc B Biol Sci. 2023;290:111–23.
Modrek B, Lee C. A genomic view of alternative splicing. Nat Genet. 2002;30(1):13–9.
Suresh S, Crease TJ, Cristescu ME, Chain FJ. Alternative splicing is highly variable among Daphnia pulex lineages in response to acute copper exposure. BMC Genom. 2020;21(1):1–14.
Brites D, McTaggart S, Morris K, Anderson J, Thomas K, Colson I, Fabbro T, Little TL, Ebert D, Pasquier LD. The Dscam homologue of the crustacean Daphnia is diversified by alternative splicing like in insects. Mol Biol Evol. 2008;25(7):1429–39.
Wang Y, Liu J, Huang BO, Xu YM, Li J, Huang LF, Lin J, Zhang J, Min QH, Yang WM, Wang XZ. Mechanism of alternative splicing and its regulation. Biomed Rep. 2015;3(2):152–8.
Pop M, Salzberg SL. Bioinformatics challenges of new sequencing technology. Trends Genet. 2008;24(3):142–9.
Ye S, Yu X, Chen H, Zhang Y, Wu Q, Tan H, Song J, Saqib HSA, Farhadi A, khwanuddin M, Ma H. Full-Length transcriptome reconstruction reveals the genetic mechanisms of eyestalk displacement and its potential implications on the interspecific hybrid crab (Scylla serrata♀× S. paramamosain♂). Biology. 2022;11(7):1026.
Liu M, Xiao F, Zhu J, Fu D, Wang Z, Xiao R. Combined PacBio Iso-Seq and Illumina RNA-Seq Analysis of the Tuta absoluta (Meyrick) Transcriptome and Cytochrome P450 Genes. Insects. 2023;14(4):363.
Marden JH. Quantitative and evolutionary biology of alternative splicing: how changing the mix of alternative transcripts affects phenotypic plasticity and reaction norms. Heredity. 2008;100(2):111–20.
Thielsch A, Brede N, Petrusek A, De Meester LUC, Schwenk K. Contribution of cyclic parthenogenesis and colonization history to population structure in Daphnia. Mol Ecol. 2009;18(8):1616–28.
R Core Team. R: A Language and Environment for Statistical Computing. 2017. https://www.R-project.org/.
MacFarland TW, Yates JM. Oneway Analysis of Variance (ANOVA). In: Using R for Biostatistics. Cham: Springer; 2021. https://doi.org/10.1007/978-3-030-62404-0_5.
Weiss LC, Heilgenberg E, Deussen L, Becker SM, Kruppert S, Tollrian R. Onset of kairomone sensitivity and the development of inducible morphological defenses in Daphnia pulex. Hydrobiologia. 2016;779:135–45.
Salmela L, Rivals E. LoRDEC: accurate and efficient long read error correction. Bioinform. 2014;30(24):3506–14.
Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinform. 2018;34(18):3094–100.
Nickel J, Schell T, Holtzem T, Thielsch A, Dennis SR, Schlick-Steiner BC, Steiner FM, Möst M, Pfenninger M, Schwenk K, Cordellier M. Hybridization dynamics and extensive introgression in the Daphnia longispina species complex: new insights from a high-quality Daphnia galeata reference genome. Genome Biol Evol. 2021;13(12):evab267.
Tseng E. cDNA cupcake. Available online at: https://github.com/Magdoll/cDNA_Cupcake.
Tseng E. Cogent. Available online at: https://github.com/Magdoll/Cogent.
Huang Y, Niu B, Gao Y, Fu L, Li W. CD-HIT Suite: a web server for clustering and comparing biological sequences. Bioinform. 2010;26(5):680–2.
Manni M, Berkeley MR, Seppey M, Simão FA, Zdobnov EM. BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, prokaryotic, and viral genomes. Mol Biol Evol. 2021;38(10):4647–54.
Pertea G, Pertea M. GFF utilities: GffRead and GffCompare. F1000Research. 2020;9:304.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, MacManes MD, Ott M, Orvis J, Pochet N, Strozzi F, Weeks N, Westerman R, William T, Dewey CN, Henschel R, LeDuc RD, Friedman N, Regev A. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Pruitt KD, Tatusova T, Maglott DR. NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res. 2007;35(suppl_1):D61–5.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: Architecture and applications. BMC Bioinform. 2009;10:421.
Buchfink B, Reuter K, Drost HG. Sensitive protein alignments at tree-of-life scale using DIAMOND. Nat Methods. 2021;18(4):366–8.
UniProt Consortium. UniProt: The universal protein knowledgebase in 2023. Nucleic Acids Res. 2023;51:D523–31.
Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, Von Mering C, Bork P. Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol Biol Evol. 2017;34(8):2115–22.
Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. Available online at: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Li B, Dewey CN. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12:e323.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57:289–300.
Ge SX, Jung D, Yao R. ShinyGO: A graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36:2628–9.
Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer. 2009. Available online at: http://had.co.nz/ggplot2/book.
Trincado JL, Entizne JC, Hysenaj G, Singh B, Skalic M, Elliott DJ, Eyras E. SUPPA2: fast, accurate, and uncertainty-aware differential splicing analysis across multiple conditions. Genome Biol. 2018;19:1–11.
Garrido-Martín D, Palumbo E, Guigó R, Breschi A. ggsashimi: Sashimi plot revised for browser-and annotation-independent splicing visualization. PLoS Comp Biol. 2018;14(8):e1006360.
Rozenberg A, Parida M, Leese F, Weiss LC, Tollrian R, Manak JR. Transcriptional profiling of predator-induced phenotypic plasticity in Daphnia pulex. Front Zool. 2015;12(1):18.
Lima SL. Nonlethal effects in the ecology of predator-prey interactions. What are the ecological effects of anti-predator decision-making? BioScience. 1998;48:25–34.
Garay-Narvaez L, Ramos-Jiliberto R. Induced defenses within food webs: the role of community trade-offs, delayed responses, and defense specificity. Ecol Complex. 2009;6(3):383–91.
Tollrian R, Harvell CD. The ecology and evolution of inducible defenses. New Jersey: Princeton university press; 1999.
Lu Y, Johnston PR, Dennis SR, Monaghan MT, John U, Spaak P, Wolinska J. Daphnia galeata responds to the exposure to an ichthyosporean gut parasite by down-regulation of immunity and lipid metabolism. BMC Genom. 2018;19:1–10.
Oka SI, Titus AS, Zablocki D, Sadoshima J. Molecular properties and regulation of NAD+ kinase (NADK). Redox Biol. 2023;59:102561.
Otte KA, Fröhlich T, Arnold GJ, Laforsch C. Proteomic analysis of Daphnia magna hints at molecular pathways involved in defensive plastic responses. BMC Genom. 2014;15:1–17.
Sopko R, Perrimon N. Receptor tyrosine kinases in Drosophila development. Cold Spring Harb Perspect Biol. 2013;5(6):a009050.
Haltom AR, Jafar-Nejad H. The multiple roles of epidermal growth factor repeat O-glycans in animal development. Glycobiology. 2015;25(10):1027–42.
Theopold U, Pintér M, Daffre S, Tryselius Y, Friedrich P, Nässel DR, Hultmark D. CalpA, a Drosophila calpain homolog specifically expressed in a small set of nerve, midgut, and blood cells. Mol Cell Biol. 1995;15(2):824–34.
Herasse M, Ono Y, Fougerousse F, Kimura EI, Stockholm D, Beley C, Montarras D, Pinset C, Sorimachi H, Suzuki K, Beckmann JS, Richard I. Expression and functional characteristics of calpain 3 isoforms generated through tissue-specific transcriptional and posttranscriptional events. Mol Cell Biol. 1999;19(6):4047–55.
Yu FX, Yin HL, Morrison-Bogorad M, Lin SC. Effects of thymosin β4 and thymosin β10 on actin structures in living cells. Cell motil cytoskelet. 1994;27(1):13–25.
Herrmann D, Hatta M, Hoffmeister-Ullerich SAH. Thypedin, the multi copy precursor for the hydra peptide pedin, is a β-thymosin repeat-like domain containing protein. Mech Dev. 2005;122(11):1183–93.
Nguyen MT, Dash R, Jeong K, Lee W. Role of Actin-Binding Proteins in Skeletal Myogenesis. Cells. 2023;12(21):2523.
Hult EF, Tobe SS, Chang BS. Molecular evolution of ultraspiracle protein (USP/RXR) in insects. PLoS ONE. 2011;6(8):e23416.
Ogura T, Tan A, Tsubota T, Nakakura T, Shiotsuki T. Identification and expression analysis of ras gene in silkworm, Bombyx mori. PLoS ONE. 2009;4(11):e8030.
Konstandi OA, Papassideri IS, Stravopodis DJ, Kenoutis CA, Hasan Z, Katsorchis T, Wever R, Margaritis LH. The enzymatic component of Drosophila melanogaster chorion is the Pxd peroxidase. Insect Biochem Mol Biol. 2005;35(9):1043–57.
Guo Z, Qin J, Zhou X, Zhang Y. Insect transcription factors: a landscape of their structures and biological functions in Drosophila and beyond. Int J Mol Sci. 2018;19(11):3691.
Yang F, Li X, Li S, Xiang J, Li F. A novel cuticle protein involved in WSSV infection to the Pacific white shrimp Litopenaeus vannamei. Dev Comp Immunol. 2020;102:103491.
Christjani M, Fink P, von Elert E. Phenotypic plasticity in three Daphnia genotypes in response to predator kairomone: evidence for an involvement of chitin deacetylases. J Exp Biol. 2016;219(11):1697–704.
An H, Do TD, Jung G, Karagozlu MZ, Kim CB. Comparative Transcriptome Analysis for Understanding Predator-Induced Polyphenism in the Water Flea Daphnia pulex. Int J Mol Sci. 2018;19:2110.
Orsini L, Brown JB, Shams Solari O, Li D, He S, Podicheti R, Stoiber MH, Spanier KI, Gilbert D, Jansen M, Rusch DB, Pfrender ME, Colbourne JK, Frilander MJ, Kvist J, Decaestecker E, De Schamphelaere KAC, De Meester L. Early transcriptional response pathways in Daphnia magna are coordinated in networks of crustacean-specific genes. Mol Ecol. 2018;27(4):886–97.
Itoh K, Nishihara S. Mucin-Type O-Glycosylation in the Drosophila Nervous System. Front Neuroanat. 2021;15:767126.
Tran DT, Ten Hagen KG. Mucin-type O-glycosylation during development. J Biol Chem. 2013;288(10):6921–9.
Zhang X, Blair D, Wolinska J, Ma X, Yang W, Hu W, Yin M. Genomic regions associated with adaptation to predation in Daphnia often include members of expanded gene families. Proc R Soc B: Biol Sci. 1955;2021(288):20210803.
Holstein TW. The evolution of the Wnt pathway. Cold Spring Harb Perspect Biol. 2012;4(7):a007922.
Richards AL, Watza D, Findley A, Alazizi A, Wen X, Pai AA, Pique-Regi R, Luca F. Environmental perturbations lead to extensive directional shifts in RNA processing. PLoS Genet. 2017;13(10):e1006995.
Gonzàlez-Porta M, Calvo M, Sammeth M, Guigó R. Estimation of alternative splicing variability in human populations. Genome Res. 2012;22(3):528–38.
Singh P, Börger C, More H, Sturmbauer C. The role of alternative splicing and differential gene expression in cichlid adaptive radiation. Genome Biol Evol. 2017;9(10):2764–81.
Wang Z, Burge CB. Splicing regulation: from a parts list of regulatory elements to an integrated splicing code. RNA. 2008;14(5):802–13.
Lin CH, Patton JG. Regulation of alternative 3’splice site selection by constitutive splicing factors. RNA. 1995;1(3):234–45.
Asselman J, De Coninck DI, Vandegehuchte MB, Jansen M, Decaestecker E, De Meester L, Vanden Bussche J, Vanhaecke L, Jassen CR, De Schamphelaere KA. Global cytosine methylation in Daphnia magna depends on genotype, environment, and their interaction. Environ Toxicol Chem. 2015;34(5):1056–61.
Acknowledgements
We acknowledged Min-Ho Mun, Jung-Il Kim, and Hyung-Eun An (Sangmyung University, Korea) for help during Daphnia breeding and the experiment.
Funding
This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) [No. NRF-2022R1A2C3002750].
Author information
Authors and Affiliations
Contributions
CBK: Funding acquisition, Resources, Conceptualization, Methodology, Writing - review & editing. TJC: Conceptualization, Resources, Data curation, Methodology, Writing original draft. SMH: Conceptualization, Resources, Data curation, Methodology. AM: Data curation, Methodology, Writing - review & editing.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Choi, TJ., Malik, A., Han, SM. et al. Differences in alternative splicing events in the adaptive strategies of two Daphnia galeata genotypes induced by fish kairomones. BMC Genomics 25, 725 (2024). https://doi.org/10.1186/s12864-024-10643-1
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12864-024-10643-1