Genetic diversity of two populations of the tufted puffin Fratercula cirrhata (Pallas, 1769)

Copyright: ©2021 The Author(s). This is an open access article distributed under the terms of the Creative Commons BY 4.0 International (Attribution) License (https://creativecommons.org/ licenses/by/4.0/legalcode), which permits the unrestricted distribution, reproduction and use of the article provided the original source and authors are credited. Genetic diversity of two populations of the tufted puffin Fratercula cirrhata (Pallas, 1769)


INTRODUCTION
The tufted puffin Fratercula cirrhata (Pallas, 1769) is a pelagic seabird that belongs to the family Alcidae, in the order Charadriiformes. This species was first described by Peter Simon Pallas in 1769 with the scientific name Alca cirrhata, although it was later renamed Lunda cirrhata (Pallas, 1811). The first to fifth editions of the Check-list of North American Birds used the genus name Lunda. Then, in the 34 th supplement (Eisenmann et al., 1982) of the fifth edition, the genus name Fratercula was adopted, and this use has continued to the current seventh edition (American Ornithologists' Union, 1998).
The genetic diversity of this species has received scant attention in the literature. Only 13 sequences for cytochrome c oxidase subunit I (COX1) are available in the Barcode of Life Data System v4 (BOLD) (Ratnasingham and Hebert, 2007) for genetic diversity analysis. Five of these are for birds from Japan (Saitoh et al., 2015), with four from Alaska (Hebert et al., 2004;Tavares and Baker, 2008), one from Canada (Kerr et al., 2007), two from Russia (Kerr et al., 2009) and one of unknown origin (Pereira and Baker, 2008). There are no clear reginal differences in the COX1 data in the BOLD system.
In this study, to improve our understanding of the genetic diversity of tufted puffins, we analyze the genetic diversity of two captive populations (originally from far eastern Russia and Alaska) using not only COX1 sequences but also the D-loop locus, which is expected to exhibit more variability than COX1 within the mitochondrial DNA (mtDNA). The D-loop locus has been effectively used for population studies in several Charadriiformes bird species (e.g., Liebers et al., 2001Liebers et al., , 2004Liebers and Helbig, 2002;Baker, 2003, 2005;Miller et al., 2010Miller et al., , 2015Sternkopf et al., 2010;Rheindt et al., 2011;Birt et al., 2012;Küpper et al., 2012;Verkuil et al., 2012;Pons et al., 2014;Trimbos et al., 2014;Pshenichnikova et al., 2015;Sonsthagen et al., 2015;Tigano et al., 2015;Wallace et al., 2015;Wang et al., 2019). In addition, we also selected one nuclear locus, Rh family B glycoprotein (RHBG), to sequence. RHBG is a member of the RH gene family, which was duplicated in the common ancestor of vertebrates (e.g., Huang and Peng, 2005;Peng and Huang, 2006;Kitano et al., 2010;Suzuki et al., 2014Suzuki et al., , 2017, during two rounds of wholegenome duplication (Dehal and Boore, 2005). A unique characteristic of the bird lineage is that it exhibits a higher rate of nonsynonymous substitutions in the RHBG genes, probably because the functional constraints on this gene are relaxed (Suzuki et al., 2017). We therefore expected to obtain more variable sequence data in the analysis of this nuclear gene.
In Japan, as of 2020, five aquariums (Aquamarine Fukushima, Kamogawa Seaworld, Osaka Aquarium Kaiyukan, Aqua World Oarai and Tokyo Sea Life Park) hold tufted puffins. The birds at four of these aquariums originally came from far eastern Russia, while those at Tokyo Sea Life Park originally came from Alaska. We believe that the genetic data provided in this study will also be helpful for the management of captive tufted puffins in aquariums.

MATERIALS AND METHODS
DNA extraction Tufted puffins from the Tokyo Sea Life Park (Kasai Rinkai Suizokuen, hereafter "Kasai") and Aqua World Oarai (hereafter "Oarai") were used. The Kasai individuals were originally from Kodiak Island and the Triplet Islands in Alaska, USA, while the Oarai individuals were originally from Kamchatka, in far eastern Russia. Individuals were identified by means of implanted microchips. The family trees of the birds were constructed based on the Japanese studbook for tufted puffins (Figs. 1 and 2). Blood collection was carried out under manual restraint, with sufficient attention given to bird welfare. A blood sample was taken from the medial metatarsal vein of each individual and stored at − 20 °C until DNA extraction. All procedures for collecting samples from birds were performed in accordance with the guidelines of the Japanese Association of Zoos and Aquariums.
Genomic DNA was extracted from the blood samples using the conventional sodium dodecyl sulfate lysis method with proteinase K and phenol-chloroform. The quality of the extracted DNA was confirmed via 0.8% agarose gel electrophoresis, and DNA concentration was measured by a Beckman DU-7500 spectrophotometer (Beckman Coulter, Brea, CA, USA).
Polymerase chain reaction and DNA sequencing PCR was performed in a 5-μl reaction volume containing 20-50 ng of template DNA, 1 × KAPA2G Fast Multiplex PCR Mix (KAPA Biosystems, Wilmington, MA, USA) and 0.2 μM each of the forward and reverse primers. Thermal cycling was performed in a Veriti 200 machine (Applied Biosystems, Waltham, MA, USA). A two-step PCR method was employed: 95 °C for 3 min, followed by 30 cycles of denaturation at 95 °C for 15 s, annealing at 68 °C for 30 s, and extension at 72 °C for 1 min.
Based on the family trees (Figs. 1 and 2), three maternal groups were recognized in Kasai and seven in Oarai. One individual was selected for the sequencing of all three loci from each maternal group in Kasai (K02, K11 and K19) and Oarai (O08, O34, O36, O39, O43, O45 and O48). Since there was no information about the parents of three individuals (K01, K12 and K30) from Kasai ( Fig. 1) and eight individuals (O01, O03, O10, O22, O28, O31, O38 and O42) from Oarai ( Fig. 2), their maternal groups could not be determined. These individuals were therefore also selected for sequencing. For O48, only the mtDNA loci (COX1 and D-loop) were sequenced, because one of its RHBG alleles should be identical to that of its father, O10. In total, six Kasai and 14 Oarai individuals were selected for the sequencing of all three loci. In addition, five individuals (K26, K44, K47, K51 and K57) from Kasai and one (O78) from Oarai were selected for the sequencing of the mtDNA loci only, because their maternal groups were unknown.
Phred/Phrap software (Ewing et al., 1998) was used for base-calling and sequence assembly, and to obtain quality scores for the assembled data. Editing of the assembled sequences was performed using Consed (Gordon et al., 1998). Essentially, both strands of the PCR products were sequenced. If only one strand of a portion of the end region was sequenced, sites with high-quality scores (Phred quality ≥ 30) were used. Two overlapping peaks for one site in a sequence were interpreted as indicating heterozygosity, and the site was designated according to the Nomenclature Committee of the International Union of Biochemistry protocol (R: A or G; Y: C or T; K: G or T; M: A or C; S: C or G; and W: A or T).
Sequence data analysis Multiple alignments were performed using MUSCLE (Edgar, 2004), implemented with MEGA7 software (Kumar et al., 2016). Haplotypes of RHBG were predicted using PHASE 2.1.1 (Stephens  , 2001). Genetic differentiation between two populations was estimated by analysis of molecular variance (AMOVA) using Arlequin 3.5 (Excoffier and Lischer, 2010). Character state-based phylogenetic networks for the haplotypes of each locus were constructed manually following Bandelt (1994) and Saitou and Yamamoto (1997). A phylogenetic tree was constructed using the neighbor-joining method (Saitou and Nei, 1987) with the distances corrected by the Kimura 2-parameter model (Kimura, 1980) in MEGA7. The bootstrap test was carried out with 1,000 replicates (Felsenstein, 1985).

RESULTS AND DISCUSSION
Nucleotide sequences of tufted puffin loci We obtained nucleotide sequences of 468 bp for COX1, 360 bp for D-loop, and 840 bp for RHBG in tufted puffins. For RHBG, a region containing part of exon 5 to part of exon 8 was sequenced. A section of intron 6 of RHBG could not be sequenced due to a poly(C) tract. All nucleotide sequence data obtained in this study were deposited in the DDBJ/EMBL/GenBank database under the accession numbers LC535135-LC535210.
We observed five haplotypes for COX1 in the aquarium birds, in addition to two that were represented only in the BOLD data (Table 1, Fig. 3). In addition to the two populations sequenced in this study, we included sequences from BOLD in the phylogenetic network. In total, we found six segregating sites in the region; all were transitional and synonymous. Haplotype I was the major haplotype, and haplotypes III, IV and V each had one mutation differentiating them from haplotype I. Haplotype II was the second largest in this study, and haplotype VII differed by one mutation from haplotype II. Haplotype VI was located midway between haplotypes I and II.
We observed seventeen D-loop haplotypes in this study (Table 1, Fig. 4). Four individuals (K11, K44, O08 and O31) each had one heterogeneous site: K11 (158Y), K44 (48Y), O08 (86Y) and O31 (164R). To confirm these heterogeneous sites, we sequenced the D-loop of K62, a grandson of K11 (Fig. 1), and of O73, a son of O08 (Fig.  2). We obtained identical sequences, with the same heterogeneous sites, from K62 and K11 and from O73 and O08 (Fig. 4). It is therefore likely that these individuals are heteroplasmic, that is, they carry more than one type of mtDNA within each cell. Although O31 and K44  also had one heterogeneous site each, this heterogeneity could not be confirmed, because O31 has no relatives in the aquarium population and K44 has no relatives on his maternal side. Berg et al. (1995) observed heteroplasmic variation in the D-loop region in four Laridae, ten Alcidae, and one Scolopacidae species. It therefore seems that Charadriiformes species tend to exhibit heteroplasmy in the D-loop region. Berg et al. (1995) observed tandem repeat heteroplasmy, whereas we observed only point mutations in this study; however, we did not sequence the heteroplasmic tandem repeat region, which is located at the 5′ end of the D-loop region. Heteroplasmic sites represented by point mutations have also been observed in the same D-loop region in Cassin's auklet (Ptychoramphus aleuticus; Wallace et al., 2015).
Two large COX1 haplotype groups (I and II) were observed for this species (Fig. 3). D-loop sequences from individuals who had COX1 haplotype II formed a cluster in the D-loop tree (Fig. 4), although the bootstrap value was low (36%). Although bootstrap values for all nodes are low, compatibility between the COX1 and D-loop sequences shows the reliability of our DNA sequencing, since mtDNA is haploid and consistent haplotype patterns are expected from anywhere in the mitochondrial genome. In total, 10 segregating sites of the D-loop sequences were observed (Fig. 5). Of these, the 125 th site had three character states (A: 1-2-3-4-7-8-9-14-15-16-17; T: 5-11-12; C: 6-10-13). There was therefore a total of 11 mutations. The changes between A and T and between A and C at the 125 th site were transversions. All other changes were transitions.
The predicted haplotypes of RHBG are shown in Supplementary Table S1. The most likely haplotype pairs for each individual were used for the subsequent analysis (Table 1, Fig. 6). In total, 18 segregating sites were observed. Sites 34, 35 and 64 were in exon 5; sites 94 and 182 were in intron 5; sites 357, 371 and 372 were in exon 6; site 441 was in intron 6; sites 599, 632, 674, 675 and 698 were in exon 7; and sites 766, 767, 820 and 904 were in intron 7. Four nonsynonymous (34, 64, 371 and 675) and seven synonymous (35,357,372,599,632,674 and 698) changes were observed. The changes at sites 94 and 182 were transversions and all others were transitions.   Table 1). Fig. 6. Haplotype network for a section of the RHBG sequence obtained for tufted puffins. The numbers on each edge are the positions of the mutations responsible for those edges; boldface denotes transversions, underlining denotes nonsynonymous changes, and italicization denotes substitutions in introns. The circle at each node (a-n) indicates the proportion of occurrences of that haplotype observed in the Kasai (black) and Oarai (gray) populations, and the size of the circle indicates the frequency of each haplotype (see Table 1).
With respect to polymorphism, there were no large differences in θ W or π values between the two populations for the two mtDNA loci (COX1 and D-loop; Table 2). The RHBG sequences from Kasai (θ W = 0.0059, π = 0.0059) showed slightly higher nucleotide diversity than those from Oarai (θ W = 0.0040, π = 0.0030). The levels of polymorphism within the D-loop sequence tended to be larger than those for COX1 or RHBG. Nucleotide diversities of D-loop of several species of Charadriiformes were observed from the lower 0.00014 of Larus cachinnans mongolicus (excluding introgressed haplotypes) (Liebers et al., 2001) to the upper 0.019 of Aethia cristatella (St.-Jonah population) (Pshenichnikova et al., 2015) (Supplementary Table S2). Therefore, the nucleotide diversity of tufted puffins (Kasai: 0.0082, Oarai: 0.0080) was considered to be moderate in Charadriiformes bird species.
No clear genetic differences between the two tufted puffin populations We observed no clear genetic differences between the Kasai (originally from Alaska) and Oarai (originally from Kamchatka, far eastern Russia) populations in the three loci we analyzed. The major haplotype of COX1, haplotype I, consists of 10 individuals (six Kasai and four BOLD) from North American birds and 12 (nine Oarai and three BOLD) from East Asian birds. The second largest haplotype, II, consists of five individuals (four Kasai and one BOLD) from North American birds and nine (five Oarai and four BOLD) from East Asian birds (Fig. 3). In the phylogenetic network of the D-loop locus (Fig. 5), the largest haplotype, '4', consists of three Kasai and three Oarai sequences, and the other major haplotype groups, including '3' (two Kasai and two Oarai), '5' (one Kasai and three Oarai) and '7' (two Kasai and two Oarai), also consist of individuals from both populations. For RHBG (Fig. 6), the largest haplotype, 'a', consists of three Kasai and 11 Oarai sequences, and the second largest, 'l', consists of two Kasai and three Oarai sequences. We also performed AMOVA within and between the two populations (Table  3). Although P-values are not statistically significant at the 1% level, low fixation indices were observed for three loci ( − 0.0623 for COX1, − 0.0128 for D-loop and 0.1072 for RHBG). Thus, it is clear that the two populations have not diverged genetically. Saitoh et al. (2015) analyzed COX1 from more than 200 bird species around the Japanese archipelago and showed that the average genetic distance within species (estimated using the Kimura 2-parameter method (Kimura, 1980)) was 0.49% (range: 0-6.13%). They showed that the genetic distance for the tufted puffin was on average 0.21%, with a maximum of 0.45%, from seven individuals (deposited in the BOLD system). Based on the merged Kasai and Oarai data in our study (27 individuals), the estimated genetic distance was on average 0.25% with a maximum of 0.65%. These values are comparable with the values found by Saitoh et al. (2015) and fit within the range of within-species genetic variability of Japanese birds.
Although tufted puffins are known to breed throughout the boreal and low Arctic, from California to Hokkaido, their range during the feeding season is uncertain. They are thought to generally feed offshore throughout the year, and often far from the continental shelf in winter, and lay one egg per pair per year (Gaston and Jones, 1998). Although adults return to breeding colonies en masse in spring (Wehle, 1980), it is not known whether all adults return to their natal colony. It is likely that pairbonded birds occupy the same burrows year after year, but information on nest-site fidelity and the duration of pair-bonds is scarce (Piatt and Kitaysky, 2002). One observational study estimated that a minimum of 29% of banded tufted puffins returned to monitored nest sites on Ugaiushak Island, Alaska (Wehle, 1980). Our molecular data suggest that the tufted puffin has not diverged genetically within its breeding range.
According to The IUCN Red List of Threatened Species (https://dx.doi.org/10.2305/IUCN.UK.2018-2.RLTS. T22694934A132582357.en), the worldwide population of tufted puffins is extremely large, and their Red List status is Least Concern. However, the current population trend is rated as "decreasing". In Japan, which represents the southern limit of the breeding range, tufted puffin numbers have declined steeply (Senzaki et al., 2020). In 1959and 1960, about 250 individuals for each year were observed on Moyururi island, Hokkaido (Fujimaki, 1961). In 1972, 111 individuals were observed on Moyururi and neighboring Yururi islands (Japan Environmental Agency, 1973), whereas Haga (1973) estimated a total of 61 individuals on Moyururi and Yururi islands. In any case, until the 1970s, tufted puffins were observed breeding not only on Moyururi and Yururi islands but also in other regions of Japan (summarized in Osa and Watanuki (2002)). More recently, however, only about 15 breeding pairs were observed, and only on Yururi and Moyururi islands (Osa and Watanuki, 2002). In 2019, unusually high mortality of this species in the eastern Bering Sea was reported (Jones et al., 2019). This may have been the result of starvation caused by a warming climate. The present study has provided genetic diversity data for three loci (COX1, D-loop and RHBG) of the tufted puffin. Although the number of populations analyzed in the study was limited, we established that tufted puffins maintain a moderate level of genetic diversity by comparisons of nucleotide diversities in the D-loop region for several species of Charadriiformes. These genetic markers will facilitate further assessment of tufted puffin population dynamics.
We thank Ms. Yoko Osuka for giving us the opportunity to conduct this study. We are grateful to Mr. Sotaro Kawakami for assisting with the blood collection at Tokyo Sea Life Park. We thank two anonymous reviewers and the editor for their valuable comments.