Pierre-Louis Stenger1,2
Chin-Long Ky1,2
Céline M. O. Reisser1,3Céline Cosseau4
Christoph Grunau4Mickaël Mege1,5Serge Planes6
Jeremie Vidal-Dupiol2*- 1IFREMER, UMR 241 Écosystèmes Insulaires Océaniens, Labex Corail, Centre du Pacifique, Tahiti, French Polynesia
- 2IHPE, Université de Montpellier, CNRS, IFREMER, Université de Perpignan Via Domitia, Montpellier, France
- 3MARBEC, Université de Montpellier, CNRS, IFREMER, IRD, Montpellier, France
- 4IHPE, Université de Montpellier, CNRS, IFREMER, Université de Perpignan Via Domitia, Perpignan, France
- 5IFREMER, PDG-RBE-SGMM-LGPMM, La Tremblade, France
- 6EPHE-UPVD-CNRS, USR 3278 CRIOBE, Labex Corail, PSL Research University, Université de Perpignan, Perpignan, France
Today, it is common knowledge that environmental factors can change the color of many animals. Studies have shown that the molecular mechanisms underlying such modifications could involve epigenetic factors. Since 2013, the pearl oyster Pinctada margaritifera var. cumingii has become a biological model for questions on color expression and variation in Mollusca. A previous study reported color plasticity in response to water depth variation, specifically a general darkening of the nacre color at greater depth. However, the molecular mechanisms behind this plasticity are still unknown. In this paper, we investigate the possible implication of epigenetic factors controlling shell color variation through a depth variation experiment associated with a DNA methylation study performed at the whole genome level with a constant genetic background. Our results revealed six genes presenting differentially methylated CpGs in response to the environmental change, among which four are linked to pigmentation processes or regulations (GART, ABCC1, MAPKAP1, and GRL101), especially those leading to darker phenotypes. Interestingly, the genes perlucin and MGAT1, both involved in the biomineralization process (deposition of aragonite and calcite crystals), also showed differential methylation, suggesting that a possible difference in the physical/spatial organization of the crystals could cause darkening (iridescence or transparency modification of the biomineral). These findings are of great interest for the pearl production industry, since wholly black pearls and their opposite, the palest pearls, command a higher value on several markets. They also open the route of epigenetic improvement as a new means for pearl production improvement.
Introduction
Pearls have captivated and amazed human civilization since 7500 BP (Charpentier et al., 2012) and have formed the basis for economic development of several tropical countries. French Polynesia started to trade pearls with Europeans at the end of the 18th century (Southgate and Lucas, 2011) and has been an industrial producer since 1964 (Le Pennec, 2010). Today, cultured pearl farming represents the second economic resource of the country, just behind tourism (Ky et al., 2019). However, since 2001, an unprecedented economic crisis has endangered the French Polynesian pearl farming sector due to the overproduction of low quality pearls (Bouzerand, 2018). To remediate to this socio-economic crisis, the policy of producing “less but better” was adopted in 2013 (Ky et al., 2013). Since then, production has focused on obtaining cultured pearls with high market value, by selecting traits such as pearl color.
Two individuals are needed to produce a pearl, a donor and a recipient oyster. The sacrificed donor is used to provide a piece of mantle (the biomineralizing graft) which is placed, together with a small marble of nacre, into the gonad of the recipient oyster. Donor oysters are selected for their inner-shell color, since the color of a cultured pearl is determined by that of the donor oyster (Ky et al., 2013). Recipient oysters are selected for their vigor. The Polynesian black-lipped pearl oyster, Pinctada margaritifera var. cumingii (Linnaeus, 1758) is the pearl oyster species showing the largest range of inner shell color (Ky et al., 2013, 2017), therefore offering a wide range of pearl colors and shades, like dark, pastel, silver, peacock, red, golden, green, blue, and even rainbow (Ky et al., 2014; Stenger et al., 2019). Moreover, the color of this bivalve has both a structural (iridescence; Liu et al., 1999) and a biological (pigments; Iwahashi and Akamatsu, 1994) basis. While the genetic determination of this color has been partly demonstrated (Ky et al., 2013, 2017), several environmental factors are also known to affect the shell color (Joubert et al., 2014; Le Pabic et al., 2016; Le Moullac et al., 2018), such as the depth at which an oyster is grown (Stenger et al., 2019). In this latter work, authors demonstrated that the transplantation of oysters from the sub-surface (−4 m) to the bottom of the lagoon (−30 m) significantly darkens the inner shell color compared with oysters maintained at the sub-surface. This induced phenotype was persistent through time (“enduring”), even after the oysters were returned to shallow water, which would seem to indicate an epigenetic control mechanism rather than a direct environmental influence acting on the darkening of the shell color (Stenger et al., 2019).
Since the first definition of epigenetic by Waddington in the late 1930s, epigenetic received many definitions (Nicoglou and Merlin, 2017). In this work we have selected the definition proposed by Russo et al. (1996), e.g., “epigenetic is the study of mitotically and/or meiotically heritable changes in gene function that cannot be explained by changes in DNA sequence” (Russo et al., 1996). Mechanistically, this memory function is based on changes in chromatin structure, such as non-coding RNA and/or histone modifications and/or DNA methylation. Here, we use the term epigenetics to describe any changes in DNA methylation that occur upon environmental cues. Coloration mediates an organism’s relationship with their environment in important ways including anti-predator defenses, social signaling, thermoregulation, or protection (Cuthill et al., 2017). Several species are known to change their coloration more than once in their lifetime in response to environmental triggers to reach an optimal phenotype in a new environment. To increase their camouflage the arctic hare Lepus arcticus (Ross, 1819), the ermine Mustela erminea (Linnaeus, 1758), and the ptarmigan Lagopus muta (Montin, 1776) changed their coat color from brown or gray in the summer to white in the winter (Zimova et al., 2018). These changes are known to be induced by temperature, photoperiod, and/or food rarefaction (Zimova et al., 2018). The corresponding changes in color expression could be brought by epigenetic processes (Hu and Barrett, 2017) like in mouses (Dolinoy, 2008). Indeed, the most striking example of color change in mammals rely on the environmentally induced differential DNA methylation of the intracisternal A-particle gene located upstream of the agouti locus (Dolinoy et al., 2006; Dolinoy, 2008). In Mollusca, a few epigenetic studies have been made, especially on color expression and variation. To date, Feng et al. (2018) provided a catalog of long non-coding RNA (lncRNA) expressed in the mantle of the Pacific oyster Crassostrea gigas (Thunberg, 1793). These authors suggested that these lncRNAs may affect the expression of pigment-related genes such as tyrosinase-like proteins, dopamine, beta-monooxygenase, chorion peroxidase, or cytochrome P450 2U1, thus leading to different shell color phenotypes. More recently, the same group (Feng et al., 2020) have studied the role of microRNAs in the regulation of the shell color of Crassostrea gigas. In this study, four miRNAs (lgi-miR-315, lgi-miR-96b, lgi-miR-317, and lgi-miR-153) were found closely associated with shell color along with the regulation of Cytochrome P450 2U1, Tyrosinase-like protein 2 and 3. The authors concluded that lgi-miR-317, its targeted mRNA encoding peroxidase, and the lncRNA TCONS_00951105 might play a key role in shell melanin synthesis.
Because color phenotype is important for the pearl market and the phenotypic plasticity of this trait is associated with the putative involvement of epigenetic mechanisms, the study of these mechanisms opens a new avenue for improvement in the pearl industry with, for example, the development of epi-markers for environmentally induced color variation testing. As a first step in exploring a possible interaction between epigenetic mechanisms and color variation, we designed a depth variation experiment to induce environmentally driven color variation. In order to disentangle the genetic factors from the epigenetic ones influencing color variation, a non-lethal sampling design was used enabling us to monitor changes in DNA methylation over time and depth within the same individuals (constant genotypes). DNA methylation was studied at the whole genome scale by whole-genome bisulfite sequencing and provided evidence for an epigenetic control of pearl oyster color variation. This approach enabled us to find any differences in DNA methylation in pearl oysters after a period at increased depth and, when this occurred, to examine whether genes related to pigmentation and/or biomineralization processes were affected by such changes. Results of this kind could allow the pearl industry to turn to more sustainable production strategies.
Materials and Methods
Biological Material and the Yo-Yo Experiment
In order to trigger an environmentally driven color change of the inner shell of P. margaritifera, we set up an in situ “yo-yo” experiment (May to August 2017) (Figure 1). Six individuals of 4 years of age (approximately 14 cm height) originating from three different families (two individuals per family Stenger et al., in press) were used. These six individuals were first maintained at 8 m depth for 1 month (May 2017). Then, three of them (1 from each family) were selected and transferred to 30 m depth (treatment) for environmental pressure while the three others were left at 8 m (control). This exposure was maintained for 1 month (June 2017). Then, the three pearl oysters that had been placed at 30 m were transferred back to 8 m depth for a final month of exposure (July 2017). During each transfer, a piece of mantle (the biomineralizing tissue responsible for the inner shell coloration) was sampled by a non-lethal method: (i) oysters were anesthetized in 20 L seawater containing 200 mL benzocaine at 120 g/L 96° ethanol under air aeration; (ii) they opened their valves under the effect of the benzocaine, a 5 mm3 fragment of the mantle was carefully sampled with tweezers and scissors; (iii) the sample was flash frozen in liquid nitrogen. Alongside the sampling of the mantle, the color of the inner shell of each individual (control and treatment) was filmed with a mini photo studio for color variation analysis. This mini standardized photo studio was composed of a tripod supporting a Nikon D3100 reflex camera equipped with a Nikon 18-55VR lens. This set up was used to film the reflection of the inner shell color on a small mirror (spatula). To assess constant exposition to light, all films were made under a blackout drape with three white LED lamps.
Figure 1. Design of the Yo-Yo experiment. Six P. margaritifera individuals were used for this experiment. Three individuals were maintained at a depth of 8 m during the whole experiment and were used as a control (C – in blue circles). The three other individuals were subjected to depth variation treatment (T – in green circles). A non-lethal sampling of mantle tissue was made every month at the time of the depth change.
Color Variation Analysis
For each sample, ten screenshots from the films were randomly captured and analyzed for color variation. Color was quantified and qualified using the R package ImaginR V.2 (Stenger, 2017) as previously described (Stenger et al., 2019). A Shapiro test (stats v3.5.0 R package) was used to assess the normal distribution of the data. The average saturation and darkness for each sample were calculated from all ten screenshots, and a Wilcoxon test (stats v3.5.0 R package based on Hollander et al., 1973 and Royston, 1995) was used to test for any difference between the groups.
DNA Extraction and BS-seq
DNA was extracted with the QIAamp DNA Mini Kit (QIAGEN®, Cat No./ID: 51306) following manufacturer’s recommendations, with the addition of two steps: (i) a RNase A treatment to remove RNA and (ii) after overnight digestion the addition of 50 μL of a saturated KCl solution (34 g KCl/100 g H2O) and a centrifugation at 14,000g for 15 min to remove the mucopolysaccharides (Sokolov, 2000). DNA quality and quantity were assessed with a NanoDropTM 2000 and 1.2% agarose gel electrophoresis. Bisulfite-conversion, library construction, and sequencing (2 × 150 bp) were performed by Genome Québec (MPS Canada) on an Illumina HiSeq X.
Bioinformatics Pipeline
Analyses were performed on the Ifremer Datarmor cluster1. Raw read quality was assessed with FastQC (Andrews, 2010; Supplementary Table 1). Reads were cleaned and adaptors removed with Trimmomatic (Bolger et al., 2014) (V. 0.36 – illuminaclip 2:30:10; leading 28, trailing 28, and minlen 40). Bismark aligner (Krueger and Andrews, 2011) (V. 0.19) was used to map reads on the draft genome of P. margaritifera (Reisser et al., 2020) using the following parameters: multi-seed length of 30 bp, 0 mismatches, default minimum alignment score function. Deduplication was done with Deduplicate_bismark (Krueger and Andrews, 2011). Bed files for methylome characterization were obtained with Bismark_methylation_extractor (Krueger and Andrews, 2011). All scripts are provided on GitHub (PLStenger/Pearl_Oyster_Colour_BS_Seq/00_scripts). Raw reads are available through the NCBI Sequence Read Archive (SRA, BioProject PRJNA663978, BioSample SAMN16191417 to SAMN16191446).
Methylation Calling, Methylome Characterization, and Differential Methylation Analysis
The R package Methylkit (Akalin et al., 2012) (V. 1.11.0) was used for methylation calling in CpG, CHH, and CHG contexts using a minimum coverage of 10, directly from BAM files with processBismarkAln. Methylkit was also used for methylation characterization in the CpG, CHH, and CHG contexts, as well as for the coverage calculation and clustering analysis (ward clustering correlation distance method). The average gene methylation was calculated with DeepTools V. 3.3.0 (Ramírez et al., 2014). The gene body methylation rate (GBMR), corresponding to the CpG methylation rate, was calculated with the map function from bedtools V. 2.27.1 (Quinlan and Hall, 2010).
The mantle’ gene expression data used for methylation/gene expression correlation analysis came from individuals studied in Stenger et al., in press (SRA BioProject PRJNA521849). Briefly, RNA sequencing was done using high quality RNA extracted from twelve individuals. Library construction and sequencing was performed (Paired–end 100-bp; Illumina® HiSeq® 4000) by Génome Québec (MPS Canada). Read quality check and trimming were done as described in the section “Bioinformatics Pipeline.” Cleaned reads were paired-mapped against the P. margaritifera draft genome with TopHat (V1.4.0) (Trapnell et al., 2012). Cufflinks (V2.2.1.0) and Cuffmerge (V2.2.1.0) were used to assemble and merge the transcriptome produced for each library, respectively (Trapnell et al., 2012). HTSeq-count (V0.6.1) (Anders et al., 2015) was used to count read-mapped per transcript. The average of the twelve count files was computed and the RPKM method was used for gene expression data normalization as previously done in Wang et al., 2014.
To identify the effect of the depth treatment, differentially methylated CpGs (DMCpGs) were identified with the getMethylDiff function of the R package Methylkit, with difference >25% and q value < 0.05 for CpG positions between depth treatment and control oysters at each sampling time.
Functional Analysis of Differentially Methylated Genes
An annotation file was obtained following the first three steps of https://github.com/enormandeau/go_enrichment, completed with PLASTX (Nguyen and Lavenier, 2009) against UniProt-Swiss-Prot and TrEMBL (e-value at 1 × 10–3). A protein domain search was then performed with InterProScan. Finally, Gene Ontology terms were assigned with Blast2GO by combining information from the two annotation files (Conesa et al., 2005).
GOATOOLS (Klopfenstein et al., 2018) was used to test for enrichment of GO terms in a selected set of genes (significantly differentially methylated, lowly and highly methylated), using a Fisher’s exact test. Histograms showing the GO terms enrichment of were generated with the ggplot2 R package (Wickham and Chang, 2019) (V. 2.2.1).
Results
Depth Variation Induces Significant Darkening
The color variation analysis with the ImaginR package detected no significant modification of hue, saturation or darkness among the controls during the 3 months of the “yo-yo” experiment (pairwise Wilcoxon tests with P values > 0.05). Among the samples subjected to the treatment, no significant change occurred for hue or saturation during the 3 months of the experiment. However, a significant difference in darkness (pairwise Wilcoxon tests with P values < 1.10e–5) was detected, with an increase in the darkness value during the month at 30 m depth. This increase was still visible and significant after the oysters were returned to 8 m depth for 1 month (pairwise Wilcoxon tests with P values < 1.10e–5; Table 1). Raw data are available for download using the following links (screenshots: https://figshare.com/articles/dataset/ImaginR_raw_data_screenshots_zip/14049983; videos: https://figshare.com/articles/dataset/ImaginR_raw_data_videos_zip/14049986).
Table 1. Darkness values obtained with the ImaginR R package analysis.
Read Processing and Read Mapping to the Reference Genome
For the control, Illumina sequencing produced averages of 105,505,318 (±1,514,514; n = 3), 103,770,226 (±5,087,254; n = 3), and 103,490,158 (±5,509,306; n = 3) raw sequence reads for sampling times 1, 2, and 3, respectively. For the depth treatment, averages of 98,184,541 (±6,186,502; n = 3), 104,727,236 (±6,141,350; n = 3), and 113,132,278 (±4,306,828; n = 3) raw sequence reads were produced for sampling times 1, 2, and 3, respectively. After cleaning and filtering, averages of 101,575,689 (±1,584,202), 100,159,677 (±5,076,306), and 100,073,964 (±5,415,409) reads were kept for the control, and 94,669,211 (±5,877,099), 101,048,705 (±6,094,775), and 109,246,931 (±4,195,331) for the depth treatment, for the three successive sampling times, respectively. Filtered reads were mapped on the reference genome with Bismark and showed similar mapping rates for all samples (∼32.1%) (Supplementary Table 2). PCR duplicates were removed and represented, on average, 0.059% (±0.002) of the total reads.
Characterization of the P. margaritifera Mantle Methylome
Methylation in the P. margaritifera mantle displayed a mosaic pattern with an enrichment in the gene showing a CpG methylation rate of 22.81% in introns, 17.32% in exons, and an average of 18.26% in gene body (Figures 2A,B). A slight enrichment in the 3 kbps upstream (12.67%) and downstream (16.58%) of genes (Figures 2A,B) is also found. The whole genome CpG methylation rate is equal to 11.53% on average.
Figure 2. (A) Pattern of methylation level in the mantle of P. margaritifera (TSS, Transcription Start Site; TES, Transcription End Site) reported on a metagene. (B) Distribution of methylated CpG in different genomic regions. (C) Distribution of methylated CpG across genes ranked according to methylation level. The left part (1) is composed of lowly methylated genes while the right part (2) is composed of highly methylated genes. (D) CpG methylation rate averaged across rank of gene expression.
The fractionation of gene body methylation rate (Figure 2C) shows a bimodal distribution characterized by two peaks. This distribution enables the classification of genes into two sets: the “lowly” methylated genes (0.00–0.19% methylation – rank 0–0.2 and marked as “1” Figure 2C) and the “highly” methylated genes (0.50–40.12% of methylation; 40.12% is the maximum methylation found in a gene in these data – rank 0.5–1 and marked as “2” Figure 2C). The enrichment analysis performed on the lowly methylated genes showed an overrepresentation of GO terms involved in reproductive process functions, such as “oocyte maturation” (GO:0001556), “primary sex determination” (GO:0007539), “male meiosis I” (GO:0007141), etc., cellular signaling, such as “negative regulation or cellular response to drug” (GO:0048523), “negative regulation of phospholipid metabolic process” (GO:0071072), etc., and “androgen receptor signaling pathway” (GO:0030521). In the highly methylated genes set, the enrichment analysis showed an overrepresentation of GO terms involved in housekeeping functions, such as “mRNA regulation” (GO:0043488), “mRNA splice site regulation” (GO:0006376), “ribosomal small unit assembly” (GO:0000028), “transcription-dependent tethering of RNA polymerase II” (GO:0000972), “DNA amplification” (GO:0006277), etc (Supplementary Figure 1).
To test for a correlation between the GBMR and the level of gene expression, the distribution of gene body methylation levels was represented according to gene expression rank (in RPKM). This distribution revealed that moderately expressed genes (100–1000 RPKM) have higher methylation levels than lowly (>100 RPKM) or highly (<1000) expressed ones in P. margaritifera (Figure 2D).
Methylation Calling and Differential Methylation Analysis
Methylation calling showed that the Cytosine methylation level was in average of 16.5% [12.0% in the CpG context, 0.9% in the CHG context, 1.0% in the CHH context, and 2.6% in another context (CN or CHN); Supplementary Figure 2]. Similarities between the methylation patterns of each sample were analyzed by a clustering approach (ward clustering correlation distance method). This showed that samples clustered first by treatment (depth vs. control), then by genotype (i.e., individuals) and, for the control only, by sampling time (Figure 3).
Figure 3. CpG methylation clustering with the distance method “correlation” and clustering method “ward” from the clusterSamples function of the R package Methylkit. C, control individual (in blue); T, depth treatment individual (in green); #, label of the individual number; S1, sampling time 1; S2, sampling time 2; S3, sampling time 3.
For the differential methylation analysis, triplicates were made according to treatment (depth and control), and nine pairwise comparisons were made in order to identify and disentangle the different effects (Figure 4). Time effect was quantified by the differential methylation analysis performed among the control individuals (C; three different genotypes), but between the sampling times (S; three samples per genotype). This led us to perform three comparisons (C-S1 vs. C-S2, C-S2 vs. C-S3, and C-S1 vs. C-S3). Overall, time effect was associated with non-redundant (not visible in other conditions) hyper- or hypomethylations of 60 and 111 CpGs, respectively. Neither hyper- nor hypomethylation were detected in CHG or CHH contexts.
Figure 4. Number of significant hypermethylated (in red) and hypomethylated (in blue) positions (q < 0.05) in the nine pairwise comparisons. Blue boxes are the control groups, green boxes are the treatment groups. T, depth treatment individuals (in green); C, control individuals (in blue); S, sampling time.
The cumulative effect of depth variation and time was quantified by differential methylation analysis performed among the individuals in the depth treatment (T; three different genotypes), but between the sampling times (S; three samples per genotype). As previously, this led us to make three pairwise comparisons (T-S1 vs. T-S2, T-S2 vs. T-S3, and T-S1 vs. T-S3). In total, 71 non-redundant hyper- and 76 non-redundant hypomethylations were identified, all in the CpG context.
The cumulative effect of depth variation and genotype was quantified by differential methylation analysis performed between treatments (depth treatment vs. control; three genotypes per condition) at the three-sampling times (T-S1 vs. C-S1, T-S2 vs. C-S2, and T-S3 vs. C-S3). These analyses revealed that: (i) in the CpG context, 2087 non-redundant cytosines were hypermethylated while 2935 non-redundant cytosines were hypomethylated; (ii) in the CHG context, 12 non-redundant cytosines were hypermethylated and 16 were hypomethylated; (iii) in the CHH context, three non-redundant cytosines were hypermethylated and three were hypomethylated.
Finally, to extract the effect of the depth variation only, we subtracted the time effect from the cumulative effect of depth and time (comparison among treatments, but between the sampling times). To do so, positions presenting differential methylation in response to the time effect were considered as non-significant when they were also differentially methylated in the cumulative effect of depth and time. None of the CpGs that were differentially methylated for the time effect were also differentially methylated for the depth variation and time effect. All the 71 hyper- and 76 hypomethylations previously identified were therefore kept for subsequent analysis.
Enrichment Analysis and Exploration of Genes With Differentially Methylated Positions
To correlate changes in DNA methylation with changes in pigmentation, we first performed an enrichment analysis and looked for biological processes and molecular functions linked to pigmentation. As a second step, we then individually screened the genes displaying DMCpGs and searched the literature for their putative involvement in pigmentation. DMCpG that were located outside of the gene body were not considered (92 DMCpGs). Genes encoding proteins of unknown function were present in the set of genes containing DMCpGs, but the lack of functional annotation prevented us from proposing a mechanism that could explain their involvement in the phenotypic changes that occurred, although we cannot exclude that they may have a role in this phenomenon.
C-S1 vs. C-S2 vs. C-S3 Comparison
The GO categories significantly enriched in control conditions in relation to time were not associated with pigmentation. They were, however, associated with a seasonal effect, shown by an enrichment in GO terms linked to growth and reproduction (Supplementary Figure 2). The methylation information for all these genes and for all comparisons are provided in Supplementary Table 3.
T-S1 (1 Month at 8 m Depth) vs. T-S2 (1 Month at 8 m Depth Followed by 1 Month at 30 m Depth)
Gene ontology (GO) terms enrichment analysis performed for biological process category showed an enrichment of several GO terms associated with pigmentation (Figure 5): “pigment metabolic process” (GO:0042440), “pteridine-containing compound metabolic process” (GO:0042558, Frost and Malacinski, 1979), and “folic acid-containing compound biosynthetic process” (GO:0009396, Katz et al., 1987). For the molecular function category, additional GO terms linked to pigmentation processes were enriched: “UDP-glycosyltransferase activity” (GO:0008194, Vajro et al., 1995), “hydroxymethyl-formyl- and related transferase activity” (GO:0016742, Moreau et al., 2012), “alpha-1,3-mannosyl-glycoprotein 2-beta-N-acetylglucosaminyltransferase activity” (GO:0003827, Zhai et al., 2016), “phosphoribosylglycinamide formyltransferase activity” (GO:0004644), and “acetylglucosaminyltransferase activity” (GO:0016262, Chakraborty et al., 1999; Shin et al., 2015).
Figure 5. Histograms of the enriched biological processes (BP) and molecular functions (MF) from the genes with differentially methylated positions in individuals from the depth treatment samplings in T-S1 vs. T-S2 (T-S1, T-S2 = depth treatment at sampling times 1 and 2, respectively).
At the gene level (Table 2), three candidate genes were identified that had hypermethylated CpGs after the month at −30 m. The first encodes a trifunctional purine biosynthetic protein adenosine-3 (GART), the second an alpha-1,3-mannosyl-glycoprotein-2-beta-n-acetylglucosaminyltransferase (MGAT1), and the third a multidrug resistance-associated protein (ABCC1). GART is known to be involved in the synthesis of purines (purine synthesis pathway) and expression disturbances of this gene can modify the production of two pigments: melanin and pheomelanin (Amsterdam et al., 2004; Ng et al., 2009). MGAT1 is involved in the glycan synthesis pathway, and is known to be essential for shell formation in the Pinctada genus (Takakura et al., 2008). Finally, the ABBC1 gene encodes an active transporter of glutathione-S-transferase (Homolya et al., 2003; Fernandes and Gattass, 2009; Rocha et al., 2014), an enzyme regulating the balance of eumelanin/pheomelanin production (Sonthalia et al., 2016).
Table 2. All significantly differentially methylated positions (q value < 0.05) for treatment comparisons.
T-S2 (1 Month at 8 m Depth Followed by 1 Month at 30 m Depth) vs. T-S3 (1 Month at 8 m, 1 Month at 30 m, and 1 Month Back at 8 m)
GO term enrichment analysis highlighted enrichment for several GO terms associated with pigmentation (Figure 6): “L-ornithine transmembrane transporter” (GO:0000064) for the biological process category; and “ornithine transport” (GO:0015822), “lysine transport” (GO:0015819), “L-amino acid transport” (GO:0015807), “basic amino acid transmembrane transporter activity” (GO:0015171), “L-amino acid transmembrane transporter activity” (GO:0015179), and “xenobiotic transporter activity” (GO:0015238) for the molecular function category.
Figure 6. Histograms of the enriched biological processes (BP) and molecular functions (MF) from the genes with differentially methylated positions in individuals from the depth treatment samplings in T-S2 vs. T-S3 (T-S2, T-S3 = depth treatment at sampling times 2 and 3, respectively).
At the gene level, 15 genes displaying DMCpGs were deemed good candidates to partly explain the phenotypic changes observed. Among the hypomethylated genes, we found those for a cationic amino acid transporter 2 (SLC7A2), a multidrug-associated protein (also present in the T-S1 vs. T-S2 comparison), and a coiled-coil domain-containing protein 79 (CCDC79). SLC7A2 is involved in the transport of arginine, lysine, and ornithine. A high number of the proteins found in the mineralized structure of P. margaritifera are known to be enriched in arginine and lysine, such as MSP-1, MSP-2, Aspein, Prismlain-14, and MSI60 (Addadi et al., 2006; Joubert et al., 2010; Gueguen et al., 2013). According to the NCBI GenPept database, CCDC79 can bind ion calcium, like the product of the MGAT1 gene (see above). Among the hypermethylated genes, we found those for a G-protein coupled receptor (GRL101), a target of rapamycin complex 2 subunit (MAPKAP1), an enolase (enolase-4), sodium/potassium/calcium exchanger 1 and a perlucin-like protein. MAPKAP1 (TOR signaling pathway) promotes dark epithelial pigmentation (Liu et al., 2017), GRL101 (rhodopsin signaling pathway) is an ortholog of the pigment dispersing factor (Tanaka et al., 2014), and enolase (glycolysis/gluconeogenesis pathway) is a biomarker of vitiligo, a human pigmentation disorder affecting melanocytes (Hamid et al., 2015). Perlucin is a well-known matrix protein found in the nacreous layer of the pearl oyster shell (Joubert et al., 2010) with a function in the biomineralization process.
We also identified several GO terms with a less important role in pigmentation, like “TOR signaling” (GO:0031929), “TORC2 signaling” (GO:0038203), “glycolytic process” (GO:0006096), and “compound eye development” (GO:0048749) for the biological process GO category; and “calcium, potassium: sodium anti-porter activity” (GO:0005432), “alkali metal ion binding” (GO:0031420), and “mannose binding” (GO:0005537) for the molecular function GO category.
T-S1 (1 Month at 8 m Depth) vs. T-S3 (1 Month at 8 m, 1 Month at 30 m, and 1 Month Back at 8 m)
As the previous enrichment analyses revealed, significant enrichment of several GO categories correlated with pigmentation were identified like in the T-S2 vs. T-S3 pairwise comparison (Figure 7): “L-ornithine transmembrane transporter” (GO:0000064) for the biological process GO category; and “ornithine transport” (GO:0015822), “lysine transport” (GO:0015819), “L-amino acid transport” (GO:0015807), “basic amino acid transmembrane transporter activity” (GO:0015171), “L-amino acid transmembrane transporter activity” (GO:0015179), and “xenobiotic transporter activity” (GO:0015238) for molecular function GO category.
Figure 7. Histograms of the enriched biological processes (BP) and molecular functions (MF) from the genes with differentially methylated positions in individuals from the depth treatment samplings in T-S1 vs. T-S3 (T-S1, T-S3 = depth treatment at sampling times 1 and 3, respectively).
At the gene level, only the GPI-anchor transamidase gene could be linked to the pigmentation process (among other pathways) since it is involved in a human disorder characterized by altered dermal pigmentation (Ng and Freeze, 2014).
We also identified several GO terms with a less important role in pigmentation, such as the “drug transmembrane transport” (GO:0006855) for the biological process GO category; and “attachment of GPI anchor to protein” (GO:0016255), “regulation of TOR signaling” (GO:0032006), and “protein glycosylation” (GO:0006486) for the molecular function GO category.
Discussion
Pearl farming, the second economic resource of French Polynesia, has been suffering a major economic crisis since 2001. To address this problem, stakeholders, together with pearl farmers and scientists, are developing an ambitious plan to reduce the volume of pearl production, but to increase pearl quality. This objective could be reached through one of the most economically interesting traits of Pinctada margaritifera: its ability to express the largest range of inner shell color (and thus pearl color) of any pearl-producing species worldwide (Ky et al., 2014; Stenger et al., 2019). Indeed, producing unique, highly valuable, pearls displaying a palette of phenotypes ranging from very dark to very pale colors constitutes an efficient way to diversify the production and appeal to different markets. Selection of donor oysters based on color phenotype was started a few years ago (Ky et al., 2013). However, the possibility of acting directly on the selected oysters to enhance their color and quality as donors is also an attractive method for improving pearl quality. For this, a better understanding of the interactions between the phenotype and gene expression correlated with environmental conditions is essential (Gavery and Roberts, 2017), and would allow the development of epigenetic marker-assisted selection or, more generally, epigenetics-assisted cultural practices.
As a first step to understanding the mechanisms behind this environmentally induced color variation, we applied an environmental forcing (culture-depth variation) known to affect pearl and inner-shell darkness (Stenger et al., 2019), and studied, under constant genotype, the DNA methylation changes induced. In addition to providing the first description of a pearl oyster methylome, our analyses identified specific methylation changes that affected candidate genes involved in the expression of shell and pearl color darkness. These genes were involved in both pigmentation (biological coloring mechanisms) and iridescence (physical coloring mechanisms based on the differential organization of biomineral crystals).
The First Methylome of the Pterioidea Super-Family Shows Similar Characteristics to Other Invertebrate Methylomes
The present study provides to our knowledge the first methylome of a member of the Pterioidea super-family. The whole-genome bisulfite sequencing (WGBS-seq) of the pearl oyster mantle tissues revealed that its methylome is of the mosaic type and similar in many ways to what is classically described in some other invertebrates. Indeed, P. margaritifera mainly displays cytosine methylation in the CpG context, as already described in Crassostrea gigas (Thunberg, 1793) by Gavery and Roberts (2013) and in Biomphalaria glabrata (Say, 1818) by Adema et al. (2017). Likewise, differential methylation was also mainly identified in the CpG context, as in Wang et al. (2014) who reported that more than 99% of DNA methylation changes were restricted to the CpG context in C. gigas. The pattern of methylation we obtained occurred essentially in the gene bodies, with some accretion upstream and downstream of the gene, as described for mosaic methylation in other mollusks (Sarda et al., 2012; Wang et al., 2014; Rondon et al., 2017). These features of a short, but densely methylated region (corresponding to the genes), interspersed by long unmethylated regions (intergenic) are characteristic of the mosaic pattern of DNA methylation (Gavery and Roberts, 2013), which shows the typical bimodal distribution generally met in invertebrates. The lowly methylated genes were involved in the reproductive process, cellular signaling, and environmentally responsive functions. The low methylation of genes involved in the response to environmental changes is something reported in many different invertebrates (Sarda et al., 2012), but the presence of functions associated with reproduction is less common. It may be explained by the tissue used to produce this methylome, the mantle, which is a tissue not involved in reproduction (see correlation between methylation levels and gene expression Figure 2D). The genes identified in highly methylated regions were essentially involved in housekeeping functions (Gavery and Roberts, 2013; Olson and Roberts, 2014; Wang et al., 2014). When comparing the gene body methylation rate with the gene expression level, we demonstrated that moderately methylated genes have higher expression levels than lowly or highly methylated ones. This result is consistent with other studies (Feng et al., 2010; Xiang et al., 2010; Zemach et al., 2010) and strengthens the emerging hypothesis that, in invertebrates, gene expression and gene body methylation functions as a negative feedback loop in which gene expression increases with gene methylation until reaching a tipping point where additional methylation decreases transcription (Dixon et al., 2018).
Environmentally Induced DNA Methylation Changes and Their Link With Pigmentation
Among the DNA methylation changes that occurred during our experiment, several occurred in genes known to be involved in pigmentation pathways. Among these, the pteridine pathway was recently identified as a key player in the expression of the yellow color phenotype in P. margaritifera (Stenger et al., in press). Indeed, different derivates of pteridine can lead to the production of sepiapterin and xanthopterins, two yellow pigments (Ng et al., 2009). The folic acid pathway is also significantly affected by DNA methylation changes. Although its involvement in molluscan pigmentation is unknown, folic acid deficiency is linked to melanosis (e.g., melanin overproduction) in mammals, which results in a black pigmentation (Sharp et al., 1980). Since the implication of melanin in pearl oyster pigmentation has previously been identified (Lemer et al., 2015), an epigenetically driven modification of gene expression in the folic-acid pathway could be associated with the darkening color phenotypes expressed in response to an increase in depth.
At the gene level, GART, a hypermethylated gene included in three enriched GO categories (Figure 8), encodes a trifunctional purine biosynthetic protein, adenosine-3. This protein is involved in the de novo purine synthesis pathway (Amsterdam et al., 2004) and is composed of three subunits (a phosphoribosylglycinamide formyltransferase, a phosphoribosylglycinamide synthetase, and a phosphoribosylaminoimidazole synthetase). Biochemically, it catalyzes steps 2, 3, and 5 of inosine monophosphate (IMP) synthesis (Amsterdam et al., 2004; Ng et al., 2009). IMP is one of the precursors initiating the pterin and the Raper-Manson pathways, two pathways leading to pigmentation in P. margaritifera (Stenger et al., in press). Ng et al. (2009) have shown that mutations in GART are associated with pigmentation defects in juvenile zebrafish Danio rerio (Buchanan-Hamilton, 1822) due to disturbances of the pterin and Rapper-Mason pathways. Wild-type zebrafish are mainly yellow with black spots, while Δ-GART juveniles are entirely black (Ng et al., 2009). We can, therefore, hypothesize that methylation changes in the GART gene may affect its expression, subsequently affecting the pterin and Rapper-Mason pathways and leading to a darkening of the shell.
Figure 8. Methylation direction changes in candidate genes linked to inner shell color darkening in the depth treatment individuals at the three successive sampling times (T-S1, T-S2, and T-S3). The first four genes are involved in pigmentation processes, and the last two in biomineralization processes. Red “+” illustrates a significant hypermethylation of at least one CpG and blue “-” illustrates a significant hypomethylation of at least one CpG. +/T-S(a) and -/T-S(b) mean that the methylation in the present sampling time is more methylated than in the T-S(a) and less than in the T-S(b). GART, trifunctional purine biosynthetic protein adenosine-3; ABCC1, multidrug resistance-associated protein 1; MAPKAP1, target of rapamycin complex 2 subunit MAPKAP1; GRL101, G-protein coupled receptor GRL101; MGAT1, alpha-1,3-mannosyl-glycoprotein-2-beta-n-acetylglucosaminyltransferase. Pictures showed are illustrative.
The gene for multidrug resistance-associated protein 1 (ABCC1) presented two hypermethylated positions after the period at 30 m, a methylation state that reverted after the return to 8 m. ABCC1 is known to mediate ATP-dependent transport of glutathione and glutathione conjugates (Homolya et al., 2003). In a previous study it was proposed that glutathione plays an essential role in the expression of the yellow and black pigments (Stenger et al., in press). Glutathione-S-transferase (GST) activity is central in regulating the production of the yellow pheomelanin and black eumelanin pigments through the Raper-Manson pathway (Sonthalia et al., 2016). Methylation changes in the ABCC1 gene could therefore promote variation in the quantity of glutathione available and modify the regulation of the production of pheomelanin and eumelanin. An overproduction of eumelanin may explain the observed darkening of the shell.
GRL101 presented a hypermethylated response to the return to 8 m. According to Tanaka et al. (2014), this gene is an ortholog of the pigment dispersing factor, a gene responsible for changes in the concentration of chromatophoral pigment in response to darkness (Rao and Riehm, 1993). In crustaceans, it was proposed that color variation due to changing light conditions was caused by the dispersion of retinal chromatophore pigments linked with the activation of GRL101 (Rao and Riehm, 1993; Auerswald et al., 2008). The methylation change of GRL101, the similarities between the environmental triggers (a decrease of light) activating GRL101 in other organisms, and the phenotypes resulting from this activation argue in favor of the involvement of GRL101 in P. margaritifera color variation.
The GPI-anchor transamidase-like gene was hypomethylated after the return to 8 m depth. According to Ng and Freeze (2014), a mutation of the GPI-anchor transamidase genes is involved in a human disorder characterized by an altered dermal pigmentation (Ng and Freeze, 2014). This was later confirmed by RNAi experiments targeting GPI-anchor transamidase transcripts and resulting in a hyper-pigmented dark swellings in the maize anthracnose fungus Colletotrichum graminicola (G.W. Wilson 1914) (Oliveira-Garcia and Deising, 2016).
Enolase-4 is another candidate gene displaying hypomethylation at 30 m. Enolases are metalloenzymes involved in glycolysis and glycogen storage. One study reported a correlation between enolase activity and pigmentation: Hamid et al. (2015) showed that patients with vitiligo (a pigmentation disorder affecting melanocytes and inhibiting pigment synthesis) synthesize antibodies directed against enolases. Since their discovery, enolases have been used as biomarkers for the diagnosis, treatment, and monitoring of vitiligo (Hamid et al., 2015). The causes of this pathology are still unclear, although both genetic and environmental factors seem to be involved (Hamid et al., 2015). In the case of the pearl oyster, the hypomethylation of the enolase gene at 30 m may be associated with a change of expression inducing a darker phenotype.
The last of the genes subject to methylation change (hypermethylation after the last period at 8 m) and displaying a functional link with pigmentation is the target of rapamycin complex 2 subunit (MAPKAP1). This gene is involved in the TORC2 and TOR signaling pathways. The activation of these signaling pathways is known to promote a dark epithelial pigmentation due to the proliferation and the migration of retinal pigmentation epithelial cells (RPE cells) (Liu et al., 2017).
Among the six genes displaying a functional link with the pigmentation process or its regulation, four (GART, ABCC1, MAPKAP1, and GRL101) were associated with the expression of darker phenotypes, while the two others were associated with pigmentation disorders. Further experiments will be necessary to confirm and define their role, such as gene expression quantification, RNAi, and, once possible, genome editing, and epigenetic engineering. Thus, our results provide the first step toward this new research field.
Biomineralization and Pearl Darkness
In addition to a darkening of the coloration (Stenger et al., 2019), previous experiments have shown that a variation in depth also affects the shape and size of the aragonite tablets of the shell of P. margaritifera (Rousseau and Rollion-Bard, 2012). Aragonite tablets are the structural unit of nacre, the CaCO3 polymorph that constitutes the inner-shell of pearl oyster and the pearl itself (Rousseau and Rollion-Bard, 2012). Variation in the organization of these aragonite tablets can induce a change in color and luster (brightness) due to a change in the physical iridescence (Liu et al., 1999; Wang et al., 2008). Although not yet demonstrated, it is suspected that pigments contributing to nacre color are constituents of the intra-lamellar silk-fibroin gel that is localized between aragonite tablets (Addadi et al., 2006). Variation in the size of these tablets could therefore lead to a variation in the quantity of pigments that can be viewed through the last biomineralized aragonite layers. Interestingly, among the genes displaying methylation changes in response to a variation in water depth, two are well-known actors of the biomineralization processes of the nacreous layer: perlucin (Joubert et al., 2010) and MGAT1 (Takakura et al., 2008).
Perlucin is a protein found in the shell organic matrix of several Mollusca, including P. margaritifera (Weiss et al., 2000; Blank et al., 2003; Marie et al., 2013; Joubert et al., 2014). Experiments with purified Mollusca perlucin have suggested its involvement in calcium carbonate precipitation by favoring nucleation, crystallization and crystal growth control (Weiss et al., 2000). A variation in its expression can therefore have a huge effect on aragonite tablet size and organization (Rousseau and Rollion-Bard, 2012) and may thus modify the iridescence and transparency properties of the top aragonite layers (Rousseau and Rollion-Bard, 2012).
MGAT1 is a gene whose product initiates carbohydrate formation and is essential for the conversion of high-mannose to hybrid and complex N-glycans. This protein is involved in the protein glycosylation pathway, which is part of Protein modification. Interestingly, Takakura et al. (2008) identified an acidic N-glycan post-translationally attached to nacrein in Pinctada fucata that allows calcium binding. Moreover, nacrein is one of the main proteins found in the nacreous part of the shell (Rousseau and Rollion-Bard, 2012). So, although the MGAT1 gene plays no role in crystal formation, a possible link affecting nacrein formation can still be found. Future proteomics studies will be necessary to better uncover the role of MGAT1 in P. margaritifera shell coloration.
Epigenetics and Pearl Culture
As previously reported, the yo-yo experiment resulted in a general darkening of the inner shell of P. margaritifera in response to an increase in depth (Stenger et al., 2019), an environmentally induced phenotype that was maintained even after a return to the control depth (8 m). Such an enduring phenotypic response could be considered as good evidence of the involvement of epigenetic control (Sutherland and Costa, 2003; Harris et al., 2012; Bräutigam et al., 2013). The maintenance of this phenotype was previously documented in this species cultured for pearl production, and can last for over 18 months (Stenger et al., 2019). In another biological model, maize, stress-induced hypermethylation of P-pr (Richards, 2006; Lukens and Zhan, 2007) was associated with a reduced pigmentation that lasted in some cases for the entire life of an individual, and could even be transmitted to the next generation (Richards, 2006; Bossdorf et al., 2008). Such effects could offer huge benefits for pearl farming. First, this long-lasting effect suggests that farmers could better control their production through dedicated conditioning of recipient and/or donor oysters. Additionally, the transgenerational effect described for maize, although not yet tested for pearl oysters, suggests that epigenetic marker-assisted selection could be envisioned. Such an approach may offer the possibility of selecting phenotypes of interest without the associated risk of eroding genetic diversity and/or the integration into natural populations of spat produced by farmed oysters (Reisser et al., 2020).
Data Availability Statement
The datasets generated for this study can be found on the NCBI (BioProject PRJNA663978). Scripts used are provided on GitHub (PLStenger/Pearl_Oyster_Colour_BS_Seq/00_scripts).
Author Contributions
JV-D, P-LS, SP, and C-LK designed the study. P-LS and MM performed the experiments. P-LS, CC, and CG analyzed the data. P-LS and JV-D drafted the manuscript. Funding was obtained by CL-K, SP, and JV-D. All authors approved the manuscript.
Funding
This study was supported by grants from the “Direction des Ressources Marines,” through the AmeliGEN project (# 10065/MEI/DRMM). P-LS’s Ph.D. was funded by AmeliGEN and the Pacific Doctoral School (EDP) (ED 469). This study was performed within the framework of the “Laboratoire d’Excellence (LABEX)” TULIP (ANR-10-LABX-41) and CORAIL (ANR-10-LABEX).
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.
Acknowledgments
The authors would like to thank the Regahiga Pearl Farm (Mangareva Island, Gambier Archipelago, French Polynesia) for providing the pearl oysters used in this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.630290/full#supplementary-material
Supplementary Figure 1 | Treemaps of biological processes obtained from the GO terms significantly enriched from the lowly methylated genes (A) (data from Figure 2C – “1” part) and the highly methylated genes (B) (data from Figure 2C – “2” part).
Supplementary Figure 2 | Treemaps of biological processes and molecular functions obtained from the GO terms significantly enriched from the genes with differentially methylated positions in control individuals (C-S1, control in time 1; C-S2, control in time 2; C-S3, control in time 3).
Supplementary Table 1 | FastQC results for raw and trimmed data.
Supplementary Table 2 | Average values by category (control/depth treatment) and sampling time for: the number of read pairs, paired-end alignments, mapping efficiency, number of duplicated alignments removed after mapping, the number of C in sequences, the number of C in percentage in CpG, CHG, CHH or in unknown contexts, and the basic statistics of the gene body methylation rate (GBM – on the 54408 genes of the reference genome).
Supplementary Table 3 | All significantly differentially methylated positions (q value < 0.05) for the controls. Scaffold name: scaffold number; Scaffold size: scaffold size; Position: position of the significantly differentially methylated positions on this scaffold; Genomic element: the genomic feature containing the DNA methylation change; q value: q value; Diff-M: the differential methylation level obtained using the Methylkit R package; Gene annotation name: Gene annotation according to the nr and swiss-prot blast top hit. Blue shading (negative values) indicates hypomethylation of the position and red shading (positive values) hypermethylation; darker shading indicates a stronger degree of hypomethylation of hypermethylation, respectively. The last column indicates in which comparison the significantly different methylated positions were found.
Footnotes
References
Addadi, L., Joester, D., Nudelman, F., and Weiner, S. (2006). Mollusk shell formation: a source of new concepts for understanding biomineralization processes. Chemistry 12, 980–987. doi: 10.1002/chem.200500980
Adema, C. M., Hillier, L. D. W., Jones, C. S., Loker, E. S., Knight, M., Minx, P., et al. (2017). Whole genome analysis of a schistosomiasis-transmitting freshwater snail. Nat. Commun. 8:15451. doi: 10.1038/ncomms15451
Akalin, A., Kormaksson, M., Li, S., Garrett-Bakelman, F. E., Figueroa, M. E., Melnick, A., et al. (2012). MethylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 13:R87. doi: 10.1186/gb-2012-13-10-R87
Amsterdam, A., Nissen, R. M., Sun, Z., Swindell, E. C., Farrington, S., and Hopkins, N. (2004). Identification of 315 genes essential for early zebrafish development. Proc. Natl. Acad. Sci. U.S.A. 101, 12792–12797. doi: 10.1073/pnas.0403929101
Anders, S., Pyl, P. T., and Huber, W. (2015). Genome analysis HTSeq — a Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169. doi: 10.1093/bioinformatics/btu638
Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data. Cambridge: Babraham Bioinforma.
Auerswald, L., Freier, U., Lopata, A., and Meyer, B. (2008). Physiological and morphological colour change in Antarctic krill, Euphausia superba: a field study in the Lazarev Sea. J. Exp. Biol. 211, 3850–3858. doi: 10.1242/jeb.024232
Blank, S., Arnoldi, M., Khoshnavaz, S., Treccani, L., Kuntz, M., Mann, K., et al. (2003). The nacre protein perlucin nucleates growth of calcium carbonate crystals. J. Microsc. 212(Pt 3), 280–291. doi: 10.1111/j.1365-2818.2003.01263.x
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Genome analysis trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bossdorf, O., Richards, C. L., and Pigliucci, M. (2008). Epigenetics for ecologists. Ecol. Lett. 11, 106–115. doi: 10.1111/j.1461-0248.2007.01130.x
Bouzerand, E. (2018). Points Forts de la Polynésie Française – Bilan – La Perle en 2016. Papeete: Institu de la Statistique de la Polynésia française.
Bräutigam, K., Vining, K. J., Lafon-Placette, C., Fossdal, C. G., Mirouze, M., Marcos, J. G., et al. (2013). Epigenetic regulation of adaptive responses of forest tree species to the environment. Ecol. Evol. 3, 399–415. doi: 10.1002/ece3.461
Chakraborty, A. K., Funasaka, Y., Ichihashi, M., Sodi, S., Bhattacharya, M., and Pawelek, J. (1999). Upregulation of mRNA for the melanocortin-1 receptor but not for melanogenic proteins in macrophage x melanoma fusion hybrids exhibiting increased melanogenic and metastatic potential. Pigment Cell Res. 12, 355–366. doi: 10.1111/j.1600-0749.1999.tb00519.x
Charpentier, V., Phillips, C. S., and Méry, S. (2012). Pearl fishing in the ancient world: 7500 BP. Arab. Archaeol. Epigr. 23, 1–6. doi: 10.1111/j.1600-0471.2011.00351.x
Conesa, A., Götz, S., García-gómez, J. M., Terol, J., Talón, M., Genómica, D., et al. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. BMC Bioinformatics 21:3674–3676. doi: 10.1093/bioinformatics/bti610
Cuthill, I. C., Allen, W. L., Arbuckle, K., Caspers, B., Chaplin, G., Hauber, M. E., et al. (2017). The biology of color. Science 357:eaan0221. doi: 10.1126/science.aan0221
Dixon, G., Liao, Y., Bay, L. K., and Matz, M. V. (2018). Role of gene body methylation in acclimatization and adaptation in a basal metazoan. Proc. Natl. Acad. Sci. U.S.A. 115, 13342–13346. doi: 10.1073/pnas.1813749115
Dolinoy, D. C. (2008). The agouti mouse model: an epigenetic biosensor for nutritional and environmental alterations on the fetal epigenome. Nutr. Rev. 66, 7–11. doi: 10.1111/j.1753-4887.2008.00056.x
Dolinoy, D. C., Weidman, J. R., Waterland, R. A., and Jirtle, R. L. (2006). Maternal genistein alters coat color and protects Avy mouse offspring from obesity by modifying the fetal epigenome. Environ. Health Perspect. 114, 567–572. doi: 10.1289/ehp.8700
Feng, D., Li, Q., Yu, H., Kong, L., and Du, S. (2018). Transcriptional profiling of long non-coding RNAs in mantle of Crassostrea gigas and their association with shell pigmentation. Sci. Rep. 8:1436. doi: 10.1038/s41598-018-19950-6
Feng, D., Li, Q., Yu, H., Liu, S., Kong, L., and Du, S. (2020). Integrated analysis of microRNA and mRNA expression profiles in Crassostrea gigas to reveal functional miRNA and miRNA-targets regulating shell pigmentation. Sci. Rep. 10:20238. doi: 10.1038/s41598-020-77181-0
Feng, S., Cokus, S. J., Zhang, X., Chen, P. Y., Bostick, M., Goll, M. G., et al. (2010). Conservation and divergence of methylation patterning in plants and animals. Proc. Natl. Acad. Sci. U.S.A. 107, 8689–8694. doi: 10.1073/pnas.1002720107
Fernandes, J., and Gattass, C. R. (2009). Topological polar surface area defines substrate transport by multidrug resistance associated protein 1 (MRP1/ABCC1). J. Med. Chem. 52, 1214–1218. doi: 10.1021/jm801389m
Frost, S. K., and Malacinski, G. M. (1979). The developmental genetics of pigment mutants in the Mexican axolotl. Dev. Genet. 1, 271–294. doi: 10.1002/dvg.1020010402
Gavery, M. R., and Roberts, S. B. (2013). Predominant intragenic methylation is associated with gene expression characteristics in a bivalve mollusc. PeerJ 2013, 1–15. doi: 10.7717/peerj.215
Gavery, M. R., and Roberts, S. B. (2017). Epigenetic considerations in aquaculture. PeerJ 5:e4147. doi: 10.7717/peerj.4147
Gueguen, Y., Montagnani, C., Joubert, C., Marie, B., Belliard, C., Tayale, A., et al. (2013). “Characterization of molecular processes involved in the pearl formation in Pinctada margaritifera for the sustainable development of pearl farming industry in French Polynesia,” in Recent Advances in Pearl Research–Proceedings of the International Symposium on Pearl Research 2011, eds S. Watabe, K. Maeyama, and H. Nagasawa (Setagaya: TERRAPUB), 183–193.
Hamid, K. M., Bitsue, Z. K., and Mirshafiey, A. (2015). Autoantibodies profile in vitiligo. J. Pigment. Disord. 2, 1–5. doi: 10.4172/2376-0427.1000187
Harris, K. D. M., Bartlett, N. J., and Lloyd, V. K. (2012). Daphnia as an emerging epigenetic model organism. Genet. Res. Int. 2012, 1–8. doi: 10.1155/2012/147892
Hollander, M., Wolfe, D. A., and Chicken, E. (1973). “The one-way layout,” in Nonparametric Statistical Methods, 3rd Edn, eds John Wiley & sons (New York, NY: John Wiley and Sons), 202–288.
Homolya, L., Váradi, A., and Sarkadi, B. (2003). Multidrug resistance-associated proteins: export pumps for conjugates with glutathione, glucuronate or sulfate. BioFactors 17, 103–114. doi: 10.1002/biof.5520170111
Hu, J., and Barrett, R. D. H. (2017). Epigenetics in natural animal populations. J. Evol. Biol. 30, 1612–1632. doi: 10.1111/jeb.13130
Iwahashi, Y., and Akamatsu, S. (1994). Porphyrin pigment in Black-Lip pearls and Its application to pearl identification. Fish. Sci. 60, 69–71. doi: 10.2331/fishsci.60.69
Joubert, C., Linard, C., Le Moullac, G., Soyez, C., Saulnier, D., Teaniniuraitemoana, V., et al. (2014). Temperature and food influence shell growth and mantle gene expression of shell matrix proteins in the pearl oyster Pinctada margaritifera. PLoS One 9:e103944. doi: 10.1371/journal.pone.0103944
Joubert, C., Piquemal, D., Marie, B., Manchon, L., Pierrat, F., Zanella-Cléon, I., et al. (2010). Transcriptome and proteome analysis of Pinctada margaritifera calcifying mantle and shell: focus on biomineralization. BMC Genomics 11:613. doi: 10.1186/1471-2164-11-613
Katz, M. L., Drea, C. M., and Robison, W. G. Jr. (1987). Dietary vitamins A and E influence retinyl ester composition and content of the retinal pigment epithelium. Biochim. Biophvs. Acta 924, 432–441.
Klopfenstein, D. V., Zhang, L., Pedersen, B. S., Ramírez, F., Vesztrocy, A. W., Naldi, A., et al. (2018). GOATOOLS: a python library for gene ontology analyses. Sci. Rep. 8:10872. doi: 10.1038/s41598-018-28948-z
Krueger, F., and Andrews, S. R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572. doi: 10.1093/bioinformatics/btr167
Ky, C. L., Broustal, F., Potin, D., and Lo, C. (2019). The pearl oyster (Pinctada margaritifera) aquaculture in French Polynesia and the indirect impact of long-distance transfers and collection-culture site combinations on pearl quality traits. Aquac. Rep. 13:100182. doi: 10.1016/j.aqrep.2019.100182
Ky, C.-L., Blay, C., Sham-Koua, M., Lo, C., and Cabral, P. (2014). Indirect improvement of pearl grade and shape in farmed Pinctada margaritifera by donor “oyster” selection for green pearls. Aquaculture 432, 154–162. doi: 10.1016/j.aquaculture.2014.05.002
Ky, C.-L., Blay, C., Sham-Koua, M., Vanaa, V., Lo, C., and Cabral, P. (2013). Family effect on cultured pearl quality in black-lipped pearl oyster Pinctada margaritifera and insights for genetic improvement. Aquat. Living Resour. 26, 133–145. doi: 10.1051/alr/2013055
Ky, C.-L., Le Pabic, L., Sham Koua, M., Nicolas, M., Seiji, N., and Devaux, D. (2017). Is pearl colour produced from Pinctada margaritifera predictable through shell phenotypes and rearing environments selections? Aquaculture 48, 1041–1057. doi: 10.1111/are.12947
Le Moullac, G., Schuck, L., Chabrier, S., Belliard, C., Lyonnard, P., Broustal, F., et al. (2018). Influence of temperature and pearl rotation on biomineralization in the pearl oyster, Pinctada margaritifera. J. Exp. Biol. 221:jeb186858. doi: 10.1242/jeb.186858
Le Pabic, L., Parrad, S., Sham Koua, M., Nakasai, S., Saulnier, D., Devaux, D., et al. (2016). Culture site dependence on pearl size realization in Pinctada margaritifera in relation to recipient oyster growth and mantle graft biomineralization gene expression using the same donor phenotype. Estuar. Coast. Shelf Sci. 182, 294–303. doi: 10.1016/j.ecss.2016.03.009
Le Pennec, M. (2010). Huître Perlière et Perle de Tahiti. Puna’auia: Université de la Polynésie française.
Lemer, S., Saulnier, D., Gueguen, Y., and Planes, S. (2015). Identification of genes associated with shell color in the black-lipped pearl oyster, Pinctada margaritifera. BMC Genomics 16:568. doi: 10.1186/s12864-015-1776-x
Liu, Y., Chen, Z., Cheng, H., Chen, J., and Qian, J. (2017). Gremlin promotes retinal pigmentation epithelial (RPE) cellproliferation, migration and VEGF production via activating VEGFR2-Akt-mTORC2 signaling. Oncotarget 8, 979–987.
Liu, Y., Shigley, J., and Hurwit, K. (1999). Iridescent color of a shell of the mollusk Pinctada margaritifera caused by diffraction. Opt. Express 4, 177–182. doi: 10.1364/OE.4.000177
Lukens, L. N., and Zhan, S. (2007). The plant genome’s methylation status and response to stress: implications for plant improvement. Curr. Opin. Plant Biol. 10, 317–322. doi: 10.1016/j.pbi.2007.04.012
Marie, B., Joubert, C., Tayalé, A., and Zanella-cléon, I. (2013). Different secretory repertoires control the biomineralization processes of prism and nacre deposition of the pearl oyster shell. Proc. Natl. Acad. Sci. U.S.A. 109, 20986–20991. doi: 10.1073/pnas.1210552109
Moreau, C., Ambrose, M. J., Turner, L., Hill, L., Noel Ellis, T. H., and Hofer, J. M. I. (2012). The b gene of pea encodes a defective flavonoid 3’,5’-hydroxylase, and confers pink flower color. Plant Physiol. 159, 759–768. doi: 10.1104/pp.112.197517
Ng, A., Uribe, R. A., Yieh, L., Nuckels, R., and Gross, J. M. (2009). Zebrafish mutations in gart and paics identify crucial roles for de novo purine synthesis in vertebrate pigmentation and ocular development. Development 136, 2601–2611. doi: 10.1242/dev.038315
Ng, B. G., and Freeze, H. H. (2014). Human genetic disorders involving glycosylphosphatidylinositol (GPI) anchors and glycosphingolipids (GSL). J. Inherit. Metab. Dis. 38, 171–178. doi: 10.1007/s10545-014-9752-1
Nguyen, V. H., and Lavenier, D. (2009). PLAST: parallel local alignment search tool for database comparison. BMC Bioinformatics 10:329. doi: 10.1186/1471-2105-10-329
Nicoglou, A., and Merlin, F. (2017). Epigenetics: a way to bridge the gap between biological fields. Stud. Hist. Philos. Sci. Part C Stud. Hist. Philos. Biol. Biomed. Sci. 66, 73–82. doi: 10.1016/j.shpsc.2017.10.002
Oliveira-Garcia, E., and Deising, H. B. (2016). The glycosylphosphatidylinositol anchor biosynthesis genes GPI12, GAA1, and GPI8 are essential for cell-wall integrity and pathogenicity of the maize anthracnose fungus Colletotrichum graminicola. Mol. Plant Microbe Interact. 29, 889–901. doi: 10.1094/MPMI-09-16-0175-R
Olson, C. E., and Roberts, S. B. (2014). Genome-wide profiling of DNA methylation and gene expression in Crassostrea gigas male gametes. Front. Physiol. 5:224. doi: 10.3389/fphys.2014.00224
Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. doi: 10.1093/bioinformatics/btq033
Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A., and Manke, T. (2014). DeepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 42, W187–W191. doi: 10.1093/nar/gku365
Reisser, C. M. O., Gendre, R., Le, Chupeau, C., Lo-Yat, A., Planes, S., et al. (2020). Population connectivity and genetic assessment of exploited and natural populations of pearl oysters within a French polynesian atoll lagoon. Genes (Basel). 11, 1–16. doi: 10.3390/genes11040426
Richards, E. J. (2006). Inherited epigenetic variation—revisiting soft inheritance. Nat. Rev. Genet. 7, 395–402. doi: 10.1038/nrg1834
Rocha, G. D. G., Oliveira, R. R., Kaplan, M. A. C., and Gattass, C. R. (2014). 3β-Acetyl tormentic acid reverts MRP1/ABCC1 mediated cancer resistance through modulation of intracellular levels of GSH and inhibition of GST activity. Eur. J. Pharmacol. 741, 140–149. doi: 10.1016/j.ejphar.2014.07.054
Rondon, R., Grunau, C., Fallet, M., Charlemagne, N., Sussarellu, R., Chaparro, C., et al. (2017). Effects of a parental exposure to diuron on Pacific oyster spat methylome. Environ. Epigenetics 3, 1–13. doi: 10.1093/eep/dvx004
Rousseau, M., and Rollion-Bard, C. (2012). Influence of the depth on the shape and thickness of nacre tablets of Pinctada margaritifera pearl oyster, and on oxygen isotopic composition. Minerals 2, 55–64. doi: 10.3390/min2010055
Royston, P. (1995). Remark AS R94: a remark on algorithm AS 181: the W-test for normality. J. R. Stat. Soc. Ser. C Appl. Stat. 44, 547–551. doi: 10.2307/2986146
Russo, V. E., Martienssen, R. A., and Riggs, A. D. (1996). Epigenetic Mechanisms of Gene Regulation. Cold Spring Harbor, NY: Cold Spring Harbor Laboratory Press.
Sarda, S., Zeng, J., Hunt, B. G., and Yi, S. V. (2012). The evolution of invertebrate gene body methylation. Mol. Biol. Evol. 29, 1907–1916. doi: 10.1093/molbev/mss062
Sharp, J. R., Insalaco, S. J., and Johnson, L. F. (1980). “Melanosis” of the duodenum associated with a gastric ulcer and folic acid deficiency. Gastroenterology 78, 366–369. doi: 10.1016/0016-5085(80)90590-9
Shin, D. H., Cho, M., Choi, M. G., Das, P. K., Lee, S. K., Choi, S. B., et al. (2015). Identification of genes that may regulate the expression of the transcription factor production of anthocyanin pigment 1 (PAP1)/MYB75 involved in Arabidopsis anthocyanin biosynthesis. Plant Cell Rep. 34, 805–815. doi: 10.1007/s00299-015-1743-7
Sokolov, E. P. (2000). An improved method for DNA isolation from mucopolysaccharide-rich molluscan tissues. J. Molluscan Stud. 66, 573–575. doi: 10.1093/mollus/66.4.573
Sonthalia, S., Daulatabad, D., and Sarkar, R. (2016). Glutathione as a skin whitening agent: facts, myths, evidence and controversies. Indian J. Dermatol. 82, 262. doi: 10.4103/0378-6323.179088
Stenger, P.-L. (2017). Package ImaginR. Available online at: https://cran.r-project.org/web/packages/ImaginR/ImaginR.pdf (accessed May 29, 2017).
Stenger, P.-L., Ky, C.-L., Reisser, C., Duboisset, J., and Dicko, H. (in press). Molecular pathways and pigments behind the colors of the pearl oyster Pinctada margaritifera var. cumingii (Linnaeus 1758). Genes.
Stenger, P.-L., Vidal-Dupiol, J., Reisser, C., Planes, S., and Ky, C. (2019). Colour plasticity in the shells and pearls of animal graft model Pinctada margaritifera through colour quantification with the HSV system. Sci. Rep. 75:20. doi: 10.1038/s41598-019-43777-4
Sutherland, J. E., and Costa, M. A. X. (2003). Epigenetics and the environment. Ann. N. Y. Acad. Sci. 983, 151–160.
Takakura, D., Norizuki, M., Ishikawa, F., and Samata, T. (2008). Isolation and characterization of the N-linked oligosaccharides in nacrein from Pinctada fucata. Mar. Biotechnol. 10, 290–296. doi: 10.1007/s10126-007-9063-8
Tanaka, Y., Suetsugu, Y., Yamamoto, K., Noda, H., and Shinoda, T. (2014). Transcriptome analysis of neuropeptides and G-protein coupled receptors (GPCRs) for neuropeptides in the brown planthopper Nilaparvata lugens. Peptides 53, 125–133. doi: 10.1016/j.peptides.2013.07.027
Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.1038/nprot.2012.016
Vajro, P., Thaler, M. M., and Blanckaert, N. (1995). Bile pigment composition and bilirubin esterification in the developing chick. Pediatr. Res. 38, 349–355. doi: 10.1203/00006450-199509000-00013
Wang, N., Lee, Y. H., and Lee, J. (2008). Recombinant perlucin nucleates the growth of calcium carbonate crystals: molecular cloning and characterization of perlucin from disk abalone, Haliotis discus discus. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 149, 354–361. doi: 10.1016/j.cbpb.2007.10.007
Wang, X., Li, Q., Lian, J., Li, L., Jin, L., Cai, H., et al. (2014). Genome-wide and single-base resolution DNA methylomes of the Pacific oyster Crassostrea gigas provide insight into the evolution of invertebrate CpG methylation. BMC Genomics 15:1119. doi: 10.1186/1471-2164-15-1119
Weiss, I. M., Kaufmann, S., Mann, K., and Fritz, M. (2000). Purification and characterization of perlucin and perlustrin, two new proteins from the shell of the mollusc Haliotis laevigata. Biochem. Biophys. Res. Commun. 267, 17–21. doi: 10.1006/bbrc.1999.1907
Xiang, H., Zhu, J., Chen, Q., Dai, F., Li, X., Li, M., et al. (2010). Single base-resolution methylome of the silkworm reveals a sparse epigenomic map. Nat. Biotechnol. 28, 516–520. doi: 10.1038/nbt.1626
Zemach, A., McDaniel, I. E., Silva, P., and Zilberman, D. (2010). Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science 328, 916–919. doi: 10.1126/science.1186366
Zhai, S., Xia, X., and He, Z. (2016). Carotenoids in staple cereals: metabolism, regulation, and genetic manipulation. Front. Plant Sci. 7:1197. doi: 10.3389/fpls.2016.01197
Keywords: pearl oyster, environmental pressure, depth, color change, pigmentation, DNA methylation, methylome characterization
Citation: Stenger P-L, Ky C-L, Reisser CMO, Cosseau C, Grunau C, Mege M, Planes S and Vidal-Dupiol J (2021) Environmentally Driven Color Variation in the Pearl Oyster Pinctada margaritifera var. cumingii (Linnaeus, 1758) Is Associated With Differential Methylation of CpGs in Pigment- and Biomineralization-Related Genes. Front. Genet. 12:630290. doi: 10.3389/fgene.2021.630290
Received: 17 November 2020; Accepted: 19 February 2021;
Published: 19 March 2021.
Edited by:
Stephanie McKay, University of Vermont, United StatesReviewed by:
Silvio Zaina, University of Guanajuato, MexicoYu Jiao, Guangdong Ocean University, China
Copyright © 2021 Stenger, Ky, Reisser, Cosseau, Grunau, Mege, Planes and Vidal-Dupiol. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jeremie Vidal-Dupiol, amVyZW1pZS52aWRhbC5kdXBpb2xAaWZyZW1lci5mcg==