Age group DNA methylation differences in lemon sharks (Negaprion brevirostris): Implications for future age estimation tools.

Age information is often non-existent for most shark populations due to a lack of measurable physiological and morphological traits that can be used to estimate age. Recently, epigenetic clocks have been found to accurately estimate age for mammals, birds, and fish. However, since these clocks rely, among other things, on the availability of reference genomes, their application is hampered in non-traditional model organisms lacking such molecular resources. The technique known as Methyl-Sensitive Amplified Polymorphism (MSAP) has emerged as a valid alternative for studying DNA methylation biomarkers when reference genome information is missing, and large numbers of samples need to be processed. Accordingly, the MSAP technique was used in the present study to characterize global DNA methylation patterns in lemon sharks from three different age groups (juveniles, subadults, and adults). The obtained results reveal that, while MSAP analyses lack enough resolution as a standalone approach to infer age in these organisms, the global DNA methylation patterns observed using this technique displayed significant differences between age groups. Overall, these results confer that DNA methylation does change with age in sharks like what has been seen for other vertebrates and that MSAP could be useful as part of an epigenetics pipeline to infer the broad range of ages found in large samples sizes.

INTRODUCTION

Overfishing constitutes the most important threat to elasmobranchs including sharks, skates, and rays (Dulvy et al., 2014), making fisheries management critical for their conservation. Consequently, the characterization of demographic information is of utmost importance to model the sustainability of populations under fishing pressure. Age information is a critical parameter of these models, which is unfortunately not available for most species due (among other reasons) to the lack of physiological and morphological traits that can be accurately used as proxies of chronological age in elasmobranchs. Particularly in the case of sharks, body length has been largely used to infer age and has often been cross‐validated with vertebral growth bands (Carlson & Goldman, 2007). Yet, recent studies have reported that the age of older sharks has been vastly underestimated given that vertebral bands do not necessarily form on a yearly basis throughout a shark's lifetime (Natanson et al., 2018; Natanson & Skomal, 2015). Not less importantly, vertebral band counts are highly invasive, require sacrificing animals, and do not work for all species (Francis et al., 2007; Huveneers et al., 2013; Natanson et al., 2018).

The emergence of methodological approaches linking molecular variation to chronological age has constituted a breakthrough, particularly in those cases where age estimation has proven difficult and/or based on invasive methods. The characterization of distinctive molecular traits is possible using small tissue samples that can be more easily obtained, in most cases without even the need to capture/release the animal. One of the first among these methods was based on the quantification of the length of telomeres, based on their progressive shortening during successive cell cycle divisions. While telomere‐based aging methods proved to be a great proxy for biological age (i.e., how old an animal looks for its age), the vast differences in the starting length observed in telomeres led numerous studies to question its applicability for determining chronological age estimation (i.e., actual age of the organism) (Nehmens et al., 2021; Remot et al., 2020; Vaiserman & Krasnienkov, 2020). Epigenetics, defined as “the study of molecules and mechanisms that can perpetuate alternative gene activity states in the context of the same DNA sequence” (Cavalli & Heard, 2019), has provided a basis to start solving that problem by establishing direct links between chronological age and particular epigenetic modifications including DNA methylation, histone variants, histone modifications, and small RNAs (López‐Otín et al., 2013), providing not only insights into how organisms age but also facilitating the estimation of age using DNA methylation (the addition of a methyl group most commonly to a cytosine in the DNA sequence) as a proxy (Horvath, 2013).

Over the last decade, several epigenetic clocks using DNA methylation to estimate chronological (actual or real age) and biological age (physiological age, how old an individual seems) have been successfully developed across diverse taxa (Beal et al., 2019; Bors et al., 2021; Chen et al., 2016; De Paoli‐Iseppi et al., 2019; Horvath & Raj, 2018; Polanowski et al., 2014; Salameh et al., 2020). Chronological epigenetic clocks specifically correlate percent DNA methylation at CpG sites (a Cytosine followed by a Guanine in the same DNA strand) to age through multiple regression analyses producing highly accurate estimations (correlation coefficient, r 2, ranging from 0.70 to 0.99). These changes in DNA methylation are linked to age‐related mechanisms and developmental milestones that often result in changes in gene expression (Gladyshev, 2016; Hannum et al., 2013; López‐Otín et al., 2013). Although largely focused on mammalian species (Barratclough et al., 2021; Beal et al., 2019; Bors et al., 2021; Fei et al., 2021; Horvath & Raj, 2018; Polanowski et al., 2014; Robeck et al., 2021; Stubbs et al., 2017; Tanabe et al., 2020), accurate clocks have been developed for other taxa including birds and fishes (Anastasiadi & Piferrer, 2020; De Paoli‐Iseppi et al., 2019; Mayne et al., 2020).

Chronological epigenetic clocks get their accuracy from the characterization of DNA methylation at particular DNA positions using high‐resolution approaches such as microarrays, reduced representation bisulfite sequencing, pyrosequencing, and qPCR (Morselli et al., 2021; Trigg et al., 2021). However, such techniques are hampered by two major limitations: the requirement of a sequenced reference genome for the species studied and the potentially high cost of conducting these analyses in large numbers of samples to be ecologically meaningful in field studies, particularly for the larger sample sizes necessary for ecologically relevant population and stock assessments. Thus, the application of less demanding and inexpensive coarse‐grained methods, still providing useful aging information at the level of age intervals, represents a compromise helping the characterization of aging structures in natural populations within realistic ecological settings. The characterization of Methyl‐Sensitive Amplified Polymorphism (MSAP) is one such method, based on the use of dual methylation‐sensitive restriction enzymes to unveil broad DNA methylation patterns found at target cut site 5’‐CCGG‐3′ (invented by Reyna‐López et al., 1997). This method has been successfully used to examine the links between environmental change and genome‐wide DNA methylation patterns involved in phenotypic responses (Beal et al., 2021; Huang et al., 2017; Morán et al., 2013; Pierron et al., 2014; Rodríguez‐Casariego et al., 2020; Suarez‐Ulloa et al., 2019; Zhao et al., 2015).

Since it is well known that gene expression patterns change during development and aging, the observation of global amounts of DNA methylation changing over an individual's lifetime comes as no surprise (Heyn et al., 2012). However, it is not yet known if such changes are large enough to be efficiently detected using MSAP as a standalone epigenetic aging method. The present study addresses that question by using MSAP in lemon shark (Negaprion brevirostris) samples belonging to three age groups (juvenile, subadult, and adult), available through the unique long‐term sampling archive from Bimini, Bahamas. This archive of samples is one of the most comprehensive datasets for sharks and is the best known age dataset to date for any shark species. The MSAP method used in the present study was previously adjusted for use in lemon sharks from populations monitored in our own research (Beal et al., 2021). The obtained results are consistent with the applicability of MSAP analyses as part of a larger epigenetic pipeline to determine chronological age, particular at earlier stages, to define age groups in large numbers of samples which can be subsequently more finely characterized using single nucleotide resolution methods.

METHODS

Samples and DNA extraction

Lemon sharks, as studied in Bimini (Bahamas) by the Bimini Shark Lab for more than 30 years, exhibit a life history of biennial parturition after reaching sexual maturity around age 12 (Brown & Gruber, 1988; Feldheim et al., 2004) as well as natal philopatry where they return back to the particular nursery site where they were born (Feldheim et al., 2014). Although age of reproductive senescence is unknown, pregnant females in the low to mid‐30s have been recorded in the Bimini area and the maximum lifespan is estimated at 37 years (Brooks et al., 2016).

Epigenetic age estimation has been successful across multiple tissue types (Horvath, 2013). Therefore, we focused on fin clips as this is the type of tissue sample most often collected during shark population monitoring efforts and therefore most relevant for managers and conservation practitioners. All lemon shark samples analyzed in this study consisted of fin clips housed in an archived collection at the Field Museum in Chicago with collection dates ranging from 1995 to 2017. Most samples were collected in Bimini, Bahamas (Bahamas Research Permits MA&MR/FIS/17B) except for three samples that were collected in Jupiter, FL (Table 1). All fin clips were preserved in Dimethyl Sulfoxide (DMSO) at room temperature for up to 2 years and then stored in −20°C freezer thereafter. Upon original collection, samples were designated as one of the following age groups/life stages following Feldheim et al. (2014). Newborn sharks had an open umbilical scar, and juvenile sharks had closed umbilical scars but did not exceed 70 cm precaudal length (PCL) (1–5 years). Subadults (6–11 years) were between 70–175 cm PCL for males and 70–185 cm for females, and adults (12+ years) were greater than 175 cm PCL for males and greater than 185 cm for females.

TABLE 1

Samples and metadata. Batch 1 samples include LS1‐LS10. Batch 2 includes all remaining samples listed

ID	Batch	Capture date	Age group for testing	Sex	Collection site
LS1	1	2/2/2017	Subadult	M	Bimini
LS2	1	12/12/2016	Adult	M	Bimini
LS3	1	3/16/2009	Subadult	F	Jupiter
LS4	1	2/22/2009	Adult	F	Jupiter
LS5	1	2/21/2009	Adult	F	Jupiter
LS6	1	1/29/2013	Juvenile	M	Bimini
LS7	1	7/16/2012	Juvenile	M	Bimini
LS8	1	5/31/2013	Juvenile	M	Bimini
LS9	1	6/12/2012	Subadult	M	Bimini
LS10	1	2/2/2015	Juvenile	F	Bimini
LS11	2	4/16/2004	Subadult	F	Bimini
LS12	2	3/25/2004	Subadult	F	Bimini
LS13	2	4/7/2004	Subadult	M	Bimini
LS14	2	11/24/1996	Subadult	M	Bimini
LS15	2	7/15/1997	Subadult	M	Bimini
LS16	2	7/16/1997	Subadult	F	Bimini
LS17	2	7/23/1997	Subadult	M	Bimini
LS20	2	3/9/1997	Subadult	M	Bimini
LS22	2	5/9/2008	Adult	F	Bimini
LS24	2	2/25/1996	Adult	F	Bimini
LS25	2	6/8/1996	Adult	M	Bimini
LS26	2	7/21/1996	Adult	F	Bimini
LS29	2	6/16/1994	Subadult	F	Bimini
A13	2	7/7/2009	Juvenile	F	Bimini
A14	2	7/8/2009	Juvenile	F	Bimini
A59	2	6/20/2002	Juvenile	M	Bimini
A60	2	6/20/2002	Juvenile	M	Bimini
A61	2	6/19/2002	Juvenile	F	Bimini
A62	2	6/20/2002	Juvenile	F	Bimini

The samples used in this study were processed in two batches. The first batch (Batch 1, LS1‐10) was analyzed as part of a proof‐of‐concept pilot study, where DNA extraction was performed at the Field Museum in Chicago using a salting‐out protocol (Sunnucks & Hales, 1996). These samples were preserved in ice and shipped to Florida International University (FIU) for MSAP analyses. Batch 2 samples were shipped as fin clip samples in DMSO, and DNA was subsequently extracted following the same salting‐out protocol used for Batch 1 samples (Beal et al., 2021). Batch 1 samples (n = 10) consisted of four juveniles, three subadults, and three adults. Batch 2 samples (n = 19) consisted of six juveniles, nine subadults, and four adults.

DNA methylation pattern analysis

DNA methylation patterns were determined using MSAP analyses. This method employs dual methylation‐sensitive enzymes that cut differentially at the same target site 5’‐CCGG‐3′ (Reyna‐López et al., 1997) creating dual fragment profiles identifying four types of DNA methylation states present (full or hypermethylation, internal C methylation, hemi‐methylation, or no methylation). MSAP analyses were conducted independently for each of the two batches of samples as detailed in Beal et al. (2021). Briefly, a parallel restriction digestion/ligation (each containing one of the methyl‐sensitive restriction enzymes MSPI or HPAII and each containing the non‐methyl sensitive rare cutting restriction enzyme EcoRI) was performed using 200 ng of DNA, followed by a pre‐selective PCR targeting the ligated ends with an additional base pair, and ending with a final selective PCR with an additional two base pairs added to the primer along with a fluorescent dye (6‐FAM) used to detect fragment length (Table 2). A total of four selective primers in two combinations were used for both batches. Fragment analysis was conducted by the FIU DNA Core on an Applied Biosystems® 3130XL Genetic Analyzer. Raw fragment data were then scored as presence (+)/ absence (−) for all fragments for each individual for each parallel reaction (MSPI and HPAII) to decipher the four types of DNA methylation status present. Presence (+)/absence (−) (MSPI/HPAII) pattern for each fragment gave methylation type (+/+ = unmethylated, +/− = hemi‐methylated, −/+ = internal C methylation, −/− = full methylation). Batch 1 samples were subjected to the MSAP protocol and fragment analysis at a separate and earlier time than Batch 2 samples. Once raw data were available for all samples, the final dataset was scored using all samples to allow for loci to be compared.

TABLE 2

Adaptors and Primers usedfor MSAP analyses. Only selective primers were fluorescently labeled (FAM) for fragment analysis

Step/combo	Oligo name	Sequence (5′ to 3′)
Restriction ligation	EcoRI Adaptor Fwd	CTCGTAGACTGCGTACC
Restriction ligation	EcoRI Adaptor Rv	AATTGGTACGCAGTCTAC
Restriction ligation	HPAII/MSPI Adaptor Fwd	CGTTCTAGACTCATC
Restriction ligation	HPAII/MSPI Adaptor Rv	GACGATGAGTCTAGAA
Preselective combo A	Pre‐EcoRI + A	GACTGCGTACCAATTCA
Preselective combo A	Pre‐MspI‐HpaII + T	GATGAGTCTAGAACGGT
Preselective combo B	Pre‐EcoRI + C	GACTGCGTACCAATTCC
Preselective combo B	Pre‐MspI‐HpaII + A	GATGAGTCTAGAACGGA
Selective Combo 1	selective MspI‐HpaII‐TTG	GATGAGTCTAGAACGGTTG‐FAM
Selective Combo 1	selective MspI‐HpaII‐TCT	GATGAGTCTAGAACGGTCT‐FAM
Selective Combo 2	selective EcoRI‐AAC	GACTGCGTACCAATTCAAC‐FAM
Selective Combo 2	selective MspI‐HpaII‐TCA	GATGAGTCTAGAACGGTCA‐FAM

Data analysis

The msap package (version 1.1.8) in R was used to categorize each locus as either unmethylated, hemi‐methylated, methylated at an internal cytosine, or fully methylated, for each sample (Pérez‐Figueroa, 2013). Although a lack of bands in both reactions could indicate full methylation or the lack of a genetic target, this case was considered as full methylation, consistent with other DNA methylation studies working with species where there is low genetic structure (Feldheim et al., 2001; Rodríguez‐Casariego et al., 2020), as has been shown in lemon sharks (Feldheim et al., 2001). Methylation‐susceptible loci (MSL) were identified as loci with a methylated state in at least 5% of the samples (Pérez‐Figueroa, 2013), and polymorphic loci were identified as those that displayed both methylation and a lack of methylation in at least two instances across samples (Herrera & Bazaga, 2010; Pérez‐Figueroa, 2013). Polymorphic MSL were the focus of all subsequent analyses of DNA methylation. For all distance‐based analyses, Gower distances based on the four possible methylation states (u, h, I, or f) were utilized (daisy function, cluster package, version 2.1.0; Gower, 1971).

Initially, all samples were analyzed together. A Non‐metric Multidimensional Scaling Analysis (NMDS; metaMDS function, vegan package version 2.5–6; Faith et al., 1987) was used to visualize patterns of DNA methylation and potential influencing factors including age group, year, and season of sampling, along with sample batch. Given the potential influence of differences in sample processing between batches, each batch was then analyzed separately. Each batch was re‐analyzed individually with the msap package to characterize methylation states and to identify polymorphic MSL within each batch. Differences in DNA methylation patterns between age groups in each batch were visualized with NMDS analyses and evaluated with permutational analysis of variance (PERMANOVA; adonis function; Anderson, 2001). The homogeneity of dispersions between age groups of each batch was assessed with the betadisper function (Anderson et al., 2006).

Given the larger sample size of Batch 2, differences in DNA methylation patterns between age groups were further assessed within this batch using a Discriminant Analysis of Principal Component (DAPC; adegenet package version 2.1.3; Jombart & Ahmed, 2011). The optimal number of principal components (PCs) to be retained (4 PCs) was determined through the xval‐Dapc function as well as an assessment of a‐scores (optim.a.score function; Jombart et al., 2010). Both discriminant functions were retained. To assess the separation between age groups, the probability of membership of each sample to each age group was evaluated. A distance‐based redundancy analysis (dbRDA; Legendre & Anderson, 1999; Oksanen et al., 2019) was performed with batch 2 samples to characterize the relationship between all the metadata available for the samples (age class, sex, year, season) and DNA methylation variation. The equation was as follows:

A permutation test was performed on the dbRDA model (anova.cca function) to assess the significance of these variables upon the model (Legendre et al., 2011).

RESULTS

Combined batch sample analysis

A preliminary MSAP analysis using Batch 1 samples (n = 10) revealed significant differences between age groups (Analysis of molecular variance (AMOVA), p < .05, Figures S1 and S2). In order to validate these findings, additional samples were incorporated into the analyses (Batch 1 + Batch 2 samples). The combined batch analysis encompassed a total of n = 29 samples (10 juveniles, 12 subadults, and seven adults) used in R analyses, specifically the msap package (Table 1). Two selective primer pairs (Table 2) were used which resulted in a total of 364 loci (133 and 231 for primers 1 and 2, respectively), out of which 346 were identified as MSL of which 76% were polymorphic and 18 as non‐methylated loci (NML) of which 33% were polymorphic using the msap package (additional information from MSAP and other R analyses can be found in the supplemental R Code file). AMOVA analyses indicated significant differences between age groups for MSL (Phi_ST = 0.108, df = 28, p = .003) but not for NML (Phi_ST = 0.005, df = 28, p = .387). The NMDS plot showed juveniles to overall cluster toward the bottom of the plot and adults toward the top with subadults having a large amount of variation; however, there was an observable difference in clustering between samples from Batch 1 and Batch 2 (Figure 1). The potential effect of differences in year of sampling and season was investigated, but neither appeared to have a significant influence in the clustering of data points (Figures S3 and S4). The only obvious difference appeared to be the fact that the samples were processed in two separate batches that were not processed by the same lab for all steps. Therefore, it was decided to investigate age group clustering within each batch effort to eliminate interference from the different sample processing and storage variables. Despite this apparent batch effect, age group differences were readily obvious and appeared robust despite possible differences between batches.

FIGURE 1

NMDS for Combined Samples (Batch 1 & 2). Both NMDS plots relay how similar/different each point is to another based on DNA methylation pattern. Plot A shows age groups for all samples. Plot B shows individual sample IDs where batch 1 (LS1‐10) clearly clusters together (towards right) and away from batch 2 samples (all other samples).

Independent batch sample analysis

The analysis of Batch 1 samples (n = 10) consisted of four juveniles, three subadults, and three adults. A total of 308 loci were identified as MSL (45% of them polymorphic) and 56 to as NML (0% polymorphic) using the msap package (Figures S5 and S6). AMOVA analyses revealed a significant difference between age groups for MSL (Phi_ST = 0.111, df = 9, p = .049). Accordingly, juveniles had less unmethylated loci compared to subadult and adult groups (Figure 2). The NMDS plot showed that juveniles cluster more closely together than individuals in subadult or adult groups (Figure 3). Similarly, PERMANOVA analyses found a significant difference between age groups (p = .05) and a marginally significant difference in dispersion between groups (p = .09). Such age group differences were similar to those observed in the combined batch sample analysis, with samples from juvenile sharks showing a consistent clustering separate from other age groups, and with subadults displaying the most variation among data points.

FIGURE 2

Percent DNA methylation per age Group for Batch 1 samples. Plot (a) depicts % loci as methylated or non‐methylated. Plot (b) depicts % loci for each of the four types of methylation status determined by MSAP analysis (none‐ non‐methylated, hemi‐hemi‐methylated, Internal C‐ internal cytosine methylated, full‐full/hypermethylated target site).

FIGURE 3

NMDS plot for batch 1. Distance between any two dots depicts the difference/similarity between two given samples

The analysis of Batch 2 samples consisted of n = 19 samples, including six juveniles, nine subadults, and four adults. This analysis revealed the presence of 341 MSL (62% polymorphic) and 23 NML (26% polymorphic). MSAP Package PCoA graph and NMDS plot for Batch 2 samples are in Figures S7 and S8. The AMOVA analysis of MSL data found a significant difference between age groups (Phi_ST = 0.308, df = 18, p < .0001), whereas AMOVA of NML did not (Phi_ST = 0.023, df = 18, p = .311). Again, juveniles displayed less unmethylated loci compared to subadults and adult groups, and a large portion of loci was found to be fully methylated in juveniles for this batch effort (Figure 4). Similar to the combined analysis and Batch 1 analysis, the NMDS plot for batch 2 showed samples from juvenile sharks clustering independently from subadults and adults (Figure 5). PERMANOVA analysis identified significant differences between age groups and dispersion between age groups (p = .001 and p = .007, respectively). To further analyze differences between age groups, a DAPC analysis was conducted on Batch 2 samples (feasible in this particular case based on its larger sample size) with discriminant function 1 accounting for the majority (~50%) of differentiation observed between juveniles and other age groups (Figure 6). Lemon shark age groups (i.e., juvenile, subadult, and adult) showed differential DNA methylation patterns (Figures 5 and 6), with the juvenile age group showing the most differentiated pattern from the other groups and with 100% accurate membership predictability compared with individuals from the other two age groups (Figure 7). These results are consistent with a progressive change in DNA methylation pattern from juvenile to adult individuals (Figure 6b). A dbRDA was conducted using all available variables (age class, season, sex, and year) to check for significant relationships with DNA methylation. The only variation found to significantly explain variance in the model was age class (p = .001). The other variable p values were season (p = .419), sex (p = .169), and year (p = .542). Age class along with these other variables was deemed to explain 60.7% of variation in the model; however, the adjusted R 2 for only significant variables (i.e., age class) was 33.3% of variation meaning that 66.7% was unexplained variation.

FIGURE 4

Percent DNA methylation per age Group for Batch 2. Plot (a) depicts % loci as methylated or non‐methylated. Plot (b) depicts % loci for each of the four types of methylation status determined by MSAP analysis (none‐non‐methylated, hemi‐hemi‐methylated, Internal C‐ internal cytosine methylated, full‐full/hypermethylated target site).

FIGURE 5

NMDS plot for batch 2. Distance between any two dots depicts the difference/similarity between two given samples

FIGURE 6

DAPC and density plots for discriminant function 1 and 2 for batch 2. Plot (a), shows spatial and temporal distribution of data points for the three age groups, x‐axis is DF1 and y‐axis is DF2. Plot (b), density plot for DF1 showing the most differentiation of juveniles from other age groups. Plot (c), density plot for DF2. In both B and C plots, we see that subadults have the most variation across both DF1 and DF2.

FIGURE 7

Probability of membership to a given age group. Samples are listed on the x‐axis as sample ID_Letter representative of age group as follows: J‐juveniles, S‐subadult, and a‐adult. Juveniles showed a 100% probability of being assigned to their correct group. Subadults had the most variation (50%–95%) in probability of correct assignment, whereas all adult samples had 60% or greater probability.

DISCUSSION

This study represents the first effort to investigate the relationship between age and global DNA methylation pattern variation in sharks. The obtained results are consistent with the presence of differences in DNA methylation patterns across different age/life stage groups of lemon sharks. Samples were provided in two different batches from different source material (extracted DNA and fin clips in DMSO), leading to a batch effect in data analysis (Figure 1). DNA methylation has been found to be stable at room temperature, and any loss of DNA methylation would be a result from DNA degradation, so it is plausible that there was a slight difference in DNA quality between the batches due to sample holding that is being reflected in the NMDS analysis. Despite this, the sensitivity of the MSAP technique to sample processing observed in the present work constitutes a valuable result of itself and further informs future studies. Additional factors potentially contributing to batch effects were also explored (e.g., sampling seasonality, sex, year sampled); however, their contribution was not significant (Figures S3 and S4; ANOVA p‐values .419, .169, and .542). These results also bear interest for studies relying on samples with different origins.

Beyond the observed batch effect, the differences by age group were visible for combined Batch (1 + 2) analysis (n = 29; Figure 1) as well as when each batch (Batch 1, n = 10, Figure 3; and Batch 2, n = 19, Figure 5) was analyzed individually. The consistency of the age group pattern regardless of the batch analyzed further supports a strong link between age and DNA methylation patterns across age groups. To evaluate the observed age group differences, further analyses were conducted using batch 2 samples only in order to remove any interference from the batch effect (Figures 6 and 7) based on the larger sample size of this batch (n = 19).

The importance of age clustering for clock development

The findings of this study provide the first evidence for the role of epigenetics during aging in elasmobranchs including sharks, a mysterious group known for the long life spans in some of its members. Previous studies on the whale shark reported that global DNA methylation patterns were similar to other vertebrate species (i.e., high levels of methylation at most CpG sites except for those at promoter regions), leading to the conclusion that this epigenetic modification plays a similar role in sharks as that characterized in other vertebrates (Peat et al., 2017). The present study further supports that conclusion, based on the similarity between the observed age‐related changes in DNA methylation patterns in lemon sharks and other vertebrates (Anastasiadi & Piferrer, 2020; Beal et al., 2019; Bors et al., 2021; De Paoli‐Iseppi et al., 2019; Fei et al., 2021; Polanowski et al., 2014; Stubbs et al., 2017). The largest age‐related change in DNA methylation patterns was found between juveniles and adults seen in the DAPC analysis (Figure 6a), with approximately 50% of the differences between groups explained by discriminant function 1 (x‐axis Figure 6a,b). Interestingly, subadult DNA methylation patterns had a large amount of variation that spanned between the juvenile and adult groups (Figure 6b). These results seem to be in line with reports describing gradual overall changes in DNA methylation throughout the genome correlated with chronological aging (Bollati et al., 2009; Ciccarone et al., 2018; Day et al., 2013; Horvath, 2013; Maegawa et al., 2010). Along the continuum of age, hypomethylation has been observed across the genome due to drift and errors in maintenance of DNA methylation on the genome (Berdyshev et al., 1967; Horvath, 2013; Martin‐Herranz et al., 2019; Shimoda et al., 2014; Wilson et al., 1987). However, this is not the case for all sites as revealed by the presence of some sites across the genome displaying age‐driven hypermethylation (Rakyan et al., 2010). Some of these changes are possibly correlated with changes in gene expression during development (Calvanese et al., 2009; López‐Otín et al., 2013; Pal & Tyler, 2016; Unnikrishnan et al., 2019). This bi‐directional change along with the specificity of individual CpG sites correlating to age could explain why subadults, as observed in this study, have a wider variation in DNA methylation points that seems to span between juveniles and adults as they transition between these two age groups. It is also possible that this same bi‐directional change is one reason that the differentiation between subadults and adults is less obvious and may require more higher resolution methods to elucidate age for these groups.

Another interesting pattern observed in this work is the larger amount of variation in DNA methylation patterns for older shark individuals. Many of these DNA methylation positions most likely are driven by the onset of epigenetic modifications linked to changes in gene expression related to aging; however, many may very well be attributed to the higher incidence of epigenetic drift later in life associated with heterogeneous environments (Jung et al., 2017), contributing to the variation observed, especially since these samples come from varying sample years. This effect is best illustrated by human twin studies revealing a variation in epigenomes driven by environmental heterogeneity (Bell & Spector, 2011; Fraga et al., 2005) but has also been captured by epigenetic clock studies (Beal et al., 2019; Robeck et al., 2021; Stubbs et al., 2017; Thompson et al., 2017).

Although the results reported in this study support a place for MSAP in aging studies, the prevalence of such random and environmentally driven DNA methylation changes during aging prevents this technique from producing accurate estimations of chronological age. On the other hand, MSAP proved useful for making general assignments of age group. Furthermore, these findings suggest that the epigenetic aging pattern observed in mammals also extends to sharks, validating the widespread role of DNA methylation in such process.

The manifold implications of the age of the shark

Future development of epigenetic clocks in sharks would benefit shark conservation similar to the benefits that have been found in being able to characterize the direct link between age, environmental stimulus, and age‐related diseases developed in humans (Bektas et al., 2018; Cavalli & Heard, 2019; Dugué et al., 2018; Zhang et al., 2017). Better understanding of this connection in sharks would allow for environmental threats impacting populations to be identified, as well as the specific age groups most vulnerable to the threats, something that is little known for most shark species. Additionally, with little to no age data existing for most shark populations, little demographic and life history information is available, and the status of populations is hard to assess (Carlson & Goldman, 2007; Goldman & Cailliet, 2004). The current methods used to age sharks (band central staining and bomb radiocarbon aging) are not applicable to all species, require sacrificing individuals, and are fraught with high variation (Campana, 2001; Francis et al., 2007; Huveneers et al., 2013; Natanson et al., 2018). Furthermore, being able to estimate shark age from tissue would greatly enhance the conservation work being done in the fin trade which has been identified as one of the top threats to endangered shark species (Walls & Dulvy, 2020). Current molecular methods only allow for the identification of species and possibly the origin population (Cardeñosa et al., 2020, 2021; Clarke et al., 2004, 2006; Fields et al., 2015, 2018) as well as investigations into the quantities of metals and other pollutants (Barcia et al., 2020; Nalluri et al., 2014) these sharks are exposed to but do not allow for the estimation of age.

Epigenetic alteration is one of the hallmarks of aging, and the study of these alterations could lead to a better understanding and ability to combat age‐related diseases (López‐Otín et al., 2013). As such, studying the specific markers that correlate with age for sharks may prove to be beneficial to understanding and learning how to treat age‐related diseases in humans. Many sharks live long lives (Christiansen et al., 2016; Natanson & Skomal, 2015; Nielsen et al., 2016), and one means to their longevity is thought to be that sharks are less susceptible to age‐related diseases. For example, sharks may be resistant to the age‐related cancers due to potential tumor‐suppressing properties found in their cartilage (although research remains inconclusive) (Oliveira et al., 2019; Walsh et al., 2013); however, there remains a lack of systematic surveys and studies to confirm the actual incidence of cancer in sharks in the first place (Ostrander et al., 2004). Findings from the white shark genome do provide evidence that if resistance to cancer does exist for sharks, it may be due to positive selection of genome stabilizing genes which would have a direct effect upon tumor suppression since many tumors form due to breaks and damage to the DNA (Marra et al., 2019). Further study into epigenetic clocks may provide more insights into the longevity many sharks have and could provide insights into human health and aging.

CONCLUSIONS

This study contributes to fill two major gaps in the study of shark aging and epigenetics: sharks do exhibit DNA methylation patterns that change with age, and the MSAP method is efficient in identifying such changes. These findings are important because they support the development of epigenetic aging tools in this group of elasmobranchs, fostering the development of additional resources including the generation of full reference genomes. In addition, this work underscores the role of MSAP as a simple, inexpensive, and widely applicable technique to perform the first round of analyses within more elaborated epigenetic pipelines able to systematically estimate chronological age in samples from diverse origins and age groups. Overall, this study provides a first step toward the development of such pipelines for the shark research community with critical implications in conservation and management (Marra et al., 2019).

AUTHOR CONTRIBUTIONS

Andria Beal: Lead Conceptualization, Lead Data curation, Formal analysis, Lead Funding acquisition, Lead Investigation, Methodology, Project administration, Lead Visualization, Lead Writing – original draft. Serena Hackerott: Lead Formal analysis, Writing – original draft, Writing – review & editing. Kevin Feldheim: Data curation, Funding acquisition, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing. Samuel Gruber: Data curation, Resources. Jose Maria Eirin‐Lopez: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Visualization, Writing – original draft‐, and Writing – review & editing.

CONFLICT OF INTEREST

All authors claim no conflict of interest.

Appendix S1

Click here for additional data file.

Figures S1–S8

Click here for additional data file.