Massive Heat Shock Protein 70 Genes Expansion and Transcriptional Signatures Uncover Hard Clam Adaptations to Heat and Hypoxia

Introduction

Heat shock proteins (HSPs) are a group of soluble intracellular proteins produced by cells in response to various environmental stresses (Welch, 1993). In the absence of stress, the content of HSPs is ~5% of the total intercellular protein content, but it can rapidly increase up to 15% or more after exposure to stress (Srivastava, 2002). Based on their molecular weights (MWs), HSPs can be divided into various families, such as small HSPs, HSP40, HSP70, HSP90, and HSP110 (Lindquist and Craig, 1988). HSP70, also known as HSPa, is an abundant heat-inducible 70 kDa HSP encoded by HSP70 and is evolutionarily conserved in terms of both structure and function (Lindquist and Craig, 1988; Daugaard et al., 2007). In general, HSP70 has three motifs: an N-terminal adenosine triphosphatase domain (ATPase; ~400 aa), a substrate-binding domain (SBD; ~180 aa), and a variable-length C-terminal domain (Nikolaidis and Nei, 2004). HSP70 first received attention for its role in the response to heat shock (Ritossa, 1962) and was also found to be critical for the folding and assembly of other cellular proteins (Gething and Sambrook, 1992). Some HSP70 genes are expressed only under stress conditions, but not all the members are induced under stress conditions; they are necessary for cell viability under normal conditions (Daugaard et al., 2007; Murphy, 2013). Due to their organelle- and tissue-specific expression, different members of the HSP70 family have distinct biological functions, although they might function redundantly in some cases (Daugaard et al., 2007; Murphy, 2013). In summary, the HSP70 family functions in buffering against extrinsic and intrinsic stress (Murphy, 2013).

The HSP70 gene family has been studied in humans (Brocchieri et al., 2008), bovines (Tripathy et al., 2021), fish (Song et al., 2016; Xu et al., 2018), the fruit fly Drosophila melanogaster (Gong and Golic, 2004), and nematodes (Heschl and Baillie, 1990; Nikolaidis and Nei, 2004; Guerin et al., 2019). Mollusca is the second largest phylum in Animalia, and as the largest group of species in the Lophotrochozoa, it is central to our understanding of the biology and evolution of this superphylum of protostomes (Zhang et al., 2012). Bivalves, which are the second class in the Mollusca, are among the oldest and evolutionarily most successful groups of invertebrates (Wang et al., 2013). They are widely distributed from intertidal areas to the deep sea and from polar to tropical regions (Li et al., 2017). Several members of HSP70 have been cloned in several bivalve species and were upregulated after exposure to heat stress and other stresses (Franzellitti and Fabbri, 2005; Liu et al., 2014; Cheng et al., 2019). As an increasing number of genomic resources, the genome-wide HSP70 gene family has been elementary studied in Pteriomorphia and Veneridae. HSP70 genes are expanded in the oysters Crassostrea gigas (Zhang et al., 2012) and C. hongkongensis (Peng et al., 2020) the Sydney rock oyster Saccostrea glomerata (Powell et al., 2018), the pearl oyster Pinctada fucata (Takeuchi et al., 2016), the scallop Patinopecten yessoensis (Cheng et al., 2016), the Manila clam Ruditapes philippinarum (Yan et al., 2019), Cyclina sinensis (Wei et al., 2020), the deep-sea vent mussel Bathymodiolus platifrons (Sun et al., 2017) and the shallow-water mussel Modiolus philippinarum (Sun et al., 2017), suggesting that these genes play important roles in adaptation to heat and other stress factors in dynamic environments with a wide variety of stress factors. As a case of molecular convergent evolution, HSP70 family members independently expanded across multiple lineages to mitigate environmental stress (Guerin et al., 2019). Phylogenetic analysis showed recent species-specific HSP70 expansions along with older mixed-species clusters dating to ancestral bivalve lineages, which was consistent with natural selection favouring the expansion (Guerin et al., 2019).

The hard clam, Mercenaria mercenaria, has a burrowing lifestyle typical of bivalves and naturally lives along the East Coast of the United States and Canada (Menzel, 2009). It has become an emerging pond culture species since being imported to China. During pond culture, warm water and low oxygen are major environmental stresses in the summer. Heat is commonly accompanied by hypoxia in the ponds. In our previous transcriptional study, most HSPa12 transcripts showed high expression in the hypoxia-stress groups, while other HSPa subfamilies (including B2, 4, 5, 8, 9 and 14) showed high expression in the heat combined hypoxia stress groups (Hu et al., 2022). These results indicated that HSP70 has a vital function in protecting cells against heat and hypoxia stress (Hu et al., 2022). However, the HSP70 gene family in this hardy species remains poorly understood.

In this study, we comprehensively analysed the chromosomal locations, duplication, structure and motif compositions of 133 HSP70 genes based on the recently completed genome sequence of the M. mercenaria. We also analysed the evolutionary relationships of HSP70 genes by performing phylogenetic and collinearity analyses. In addition, we studied the expression of HSP70 genes under thermal and hypoxia stress. This study provides a comprehensive understanding of the HSP70 gene family in the M. mercenaria and lays a significant foundation for further studies on the functional characteristics of MmHSP70 genes during exposure to acute thermal and hypoxia stress.

Materials and Methods

Identification and Sequence Analysis of HSP70 Genes

The M. mercenaria genome was downloaded from the National Center for Biotechnology Information (BioProject number: PRJNA596049). Gene family identification methods include based on sequence homology (Blast) and conserved domain (HMMER). In this study, we used combined methods to identify M. mercenaria HSP70 gene family. Firstly, HSP70 protein sequences from the oyster C. gigas, which were downloaded from UniProt (https://www.uniprot.org/), were treated as a query database to perform Basic Local Alignment Search Tool algorithm program (BLASTP) search with an E-value ≤ e-5 and identity ≥ 30% against the M. mercenaria genome. These candidate genes were believed to have sequence similarities with C. gigas HSP70. Then, the candidate genes were filtered by conserved domain. HMMER is a common domain prediction software. Cut_tc is the strictest option controlling model-specific thresholding. The hidden Markov model (HMM) profile of the HSP70 domain (PF00072) was downloaded from the Pfam protein family database (https://pfam.xfam.org/). The putative M. mercenaria HSP70 genes were filtered using HMMER with cut_tc algorithm. Finally, protein sequences with fewer than 300 amino acids were excluded from further analyses. The putative MmHSP70 genes functionally annotated as HSP70 in Swiss-prot database were manually selected as the final identified M. mercenaria HSP70. The ExPASy ProtParam tool (https://web.expasy.org/protparam/) was used to analyse molecular weights (MWs).

To study the evolution of HSP70 in Veneridae, another clam, C. sinensis, was used for comparison. Whole-genome proteins of C. sinensis were downloaded from CNGBdb (https://db.cngb.org/, CNA0003280) (Wei et al., 2020). The CsHSP70 genes were identified as described above.

Chromosomal Distribution, Gene Structure and Conserved Motif Characterization

HSP70 genes were mapped to the chromosomes and contigs according to the Generic Feature Format (GFF) file, and the results were visualized by gene location visualization of the GTF/GFF file in TBtools (Chen et al., 2020). Those genes were named according to their position on the chromosome. MCScanX was applied to search for gene duplicate types in the HSP70 family (Wang et al., 2012). The whole-genome protein sequences were compared in pairs by Diamond software with the parameters max-target-seqs 5 and evalue 1e-10 (Buchfink et al., 2015). Calculator 2.0 software was employed to calculate the Ka and Ks values of tandem HSP70s (Wang et al., 2010). HSP70 gene structure was investigated based on the coding sequence (exon), untranslated region, and intron data from the GFF file and visualized by TBtools (Chen et al., 2020). Conserved motif analysis of the HSP70 proteins was conducted by MEME (5.4.1), and the parameters were set to anr mode, an optimum mode width of 6 to 200 and a maximum number of motifs of 10 (Bailey et al., 2015).

Sequence Alignment and Phylogenetic Analyses

We perform two alignment process based on protein sequences. The first one was based on 133 M. mercenaria HSP70 protein sequences, the second one was based on 133 M. mercenaria HSP70 protein sequences and 60 C. sinensis HSP70 protein sequences. Multiple sequence alignment of HSP70 protein sequences was carried out with the ClustalW algorithm in MEGA 7.0 software (Kumar et al., 2016). Two phylogenetic trees (one was M. mercenaria HSP70 tree, the other was M. mercenaria and C. sinensis HSP70 tree) were constructed using the neighbour-joining (NJ) method and 1000 bootstrap replicates. Gap data Treatment was pairwise deletion.

Synteny Analysis

Diamond software was used to perform two-way BLASTP between the M. mercenaria whole-genome protein sequences and those of C. sinensis (Buchfink et al., 2015). In details, M. mercenaria whole-genome protein sequence file was firstly blastp with C. sinensis whole-genome protein sequence file. Then, C. sinensis whole-genome protein sequence file was blastp with M. mercenaria whole-genome protein sequence file. The blastp parameters were as following: max-target-seqs was 5 and evalue was 1e-10. The two blastp files were merged into one blast file named Mme_Csi.Blast. The M. mercenaria and C. sinensis whole genome GFF files were also merged into one GFF file named Mme_Csi.gff. MCScanX was used to search for duplicate genes between the M. mercenaria and C. sinensis (Wang et al., 2012). The Mme_Csi.Balst file and Mme_Csi.gff file were used as the input to MCScanX software. The results were visualized by JCVI. The syntenic HSP70 gene pairs were highlight in red lines.

Tissue Expression of MmHSP70 Genes

The healthy adult hard clam tissue transcriptome libraries were constructed in our previous study including testis, ovary, mantle, gill, foot, intestine, liver, stomach, adductor muscle, and hemolymph (Song et al., 2021). The raw data have been deposited in the SRA of NCBI with the accession number PRJNA596049. Raw data were filtered using in-house Perl scripts to remove adapter reads, poly-N sequences and low-quality reads. After quality control, clean reads were mapped to the hard clam genome by HISAT v 2.0.4 (default parameters). HTSeq v0.6.1 (union model) was used to calculate the numbers of reads mapped to each gene. The per kilobase per million mapped reads (FPKM) value of each gene was calculated. Transcript abundance was evaluated using FPKM values. A heatmap was generated using OmicShare tools, a free online platform for data analysis (http://www.omicshare.com/tools).

Expression Analysis of MmHSP70 Genes Under Heat and Hypoxia Stress

Transcriptome libraries of M. mercenaria exposed to heat, hypoxia, and heat combined with hypoxia were constructed in our previous study (Hu et al., 2022). In brief, heat, hypoxia and combined stress challenge experiments were carried out using a novel hypoxia simulation device (Li et al., 2019). The M. mercenaria were abruptly exposed to the experimental temperature and dissolved oxygen conditions and remained in these experimental conditions for 3 days. Six groups were established for the six treatments in the experiment: 20°C, 6 mg/L DO (control, C_6); 20°C, 2 mg/L DO (moderate hypoxia, C_2); 20°C, 0.2 mg/L DO (severe hypoxia, C_02); 35°C, 6 mg/L DO (heat, H_6); 35°C, 2 mg/L DO (heat plus moderate hypoxia, H_2); and 35°C, 0.2 mg/L DO (heat plus severe hypoxia, H_02). All M. mercenaria under all conditions were under stress for 3 days. Gills were selected as the target tissue for RNA-seq. RNA-seq data (PRJNA764366 and PRJNA764372) were used to examine MmHSP70 gene expression profiles in response to heat and hypoxia stress in the M. mercenaria (Hu et al., 2022). The similar methods were used to evaluate transcript abundance as above.

Results

Identification and Sequence Analysis of HSP70 Genes

Compared to those of previously studied bivalves, including the oyster C. gigas (Zhang et al., 2012) and the scallops P. yessoensis (Cheng et al., 2016) and C. farreri (Hu et al., 2019), the genome of the M. mercenaria had more HSP70 genes. In the present study, 133 MmHSP70 genes were identified in the genome of M. mercenaria by BLASTP and HMM methods, which gives the M. mercenaria the largest HSP70 gene repertoire in bivalve, to the best of our knowledge. The M. mercenaria MmHSP70 genes included 4 HSPa B2 genes, 1 HSPa4 gene, 1 HSPa5 gene, 1 HSPa8 gene, 2 HSPa9 genes, 122 HSPa12 genes, 1 HSPa14 gene and 1 HSPa17 gene. The expanded HSP70 genes mainly included 122 (91.73%) HSPa12 genes, in accordance with findings in the oyster C. gigas (Zhang et al., 2012) and the scallops P. yessoensis (Cheng et al., 2016) and C. farreri (Hu et al., 2019). Information including genome location, amino acid length, annotation and MW is summarized in Supplementary File 1. The encoded proteins ranged from 303 aa (MmHSPa12_90) to 1387 aa (MmHSPa12_65). The predicted MWs of the HSP70 proteins varied from 33.78 kDa (MmHSPa12_90) to 156.41 kDa (MmHSPa12_65). Moreover, 60 CsHSP70 genes were identified in the C. sinensis genome, 47 (78.33%) of which were HSPa12 genes.

Chromosomal Distribution and Duplication of HSP70 Genes

A total of 133 MmHSP70 genes were unevenly distributed among 19 M. mercenaria chromosomes and 4 scaffolds (Figure 1). Among the 133 HSP70 genes, 41 (33.08%) were located on Chr 7. There were fewer than 11 HSP70 genes on other chromosomes and contigs. Only one HSP70 gene was observed on Chr 6, Chr 10 and Chr 19. Similarly, 60 CsHSP70 genes were distributed on 14 chromosomes, and 28 genes were located on Chr 14.

FIGURE 1

Figure 1 Chromosomal distribution of Mercenaria mercenaria HSP70 genes. The HSPa12 genes were highlight in red.

Forty-five and thirty-five MmHSP70 genes were tandemly duplicated and proximally duplicated, respectively, which accounted for 60.15% of the HSP70 gene family. Twenty-six and seven CsHSP70 genes were tandemly duplicated and proximally duplicated, respectively. In total, 17 tandemly duplicated HSP70 gene pairs were found. There were two groups of three tandemly duplicated gene pairs on Chr 2 and Chr 4. Genes in the same tandemly duplicated pairs had short physical distances. The Ka and Ks values of tandemly duplicated HSP70 gene pairs were calculated by Calculator 2.0 with the MA method (Wang et al., 2010) to determine if they experienced various selection pressures (Song et al., 2021). The Ka/Ks values of all tandem duplication HSP70 gene pairs were significantly less than 1 (ranging from 0.0548205 to 0.417659, P < 0.05, Table 1). All tandem pair HSP70 genes had multiple exons. MmHSPa B2_1, MmHSPa B2_2 and MmHSPa B2_3 had similar gene structures (Figure 2). A total of 10 conserved motifs were detected in tandem pair HSP70 genes, whose length ranged from 29 to 200 aa. Highly similar conserved motifs were found in the same pairs of tandem duplicated genes (Figure 3). Moreover, the tandemly duplicated gene pairs on Chr 7 also had a similar motif pattern (Figure 3).

TABLE 1

Table 1 Selection pressure on the tandem duplication HSP70 gene pairs.

FIGURE 2

Figure 2 Gene structure of M. mercenaria tandem duplication MmHSP70 gene pairs. Green boxes, black lines and yellow boxes represent exons, introns and untranslated region (UTR), respectively. The tandemly duplicated gene pairs were shown by adjacent two or three genes with the same colours.

FIGURE 3

Figure 3 Distribution of conserved motifs in tandem duplication MmHSP70 gene pairs. The tandemly duplicated gene pairs were shown by adjacent two or three genes with the same colours.

Phylogenetic Analysis of HSP70 Genes

To study the evolutionary relationships of HSP70, two NJ phylogenetic trees were constructed based on the protein sequences. In Figure 4, the tree showed that the M. mercenaria HSP70 gene family mainly had two clusters, including the HSPa12 subfamily and other subfamilies. One ancestral branch consisted of 11 HSP70 genes, including HSPa4, HSPa5, HSPa8, HSPa9, HSPa14, HSPa17 and HSPa B2. The M. mercenaria HSPa B2s were first clustered together. The tandemly duplicated HSPa9 gene pairs were also clustered together. The other branch consisted of 122 HSPa12 subfamily members (Figure 4). Another NJ phylogenetic tree of HSP70 protein sequences from the M. mercenaria and C. sinensis was also constructed (Figure 5). These two trees had similar topological structures. The ancestral branch consisted of 24 HSP70 genes, including HSPa4, HSPa5, HSPa8, HSPa9, HSPa14, HSPa17 and HSPa B2. The clam HSPa4, HSPa5, HSPa8, HSPa9, HSPa14, HSPa17 and HSPa B2 were clustered together. The expanded HSPa12 branch consist of 122 M. mercenaria HSPa12 genes and 47 C. sinensis HSPa12 genes. Species-specific HSPa12 expansions could also be observed.

FIGURE 4

Figure 4 Phylogenetic analysis of 133 M. mercenaria HSP70 protein sequences. The phylogenetic tree was created with MEGA 7.0 software using the Neighbor-Joining method. The HSPa12s were highlight in red. Support value was assessed using bootstrapping (100 pseudoreplicates), if value was less than 50, the data was not shown.

FIGURE 5

Figure 5 Phylogenetic analysis of 193 clam HSP70 protein sequences (133 M. mercenaria HSP70 and 60 Cyclina sinensis HSP70 protein sequences). M. mercenaria HSP70s and C. sinensis HSP70s were marked with red triangles and blue squares, respectively. The HSPa12 branches were highlight in red. Support value was assessed using bootstrapping (100 pseudoreplicates), if value was less than 50, the data was not shown.

Synteny Analysis of HSP70 Genes Between the M. Mercenaria and C. Sinensis

We performed synteny analysis of the HSP70 genes between the M. mercenaria and C. sinensis. The HSP70 genes in the M. mercenaria were homologous to genes in C. sinensis, and syntenic conservation was observed in C. sinensis (22 pairs of orthologous gene pairs distributed on several chromosomes). Interestingly, a high-density collinear gene block on Chr 7 could be clearly seen (Figure 6). Most MmHSP70 genes located on Chr 7 were HSPa12 genes. These results indicated that the M. mercenaria HSP70 gene family was conserved and that the HSPa12 genes expanded on Chr 7 of the M. mercenaria genome might have evolved from those of a common ancestor with C. sinensis.

FIGURE 6

Figure 6 Synteny analyses between the HSP70 genes of the M. mercenaria and C. sinensis. Grey lines in the background indicate collinear blocks in the M. mercenaria and C. sinensis genomes, while red lines highlight syntenic HSP70 gene pairs.

Expression Patterns of MmHSP70 Genes in Different Tissues

RNA-seq data of ten tissues, namely, the testis, ovary, mantle, gill, foot, intestine, hepatopancreas, stomach, adductor muscle, and hemolymph, from healthy adult M. mercenaria were used to characterize the expression profiles of MmHSP70s (Song et al., 2021). The expression data was shown in Supplementary File 2. In general, the average expression of HSP70 genes was highest in the hemolymph and gill and lowest in the foot and adductor muscle. Based on FPKM values, a heatmap of HSP70 genes in various tissues was created (Figure 7). Interestingly, the expression pattern of HSP70s in the gill was different from that in the other tested tissues. Some HSP70 genes showed highly tissue-specific expression. For example, MmHSPa12_67, MmHSPa_81 et al. had high expression in the gill. MmHSPa12_10, MmHSPa12_105 et al. showed high expression in the hemolymph. MmHSPa12_64, MmHSPa12_54 et al. showed high expression in the hepatopancreas. Most HSP70 genes were expressed in at least one tissue.

FIGURE 7

Figure 7 Heatmap analysis of MmHSP70 in ten tissues of adult M. mercenaria based on the FPKM. The name of tissues was abbreviated as Te, testis; Ov, ovary; Ma, mantle; Gi, gill; Fo, foot; In, intestine; Hep, hepatopancreas; St, stomach; Ad, adductor muscle; and Hem, hemolymph. The colour scale represented Z-score.

Expression Profiles of Tandem Duplication MmHSP70 Genes in Response to Heat and Hypoxia Stress

The tightly linked gene modules formed by tandem duplication play an important role in adaptation to environmental stress in bivalves (Zhang et al., 2012; Song et al., 2021; Li et al., 2021). To provide insight into the functions of tandem duplication HSP70 genes in the M. mercenaria, we determined the expression profiles of such genes by using RNA-seq data. The expression data was shown in Supplementary File 3. The tandem duplication HSP70 gene pairs seemed to have similar expression patterns (Figure 8). Most of the gene pairs had high expression under hypoxic stress (C_2, C_02). The gene pairs consisting of MmHSPa B2_1, MmHSPa B2_2 and MmHSPa B2_3 had high expression under heat plus extreme hypoxia stress (H_02). Moreover, the gene pair consisting of MmHSPa12_92 and MmHSPa12_93 also had high expression in the H_02 group.

FIGURE 8

Figure 8 In silico gene expression analysis of M. mercenaria tandem duplication MmHSP70 genes under heat, hypoxia and combined stress. The colour scale represented Z-score. C_6: 20°C, 6 mg/L DO, control; C_2: 20°C, 2 mg/L DO, moderate hypoxia stress; C_02: 20°C, 0.2 mg/L DO, severe hypoxia stress; H_6: 35°C, 6 mg/L DO, heat stress; H_2: 35°C, 2 mg/L DO, heat plus moderate hypoxia stress; H_02: 35°C, 0.2 mg/L DO, heat plus severe hypoxia stress. The tandemly duplicated gene pairs were shown by adjacent two or three genes with the same colours.

Discussion

The HSP70 gene family functions in a wide range of housekeeping and stress-related activities (Rosenzweig et al., 2019). The function of HSP70 family members has been investigated in some bivalves (Zhang et al., 2012; Cheng et al., 2016; Hu et al., 2019), but to date, no detailed analysis of the HSP70 family has been reported in the M. mercenaria. In this study, we performed exhaustive research on the M. mercenaria genome, conducted phylogenetic and synteny analyses with selected bivalve species, and determined expression patterns by examining RNA-seq datasets for M. mercenaria under heat and hypoxia stress. This study provided a comprehensive understanding of the HSP70 gene family in the M. mercenaria and laid a significant foundation for further studies on the functional characteristics of MmHSP70 genes during exposure to heat and hypoxia stress.

HSP70 genes are widespread in all domains of life, and the copy number varies among species. Compared with those in other representative animals, HSP70 genes exhibit expansion in bivalves (Guerin et al., 2019; Hu et al., 2019). In this study, a total of 133 and 60 HSP70 genes were identified in the M. mercenaria and C. sinensis genomes, respectively. Compared with other bivalve species, for example, C. gigas (88) (Zhang et al., 2012), P. yessoensis (61) (Cheng et al., 2016) and C. farreri (65) (Hu et al., 2019), the M. mercenaria had the largest HSP70 gene repertoire, to the best of our knowledge. The HSP70 gene copy number varied greatly among bivalves, which might be a reflection of regulatory variation affecting physiological differences in response to fluctuations in environmental factors, especially temperature (Zhang et al., 2012). There were 8 subfamilies of M. mercenaria HSP70 genes. Among these, 122 HSPa12 genes accounted for the vast majority of HSP70 members, which was consistent with the results of previous reports (Zhang et al., 2012; Cheng et al., 2016; Hu et al., 2019). In both previous studies and our study, HSPa1, HSPa2, HSPa6 and HSPa7 were absent in the bivalve genome, suggesting species-specific loss (Fu et al., 2021). Notably, we identified the HSPa17 gene in the M. mercenaria and C. sinensis genomes, which was not reported in the oyster (Zhang et al., 2012) and scallop genomes (Cheng et al., 2016; Hu et al., 2019). The HSPa17 gene was also identified in the human (Brocchieri et al., 2008), bovine (Tripathy et al., 2021), channel catfish (Song et al., 2016), large yellow croaker (Xu et al., 2018), rainbow trout (Ma and Luo, 2020), and ascidian (Fujikawa et al., 2010) genomes. However, we did not identify HSPa13 in the M. mercenaria genome. MmHSP70 genes were unevenly located on 19 chromosomes and 4 scaffolds, and 41 were located on Chr 7 (Figure 1). In addition, 46.67% of CsHSP70 genes were located on Chr 14. It might be a common phenomenon for members of large gene families to be unequally distributed on chromosomes. For example, 19 (a total of 66) transient receptor potential channel genes were located on Chr 2 in an oyster (Fu et al., 2021), and 59 (a total of 159) inhibitors of apoptosis genes were densely tandemly linked on Chr 5 in the M. mercenaria (Song et al., 2021).

Gene duplication can provide materials for evolutionary novelty (Zhang, 2003; Roth et al., 2007; Song et al., 2021). Tandem duplication is an important driving force of gene duplication (Chen et al., 2013; Long et al., 2013). After tandem duplication, dose-sensitive genes are lost through selective sweeps to ensure normal organismal function, and environmental stress-related genes tend to be amplified. The HSPa2 gene has undergone tandem duplication in several teleost fishes (Evgen'Ev et al., 2014; Metzger et al., 2016). Tandem duplication was responsible for the expansion of the oyster and scallop HSP70 gene families (Zhang et al., 2012; Cheng et al., 2016; Hu et al., 2019). Forty-five MmHSP70 genes were tandemly duplicated in the M. mercenaria genome, while twenty-six CsHSP70 genes were tandemly duplicated in the C. sinensis genome. Differences in the number of tandemly duplicated HSP70 genes between species may be explained by tandemly duplicated genes having independent origins or there having been some genomic rearrangement in this region (Metzger et al., 2016). In our study, the same pairs of tandemly duplicated HSP70 genes seemed to have highly similar conserved motifs (Figure 3). Moreover, all HSPa B2 tandem duplication pair genes (MmHSPa B2_1, MmHSPa B2_2 and MmHSPa B2_3) had the same number of exons (Figure 2). In the early phase of evolution, the function of duplicated genes is retained through purifying selection (Kondrashov et al., 2002). The Ka/Ks values of all tandemly duplicated HSP70 gene pairs were significantly less than 1 (ranging from 0.0548205 to 0.417659, P < 0.05, Table 1). Moreover, the tandemly duplicated HSP70 gene pairs seemed to have similar expression patterns under heat and hypoxia stress (Figure 8). These results indicated that duplicated MmHSP70 genes had experienced different levels of purifying selection and that purifying selection could be an important source of sequence and functional constraints (Song et al., 2021).

In the present study, MmHSP70 genes were expressed in a tissue-specific and stress-specific pattern. Specifically, most MmHSP70 genes were expressed in at least one tissue, while some MmHSP70s showed tissue-specific expression. Some HSP70 genes were highly expressed in gills in the Manila clam (Liu et al., 2015; Nie et al., 2017) and scallops (Cheng et al., 2019). In our present study, the expression pattern of HSP70 in the gill was different from that in the other tested tissues, and some MmHSP70 genes had high expression in the gill. Bivalve gills are sensitive to environmental stress, and high HSP70 expression promotes the regulation of the environmental stress response (Cheng et al., 2019). In marine ecosystems, bivalves always encounter stresses such as high temperature, low osmotic pressure, low dissolved oxygen, heavy metals, toxic dinoflagellates and bacterial invasion, and numerous studies have shown that HSP70 genes play a vital role in maintaining cellular homeostasis to defend against abiotic and biotic stresses (Piano et al., 2002; Franzellitti and Fabbri, 2005; Liu et al., 2014; Cheng et al., 2016; Nie et al., 2017; Nie et al., 2018; Cheng et al., 2019; Clark et al., 2021; Hu et al., 2022). In our previous study, many HSP70 genes were differentially expressed under heat, hypoxia and combined stress (Hu et al., 2022). Tandemly duplicated gene pairs played vital roles in stress adaptation in bivalves (Zhang et al., 2012; Li et al., 2021; Song et al., 2021). In the present study, we determined the expression profiles of tandemly duplicated HSP70 genes by analysing RNA-seq data. As shown in Figure 8, the tandemly duplicated HSP70 gene pairs seemed to have similar expression patterns. MmHsp70s exhibited diverse expression patterns in the gill when bivalves were exposed to heat, hypoxia and combined stress. This could be partly explained by HSP70 subfunctionalization after gene duplication (Ramsøe et al., 2020). Most tandemly duplicated HSP70 gene pairs had higher expression under hypoxic stress. Interestingly, the tandemly duplicated HSPa B2 gene pairs (MmHSPa B2_1, MmHSPa B2_2 and MmHSPa B2_3) showed significantly high expression under heat plus severe hypoxia stress. Above all, the considerable expansion of HSP70 genes may enhance transcriptional complexity and play an important role in adaptation to diverse temperatures and dissolved oxygen conditions in the M. mercenaria.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA596049, https://www.ncbi.nlm.nih.gov/, PRJNA764366, https://www.ncbi.nlm.nih.gov/, PRJNA764372.

Author Contributions

ZH contributed to Data curation, Software, Formal analysis, Writing - original draft and Writing - review and editing. HS contributed to Data curation, Formal analysis, Writing - review and editing. JF contributed to Formal analysis, Writing - review and editing. CZ contributed to Data curation, Formal analysis. M-JY contributed to Formal analysis, Software. PS contributed to Formal analysis, Software. Z-LY contributed to Visualization. Y-RL contributed to Resources. Y-JG contributed to Resources. H-ZL contributed to Resources. TZ contributed to Supervision, Funding acquisition, Writing - review and editing. All authors read and approved the final manuscript.

Funding

This research was supported by the China National Key R & D Program (2019YFD0900800), Tianjin science and technology commission project (20YFZCSN00240), Tianjin Agricultural Committee Project (202103010), the China Agriculture Research System of MOF and MARA, Shandong Province Agricultural Major Applied Technology Innovation Project (Grant No. SF1405303301), Taishan Industrial Leading Talents Project (Grant No. LJNY201704, Recipient: Tao Zhang), the ‘Double Hundred’ Blue Industry Leader Team of Yantai (Recipient: Tao Zhang), and the Creative Team Project of the Laboratory for Marine Ecology and Environmental Science, Qingdao Pilot National Laboratory for Marine Science and Technology (Grant No. LMEES-CTSP-2018-1). The funders had no role in the study design, data collection or analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest

Author H-ZL was employed by Shandong Fu Han Ocean Sci-Tech Co., Ltd.

The remaining 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.

Acknowledgments

The authors thanked Novogene Bioinformatics Institute for bioinformatics analysis. The authors also thanked the AJE (https://www.aje.cn/) for editing the language of the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.898669/full#supplementary-material

References

Bailey T. L., Johnson J., Grant C. E., Noble W. S.. (2015). The MEME Suite. Nucleic Acids Res. 43 (W1), W39–W49. doi: 10.1093/nar/gkv416

PubMed Abstract | CrossRef Full Text | Google Scholar

Brocchieri L., De Macario E. C., Macario A. J. (2008). Hsp70 Genes in the Human Genome: Conservation and Differentiation Patterns Predict a Wide Array of Overlapping and Specialized Functions. BMC Evol. Biol. 8 (1), 19–19. doi: 10.1186/1471-2148-8-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchfink B., Xie C., Huson D. H. (2015). Fast and Sensitive Protein Alignment Using DIAMOND. Nat. Methods 12 (1), 59–60. doi: 10.1038/nmeth.3176

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen C., Chen H., Zhang Y., Thomas H. R., Frank M. H., He Y, et al. (2020). TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Mol. Plant 13 (8), 1194–1202. doi: 10.1016/j.molp.2020.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng D., Liu H., Zhang H., Soon T. K., Ye T., Li S., et al. (2019). Differential Expressions of HSP70 Gene Between Golden and Brown Noble Scallops Chlamys Nobilis Under Heat Stress and Bacterial Challenge. Fish Shellfish Immunol., 94 , 924–933. doi: 10.1016/j.fsi.2019.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng J., Xun X., Kong Y., Wang S., Yang Z., Li Y., et al. (2016). Hsp70 Gene Expansions in the Scallop Patinopecten Yessoensis and Their Expression Regulation After Exposure to the Toxic Dinoflagellate Alexandrium Catenella. Fish Shellfish Immunol., 58, 266–273. doi: 10.1016/j.fsi.2016.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen S., Krinsky B. H., Long M. (2013). New Genes as Drivers of Phenotypic Evolution. Nat. Rev. Genet. 14 (9), 645–660. doi: 10.1038/nrg3521

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark M. S., Peck L. S., Thyrring J. (2021). Resilience in Greenland Intertidal Mytilus: The Hidden Stress Defense. Sci. Total Environ. 767, 144366. doi: 10.1016/j.scitotenv.2020.144366

PubMed Abstract | CrossRef Full Text | Google Scholar

Daugaard M., Rohde M., Jaattela M. (2007). The Heat Shock Protein 70 Family: Highly Homologous Proteins With Overlapping and Distinct Functions. FEBS Lett. 581 (19), 3702–3710. doi: 10.1016/j.febslet.2007.05.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Evgen'Ev M. B., Garbuz D. G., Zatsepina O. G. (2014). Heat Shock Proteins and Whole Body Adaptation to Extreme Environments (Switzerland: Springer).

Google Scholar

Franzellitti S., Fabbri E. (2005). Differential HSP70 Gene Expression in the Mediterranean Mussel Exposed to Various Stressors. Biochem. Biophys. Res. Commun. 336 (4), 1157–1163. doi: 10.1016/j.bbrc.2005.08.244

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu H., Jiao Z., Li Y., Tian J., Ren L., Zhang F., et al. (2021). Transient Receptor Potential (TRP) Channels in the Pacific Oyster (Crassostrea Gigas): Genome-Wide Identification and Expression Profiling After Heat Stress Between C. Gigas and C. Angulata. Int. J. Mol. Sci. 22 (6), 3222. doi: 10.3390/ijms22063222

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujikawa T., Munakata T., Kondo S. I., Satoh N., Wada S.. (2010). Stress Response in the Ascidian Ciona Intestinalis: Transcriptional Profiling of Genes for the Heat Shock Protein 70 Chaperone System Under Heat Stress and Endoplasmic Reticulum Stress. Cell Stress Chaperones 15 (2), 193–204. doi: 10.1007/s12192-009-0133-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gething M., Sambrook J. (1992). Protein Folding in the Cell. Nature 355 (6355), 33–45. doi: 10.1038/355033a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong W. J., Golic K. G. (2004). Genomic Deletions of the Drosophila Melanogaster Hsp70 Genes. Genetics 168 (3), 1467–1476. doi: 10.1534/genetics.104.030874

PubMed Abstract | CrossRef Full Text | Google Scholar

Guerin M. N., Weinstein D. J., Bracht J. R. (2019). Stress Adapted Mollusca and Nematoda Exhibit Convergently Expanded Hsp70 and AIG1 Gene Families. J. Mol. Evol. 87 (9), 289–297. doi: 10.1007/s00239-019-09900-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Heschl M. F., Baillie D. L. (1990). The Hsp70 Multigene Family of Caenorhabditis Elegans. Comp. Biochem. Physiol. B 96 (4), 633–637. doi: 10.1016/0305-0491(90)90206-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu Z., Feng J., Song H., Zhou C., Yu Z. L., Yang M. J., et al. (2022). Mechanisms of Heat and Hypoxia Defense in Hard Clam: Insights From Transcriptome Analysis. Aquaculture, 549, 737792. doi: 10.1016/j.aquaculture.2021.737792

CrossRef Full Text | Google Scholar

Hu B., Li M., Yu X., Xun X., Lu W., Li X., et al. (2019). Diverse Expression Regulation of Hsp70 Genes in Scallops After Exposure to Toxic Alexandrium Dinoflagellates. Chemosphere 234, 62–69. doi: 10.1016/j.chemosphere.2019.06.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Kondrashov F. A., Rogozin I. B., Wolf Y. I., Koonin E. V.. (2002). Selection in the Evolution of Gene Duplications. Genome Biol. 3 (2), 1–9. doi: 10.1186/gb-2002-3-2-research0008

CrossRef Full Text | Google Scholar

Kumar S., Stecher G., Tamura K. (2016). MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Mol. Biol. Evol. 33 (7), 1870–1874. doi: 10.1093/molbev/msw054

PubMed Abstract | CrossRef Full Text | Google Scholar

Li A., Dai H., Guo X., Zhang Z., Zhang K., Wang C., et al. (2021). Genome of the Estuarine Oyster Provides Insights Into Climate Impact and Adaptive Plasticity. Commun. Biol. 4 (1), 1–12. doi: 10.1038/s42003-021-02823-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindquist S., Craig E. A. (1988). The Heat-Shock Proteins. Annu. Rev. Genet. 22 (1), 631–677. doi: 10.1146/annurev.ge.22.120188.003215

PubMed Abstract | CrossRef Full Text | Google Scholar

Li Y., Sun X., Hu X., Xun X., Zhang J., Guo X., et al. (2017). Scallop Genome Reveals Molecular Adaptations to Semi-Sessile Life and Neurotoxins. Nat. Commun. 8 (1), 1–11. doi: 10.1038/s41467-017-01927-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Li Q., Sun S., Zhang F., Wang M., Li M., et al. (2019). Effects of Hypoxia on Survival, Behavior, Metabolism and Cellular Damage of Manila Clam (Ruditapes Philippinarum). PloS One 14 (4), e0215158. doi: 10.1371/journal.pone.0215158

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu H., He J., Chi C. F., Shao J. (2014). Differential HSP70 Expression in Mytilus Coruscus Under Various Stressors. Gene 543 (1), 166–173. doi: 10.1016/j.gene.2014.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu T., Pan L., Cai Y., Miao J. (2015). Molecular Cloning and Sequence Analysis of Heat Shock Proteins 70 (HSP70) and 90 (HSP90) and Their Expression Analysis When Exposed to Benzo (a) Pyrene in the Clam Ruditapes Philippinarum. Gene 555 (2), 108–118. doi: 10.1016/j.gene.2014.10.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Long M., VanKuren N. W., Chen S., Vibranovski M. D.. (2013). New Gene Evolution: Little did We Know. Annu. Rev. Genet. 47, 307–333. doi: 10.1146/annurev-genet-111212-133301

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma F., Luo L. (2020). Genome-Wide Identification of Hsp70/110 Genes in Rainbow Trout and Their Regulated Expression in Response to Heat Stress. PeerJ 8, e10022. doi: 10.7717/peerj.10022

PubMed Abstract | CrossRef Full Text | Google Scholar

Menzel R. W. (2009). Quahog Clams and Their Possible Mariculture. Proc. Annu. Workshop - World Mariculture Soc. 2 (1-4), 21–36. doi: 10.1111/j.1749-7345.1971.tb00029.x

CrossRef Full Text | Google Scholar

Metzger D. C., Hemmer-Hansen J., Schulte P. M. (2016). Conserved Structure and Expression of Hsp70 Paralogs in Teleost Fishes. Comp. Biochem. Physiol. Part D: Genomics Proteomics 18, 10–20. doi: 10.1016/j.cbd.2016.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy M. E. (2013). The HSP70 Family and Cancer. Carcinogenesis 34 (6), 1181–1188. doi: 10.1093/carcin/bgt111

PubMed Abstract | CrossRef Full Text | Google Scholar

Nie H., Liu L., Huo Z., Chen P., Ding J., Yang F., et al. (2017). The HSP70 Gene Expression Responses to Thermal and Salinity Stress in Wild and Cultivated Manila Clam Ruditapes Philippinarum. Aquaculture 470, 149–156. doi: 10.1016/j.aquaculture.2016.12.016

CrossRef Full Text | Google Scholar

Nie H., Liu L., Wang H., Huo Z., Yan X. (2018). Stress Levels Over Time in Ruditapes Philippinarum: The Effects of Hypoxia and Cold Stress on Hsp70 Gene Expression. Aquacul. Rep. 12, 1–4. doi: 10.1016/j.aqrep.2018.08.003

CrossRef Full Text | Google Scholar

Nikolaidis N., Nei M. (2004). Concerted and Nonconcerted Evolution of the Hsp70 Gene Superfamily in Two Sibling Species of Nematodes. Mol. Biol. Evol. 21 (3), 498–505. doi: 10.1093/molbev/msh041

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng J., Li Q., Xu L., Wei P., He P., Zhang X., et al. (2020). Chromosome-Level Analysis of Crassostrea Hongkongensis Genome Reveals Extensive Duplication of Immune-Related Genes in Bivalves. Mol. Ecol. Resources 00, 1–15. doi: 10.1111/1755-0998.13157

CrossRef Full Text | Google Scholar

Piano A., Asirelli C., Caselli F., Fabbri E. (2002). Hsp70 Expression in Thermally Stressed Ostrea Edulis, a Commercially Important Oyster in Europe. Cell Stress Chaperones 7 (3), 250. doi: 10.1379/1466-1268(2002)007<0250:HEITSO>2.0.CO;2

PubMed Abstract | CrossRef Full Text | Google Scholar

Powell D., Subramanian S., Suwansaard S., Zhao M., O’Connor W., Raftos D., et al. (2018). The Genome of the Oyster Saccostrea Offers Insight Into the Environmental Resilience of Bivalves. DNA Res. 25 (6), 655–665. doi: 10.1093/dnares/dsy032

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramsøe A., Clark M. S., Sleight V. A. (2020). Gene Network Analyses Support Subfunctionalization Hypothesis for Duplicated Hsp70 Genes in the Antarctic Clam. Cell Stress Chaperones 25, 1111–1116. doi: 10.1007/s12192-020-01118-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritossa F. (1962). A New Puffing Pattern Induced by Temperature Shock and DNP in Drosophila. Experientia 18 (12), 571–573. doi: 10.1007/BF02172188

CrossRef Full Text | Google Scholar

Rosenzweig R., Nillegoda N. B., Mayer M. P., Bukau B. (2019). The Hsp70 Chaperone Network. Nat. Rev. Mol. Cell Biol. 20 (11), 665–680. doi: 10.1038/s41580-019-0133-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth C., Rastogi S., Arvestad L., Dittmar K., Light S., Ekman D. (2007). Evolution After Gene Duplication: Models, Mechanisms, Sequences, Systems, and Organisms. J. Exp. Zool. Part B: Mol. Dev. Evol. 308 (1), 58–73. doi: 10.1002/jez.b.21124

CrossRef Full Text | Google Scholar

Song H., Guo X., Sun L., Wang Q., Han F., Wang H., et al. (2021). The Hard Clam Genome Reveals Massive Expansion and Diversification of Inhibitors of Apoptosis in Bivalvia. BMC Biol. 19 (1), 1–20. doi: 10.1186/s12915-020-00943-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Song L., Li C., Xie Y., Liu S., Zhang J., Yao J., et al. (2016). Genome-Wide Identification of Hsp70 Genes in Channel Catfish and Their Regulated Expression After Bacterial Infection. Fish Shellfish Immunol. 49, 154–162. doi: 10.1016/j.fsi.2015.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Srivastava P. K. (2002). Roles of Heat-Shock Proteins in Innate and Adaptive Immunity. Nat. Rev. Immunol. 2 (3), 185–194. doi: 10.1038/nri749

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun J., Zhang Y., Xu T., Zhang Y., Mu H., Zhang Y., et al. (2017). Adaptation to Deep-Sea Chemosynthetic Environments as Revealed by Mussel Genomes. Nat. Ecol. Evol. 1 (5), 1–7. doi: 10.1038/s41559-017-0121

PubMed Abstract | CrossRef Full Text | Google Scholar

Takeuchi T., Koyanagi R., Gyoja F., Kanda M., Hisata K., Fujie M., et al. (2016). Bivalve-Specific Gene Expansion in the Pearl Oyster Genome: Implications of Adaptation to a Sessile Lifestyle. Zool. Lett. 2 (1), 1–13. doi: 10.1186/s40851-016-0039-2

CrossRef Full Text | Google Scholar

Tripathy K., Sodhi M., Kataria R. S., Chopra M., Mukesh M.. (2021). In Silico Analysis of HSP70 Gene Family in Bovine Genome. Biochem. Genet. 59 (1), 134–158. doi: 10.1007/s10528-020-09994-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang S., Hou R., Bao Z., Du H., He Y., Su H., et al. (2013). Transcriptome Sequencing of Zhikong Scallop (Chlamys Farreri) and Comparative Transcriptomic Analysis With Yesso Scallop (Patinopecten Yessoensis). PloS One 8 (5), e63927. doi: 10.1371/journal.pone.0063927

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang Y., Tang H., DeBarry J. D., Tan X., Li J., Wang X., et al. (2012). MCScanX: A Toolkit for Detection and Evolutionary Analysis of Gene Synteny and Collinearity. Nucleic Acids Res. 40 (7), e49. doi: 10.1093/nar/gkr1293

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang D., Zhang Y., Zhang Z., Zhu J., Yu J., et al. (2010). KaKs_Calculator 2.0: A Toolkit Incorporating Gamma-Series Methods and Sliding Window Strategies. Genomics Proteomics Bioinf. 8 (1), 77–80. doi: 10.1016/S1672-0229(10)60008-3

CrossRef Full Text | Google Scholar

Wei M., Ge H., Shao C., Yan X., Nie H., Duan H., et al. (2020). Chromosome-Level Clam Genome Helps Elucidate the Molecular Basis of Adaptation to a Buried Lifestyle. iScience 23, 101148. doi: 10.1016/j.isci.2020.101148

PubMed Abstract | CrossRef Full Text | Google Scholar

Welch W. J. (1993). How Cells Respond to Stress. Sci. Am. 268 (5), 56–56. doi: 10.1038/scientificamerican0593-56

CrossRef Full Text | Google Scholar

Xu K., Xu H., Han Z. (2018). Genome-Wide Identification of Hsp70 Genes in the Large Yellow Croaker (Larimichthys Crocea) and Their Regulated Expression Under Cold and Heat Stress. Genes 9 (12), 590. doi: 10.3390/genes9120590

CrossRef Full Text | Google Scholar

Yan X., Nie H., Huo Z., Ding J., Li Z., Yan L., et al. (2019). Clam Genome Sequence Clarifies the Molecular Basis of Its Benthic Adaptation and Extraordinary Shell Color Diversity. iScience, 19,1225–1237. doi: 10.1016/j.isci.2019.08.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang J. (2003). Evolution by Gene Duplication: An Update. Trends Ecol. Evol. 18 (6), 292–298. doi: 10.1016/S0169-5347(03)00033-8

CrossRef Full Text | Google Scholar

Zhang G., Fang X., Guo X., Li L., Luo R., Xu F., et al. (2012). The Oyster Genome Reveals Stress Adaptation and Complexity of Shell Formation. Nature 490 (7418), 49–54. doi: 10.1038/nature11413

PubMed Abstract | CrossRef Full Text | Google Scholar