Insights into adaptive divergence of Japanese mantis shrimpOratosquilla oratoriainferred from comparative analysis of full-length transcriptomes

1 Introduction

In marine ecosystems, oceanographic heterogeneity (e.g., current interfaces, habitat transition zones, ecological gradients) can impose spatially varying selective pressures on populations or species distributed along the gradients. Such heterogeneous seascapes in coastal waters have been demonstrated to strongly influence adaptive divergence, particularly for species with large population sizes where selection is expected to be highly efficient (年代andoval-Castillo et al., 2018;Coscia et al., 2020;Pratt et al., 2022）.The Japanese mantis shrimpOratosquilla oratoria(De Haan, 1844) is benthic, neritic and burrowing shrimp with a wide distribution in the Northwestern Pacific (NWP) (Manning, 1971）.Due to its high productivity as well as excellent meat quality,O. oratoriais commercially exploited in coastal waters throughout Japan, Korea and China (Kodama et al., 2006;Ahyong, 2012）.This species occurs within a highly heterogeneous seascape that spans coastal bioregions with strong gradients in sea surface temperature (SST) and salinity (Figure 1）.More specifically, SST is strongly influenced by the complex current systems with different temperature characteristics, including the cold Oyashio Current (OC), China Costal Current (CCC), a low-temperature Yellow Sea cold water mass and the warm Kuroshio Current as well as its branches (Figure 1A,Liu, 2013）.A biogeographic boundary lies in line with the Yangtze River Estuary where the annual SST is above 20°C on the southern side of the estuary, and rapidly decreases to the value less than 15°C on the northern side (Figure 1B,Johnson and Boyer, 2015）.Taxa with narrow thermal tolerance ranges are restricted to one side of the boundary as a result of this steep temperature gradient across the Yangtze River Estuary (Liu, 2013）.In addition, a strong salinity gradient is generated by the collision of the Changjiang diluted water (CDW) with several coastal currents (e.g., CCC and Subei Coastal Current) (Figure 1C,年代u and Yuan, 2005）.Therefore, the wide distribution range across heterogeneous seascapes makesO. oratoriaan ideal system to study how ecologically based natural selection has driven adaptive divergence of marine coastal species.

年代everal recent researches have revealed a remarkable pattern of genetic differentiation in Japanese mantis shrimp despite high dispersal potential during its planktonic larval stage (4–6 weeks,Hamano and Matsuura, 1987）.Our previous analyses have revealed two highly divergent lineages inO. oratoriawith clear allopatric distribution by integrating mitochondrial and nuclear evidence (Cheng and Sha, 2017;Cheng et al., 2020）.The distribution range of the twoO. oratorialineages exhibits a remarkable latitudinal cline: the northern lineage is restricted to cold waters with temperate affinities, while the southern lineage inhabits the subtropical and tropical regions influenced by the Kuroshio Current and its branches. The deep genetic divergence, which is so far unlinked to any apparent morphological variation, substantiates the existence of cryptic speciation in Japanese mantis shrimp (Cheng and Sha, 2017）.In the coast of China, the sharp genetic break ofO. oratoriais in line with the Yangtze River Estuary (Du et al., 2016;Cheng and Sha, 2017;Cheng et al., 2020）.Nevertheless, these studies are constrained by small numbers of markers that precludes an accurate description of genome-wide DNA sequence divergence. Meanwhile, we know little about how seascape heterogeneity in the NWP influences adaptative evolution ofO. oratoriaas well as the roles of genetic variation and plasticity in local adaptation to environmental gradients. The polygenic basis of traits governing physiological tolerance puts forward challenges for linking genetic variation to ecological differences among populations adapting to different environments (年代torz and Wheat, 2010;Rockman, 2012）.Thus, the genetic underpinnings of diversification and adaptation in Japanese mantis shrimp remain unclear due to the lack of genomic information.

Transcriptome represents a sample of the spatiotemporally expressed genome that can be generated in a short time at a reduced cost compared to whole genome sequencing, and can be used as an alternative to genomic approaches for non-model organisms (Morozova and Marra, 2008）.Beyond uncovering the molecular bases of physiological responses to a particular environmental challenge (Gracey and Cossins, 2003), transcriptomics has also been used in a comparative framework to reveal numerous aspects of ecological speciation and adaptation (e.g.,Zhang et al., 2015a;Cheng et al., 2019）.As a low-cost next-generation sequencing technology, RNA sequencing (RNA-seq) has become a versatile tool for studying transcriptomics. However, RNA-seq is incapable of providing full-length (FL) transcripts due to its inherent length limitations, which presents major challenges for further structural and functional genomic studies (年代teijger et al., 2013;Tilgner et al., 2013）.Another limitation is that the presence of different isoforms and differences in transcript abundance has greatly hampered transcriptome assembly from short reads. The third-generation sequencing technology seems to offer an opportunity to overcome the shortcomings of short-read sequences by enabling the generation of kilobase-sized sequencing reads without assembly (年代haron et al., 2013）.To date, an increasing number of FL transcriptome analyses have been conducted in a variety of marine organisms (e.g.,Yi et al., 2018;Pootakham et al., 2020;Wang et al., 2022), which greatly facilitate the transcriptome research in the ocean life.

In the present study, we deep-sequenced the FL transcriptomes of northern and southernO. oratoriaby using PacBio single molecule real-time (SMRT) sequencing technology. Based on the obtained FL transcriptomes, we first performed transcription factor prediction, long noncoding RNA prediction and transcript functional annotation. A comparative transcriptomic analysis was then carried out to identify orthologous genes and to evaluate the transcriptome-wide genetic divergence between northern and southernO. oratoria。With these datasets, we aimed to trace signatures of positive selection inO. oratoriaduring adaptation to heterogeneous environments and identify environmentally responsive genes that ecologically based natural selection may have acted on. Finally, we incorporated expression profiles ofO. oratoriaunder heat stress to investigate the role of evolving plasticity in local adaptation to environmental gradients. As such, this study provides abundant transcriptome resources and new insights into the divergence and adaptation of Japanese mantis shrimp, and also provides ample evidence for adaptive evolution of marine coastal species in response to heterogeneous seascapes.

2 Materials and methods

2.1 Sample collection and RNA preparation

The adultO. oratoriaindividuals were collected from the coastal waters of Qingdao (QD) and Xiamen (XM), China, which correspond to northern and southernO. oratoria, respectively (Figure 1A）.All individuals were acclimated over 24 h at 20°C. One individual at each sampling site was randomly selected and the gill, hepatopancreas, digestive tract, nervous chain, and abdominal muscle were immediately dissected out, frozen and stored in liquid nitrogen. Subsequently, all ten samples were subjected to RNA extraction using the TRIzol kit (Invitrogen, USA) according to the manufacturer’s instructions, and RNA degradation and contamination were detected with 1% agarose gels. The integrity and purity of RNA were assessed using NanoDrop 2000 (Thermo Scientific, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Only qualified RNA samples were used for cDNA library constructions.

2.2 Library construction and sequencing

RNA样品从每个在五个不同的组织dividual were mixed in equal amount to construct the PacBio sequencing library. Specifically, mRNA was purified from mixed total RNA using poly (T) oligo-attached magnetic beads and reverse-transcribed into FL cDNA using the SMARTer PCR cDNA Synthesis Kit (Clontech, USA). After PCR amplification, quality control and purification, the BluePippin Size Selection system (Sage Science, USA) was used for size selection of the FL cDNA and for producing libraries of differently sized cDNA. The screened cDNA was re-amplified by PCR, end repaired, connected to the SMRT dumbbell-type connector, and exonuclease digested. The cDNA products were then subjected to construction of SMRTbell Template libraries using SMRTBell Template Prep Kit. The concentration and quality of the libraries was assessed using Agilent 2100 Bioanalyzer and Qubit 2.0 (Life Technologies, USA). Finally, qualified libraries were sequenced on PacBio Sequel platform (Pacific Biosciences, USA).

The remaining RNA of each tissue was used for Illumina sequencing library construction. Five separate Illumina libraries were constructed for northern and southernO. oratoriausing the protocol of NEBNext UltraTM RNA Library Prep Kit (NEB, USA), respectively. Briefly, poly (A) mRNA was purified from total RNA using Oligo (dT) magnetic beads and then broken into short fragments to synthesize first strand cDNA using random hexamer primer. Second-strand cDNA was then synthesized using DNA polymerase I and RNaseH. The cDNA was subjected to end-repair, phosphorylation, 3’ adenylation, and ligation to sequencing adaptors. Afterwards, cDNA libraries were generated by PCR amplification, and the library preparations were paired-end sequenced at 150 bp on an Illumina HiSeq X Ten platform. Clean Illumina reads were produced after removing adaptor sequences, reads containing ploy-N (with the ratio of ‘N’ to be more than 10%) and low quality reads (with quality score less than 5), and then all clean reads from the same site were merged together forde novoassembled by using Trinity v2.5.1 (Grabherr et al., 2011) with default parameters.

2.3 Pacbio read error correction

The SMRT Link 5.1 pipeline was used for PacBio data processing according to the official protocol. Circular consensus sequences (CCSs) were extracted from raw reads with the following parameters: max_drop_fraction, 0.8; max_length, 18000; min_length, 200; min_predicted_accuracy, 0.8; min_passes, 1; min_zscore, -999. After discarding CCS reads with length shorter than 50 bp, the remaining CCS reads were classified into full length non-chimeric (FLNC) or non-full-length (NFL) transcripts according to whether 5’/3’ cDNA primers and poly(A) tail were simultaneously observed. The iterative clustering for error correction (ICE) algorithm was used to obtain consensus isoforms by clustering FLNC sequences, and Arrow software (https://github.com/PacificBiosciences/GenomicConsensus) was used to refine consensus isoforms using the NFL reads to produce polished consensus sequences. All polished consensus sequences were corrected in aid of Illumina RNA-seq data using LoRDEC (年代almela and Eric, 2014）.Furthermore, redundant sequences were removed from the transcriptome isoform sequences to obtain unigenes by using CD-HIT (Fu et al., 2012) with an amino-acid sequence identity threshold of 95%. Finally, the completeness of two FL transcriptomes was assessed by examining the coverage of the 1066 conserved core genes of Arthropoda (https://busco.ezlab.org/) using BUSCO 3 (年代imão et al., 2015）.

2.4 Coding sequences, transcription factor and long noncoding RNA prediction

All the isoforms were used to predict the coding sequences (CDS) and protein sequences using ANGEL software (年代himizu et al., 2006）.We used the confidence protein sequences ofO. oratoriaor closely related species for ANGLE training, and then run the ANGLE predictions for given sequences. The 5’UTR and 3’UTR (untranslated regions) sequences were predicted based on the CDS, start and stop codons. Transcripts containing the 5’ and 3’UTRs and complete CDSs were defined as FL transcripts. Transcription factors (TFs) were predicted using AnimalTFDB v2.0 (Animal Transcription Factor Database) (Zhang et al., 2015b）.BecauseO. oratoriawas not included in the database, hmmsearch was used to identify TFs based on the Pfam search results of the TF family.

Long noncoding RNAs (lncRNAs) are defined as RNAs which are at least 200 nucleotides in length and lack protein-coding capacity (Rinn and Chang, 2012）.On the basis of this, lncRNA ofO. oratoriatranscriptome was first predicted by screening the coding potential of transcripts using PLEK (Li et al., 2014), Coding-Non-Coding-Index (CNCI) (年代un et al., 2013) and Coding Potential Calculator (CPC) (Kong et al., 2007）.The transcript sequences predicted using PLEK, CNCI and CPC tools were further used to search against the Pfam-A and Pfam-B databases using hmmscan. Transcripts with predicted coding potential according to the results of all four methods were filtered out, and those without coding potential constituted the candidate set of lncRNAs.

2.5 Gene function annotation

The non-redundant transcript sequences were mapped to public databases to obtain the annotation information. Transcripts were compared against NCBI NT (non-redundant nucleotide sequences) using BLAST v2.2.31 (Altschul et al., 1997) with cut-off E-value < 1E-5, NCBI NR (non-redundant protein sequences), Swiss-Prot (http://www.ebi.ac.uk/uniprot/), KOG (euKaryotic Ortholog Group), KEGG (Kyoto Encyclopedia of Genes and Genomes) using diamond v0.8.36 (Li et al., 2002) with cut-off E-value < 1E-5. According to the annotation results of the NR database, the Blast2GO v2.5 software (Conesa et al., 2005) was used to perform Gene Ontology (GO) annotation and classification. KEGG classification was performed using KASS (vr140224) (Moriya et al., 2007) and the KEGG Automatic Annotation Server. HMMER v3.1 (Nastou et al., 2016)被用来比较tran的氨基酸序列scripts against the Protein family (Pfam) database for Pfam annotation.

2.6 Putative orthologs identification and sequence divergence analyses

Orthologous groups were constructed from the BLASTP results using OrthoMCL (Li et al., 2003) based on the Markov Cluster algorithm (mcl) with default settings. The 5’UTR, coding and 3’UTR regions were separately extracted from each pair of orthologs. The CDS and UTR regions were aligned separately to each other using MUSCLE (Edgar, 2004) with default settings. For the CDS region, pair-wise alignments were performed for the putative orthologous pairs based on protein sequences, and back-translated to DNA sequences for subsequent analysis. The phylogenetic tree was conducted on the basis of amino acid sequence alignment using the maximum likelihood method as implemented in PhyML v3.0 (Guindon et al., 2010）.The node reliability was calculated from 1000 bootstrap replications. According to the method ofYe et al. (2014), sequence divergences between the homologous regions of each gene pair were calculated by dividing the number of substitutions by the number of base pairs compared.

2.7 Test for positive selection

The ratio of non-synonymous to synonymous nucleotide substitutions provides information about the evolutionary forces operating on a gene (Biswas and Akey, 2006）.We estimated non-synonymous substitutions per nonsynonymous site (Ka) and synonymous substitutions per synonymous site (Ks) among putatively orthologous coding regions using a maximum-likelihood method implemented in the codeml program in the PAML package (Yang, 2007) with one-ratio model (model = 0). In this study, Ka/Ks > 1 was considered to indicate that genes were under positive selection and Ka/Ks < 0.1 was regarded as a signature of purifying selections. Positively selected genes (PSGs) with Ka/Ka > 1 were characterized with GO enrichment analysis implemented in the GOseq R packages (Young et al., 2010) based on the Wallenius noncentral hypergeometric distribution. Go terms with aPvalue < 0.05 were considered significantly enriched. The R bioconductor package topGO (Alexa and Rahnenführer, 2010) was used to present the directed acycline graph of the enriched GO category.

2.8 Expression of PSGs response to thermal stress

To infer the potential role of positive selection inO. oratoriaadaptation, we explored the expression level of PSGs in two representative geographic populations (Qingdao in the Yellow Sea and Zhoushan in the East China Sea) in response to the same heat stress (20°C–28°C,Lou et al., 2019a）.Taking advantage of range-wide sampling, our previous analyses of mitochondrial and nuclear sequences indicated genetic homogeneity amongO. oratoriapopulations from the East and South China Seas (Cheng and Sha, 2017;Cheng et al., 2020）.These two heat-stressed populations reported inLou et al. (2019a)correspond to northern and southernO. oratoria在这项研究中,分别。加热方案啊f the two populations were similar, with detailed description inLou et al. (2019a)。Briefly, the control condition was kept at a water temperature of 20°C after acclimation, whereas water temperature for the heat stress condition was raised at 1°C per day to 28°C and then held constant for 24h. All raw reads in the FASTQ format (BioProject accession: PRJNA475657) were filtered by removing the reads with sequencing adaptors, unknown nucleotides (N ratio > 10%) and low quality (quality scores ≤ 5). Clean reads of each sample were then mapped to the orthologous assembly of each respective group by using RSEM (Li and Dewey, 2011）.The expression level was determined by averaging the expression values of three biological replications for each experimental condition. The differentially expressed orthologous genes were identified using the DESeq R package v1.18.0 (Anders and Huber, 2010) with the filtering thresholds of |log2FC| ≥ 1 andP≤ 0.05. Hierarchical clustering was performed using the heatmap.2 function in the gplots R package (RDevelopment Core Team, 2022）.

3 Results

3.1 Full-length transcriptome sequencing

Bases on the PacBio SMRT sequencing technology, a total of 116.92 Gb and 110.32 Gb subreads were obtained, with N50 of 3,465 bp and 2,350 bp for the northern and southernO. oratoria, respectively. As shown inTable 1andFigure S1, these subreads yielded 1,192,616 CCSs with a mean length of 3,714 bp for the northernO. oratoriaand 1,312,861 CCSs with a mean length of 2,390 bp for the southernO. oratoria。According to the presence or absence of 5’/3’ cDNA primers and poly(A) tail, 1,022,216 (85.71%) FLNC reads and 162,789 NFL reads were further identified from the CCS reads of the northernO. oratoria。For the southernO. oratoria988172 (75.27%) FLNC读取和311339年NFL阅读s were generated from the CCS reads. After isoform-level clustering based on the ICE algorithm and polishing based on the Arrow algorithm, a total of 46,871 and 55,408 polished FL consensus isoforms with an average length of 3,439 bp and 2,275 bp were generated from the FLNC reads of northern and southernO. oratoria, respectively. After removing redundancy, the consensus isoforms were finally clustered into a total of 18,630 unigenes with N50 value of 3,852 bp for the northernO. oratoria, and 21,843 unigenes with N50 value of 3,061 bp for the southernO. oratoria(Table 1）.Of the 1066 conserved core genes of Arthropoda, 682 genes (64.0%) and 840 genes (78.8%) were mapped to the FL transcriptome of northern and southernO. oratoria, respectively (Table S1), indicating the sequencing data in this study were relatively complete and suitable for downstream analyses.

3.2 Comparison between PacBio and Illumina unigenes

The transcriptomes of five tissues of northern and southernO. oratoriawere separately sequenced using the Illumina platform (Table S2）.After trimming and filtering, clean reads were assembled into 60,802 unigenes with an average length of 1,212 bp and N50 value of 2,140 bp for the northernO. oratoria, and 117,403 unigenes with an average length of 997 bp and N50 value of 1,563 bp for the southernO. oratoria, which were obviously shorter than unigenes in both FL transcriptomes. Most of PacBio unigenes of northern and southernO. oratoriahad lengths > 2,000 bp, accounting for 87.1% and 61.0% of the total number, while most Illumina unigenes of northern and southernO. oratoriahad lengths < 2,000 bp, accounting for 83.5% and 88.7%, respectively (Figure 2）.The similarity between PacBio and Illumina unigenes were further conducted using BLAST v2.2.31 with the parameter set to 1E-5. The results showed that only 9.78% (5944) and 7.91% (9288) of Illumina unigenes had high similarity to 83.0% (15394) and 86.1% (18790) of PacBio unigenes in northern and southernO. oratoria, respectively.

3.3 Full-length transcriptome structure

3.3.1 Coding sequences

Through the ANGEL predictions, a total of 18,716 CDSs were generated for the northernO. oratoria, of which 10,351 containing the start and stop codons were defined as complete open reading frames (ORFs). For southernO. oratoria, a total of 21,275 CDSs were predicted, including 10,286 complete ORFs. The length distributions of CDSs are shown inFigure S2。The CDS length of northernO. oratoriaranged from 49 bp to 3,590 bp with an average length of 477 bp, and the CDS length of southernO. oratoriaranged from 49 bp to 3,943 bp with an average length of 614 bp.

3.3.2 Transcription factor

In total, 1,549 putative TFs from 29 TF gene families were identified in the northernO. oratoria, and 944 putative TFs from 29 TF gene families were identified in the southernO. oratoria(Figure S3A）.In both FL transcriptomes, Zinc finger C2H2 (zf-C2H2) was the most abundant TF gene family, followed by the family of BTB domain and Zinc Finger-containing (ZBTB) transcription factors. A GO enrichment analysis was further conducted to determine the potential functions of genes in which TFs have been determined. For both FL transcriptomes, binding (GO:0005488), metal ion binding (GO:0046872) and cation binding (GO:0043169) were the most enriched GO terms (Figure S3B）.

3.3.3 Long noncoding RNA

Four prediction approaches were applied to identify lncRNAs, and a total of 6,579 and 9,048 were predicted in the FL transcripts of northern and southernO. oratoria, respectively. A Venn diagram was generated to visualize the respective contribution of different approaches to lncRNA prediction, and 565 and 1,624 were shared among the four approaches for northern and southernO. oratoria, respectively (Figure S4）.The linkages between the identified lncRNAs and CDSs were further checked. The lncRNAs were only related to 6,415 (34.3%) CDSs in the northernO. oratoriaand 8614 (40.5%) in the southernO. oratoria, suggesting that lncRNAs were not fully determined.

3.4 Functional annotation ofO. oratoriafull-length transcriptome

To obtain more comprehensive genetic information ofO. oratoria, we combined the FLNC reads of the two FL transcriptomes to produce 123,824 polished FL consensus isoforms with N50 of 3,250 bp, and 42,735 unigenes with N50 of 3,472 bp after removing redundant sequences (Table S3）.A total of 33,741 (80.0%) unigenes had at least one significant hit in the Nr, Nt, SwissProt, KOG, GO, KEGG or Pfam database, and a total of 8,517 (19.9%) unigenes were annotated in all databases (Figure S5A）.通过查询Nr数据库,我们found that the largest number of unigenes (11,309, 33.5%) were aligned toHyalella azteca, followed byZootermopsis nevadensis(1,616, 4.8%) andLimulus polyphemus(1,077, 3.2%) (Figure S5B）.In KOG analysis, the main function classifications were found to be ‘General function prediction only’, ‘Signal transduction mechanisms’, ‘Cytoskeleton’ and ‘Posttranslational modification, protein turnover, chaperones’ (Figure S5C）.In addition, 31,883 unigenes were mapped to 357 KEGG pathways and clustered significantly in ‘signal transduction’ and ‘transport and catabolism’ (Figure S5D）.We further used the GO classification system to functionally categorize the unigenes based on Nr annotations. In total, 24,028 (56.2%) unigenes were successful annotated to three GO categories, and the largest category was biological process (49,789, 43.6%) followed by cellular component (34,354, 30.1%) and molecular function (29,974, 26.3%) (Figure S5E）.Notably, in the biological process category, we found that 2,866 genes were involved in response to stimulus, 212 genes were related to immune system process, and 464 genes were associated with reproduction. In the molecular function category,O. oratoriahad 48 genes related to antioxidant activity. These annotation and classification will aid understanding of the gene function inO. oratoria。

3.5 Sequence divergence between northern and southernO. oratoriaorthologous genes

A total of 2,182 orthologous gene pairs were identified between northern and southernO. oratoria,系统发生的树了other two arthropodsDaphnia pulexandPenaeus vannamei(Figure 3A）.The UTR regions of each unigene pair was identified based on the predicted coding region. Among the 2182 pairs of orthologs, 884 pairs contain 5’UTR, coding and 3’UTR regions. Differences between coding regions of northern and southern orthologous genes occur at 1.44% of the positions, while the overall difference of 5’UTR and 3’UTR between northern and southernO. oratoriais 2.79% and 1.46%, respectively (Table S4）.

3.6 Identification of genes under positive selection

The ratio of Ka/Ks is an indicator of selection acting on a protein-coding gene. Among the 2,182 pairs of orthologs identified, both Ka and Ks could be calculated for 1,509 (69.2%) orthologs. For the remaining orthologs, we could only calculate either Ka (84 orthologs, 3.8%) or Ks (477 orthologs, 21.9%), or the orthologs were identical (112, 5.1%). For the orthologous pairs for which Ka/Ks ratio could be calculated, the mean values of Ka, Ks, and Ka/Ks were 0.452, 0.013, and 0.951, respectively. Among them, 98 orthologs (6.5%) had a Ka/Ks > 1, indicating these genes have evolved under positive selection (Figure 3B,Table S5）.Except the orthologous pairs with high Ka/Ks, 361 pairs (24.0%) had a Ka/Ks < 0.1, suggesting that these genes have evolved under high selective constraint (Figure 3B,Table S5）.

GO enrichment analysis of 98 genes under positive selection showed that 50 GO terms covering 23 genes were significantly over-presented (P< 0.05) (Table S6）.The top enriched gene groups were mainly related to cellular ion homeostasis (GO:0048878, GO:0050801, GO:0055082) and regulation of cell communication (GO:0010646) (Figure 3C）.Notably, GO terms (GO:0010628, GO:0045893) associated with transcription regulation were also significantly enriched. Although not enriched in the GO analysis, the following genes under positive selection may also play important roles in local adaptation inO. oratoria(Table 2): two gene involved in cytoskeletal organization (Rho guanine nucleotide exchange factor 7, ARHGEF7;cell adhesion molecule 2, CADM2), six genes involved in stress response (glutathione peroxidase 3, GPx3;ankyrin repeat and zinc finger domain-containing protein 1, ANKZF1; two copies offerritinandpoly [ADP-ribose] polymerase 12, PARP12), five genes involved in innate immune response (L-type lectin, LTL;fibrinogen-related protein 2, FREP2;prophenoloxidase activating factor, PPAF; two copies ofantilipopolysaccharide factor, ALFs), five genes involved in genetic information processing (reverse transcriptase;ribonuclease Z;zinc finger protein 271, ZNF271; two copies ofzinc finger protein 791, ZNF791).

3.7 Expression analysis of PSGs

The publicly available RNA-seq data of two ecologically divergentO. oratoriapopulations (corresponding to northern and southernO. oratoriain this study) under thermal stress (Lou et al., 2019a) was used for expression analysis of PSGs. We firstly compared gene expression between genes under positive selection and those under neutral and purifying selection. Interestingly, relatively high levels of gene expression were consistently found in PSGs in comparison with genes under neutral and purifying selection regardless of experimental condition (Figure 4）.Then, we focused on the expression profiles of PSGs in the northern and southernO. oratoriain response to the same heat stress (20°C–28°C). The expression heatmap showed remarkably different expression patterns of 98 PSGs (Figure 5A）.We found specific clusters of highly expressed PSGs in the northern and southernO. oratoriaexposed to heat stress. When considered separately for population-specific cluster, there was a significantly higher number of up-regulated PSGs under heat stress in the northernO. oratoria(28/98) than its southern counterpart (4/98). The highly expressed PSGs involved in “response to stimulus” (GO:0051716), “signaling” (GO:0035556) and “biological regulation” (GO:0008199) were only detected in heat-stressed northernO. oratoria(Figure 5B）.However, only two proteins with known function were abundant in heat-stressed southernO. oratoria, a ribosome-binding protein and a CD63 antigen.

4 Discussion

4.1 Long read reference reconstruction of the full-length transcripts

This is the first FL transcriptome study on the mantis shrimp. Deep PacBio sequencing is currently one of the best methods for retrieving nucleic acid sequences from species with ultra-large and complex genomes. In this study, the FL transcriptome ofO. oratoriahas been deep-sequenced using the PacBio Sequel platform with a pooled RNA sample from different tissues to capture as many FL transcript isoforms as possible. Compared with thede novoassembled unigenes obtained from RNA-seq, the length is greatly longer in the non-assembled unigenes generated by PacBio SMRT sequencing. The annotation rate (80.0%) of PacBio unigenes is much higher than that reported in previous next-generation transcriptome studies (20.4%-33.7%,Yan et al., 2018;Lou et al., 2019b;Lou et al., 2019c）.Although fewer unigenes were obtained by PacBio sequencing than by RNA-seq, a higher proportion of PacBio unigenes have been annotated than Illumina unigenes, indicating high-efficiency of PacBio SMRT sequencing in recovering FL transcripts. Thus, the FL transcripts produced in this study would serve as important basis for gene discovery, genome assembly and annotation ofO. oratoria, and contribute towards better understanding of adaptive evolution in this species.

TFs play key roles in the transcriptional regulation of functional genes by sequence-specific binding to regulatory regions throughout the genome (Fulton et al., 2009）.In this study, we have predicted 1,549 TFs in the northernO. oratoriaand 944 TFs in the southernO. oratoria, laying a foundation for the next step in studying the potential regulatory roles of TFs inO. oratoria。Among them, zf-C2H2 and ZBTB are the most abundant TF families. C2H2 zinc-finger proteins are one of the most widespread transcription factor families in eukaryotes (Kubo et al., 1998), and participate in a wide range of biological processes, such as growth, development, and stress response (Liu et al., 2015）.ZBTB TFs are generally considered transcriptional repressors, involving in a variety of developmental and cellular processes (Ren et al., 2019）.This TF family has been shown to involve in developmental regulation of regenerative potential inDrosophila(Narbonne-Reveau and Maurange, 2019）.

LncRNAs are a group of RNA molecules with the structure similar to mRNA and the length longer than 200 nucleotides (Rinn and Chang, 2012）.Despite the lack of protein-coding capacity, LncRNAs can exert powerful biological functions by regulating gene expression at epigenetic, transcriptional and post-transcriptional levels (Kapranov et al., 2007;Cao, 2014）.To date, the LncRNAs ofO. oratoriahas not been reported yet. We have identified 565 and 1,624 LncRNAs for northern and southernO. oratoria, respectively. Considering the important roles of lncRNAs in various biological processes regulation, our data provide a reference resource for further investigation on the underlying mechanism of lncRNA-related regulation inO. oratoria。

4.2 Environmental heterogeneity driving adaptive divergence inO. oratoria

For theO. oratoriacomplex, two cryptic species have been delineated based on a number of phylogenetic analysis (Cheng and Sha, 2017;Cheng et al., 2020）.However, the genetic factors driving the evolution ofO. oratoriaare almost unknown due to the fact that few molecular data are available for inferring the evolutionary history and genetic divergence of this species. In this study, we have identified 2,182 orthologous gene pairs between northern and southernO. oratoriaand determined the average sequence divergence to be 1.44% for the coding regions. This is comparable to the average divergence of two invasive whitefly species Asia II 3 and MEAM1 transcriptomes (1.73% inWang et al., 2012), and much higher than the level of sequence divergence that was found between human and chimpanzee (0.45% inHellmann et al., 2003) as well as between the two other whitefly cryptic species MEAM1 and MED (0.83% inWang et al., 2011）.When it comes to non-coding regions, the genetic divergence is more obvious for 5’UTR (2.79%) and 3’UTR (1.46%) regions, for which both ratios are higher than the reported mean 1.12% divergence of 5’UTR and 0.86% divergence of 3’UTR regions between human and chimpanzee (Hellmann et al., 2003）.The relatively higher level of similarity at the coding regions compared to non-coding regions could be attributed to the presence of functional elements that are subject to purifying selection (Wang et al., 2011）.Collectively, the substantial sequence divergence at the transcriptomic level in combination with previous phylogenetic results concur in suggesting that northern and southernO. oratoriaconstitute different species.

The heterogeneous seascapes in the NWP, such as temperature gradient governed by the oceanic currents system and salinity gradient surrounding the Yangtze River mouth, can impose spatially varying selective pressures on local populations distributed along the gradients. As expected, our analyses have identified a total of 98 orthologs under positive selection (Ka/Ks > 1), suggesting that those genes might play important roles in the speciation and adaptive evolution ofO. oratoria。The candidate genes have been found to be involved in many biological processes, including stress response, immunity, and cytoskeletal organization. It should be noted that compared to a battery of thermally responsive genes, the relatively low number of salinity-tolerance associated genes identified here may be due to limited genomic information provided by transcriptome sequencing or relatively weak selective force generated by salinity gradient, or a combination of these. As discussed below, the functional inferences based on candidate genes suggest that seascape heterogeneity in the NWP, in particular differences in temperature, has facilitated adaptive evolution ofO. oratoria。

It has been clearly established that any intense stress is usually accompanied by enhanced reactive oxygen species (ROS) generation, causing oxidative damage (Lushchak, 2011）.The role of changes in temperature and salinity in induction of oxidative stress is highlighted in aquatic organisms (Liu et al., 2007;Madeira et al., 2013）.In response to the production of ROS, these organisms have developed the cellular scavenging system to remove these radicals involving various enzymes. GPx3 is a secreted plasma protein that scavenges ROS in the extracellular compartment, and thereby protects cells against oxidative damage (Jin et al., 2011）.ANKZF1 involves in the cellular response to hydrogen peroxide and plays a role in the maintenance of mitochondrial integrity under conditions of cellular stress (van Haaften-Visser et al., 2017）.Ferritin is a highly conserved iron storage protein that can concentrate cellular iron to prevent the harmful ROS generation and reduce oxidative damage (Harrison and Arosio, 1996;McCord, 1996）.Previous findings have suggested the roles of ferritin in salinity stress adaptation in the swimming crabPortunus trituberculatus(Huang and Xu, 2016）.A ROS signature in cells can also trigger the activation of DNA repair (Mittler et al., 2011）.RAPA proteins are involved in the detection of DNA damage and initiation of DNA repair by binding to damaged parts of DNA (Druzhyna et al., 2000）.Therefore, selection on GPx3, ANKZF1, ferritin and RAPA12 may be associated with adaptation ofO. oratoriain response to oxidative damage produced by temperature and salinity variation.

Water temperature is probably the most important environmental variable because it directly affects growth, metabolism and survival of marine organisms (Chen et al., 1995;Hennig and Andreatta, 1998）.In the crabCarcinus maenas,Chisholm and Smith (1994)found the impact of the seasonal changes in water temperature on the antibacterial activity of haemocytes. Notably, representative genes involved in innate immune were found to be positively selected in warm- and cold-tolerantO. oratoria。Enzymes encoded by these genes include LTL that functions as pattern recognition receptors in shrimp immunity (Xu et al., 2014), FREP2 which is also known as fibrinogen-like domain immunolectins and plays a role in defense and protection against infection (Hanington and Zhang, 2011), PPAF involving in the prophenoloxidase activation pathway of innate immune system in invertebrates (Wang et al., 2015), and ALF which is a type of antimicrobial peptides with a vital role in crustacean antimicrobial defense (Li and Li, 2020）.Previous common-garden experiments also suggested that differentially expressed genes after heat treatment inO. oratoriawere enriched in pathways associated with immune (Lou et al., 2019a）.Consequently, we hypothesize that these functional genes associated with immunity seem to evolve rapidly inO. oratoriato increase its capacity to manage thermal stresses. Cytoskeletal reorganization has been reported to undergo pronounced transformation under thermal stress (Richter et al., 2010）.Herein, evolutionary forces acting on cytoskeletal organization-related genes (ARHGEF7 and CADM2) provided additional evidence reflective of thermal selection inO. oratoria。ARHGEF7可以调节ρGTPase活性和解放军ys a role in regulating cytoskeletal and cell adhesion dynamics (Cheng et al., 2021), while CADM2 is a cell adhesion molecule that functions in cell recognition, adhesion, migration and differentiation (Edelman, 1983）.年代election on genes involved in cytoskeletal organization has also been implicated in the warm-adapted marine musselMytilus galloprovincialiscompared with its cold-adapted congeners (Popovic and Riginos, 2020）.Collectively, these findings suggest that adaptation was driving the evolution of key genes associate with immune and cytoskeletal organization, potentially indicating their importance forO. oratoria宽容的栖息地的温度变化。进一步validation of these candidates is necessary to determine their functional association with local thermal adaptation ofO. oratoria。

4.3 Selection towards plasticity facilitating thermal adaptation ofO. oratoria

Genetic variation and phenotypic plasticity both play important roles in adaptive evolution (Davis et al., 2005;Kelly, 2019）.In particular, plasticity in gene expression is favored by organisms under dynamic environments (李et al ., 2018）.Closely related species and genetically differentiated populations may exhibit differences in the patterns and extents of plastic responses, indicating there is genetic variation for plastic responses (Gao et al., 2008）.This means that phenotypic plasticity can evolve in response to natural selection, which in turn suggests the existence of evolving plasticity (Pigliucci, 2005）.As observed in the polar diatomFragilariopsis cylindrus(Mock et al., 2017), we found PSGs inO. oratoriahad significantly higher expression levels than those genes under neutral or purifying selection, indicative of a role of natural selection in driving gene expression.

Comparing gene expression responses to environmental conditions among locally adapted populations offers an opportunity to test for among-population variation in plasticity (DeBiasse and Kelly, 2016）.The northern and southernO. oratoriaare naturally distributed in different climate regions in the NWP and may vary in their sensitivity to thermal stress. As expected, they exhibit differences in the patterns of expression plastic response to the same heat stress (Figure 5）.The northernO. oratoria适应寒冷气候differenti更大al gene expression in PSGs between control and heat-stressed samples than its southern counterpart, suggesting more sensitive of northernO. oratoriato heat stress. The higher expression plasticity in the heat-sensitive populations has also been observed in copepodTigriopus californicus(年代choville et al., 2012), coralPorites astreoides(Kenkel et al., 2013) and snailChlorostoma funebralis(Gleason and Burton, 2015）.The observation of population-specific clusters of differentially expressed PSGs implies that natural selection may act on transcriptional plasticity to facilitate the evolution of lineage specific tolerance to heat stress inO. oratoria。年代imilarly, selection towards expression plasticity has been demonstrated in other aquatic animals, such as oyster populations from different environmental gradients (李et al ., 2018;Li et al., 2021）.Along with previous studies, our results point to the potential importance of evolving plasticity in adaptation of marine organisms to heterogeneous seascapes.

5 Conclusion

In this study, we have demonstrated the advantage of PacBio sequencing to rapidly reconstruct FL transcripts compared to RNA-seq. By comparing transcriptome sequence divergence, our results support previous proposition that northern and southernO. oratoriaconstitute different species. The abundant PSGs related to environmental responsiveness are indicative of a strong selective force thatO. oratoriamight undergo during the adaptation process to heterogeneous seascapes in the NWP. In addition, we found genes underwent positive selection also exhibit divergence in expression plastic response to heat stress, suggesting that natural selection may act on the expression plasticity to facilitateO. oratoriaadaptation. As such, our study strengthens the view that genetic variation and plasticity should be taken into account when attempting to understand the evolution of marine species distributed along the environmentally heterogeneous coast. In light of the important role of evolving plasticity in species’ responses to climate change, mechanistic studies are warranted to investigate the underlying processes that may mediate different adaptive potential, which is critical forO. oratoriain the context of a changing climate.

Data availability statement

Raw sequencing data analyzed for this study have been submitted to NCBI Sequence Read Archive (SRA) under the BioProject accession number PRJNA853060.

Author contributions

JC and ZS conceived and designed the project. JC and YL collected the samples. JC conducted laboratory work. JC and LZ performed bioinformatics and statistical analyses. JC authored drafts of the manuscript. MH and LZ reviewed drafts of the paper. All authors contributed to the manuscript and approved it for submission and publication.

Funding

This work was supported by the National Natural Science Foundation of China (No. 31872569), the Open Fund of CAS Key Laboratory of Experimental Marine Biology, Institute of Oceanology, Chinese Academy of Sciences (No. KF2021N002), the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDB42000000) and the National Key R&D Program of China (No. 2018YFD0900804).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

年代upplementary material

The Supplementary Material for this article can be found online at://www.gosselinpr.com/articles/10.3389/fmars.2022.975686/full#supplementary-material

年代upplementary Figure 1 |Length distributions of full-length transcriptomes of the northern and southernO. oratoria。(A)Distribution of the number and length of CCS reads.(B)Distribution of the number and length of FLNC sequences.(C)Distribution of the number and length of consensus isoforms.

年代upplementary Figure 2 |Length distribution of CDSs of the northern(A)and southern(B)O. oratoria。的长度预测沿着x - cd被绘制axis, while number of CDS transcripts was plotted along the left y-axis. The yellow line that represents the percentage of CDS length was plotted along the right y-axis.

年代upplementary Figure 3 |Transcription factors identified in the northern and southernO. oratoria。(A)Classification of the detected TF families. Different types of transcript family were plotted along the x-axis, while the number of transcription factors were plotted along the y-axis.(B)GO enrichment reveals the functions of genes that have been determined to be TFs.

年代upplementary Figure 4 |Venn diagram of the number of lncRNAs predicted in the northern(A)and southern(B)O. oratoria。

年代upplementary Figure 5 |Functional annotation ofO. oratoriafull-length transcriptome.(A)Classification of unigene annotation in all databases, including Nr, SwissProt, KEGG, KOG, GO, Nt and Pfam.(B)年代pecies distribution of highest scoring blastp match in Nr database.(C)KOG function classification of allO. oratoriaunigenes.(D)KEGG pathway assignment ofO. oratoriaunigenes.(E)Distribution of GO terms for all annotated unigenes in biological process, cellular component and molecular function categories.

References