Intraspecific variation in the mitochondrial genome among local populations of Medaka Oryzias latipes

Gene 457 (2010) 13–24

Contents lists available at ScienceDirect

Gene j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / g e n e

Intraspecific variation in the mitochondrial genome among local populations of Medaka Oryzias latipes Makoto Hirayama a,1, Takahiko Mukai a,2, Masaki Miya b, Yasuhiko Murata a, Yoshio Sekiya a, Toshikazu Yamashita a, Mutsumi Nishida c, Shugo Watabe d, Shoji Oda a, Hiroshi Mitani a,⁎ a

Department of Integrated Biosciences, Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8562, Japan Department of Zoology, Natural History Museum and Institute, 955-2 Aoba-cho, Chuo, Chiba 260-8682, Japan Ocean Research Institute, The University of Tokyo, 1-15-1 Minamidai, Nakano, Tokyo 164-8639, Japan d Laboratory of Aquatic Molecular Biology and Biotechnology, Graduate School of Agricultural and Life Sciences, The University of Tokyo, 1-1-1 Yayoi, Bunkyo, Tokyo 113-8657, Japan b c

a r t i c l e

i n f o

Article history: Received 30 September 2009 Received in revised form 19 February 2010 Accepted 23 February 2010 Available online 26 February 2010 Received by M. Schartl Keywords: Mitochondrial DNA Control region Tandem repeated sequence Temperature adaptation

a b s t r a c t The draft genome data of Medaka Oryzias latipes shows that it has distinct intraspecific genetic variation. To survey the genetic variations contributing to environmental adaptation, we focused on the mitochondrial DNA (mtDNA). The complete mtDNA sequences of Medaka were compared among 8 local population stocks and 4 inbred strains established from genetically divergent groups. Inbred strain HSOK, derived from the Eastern Korean group of Medaka, has a mitochondrial gene order that was distinct from other Medaka groups. Phylogenetic trees based on the mitochondrial genome sequences indicated that the mitogenome from the Shanghai stock (China) and HSOK strain were highly diverged from Japanese Medaka, and that the Japanese Medaka mitogenome was diverged into two groups; this result was fully consistent with those of the previous study using mtDNA-encode gene sequences. Among tRNA genes, the most divergent was the tRNAThr gene as reported in humans previously. The number of tandemly repeated 11 nucleotide units in the Medaka mtDNA control region (CR) varied greatly among local populations. The number of repeats was more variable in the Northern Japanese group (10–34) than in the Southern group (7–12), while two other Oryzias species, inhabiting tropical regions, had no repeats. A comprehensive comparison between the number of repeat units and meteorological data indicated that the number of repeats correlated to the index data of a cold environment and seasonal climatic change. In cold (5 °C) acclimated fish, the mRNA levels varied among mitochondria coding genes. mRNA of the cytochrome oxidase subunit I gene in some local stocks was induced by cold temperature and seemed to be correlated with the number of repeated sequences in the CR. This study revealed that the repeated sequences in the mtDNA CR might function for mtDNA gene expression and that the number of tandem repeats in Medaka mtDNA is likely related to adaptation to a harsh habitat. © 2010 Elsevier B.V. All rights reserved.

1. Introduction

Abbreviations: ATPX, ATPase subunit X; COX, cytochrome oxidase subunit X; CR, control region; CSB, conserved sequence block; Cyt b, cytochrome b; dN, nonsynonymous substitutions per nonsynonymous site; dS, synonymous substitutions per synonymous site; K2P, Kimura's two parameter; NBRP, National Biological Resource Project; NDX, NADH dehydrogenase subunit X; OH, origin of heavy-strand replication; OL, origin of light-strand replication; PH, heavy-strand promoter; PL, light-strand promoter; qPCR, quantitative polymerase chain reaction; qRT-PCR, quantitative reverse transcription-polymerase chain reaction; SNP, single nucleotide polymorphism; TA, air temperature; TAS, termination-associated sequences. ⁎ Corresponding author. Tel.: +81 4 7136 3670; fax: +81 4 7136 3669. E-mail address: [email protected] (H. Mitani). 1 Present address: Graduate School of Biosphere Science, Hiroshima University, 1-4-4 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8528, Japan. 2 Present address: Faculty of Regional Studies, Gifu University, 1-1 Yanagido, Gifu 501-1193, Japan. 0378-1119/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.gene.2010.02.012

Inter- and intraspecific mitochondrial DNA (mtDNA) variation is a very useful marker to study the evolution of animals, especially in taxonomy, systematics, ecology and population biology (Avise, 1994; Avise et al., 1987; Moritz et al., 1987). mtDNA analyses based on partial nucleotide sequences usually assume that most of the variations conform to a neutral model of molecular evolution, but some experimental studies and human mitogenome data have suggested that selective forces might act upon mtDNA (Ruiz-Pesini et al., 2004; William et al., 1995). Vertebrate mtDNA codes for 13 proteins, 2 rRNA and 22 tRNA genes and includes one non-coding control region (CR). Mitochondria serve a critical function in the maintenance of cellular energy stores, thermogenesis and apoptosis (Kelly and Scarpulla, 2004). Human mtDNA variants seem to be associated with some adaptive functions such as longevity (De Benedictis et al., 1999; Ross et al., 2001; Tanaka et al., 1998) and climatic selection (Mishmar et al., 2003).

14

M. Hirayama et al. / Gene 457 (2010) 13–24

Temperature adaptation mechanisms through mtDNA-encoded genes have been analyzed in ectotherms (Doiron et al., 2002) and endotherms (Hittel and Storey, 2002; Mishmar et al., 2003). In fish, during phenotypic cold acclimation the compensatory mechanisms in mitochondrial volume density and properties are different among species (Guderley, 2004). Moreover, adjustments in the mitochondrial properties and capacities seem to be crucial in acclimating to seasonal cold and in cold adaptation of fish (Itoi et al., 2003; Lucassen et al., 2003). If selective pressures affect the survival rate of mtDNA lineages, intra-species mitogenomic differentiation of wild animals are expected to exist. Medaka (Oryzias latipes) is a small freshwater fish distributed in East Asia that shows apparent genetic differentiation among local populations, revealed by the phylogenetic analysis using allozymes (Sakaizumi et al., 1983) and the partial mtDNA sequences (Matsuda et al., 1997; Takehana et al., 2003, 2004). The populations of Medaka have been classified into 4 genetically divergent groups: the Northern Japanese, the Southern Japanese, the Eastern Korean and the China-Western Korean groups. The Northern Japanese and the Southern Japanese groups are estimated to have diverged approximately 4–18 million years ago, and the genome-wide single nucleotide polymorphism (SNP) rate (3.42%) between the two inbred strains, HNI (Northern Japanese group) and Hd-rR (Southern Japanese group), is the highest SNP rate observed in any vertebrate species reported (Kasahara et al., 2007; Setiamarga et al., 2009). Despite the accumulated genetic variation, however, these groups can mate and produce healthy, fertile offspring. The geographical distribution area of Medaka is wide, ranging from the subarctic (ex. northernmost Japan) to the subtropical zones (southernmost Japan and China). Several studies have suggested that climate may create genetic difference among populations of this species. Iwamatsu (1997) reported that the range of the highest survival temperature in Medaka and related species is narrow (from 41 to 42 °C), whereas the lowest survival temperature is considerably broader (from 0 to 10 °C). Hirayama et al. (2006) demonstrated that at low temperatures, inter-group variations in the cell proliferation rate also exist. It was shown that cell lines derived from the Northern Japanese and Eastern Korean groups could proliferate at lower temperatures than those from the Southern Japanese group and the Medaka-related species O. celebensis, which inhabits a tropical zone. In addition, the inter- and intra-population variations in thermal reaction norms for the growth rate of the Northern Japanese group at different temperatures were determined, including the individual growth rates during the juvenile stage within 12 latitudinal populations of the Northern Japanese group in a temperature controlled (28 °C) environment (Yamahira et al., 2007; Yamahira and Takeshi, 2008). The results suggested that natural selection in high latitudes favors individuals that grow faster within a shorter growing season to individuals that have longer growing seasons. Moreover, Matsumoto et al. (2009) indicated that some of the SNPs for several genes with apparent orthology between Medaka and humans show signals of positive selection in Medaka. Thus, Medaka (Oryzias species) is an ideal model animal to how higher latitude poikilotherms have acquired a greater capacity for temperature acclimation at a genetic level. In this study, we compared 12 complete and nearly complete mtDNA sequences of Medaka from different local populations. The molecular evolutionary rate, phylogenetic utility, and nonsynonymous/synonymous substitution rates of each mitochondrial gene were assessed. We also compared CR of 36 Medakas, including 30 local populations, 4 inbred strains and 2 other Oryzias species, to find any features of environmental adaptation in Medaka populations.

in Fig. 1 and Table 1. The local stocks, which were collected from local sites have been maintained for many generations as closed colonies in Graduate School of Frontier Sciences, the University of Tokyo (Mitani et al., 2006; Shima et al., 1985). O. curvinotus and O. luzonensis were obtained from the National Biological Resource Project (NBRP). A portion of the epaxial musculature was excised from fresh specimens and immediately preserved in 99.5% ethanol. In addition to these samples, we used an already published mtDNA sequence of Medaka (DDBJ accession number AP004421, Miya et al., 2003), which was based on specimens purchased from an aquarium retailer in Tokyo, Japan. 2.2. Cold temperature acclimation of fish Adult fish of 4 local population stocks consisting of Hirosaki and Kaga stocks from the Northern Japanese group, and Ichinoseki and Nago stocks from the Southern Japanese group were acclimated at 25 °C for a month. For cooling, four fish from each local stock were subjected to a stepped cooling regime of 1 °C/h, with a maximum decrease in temperature of 3 °C/day over a 7 day period, until a final temperature of 5 °C was reached. The fish were then maintained at 5 °C for 3 weeks. Four control fish from each stock were maintained at 25 °C. Fish were fed ad libitum with brine shrimp (Artemia sp.) and commercial pellet (crushed TetraMin) and kept under the artificial conditions of 14 h:10 h light:dark. Epaxial musculature from each individual acclimated at 5 °C and 25 °C was sampled and immediately frozen with liquid nitrogen. The muscle samples were stored at −80 °C until used for DNA and RNA extractions. 2.3. DNA and RNA extractions and first strand cDNA synthesis Total genomic DNA was extracted from the muscle samples using the DNeasy Tissue Kit (Qiagen, Tokyo, Japan) following the manufacturer's protocol. Medaka musculatures from the fish acclimated to 5 and 25 °C were lysed with ISOGEN (Nippon Gene, Tokyo, Japan) and total RNA was extracted from the lysates following the manufacturer's instructions. The quality of the extracted RNA was verified by electrophoresis using denatured agarose gels (Skrypina et al., 2003). The RNA was then subjected to genome DNA degradation and RNA cleanup using an RNase-free DNase kit and an RNeasy MinElute Cleanup kit (Qiagen), respectively, according to the manufacturer's procedures. The synthesis of first strand cDNA was performed using a ReverTra Aceα-system (Toyobo, Osaka, Japan) following the manufacturer's instructions. 2.4. Long polymerase chain reaction for amplification of entire mitochondrial DNA Using a long PCR technique, two universal fish PCR primer pairs (L2508-16S and H1065-12S, and L12321-Leu and S-LA-16SH; for locations and sequences of these primers, see Miya and Nishida, 2000) were used to amplify the entire mitochondrial genome in a single reaction. Long PCR for 11 stocks and strains was performed using TaKaRa LA Taq polymerase (TaKaRa, Shiga, Japan) following the manufacturer's protocol. The long PCR products were diluted with deionized distilled water (1:10–20) for subsequent use as PCR templates. 2.5. Sequencing of mitochondrial DNA

2. Materials and methods 2.1. Medaka fish The local Medaka (O. latipes, Japanese killifish) populations, inbred strains and Medaka-related species used in this study are summarized

A total of 110 universal fish PCR primers (identical to those used in Miya et al., 2003) were used in various combinations to amplify contiguous, overlapping segments of the entire mitochondrial genome. Twenty-two Medaka-specific primers were newly designed in cases where appropriate universal fish primers were unavailable

M. Hirayama et al. / Gene 457 (2010) 13–24

15

Fig. 1. Collection sites of isolated wild stocks of local populations of Medaka (Oryzias latipes). Circles, Northern Japanese group; rectangles, Southern Japanese group; triangle, Eastern Korean group; inverted triangle, China-Western Korean group; shaded symbols, local stocks and inbred strains whose whole mtDNAs were sequenced in this study; 1, Kamikita; 2, Hirosaki; 3, Odate; 4, Yokote; 5, Tsuruoka; 6, Yamagata; 7, Ryotsu; 8, Niigata and HNI inbred strain; 9, Kosugi; 10, Kaga and KAGA inbred strain; 11, Sabae; 12, Obama; 13, Miyazu; 14, Maizuru; 15, Kinosaki; 16, Toyooka; 17, Ichinoseki; 18, Kesennuma; 19, Iwaki; 20, Tokyo; 21, Amino; 22, Kumihama; 23, Hamasaka; 24, Kasumi; 25, Oku; 26, Hagi; 27, Matsuyama; 28, Nago; 29, HSOK inbred strain; 30, Shanghai. The location IDs are in common with Table 1. The original localities of Hd-rR and AP004421 (Miya et al., 2003) were unknown. O. curvinotus and O. luzonensis which inhabit tropical region were provided by the National Biological Resource Project in Shinshu University.

(see Supplementary Table 1). TaKaRa ExTaq DNA polymerase (TaKaRa) was used for PCR. Thermal cycling conditions consisted of an initial step at 95 °C for 3 min, followed by 30 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and polymerization at 70 °C for 1 min, with a final extension step at 70 °C for 5 min. The amplified products were treated by ExoSAP-IT kit (GE Healthcare BioSciences KK, Tokyo, Japan). Purified PCR products were subjected to direct cycle sequencing with the BigDye Terminator ver. 3.1 Cycle Sequencing Kit (Applied Biosystems, Tokyo, Japan) using the same primers as those for PCR. All labeling reactions were performed according to the manufacturer's instructions. Labeled fragments were analyzed on an ABI Prism 3100 Genetic Analyzer (Applied Biosystems). All DNA sequence electropherograms were carefully examined for the presence or absence of mitochondrial pseudogenes in the nuclear genome (Mindell et al., 1999). The CR of 22 local stocks consisting of 13 and 9 populations from the Northern and Southern Japanese groups, respectively, and the two Medaka-related species O. curvinotus and O. luzonensis were sequenced. The DNA sequencing methods were identical to those mentioned above, and extracted total DNA was used directly as a template for PCR. The primer pair (L15770 and H16484) was used for the PCR of the Northern Japanese Medaka and the pair (L15765 and H16484) was used for all other fish (see Supplementary Table 1). 2.6. Phylogenetic analyses of whole mitochondrial genome sequences Mitogenome sequences were subjected to multiple alignment using MAFFT ver. 6.707b (Katoh and Toh, 2008). The CR sequences were

entirely excluded from the aligned sequences because of ambiguous positional homology and the resulting dataset comprised of 15,723 positions. Unambiguously aligned sequences were subjected to maximumlikelihood (ML) analysis using RAxML ver. 7.2.0 (Stamatakis, 2006) with the two continental populations (HSOK and Shanghai) used as collective outgroups to root the tree. A general time reversible (GTR) model with sites following a discrete gamma distribution (GTR + Γ) was used and a rapid bootstrap (BS) analysis was conducted with 1000 replications (–f a option). This option performs BS analysis using GTRCAT, which is a GTR approximation with optimization of individual per site substitution rates, and classification of those individual rates into a certain number of rate categories. After implementing the BS analysis, the program uses every fifth BS tree as a starting point to search for the ML tree using the GTR + Γ model of sequence evolution to obtain more stable likelihood values. Finally, BS probabilities for the best-scoring ML tree were calculated and written on the tree file. Uncorrected p-distances and genetic distances were calculated for the sequences of 13 protein-coding genes, two rRNA and 22 tRNA genes, and for the combined data set of 22 tRNA genes, consisting of very short sequences (65–72 bp), and putative amino acid sequences of 13 protein-coding genes based on Kimura's two parameter (K2P) model (Kimura, 1980) by using PAUP*4.0b10 (Swofford, 2002). To assess the natural selection on each lineage of each gene, nonsynonymous substitutions per nonsynonymous site (dN) and synonymous substitutions per synonymous site (dS) were estimated by PAML ver. 3.13a (Yang, 1997) using the phylogenetic tree constructed based on the total mtDNA sequences.

16

M. Hirayama et al. / Gene 457 (2010) 13–24

Table 1 List of Medaka populations used in this study, indicating the location, number of tandem repeats and latitude. Group

Location (prefecture or country)

Location IDa

Wild stock (species)

No. of repeats

Latitude

Northern Japan

Aomori

Okayama Yamaguchi Ehime Okinawa – – Korea China

1 2 3 4 5 6 7 8 8 9 10 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 – – 29 30

Kamikita Hirosaki Odate Yokote Tsuruoka Yamagata Ryotsu Niigata HNI (inbred) Kosugi Kaga KAGA (inbred) Sabae Obama Miyazu Maizuru Kinosaki Toyooka Ichinoseki Kesennuma Iwaki Tokyo Amino Kumihama Hamasaka Kasumi Oku Hagi Matsuyama Nago Hd-rR (inbred) AP002241 HSOK (inbred) Shanghai

16 18 34 28 21 28 17 15 13 16 12 13 16 14 19 19 10 12 12 9 9 9 7 7 7 10 10 8 9 7 8 8 7 11

40.49 40.37 40.17 39.19 38.44 38.15 38.04 37.55 37.55 36.45 36.16 36.16 35.55 35.29 35.33 35.27 35.38 35.32 38.56 38.55 37.03 35.44 35.43 35.35 35.37 35.38 34.40 34.25 33.51 26.36 – – 38.14 31.12

China Philippines

– –

O. curvinotus O. luzonensis

Akita Yamagata Niigata

Toyama Ishikawa Fukui Kyoto Hyogo Southern Japan

Iwate Miyagi Fukushima Tokyo Kyoto Hyogo

Eastern Korea China-Western Korea Tropical region a

0 0

– –

when the amount of precipitation was N1 mm, N10 mm, N30 mm, N50 mm, N70 mm or N100 mm. 6) Adaptation to continuous low temperature; amount of snowfall (also in adaptation to immediate temperature change), maximum amount of snowfall per day (also in adaptation to immediate temperature change), maximum amount of snow cover, annual number of days when the amount of snow cover was N5 cm, N10 cm, N20 cm, N50 cm or N100 cm. 2.8. Quantitative real-time PCR analysis To elucidate whether any correlations existed between the number of tandem repeats and the amount of mtDNA or expression levels of mtDNA genes, real-time quantitative PCR (qPCR) and reverse transcription (RT)-PCR (qRT-PCR) using SYBR Premix Ex Taq (TaKaRa) were performed for adult Medaka fish kept at 5 and 25 °C with a Smart Cycler II system (Cepheid, Sunnyvale, CA) following the manufacturer's instructions. Specificity of the primers for the genes to be amplified was examined using melting curves following the manufacturer's instructions. Primer pairs for the amplification of DNA and cDNA coding NADH dehydrogenase 1 (ND1), cytochrome oxidase subunit I (COI) and ATPase subunits 6 (ATP6), which are components of the respiratory chain complex I, IV and V, respectively, encoded by mtDNA and the βactin gene, encoded by nuclear DNA, were used (Supplementary Table 1). The relative amount of mtDNA to nuclear DNA was quantified by normalizing the values of mtDNA-encoded genes (ND1, COI, and ATP6) with those of nuclear DNA-encoded gene (β-actin) in quantitative PCR using total genomic DNA as template. The relative mRNA levels of mtDNA-encoded genes were quantified by normalizing the values of ND1, COI and ATP6 with those of β-actin in quantitative RT-PCR using first strand cDNAs as templates. Student's t-test was employed for statistical comparison between the levels at 5 and 25 °C. 3. Results 3.1. Mitochondrial genome organization

See Fig. 1 for the sampling location.

2.7. Meteorological data To examine the effect of geographical variation and climate difference on the number of the tandem repeats in the mtDNA CR of the Japanese Medaka populations, the meteorological data of each locality from 1976 to 2007 (cited from the database of the Japan Meteorological Agency (http://www.jma.go.jp/)) were averaged and compared (Supplementary Table 2). When the data from the original location was incomplete or absent, we used the data from the nearest city in its place. The following meteorological data were used as reference indices: 1) General adaptation to temperature; average air temperature (TA) 2) Adaptation to high temperature; mean highest TA, maximum highest TA, annual number of days of TA N25 °C, annual number of days of lowest TA N25 °C, annual number of days of highest TA N25 °C, N30 °C, or N35 °C 3) Adaptation to low temperature; mean lowest TA, minimum lowest TA, annual number of days of average TA b0 °C, annual number of days of lowest TA b0 °C, annual number of days of highest TA b0 °C 4) Adaptation to seasonal temperature change; difference between mean highest and lowest TAs, difference between maximum highest and minimum lowest TAs 5) Adaptation to immediate temperature change; amount of precipitation, amount of snowfall, maximum amount of precipitation per day, maximum amount of snowfall per day, annual number of days

The newly determined complete light (L)-strand nucleotide sequences of the Niigata stock, HNI, KAGA and Hd-rR strains, and the nearly complete sequences (except for a portion of the CR) of the other 6 stocks and the HSOK strain have been registered in DDBJ/EMBL/ GenBank under the accession nos. AP008938–AP008948. The genomic content of 11 of 12 Medaka mitogenomes including AP004421, contained two rRNA, 22 tRNA, and 13 protein-coding genes and the putative CR (Fig. 2). The inbred strain HSOK, established from the East Korean group, had an additional two tRNA-like, one ND4L-like, and two unidentified sequences between the COIII and ND4L genes. 3.2. Sequence divergence and phylogenetic relationships in mitochondrial genes The sequence divergence and phylogenetic relationships of Medaka mtDNA were estimated by the ML method using total mtDNA sequences but excluded the CR and pseudogenes (Fig. 3). The mitogenome of the Shanghai stock (China) and HSOK strain (Eastern Korea) were highly divergent from that of the Japanese Medaka, and also from each other. The Japanese Medaka mitogenome had diverged into two groups: the Northern Japanese group (Hirosaki, HNI, Niigata, KAGA and Toyooka) and the Southern Japanese group (Amino, Hd-rR, Matsuyama, AP004421 and Nago). The most similar sequences were shared between the HNI strain and Niigata stock, which is not surprising as the inbred strain HNI was established from the local stock of Niigata (Hyodo-Taguchi and Sakaizumi, 1993). However, there were four single nucleotide changes in the COI, cytochrome b (Cyt b), ND4 and tRNAGlu genes.

M. Hirayama et al. / Gene 457 (2010) 13–24

17

Fig. 2. Gene organization of Medaka mitochondrial genome (excluding HSOK mtDNA) (A) and a hypothetical model of gene rearrangement of the COIII–ND4 region in HSOK mtDNA (B, C). All protein-coding genes are encoded by the heavy (H)-strand with the exception of ND6, which is coded by the light (L)-strand. Transfer RNA genes are designated by singleletter amino acid codes and those encoded by the H- and L-strands are shown outside and inside the gene map, respectively. 12S and 16S indicate genes of the 12S and 16S ribosomal RNAs; ND1-6 and ND4L, NADH dehydrogenase subunits 1–6 and 4L; COI–III, cytochrome oxidase subunits I–III; ATP6 and 8, ATPase subunits 6 and 8; Cyt b, cytochrome b; CR, control region. OH, origin of H-strand replication; OL, origin of L-strand replication.

3.3. Molecular evolutionary rate of mitochondrial genes and control region Of the 13 protein-coding genes, COI and ND5 genes showed length variation (Table 2). In the COI gene, three of the Northern Japanese groups (Hirosaki, HNI and Niigata) shared the same length polymor-

Fig. 3. Maximum-likelihood tree derived from RAxML analysis of unambiguously aligned mitogenome sequences from the 12 Medaka populations (8 local and 4 inbred) (15,723 positions excluding the control region). Numerals beside the internal branches are bootstrap probabilities based on 1000 replicates. Scale indicates substitution per site.

phism, which was due to a single nucleotide deletion at the 1537th position. This deletion caused a frame-shift in the original stop codon and extended the gene length from 1557 to 1566 bp. In the ND5 gene, the HSOK strain lacks a single codon, encoding for the 611th amino acid residue (1831–1833th nucleotide position). The magnitude of sequence difference for each gene also varied. The maximum values of the uncorrected pairwise sequence differences were 11.7 (COII) to 24.9% (ND3) in 13 protein-coding genes, and 5.3, 7.7 and 6.0% in 12S rRNA, 16S rRNA and combined tRNA gene sequences, respectively (Table 2). In comparing the pairwise genetic distances (K2P model) of each gene to that of the total mitogenomic sequences (excluding CR), the relative evolutionary rates were 0.85 ± 0.13 (average ± S.D., COII) to 1.54± 0.17 (ND2) in 13 protein-coding genes (Table 2). The ATPase subunit and cytochrome oxidase subunit genes showed very similar evolutionary rates (0.85–1.04), while the NADH dehydrogenase subunit genes evolved faster (1.11–1.54). The relative evolutionary rates of rRNA and tRNA genes were slower, 0.35 ± 0.11 (12S), 0.52 ± 0.08 (16S), and 0.37 ± 0.08 (combined total tRNA). Thus, the relative evolutionary rate of the fastest gene (ND2 or ND5) was approximately 4.4 times faster than the slowest gene (12S rRNA). Concerning the individual gene sequences, we attempted to calculate the dN and dS based on the tree topology of total mtDNA sequences, because if adaptive selection occurred in a lineage, dN would exceed the dS in the branch of the tree

18

M. Hirayama et al. / Gene 457 (2010) 13–24

Table 2 Sequence variation of mitochondrial genome-encoded genes. Gene

Length

Dmax (%)

Relative rate

dN/dS

ATPase 6 ATPase 8 COI COII COIII Cyt b ND1 ND2 ND3 ND4L ND4 ND5 ND6 12S rRNA 16S rRNA Total tRNAa

683 168 1557–1566 691 784 1141 975 1051 349 297 1381 1836–1839 522 941–944 1670–1674 1549–1551

18.0 14.9 15.1 11.7 15.1 14.5 20.8 21.6 24.9 16.2 18.1 19.3 18.4 5.3 7.7 6.0

1.02 ± 0.20 0.87 ± 0.20 1.04 ± 0.08 0.85 ± 0.13 0.97 ± 0.13 1.03 ± 0.12 1.42 ± 0.19 1.54 ± 0.17 1.27 ± 0.34 1.11 ± 0.56 1.29 ± 0.10 1.53 ± 0.14 1.30 ± 0.15 0.35 ± 0.11 0.52 ± 0.08 0.37 ± 0.08

0.04 ± 0.05 0.18 ± 0.34 0.01 ± 0.03 0.07 ± 0.22 0.03 ± 0.06 0.08 ± 0.24 0.07 ± 0.10 0.12 ± 0.13 0.09 ± 0.16 0.07 ± 0.09 0.10 ± 0.21 0.15 ± 0.10 0.09 ± 0.14 – – –

Dmax, the maximum value of the pairwised sequence differences; relative rate, mean ±S.D. of the relative molecular evolutionary rates between the gene and total mtDNA sequences; dN/dS, mean± S.D. of dN per dS of all branches in a tree. Relative rates were using the pairwised genetic distances (Kimura's two parameter model), and each pairwised value of the gene was divided by each pairwised value of total mtDNA sequences. COI–III, cytochrome oxidase subunits I–III; Cyt b, cytochrome b; dN, nonsynonymous substitutions per nonsynonymous site; dS, synonymous substitutions per synonymous site; ND1–6, and ND4L, NADH dehydrogenase subunits 1–6 and 4L. a 22 tRNA gene sequences were combined into a single data set.

(Yang, 1997). Therefore we estimated dN and dS values on the 273 branches (21 branches in a tree by 13 gene trees). The results indicated that dN exceeded dS in 11 branches of 8 genes (ATP6, ATP8, COI, COII, Cyt b, ND3, ND4 and ND6), but in most of the cases were there was only one nonsynonymous, and no synonymous substitutions (9 of 11 cases). Two cases evolved two nonsynonymous per one synonymous substitution, with likelihood estimations of dN/dS= 1.26 and 1.03. Of the 11 branches, only two were internal branches. One branch was between the node of the Toyooka stock and the node of the KAGA strain in COI, and the other was between the node of the Hirosaki stock and the node of the HNI strain and Niigata stock in the Cyt b tree. Most of the dN/dS ratios in each gene tree were not high, with an average ratio between 0.01 and 0.18 (Table 2). Of the 22 mtDNA-encoded tRNA genes, both conserved and variable genes were observed (Table 3). The maximum sequence variation observed in the tRNA genes, excluding tRNAAla which was Table 3 Sequence variation of mitochondrial tRNA genes. Gene

Length

Dmax (%)

tRNAAla tRNAArg tRNAAsn tRNACys tRNAAsp tRNAGln tRNAGlu tRNAGly tRNAHis tRNAIle tRNALeu (CUN) tRNALeu (UUR) tRNALys tRNAMet tRNAPhe tRNAPro tRNASer (AGY) tRNASer (UCN) tRNAThr tRNATrp tRNATyr tRNAVal

69 69 73 65–66 69–70 71 68 69–70 71 68 73 74 73 69–70 69 70 69 71 73 72 70 72

0 5.8 4.1 12.3 12.9 7.0 10.3 5.7 11.3 7.4 5.5 5.4 8.2 4.3 5.8 4.3 11.6 7.0 13.7 9.7 8.6 4.2

Dmax, the maximum value of the pairwised sequence differences.

perfectly conserved, was 4.1–13.7%. The most divergent gene encoded for tRNAThr, with 10 of 73 nucleotides differing between the HSOK and KAGA strains or Toyooka stock (Fig. 4). The complete CR sequence of HNI of the Northern Japanese group contained a putative termination-associated sequence (TAS) (338–353 bp from the start position of CR) and 3 conserved sequence blocks (CSBs) (CSB1, 797–816 bp; CSB2, 961–976 bp; CSB3, 1005– 1023 bp), similar to typical vertebrate mitogenomes (Fig. 5). As the several samples (Hirosaki, Toyooka, Amino, Matsuyama, Nago, Shanghai and HSOK) had an undetermined region in the posterior CR, we only compared the anterior segment of CR. In the anterior half of the CR, all of the mitogenomes contained tandem repeats with very similar flanking regions. When we compared the flanking region (408 bp) following the repeated sequence, the sequence variation ranged from 10.6 to 13.0% between the Japanese Medaka and HSOK strain or Shanghai stock, 7.6 to 9.4% between the Northern and the Southern Japanese groups, and only 0 to 2.3% and 0.2 to 3.0% within the mtDNA of the Northern and the Southern Japanese groups, respectively. These sequence variations were lower than those observed in the protein genes (see Table 2). 3.4. Variable number of tandem repeats in the mtDNA control region The most repeated motif in the Northern Japanese group and Toyooka stock was GCGTTGCATGT, and in the Southern Japanese group, HSOK strain and Shanghai stock was GCATGCGCGTT (indeed, the first and last 5 nucleotide sequences conformed to the 6–10th and 1–5th of the most repeated motif in the Northern Japanese group, respectively). Most of the tandem repeats consisted of perfect repeats, however, the last motif in the Northern Japanese group differed by one nucleotide (GCGTTGCATGA) (Fig. 5). For O. curvinotus and O. luzonensis, the closest relatives of Medaka (Takehana et al., 2005), the CR sequences did not contain the repeated sequences. A total of 34 local stocks and inbred strains, including data from 12 mitogenomes, showed extensive variability in the number of tandem repeats (Table 1). Twelve of the 18 stocks in the Northern Japanese group had perfectly conserved repeats, while the sequence of Kinosaki stock consisted of GCGTTGCATGT and ACGTTGCATGT, those of Odate, Yokote, Tsuruoka, Obama and Miyazu stocks consisted of GCGTTGCATGT and GCATTG CATGT. The sample from Yamagata stock had a sequence consisting of one GCGTTGCATGT, twenty-three GCATTGCATGT, and three GCATCG CATGT. The Southern Japanese group, HSOK strain and Shanghai stock had a conserved repeat unit of 11 nucleotides (GCATGCGCGTT), which was repeated from 7 to 12 times. Interestingly, there was a correlation between the number of tandem repeats and latitude, with the repeat number tending to increase in higher latitudes (R2 = 0.3348, P b 0.05, in the Northern Japanese group; R2 = 0.236, P = 0.078, in the Southern Japanese group; R2 = 0.3239, P b 0.005, in composite data of two groups) (data not shown). Moreover, the specific meteorological data provided by the Japan Meteorological Agency clearly correlated with the number of tandem repeats in the mtDNA (see Table 4, Fig. 6, Supplementary Table 2). Large values of multiple correlation coefficients were shown between the number of tandem repeats in the mtDNA and the following meteorological data: in the Northern Japanese group, index data of adaptation to low temperature including mean lowest TA, minimum lowest TA, annual number of days of average TA b0 °C, annual number of days of lowest TA b0 °C, and index data of adaptation to seasonal temperature change including the difference between maximum highest and minimum lowest TAs; in the Southern Japanese group, index data of adaptation to immediate temperature change including amount of precipitation, annual number of days with precipitation N10 mm, annual number of days with precipitation N30 mm; and in a composite of the two groups, index data of adaptation to low temperature including minimum lowest TA, annual number of days of average TA b0 °C,

M. Hirayama et al. / Gene 457 (2010) 13–24

19

annual number of days of lowest TA b0 °C, annual number of days of highest TA b0 °C, index data of adaptation to seasonal temperature changes including the difference between maximum highest and minimum lowest TAs, index data of adaptation to continuous low temperature including amount of snowfall (this also meant index data of adaptation to immediate temperature change), and annual number of days with snow cover N5 cm, N10 cm, N20 cm, and N50 cm. The correlations were all significant (P b 0.01, P b 0.05 and P b 0.0005 for the Northern Japanese group, Southern Japanese group, and composite of two groups, respectively). There was no comparison which showed a statistically significant correlation between the number of tandem repeats and the index data of adaptation to high temperature. 3.5. Relative DNA and mRNA levels of mtDNA-coding genes in 5 °C-acclimated fish

Fig. 4. Polymorphisms in Medaka tRNAThr between the Northern and Southern Japanese groups. The nucleotide position, and both high and low frequency (N80% conserved) single nucleotide polymorphisms between the Northern and Southern Japanese groups are indicated.

We compared the relative DNA and mRNA expression levels of mtDNA genes (ND1, COI and ATP6) in 5 °C- and 25 °C-acclimated fish from 4 Japanese stocks with different numbers of tandem repeats in the CR. From the qPCR analysis, almost identical results were obtained for the DNA levels of all 3 mtDNA-encoded genes and were not significantly different among the 4 local stocks, regardless of the temperature shift (Fig. 7A). These results demonstrated that the cellular level of mtDNA did not change significantly in response to shifts in temperature. On the other hand, results from the qRT-PCR analysis indicated that the mRNA levels of tested genes varied among the stocks. Generally, the levels in the fish acclimated at 5 °C tended to be higher than those at 25 °C (Fig. 7B). Although the standard deviations of the mRNA levels in 5 °C-acclimated fish were relatively large for each local stock, it seemed to have a correlation between the number of tandem repeats, especially in COI, and the mRNA level. It should be noted that the mRNA levels of all 3 genes in Nago, which

Fig. 5. The schematic diagram and complete sequence of the mtDNA control region of the Medaka HNI inbred strain. (A) Schematic diagram of HNI mitochondrial DNA around the control region. tRNApro is encoded by the L-strand. The number of tandem repeats varies according to local stocks and strains. Positions of OH, PH, and PL (Falkenberg et al. 2007) are indicated. (B) The full nucleotide sequence of Medaka mitochondrial control region (15,641–16,745th nucleotide position in HNI mitogenome). Tandem repeated sequences are boxed with broken line. A TAS and CSBs are labeled with bold and inversion, respectively. 12S, 12S rRNA; CSB, conserved sequence block; Cyt b, cytochrome b; OH, origin of heavystrand replication; PH, heavy-strand promoter; PL, light-strand promoter; TAS, termination-associated sequence.

20

M. Hirayama et al. / Gene 457 (2010) 13–24

Table 4 Multiple correlation coefficient between the number of tandem repeats in mtDNA and the meteorological data. Index of General environmental temperature adaptation Group

Mean average TA (°C)

Northern Japan 0.393 Southern Japan 0.232 Composite data 0.328

High temperature

Seasonal temperature change

Air temperature (TA) Mean Max. highest highest Average TA TA (°C) TA (°C) N 25 °C

Lowest TA N 25 °C

Highest TA N 25 °C

Highest TA N30 °C

Highest TA N35 °C

0.269 0.229 0.290

0.314 0.146 0.129

0.072 0.190 0.096

0.115 0.169 0.105

0.061 0.068 0.001

0.089 0.002 0.000

Annual number of days

0.295 0.171 0.193

(Mean highest TA) − (mean lowest TA) (°C)

(Max. highest TA) − (Min. lowest TA) (°C)

0.136 0.081 0.093

0.646 0.279 0.411

Index of environmental adaptation

Low temperature

Continuous low temperature

TA

Snowfall

Group

Mean lowest TA (°C)

Min. lowest TA (°C)

Annual number of days Average TA b0 °C

Lowest TA b 0 °C

Highest TA b 0 °C

0.460 0.241 0.333

0.607 0.315 0.414

0.464 0.334 0.463

0.512 0.289 0.411

0.358 0.322 0.418

Northern Japan Southern Japan Composite data

Index of environmental adaptation

Continuous low temperature/ immediate temperature change

Immediate temperature change

Group

Snowfall

Precipitation (PRCP)

Northern Japan Southern Japan Composite data

Max. amount of snow cover (cm)

Annual number of days of amount of snow cover N 5 cm

N 10 cm

N 20 cm

N50 cm

N 100 cm

0.208 0.075 0.380

0.327 0.022 0.507

0.336 0.065 0.507

0.336 0.109 0.506

0.330 0.161 0.423

0.206 0.073 0.277

Amount of snowfall (cm)

Max. snowfall per day (cm)

Amount of PRCP (mm)

Max. amount of PRCP per day (mm)

Annual number of days of amount of PRCP N 1 mm

N 10 mm

N 30 mm

N 50 mm

N70 mm

N100 mm

0.216 0.012 0.416

0.006 0.052 0.126

0.241 0.433 0.028

0.347 0.124 0.313

0.000 0.178 0.156

0.257 0.478 0.010

0.269 0.610 0.160

0.245 0.149 0.263

0.294 0.186 0.277

0.251 0.230 0.189

The meteorological data were cited from the database of Japan Meteorological Agency (http://www.jma.go.jp/). The data from 32 years spanning 1976–2007 were averaged.

contains the lowest repeat number among the stocks, were not significantly changed. 4. Discussion 4.1. Intraspecific variation in the mitochondrial gene order The mitochondrion plays a critical role in the maintenance of cellular energy stores, thermogenesis, and apoptosis. In this study, the complete mitogenome sequences of 8 Medaka stocks and 4 inbred strains were compared. A local duplication was found between the ND3 and ND4 genes in the mtDNA of the HSOK strain, but the duplicated genes had lost all function. The arrangement of the pseudogenes between the COIII and ND4 gene region can be explained by a single duplication event of the segment tRNAGly–ND3–tRNAArg– ND4L–ND4, followed by accumulated mutations (Fig. 2). Gene rearrangements reported for vertebrate mtDNAs are most typically, but not exclusively, local rearrangements among clustering tRNA and protein genes, but can also involve tandem duplication of a continuous mtDNA segment (Amer and Kumazawa, 2007). To the best of our knowledge, our results represent the first report of intraspecific variation in the mitochondrial gene order. 4.2. Phylogenetic relationship of Medaka populations and evolutionary rates of the mtDNA-coding genes Natural selection acting on mtDNA has been argued for many organisms (Elson et al., 2004; Moilanen and Majamaa, 2003). Phylogenetic analysis using the DNA and amino acid sequences of the proteincoding genes showed the same topology with each other (data not

shown). The phylogenetic relationship of Medaka mtDNA was previously reported by using Cyt b gene sequences from over 100 individuals (Takehana et al., 2003, 2004), and the results were fully consistent with those of the present study (Fig. 3). In addition, most of the dN/dS ratios calculated for each gene tree reported here were not high, and averaged 0.01–0.18 (Table 2). These results indicate that most of the nucleotide substitutions in the mitochondrial protein-coding genes may be nondeleterious and non-adaptive ones, but a few mutations may be retained by adaptive selection. The structural and functional properties of proteins, such as the conservation of ligand binding and catalytic rate, can be affected by minor changes in the sequence, and help species adapt to different temperatures (Somero, 1995, 2003). The relative evolutionary rates of rRNA and tRNA genes were slower than the protein-coding genes (Table 2). Among Medaka tRNA genes, the most divergent was the tRNAThr gene, as is also found in human mtDNA. Kivisild et al. (2006) analyzed 277 human mtDNA sequences and found that the mitochondria-encoded tRNAThr varied significantly more than any other tRNA gene. We found that 15 out of 73 nucleotides were varied among 12 Medaka tRNAThr sequences (Fig. 4). Although there are also polymorphisms in 15 positions out of 66 nucleotides among 277 human tRNAThr sequences, there is no significant relationship between the position of the polymorphism in these two organisms. In addition, although human tRNAThr has three hot spot positions, none of these positions in Medaka tRNAThr was variable. Positions 15,539 and 15,559, where 2 out of 6 nucleotides differences occur between the two Japanese groups, are located in the anticodon stem and the base of the TψC loop. It is probable that these are important sites for the proper functioning of tRNAThr because pathological mutations have also been found in the corresponding positions of other tRNAs (Kivisild et al., 2006).

M. Hirayama et al. / Gene 457 (2010) 13–24

21

Fig. 6. The correlation between the number of tandem repeats in the mitogenome and meteorological data. Two typical correlations in each group were indicated. The correlation between the number of tandem repeats of the Northern Japanese Medaka and (A) minimum lowest temperature or (B) difference between the two extreme temperatures in each collection site. The correlation between the number of tandem repeats of the Southern Japanese Medaka and (C) annual number of days with precipitation N 10 mm, or (D) annual number of days with precipitation N 30 mm in each collection site. The correlation between the number of tandem repeats for composite data of the Northern and Southern Japanese groups and (E) difference between the two extreme temperatures, or (F) annual number of days with snow cover N10 cm in each collection site.

4.3. Structural properties of the control region in Medaka mtDNA The mtDNA CR is known to contain regulatory elements and encompasses the heavy (H)- and L-strand promoters (PH, PL), the origin of H-strand replication (OH), three CSBs (CSB1, 2, 3), and TAS (Coskun et al., 2003). In vertebrate mtDNA, transcription is initiated bidirectionally from two promoters, PH and PL (Kelly and Scarpulla, 2004; Shadel and Clayton, 1997). Both promoters are associated with binding sites for the nuclear DNA-coded transcription factor Tfam. We

found specific tandem repeats in the mtDNA CR of Medaka (Fig. 5) which have also been reported in various species of invertebrates and vertebrates, including endotherms and ectotherms (Lunt et al., 1998) and several other fish species (Faber and Stepien, 1997; Hoarau et al., 2002; Ludwig et al., 2000; Nesbo et al., 1998). In Medaka, the tandem repeated sequences were located in the peripheral domain, which is approximately 60 bp from tRNAPro (Fig. 5). Located approximately 140 bp from the tandem repeat sequences is a putative TAS sequence, which is a short DNA element conserved in vertebrates and located at

22

M. Hirayama et al. / Gene 457 (2010) 13–24

Fig. 7. Relative DNA (A) and mRNA expression levels (B) of 3 mtDNA genes from 4 local stocks acclimated to 5 °C. Relative levels of ND1, COI and ATP6 are represented. Relative DNA and mRNA levels of each gene were normalized with those of β-actin and also with the 25 °C-acclimated fish. Asterisks indicate significant differences (P b 0.05) between the 5 °C and 25 °C groups for each stock as determined by Student's t-test. ATP6, ATPase subunit 6; COI, cytochrome oxidase subunit I; ND1, NADH dehydrogenase subunit 1.

the 3′ end of nascent D-loop H-strands (Doda et al., 1981). TAS sequences have been proposed as a major regulation point of mtDNA replication involving a 48 kDa protein of unknown identity which binds the TAS sequence and promotes D-loop formation in bovine mitochondria (Brown and Clayton, 2002; Madsen et al., 1993; Roberti et al., 1998). The sequence variations in the flanking region (408 bp) following the repeated sequences in the CR among the 12 Medakas were lower than those of the protein-coding genes. This suggests that the region is functionally important as it contains TAS sequences and other putative regulatory elements. The T414G mutation in the human mtDNA CR occurs in close proximity to the Tfam binding site of the PL and thus could perturb both L-strand transcription and, because the L-strand acts as its primer, H-strand replication (Michikawa et al., 1999). Nuclear DNA-encoded regulatory factors may bind to CR sequences, even in the region of tandem repeated sequences. We also found much diversity in the number of tandem repeated sequences in the mtDNA CR among local stocks when partial CR sequences of a total of 36 stocks, strains and related species were compared. Particularly in the Northern Japanese group, the conserved 11 nucleotides sequence was repeated from 10 to 34 times (Table 1). Since a high mutation rate in mtDNA CR including the tandem repeated sequence variances is observed in various organisms (Berg et al., 1995; Fumagalli et al., 1996; Lunt et al., 1998), we amplified the region including tandem repeated sequences by PCR using 8 fish from the Hirosaki stock to estimate the variance of number of tandem repeats within a local stock (Supplementary Fig. 2). The major amplified products in each fish showed similar size, although minor amplified products which indicated the

heteroplasmy of mtDNA in Medaka were found. The local stocks used in this study are originally collected from local-wild schools. They are maintained as closed colonies more than 25 years and they might have reduced diversity in each generation. However, Katsumura et al. (2009) showed that the grid-based sampling using our local stocks gives more abundant variations in both CR and Cyt b gene than the deme-based sampling. Their phylogenetic network in Northern Japanese population based on CR sequences including the 5′ seven repetitive sequences and excluding in/del polymorphisms is not related to the distribution of number of repeated sequences in CR. Takehana et al. (2003) also reported that a certain subclade in the Northern Japanese group is characterized by limited differentiation, as evidenced by the short branch on the neighbor joining tree analyzed by PCR-restriction fragment length polymorphism analysis of the mitochondrial Cyt b gene. Takehana et al. (2003) also surmised that this decline in genetic variation can be attributed to a bottleneck effect, and that the subclade recently expanded its distribution area. In this study, however, the number of tandem repeats in the mtDNA varied highly in these local stocks. This implies that the tandem repeats rapidly developed in this group, possibly by positive selection to adapt to a certain environment. 4.4. Correlation between the number of tandem repeated sequences in the mtDNA control region of Medaka stocks and the meteorological data in those habitats For both the Northern and Southern Japanese groups, the tandem repeat number clearly increased in cooler environments in higher

M. Hirayama et al. / Gene 457 (2010) 13–24

latitudes. As the climatic conditions in the Japanese Archipelago are affected by the geographical features (mountain shape), monsoon, and oceanic currents, we chose the index data of climatic differences among local areas using a database provided by the Japan Meteorological Agency. We comprehensively compared the number of tandem repeats in the mtDNA with the meteorological data, assuming that (1) water temperature would be parallel to the atmospheric temperature (Caissie et al., 2001; Langan et al., 2001), (2) rainfall and snowfall could change the environmental water temperature directly and immediately and (3) snow cover could also affect water temperature. To date, most studies have exclusively focused on the phylogenetic and evolutionary systems of those sequences. Although mtDNA length differences in endotherms and ectotherms, which could be derived from differences in the tandem repeats, were reviewed by Rand (1993), our study is the first to examine the functional relationship of the tandem repeats to environmental adaptation. The repeat number was negatively correlated with the minimum lowest air temperature in the original habitats (R 2 = 0.414, P b 0.0005). Interestingly, the meteorological data were more closely correlated with the number of tandem repeats in the Northern Japanese group (R2 = 0.607, P b 0.0005) than with those in the Southern one (R2 = 0.315, P b 0.05) (Table 4). This implies that the mtDNA of the Northern Japanese group has been affected by the cooler environment. Although the function and evolutionary significance of the repeated sequences were unknown, the repeated sequence has a secondary structure with negative associated energies (Faber and Stepien, 1997; Supplementary Fig. 1). If the structure, or repeated sequence itself, was a recognition site for enzymes of replication or transcription of the mitogenome, a larger structure could be more easily recognized by enzymes which are depressed in cooler conditions. The higher latitude Medaka have larger tandem repeats in the mtDNA CR which may imply that variations in the growth rate (Yamahira et al., 2007; Yamahira and Takeshi, 2008) and in the repeat numbers could be also correlated. 4.5. Relationship between the number of tandem repeated sequences in mtDNA and the replication and transcription levels of mtDNA in the low temperature-acclimated Medaka The index data of adaptation to continuous low temperature in the composite data of the two groups showed a strong correlation with the number of tandem repeats (R2 =0.507, Pb 0.0005). We examined the DNA and mRNA levels of 3 mtDNA-encoded genes in the 5 °C-acclimated fish from 4 local stocks (Fig. 7). We found inter-populational variations in the mRNA levels of these 3 genes at 5 °C, and interestingly, the mRNA levels of COI seem to correlate with repeat number. Battersby and Moyes (1998) reported that although no significant differences were detected between the two acclimation temperatures (either 4 or 18 °C for 2 months) in the abundance of steady-state mitochondrial mRNA (COI, ATP6), rRNA (16S), or DNA copy number in rainbow trout Oncorhynchus mykiss, mitochondrial enzyme activities in both red and white muscle were higher. Oleksiak et al. (2005) also suggested that subtle variation of gene expression can affect the physiological performance of organisms. The differences in the mRNA expression response of mtDNA genes in coldacclimated Medaka, in which the number of the tandem repeats varies, can differentiate their ability to adapt to cold temperature environment. Moreover, the mRNA levels of all 3 mtDNA-coding genes in the 5 °Cacclimated fish from the Ichinoseki stock (12 repeats), which belongs to the Southern Japanese group and inhabits high latitude, were significantly higher than the levels in the 25 °C-acclimated fish (Fig. 7B). It implies that nuclear genome-encoded genes, such as Tfam as described above, may affect the transcription of mtDNA-encoded genes and that the coevolution of mitochondrial and nuclear genomes may have occurred. To investigate the functional relationship of the number of tandem repeats in the CR to temperature adaptation in Medaka, crossbreeding between the Northern Japanese (Odate stock) and Southern Japanese Medaka (Hd-rR strain),

23

which have different mitochondrion, and quantification of their expression levels in mtDNA-coding gene (COI) were examined using 16 °C- and 25 °C-incubated embryos of F0 (Odate and Hd-rR) and their F1 hybrids (Odate/Hd-rR (female/male) and Hd-rR/Odate) (Supplementary Fig. 3). The incubation of Medaka embryos at 16 °C does not affect their survival and hatchability. However, tandem repeats number-dependent manner in the expression pattern of COI were not found even in F1 hybrids, which had the different number of tandem repeats in mtDNA (Supplementary Fig. 3A). Comparing to the results in the adult fish adapted to 5 °C (Fig. 7), the tandem repeats could contribute to the mtDNA expression levels in a stricter environment. To clarify this issue, further investigations will be required. In conclusion, this study revealed that intra-species mitogenomic variations in Medaka and that the repeated sequences in the mtDNA CR might function for mtDNA gene expression. Moreover, the number of tandem repeats in the Medaka CR is possibly related to adaptation to a certain habitat, especially a strict cold environment. Acknowledgments We would like to express our sincere thanks to Dr. G.N. Somero, Stanford University, USA, for critical reading of the manuscript. Thanks are also due to Drs. Y. Ishikawa (National Institute of Radiological Sciences) and N. Shibata (Shinshu University) for kindly providing inbred Medaka strains and foreign species, respectively. We are grateful to Mr. Y. Hashiguchi (Ocean Research Institute, The University of Tokyo) for providing valuable technical advice concerning the computational analysis. This work was supported in part by a Grant-in-Aid for Scientific Research on the Priority Areas “Study of Medaka as a Model for Organization and Evolution of the Nuclear Genome” (no. 813) and “Comparative Genomics” (no. 015) from the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT) to H.M. M.H. was supported by Research Fellowships of the Japan Society for the Promotion of Science (JSPS) for Young Scientist. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.gene.2010.02.012. References Amer, S.A., Kumazawa, Y., 2007. The mitochondrial genome of the lizard Calotes versicolor and a novel gene inversion in South Asian draconine agamids. Mol. Biol. Evol. 24, 1330–1339. Avise, J.C., 1994. Molecular Markers, Natural History and Evolution. Chapman & Hall, New York. Avise, J.C., et al., 1987. Intraspecific phylogeography: the mitochondrial DNA bridge between population genetics and systematics. Ann. Rev. Ecol. Syst. 18, 489–522. Battersby, B.J., Moyes, C.D., 1998. Influence of acclimation temperature on mitochondrial DNA, RNA, and enzymes in skeletal muscle. Am. J. Physiol. 275, R905–R912. Berg, T., Moum, T., Johansen, S., 1995. Variable numbers of simple tandem repeats make birds of the order Ciconiiformes heteroplasmic in their mitochondrial genomes. Curr. Genet. 27, 257–262. Brown, T.A., Clayton, D.A., 2002. Release of replication termination controls mitochondrial DNA copy number after depletion with 2′, 3′-dideoxycytidine. Nucleic Acids Res. 30, 2004–2010. Caissie, D., El-Jabi, N., Satish, M.G., 2001. Modelling of maximum daily water temperatures in a small stream using air temperatures. J. Hydrol. 251, 14–28. Coskun, P.E., Ruiz-Pesini, E., Wallace, D.C., 2003. Control region mtDNA variants: longevity, climatic adaptation, and a forensic conundrum. Proc. Natl. Acad. Sci. USA 100, 2174–2176. De Benedictis, G., et al., 1999. Mitochondrial DNA inherited variants are associated with successful aging and longevity in humans. FASEB J. 13, 1532–1536. Doda, J.N., Wright, C.T., Clayton, D.A., 1981. Elongation of displacement-loop strands in human and mouse mitochondrial DNA is arrested near specific template sequences. Proc. Natl. Acad. Sci. USA 78, 6116–6120. Doiron, S., Bernatchez, L., Blier, P.U., 2002. A comparative mitogenomic analysis of the potential adaptive value of arctic charr mtDNA introgression in brook charr populations (Salvelinus fontinalis Mitchill). Mol. Biol. Evol. 19, 1902–1909.

24

M. Hirayama et al. / Gene 457 (2010) 13–24

Elson, J.L., Turnbull, D.M., Howell, N., 2004. Comparative genomics and the evolution of human mitochondrial DNA: assessing the effects of selection. Am. J. Hum. Genet. 74, 229–238. Faber, J.E., Stepien, C.A., 1997. The utility of mitochondrial DNA control region sequences for analyzing phylogenetic relationships among populations, species, and genera of the percidae. In: Kocher, T.D., Stepien, C.A. (Eds.), Molecular Systematics of Fishes. Academic Press, San Diego, pp. 129–143. Falkenberg, M., Larsson, N.G., Gustafsson, C.M., 2007. DNA replication and transcription in mammalian mitochondria. Annu. Rev. Biochem. 76, 679–699. Fumagalli, L., Taberlet, P., Favre, L., Hausser, J., 1996. Origin and evolution of homologous repeated sequences in the mitochondrial DNA control region of shrews. Mol. Biol. Evol. 13, 31–46. Guderley, H., 2004. Metabolic responses to low temperature in fish muscle. Biol. Rev. 79, 409–427. Hirayama, M., Mitani, H., Watabe, S., 2006. Temperature-dependent growth rates and gene expression patterns of various medaka Oryzias latipes cell lines derived from different populations. J. Comp. Physiol. B 176, 311–320. Hittel, D.S., Storey, K.B., 2002. Differential expression of mitochondria-encoded genes in a hibernating mammal. J. Exp. Biol. 205, 1625–1631. Hoarau, G., Holla, S., Lescasse, R., Stam, W.T., Olsen, J.L., 2002. Heteroplasmy and evidence for recombination in the mitochondrial control region of the flatfish Platichthys flesus. Mol. Biol. Evol. 19, 2261–2264. Hyodo-Taguchi, Y., Sakaizumi, M., 1993. List of inbred strains of the medaka, Oryzias latipes, maintained in the Division of Biology, National Institute of Radiological Sciences. MEDAKA 5, 29–30. Itoi, S., Kinoshita, S., Kikuchi, K., Watabe, S., 2003. Changes of carp FoF1-ATPase in association with temperature acclimation. Am. J. Physiol. Regul. Integr. Comp. Physiol. 284, R153–R163. Iwamatsu, T., 1997. The Integrated Book for the Biology of the Medaka. Daigaku Kyoiku Shuppan, Okayama. Kasahara, M., et al., 2007. The medaka draft genome and insights into vertebrate genome evolution. Nature 447, 714–719. Katoh, K., Toh, H., 2008. Recent developments in the MAFFT multiple sequence alignment program. Briefings Bioinformat. 9, 286–298. Katsumura, T., et al., 2009. Genetic differentiation among local populations of medaka fish (Oryzias latipes) evaluated through grid- and deme-based sampling. Gene 443, 170–177. Kelly, D.P., Scarpulla, R.C., 2004. Transcriptional regulatory circuits controlling mitochondrial biogenesis and function. Genes Dev. 18, 357–368. Kimura, M., 1980. A simple method for estimating evolutionary rates of base substitution through comparative studies of nucleotide sequences. J. Mol. Evol. 16, 111–120. Kivisild, T., et al., 2006. The role of selection in the evolution of human mitochondrial genomes. Genetics 172, 373–387. Langan, S.J., Johnston, L., Donaghy, M.J., Youngson, A.F., Hay, D.W., Soulsby, C., 2001. Variation in river water temperatures in an upland stream over a 30-year period. Sci. Total Environ. 265, 195–207. Lucassen, M., Schmidt, A., Eckerle, L.G., Pörtner, H.-O., 2003. Mitochondrial proliferation in the permanent vs. temporary cold: enzyme activities and mRNA levels in Antarctic and temperature zoarcid fish. Am. J. Physiol. Regul. Integr. Comp. Physiol. 285, R1410–R1420. Ludwig, A., May, B., Debus, L., Jenneckens, I., 2000. Heteroplasmy in the mtDNA control region of sturgeon (Acipenser, Huso, and Scaphirhynchus). Genetics 156, 1933–1947. Lunt, D.H., Whipple, L.E., Hyman, B.C., 1998. Mitochondrial DNA variable number tandem repeats (VNTRs): utility and problems in molecular ecology. Mol. Ecol. 7, 1441–1455. Madsen, C.S., Ghivizzani, S.C., Hauswirth, W.W., 1993. Protein binding to a single termination-associated sequence in the mitochondrial DNA D-loop region. Mol. Cell. Biol. 13, 2162–2171. Matsuda, M., Yonekawa, H., Hamaguchi, S., Sakaizumi, M., 1997. Geographic variation and diversity in the mitochondrial DNA of the medaka, Oryzias latipes, as determined by restriction endonuclease analysis. Zool. Sci. 14, 517–526. Matsumoto, Y., et al., 2009. Medaka: a promising model animal for comparative population genomics. BMC Res. Notes 2, 88. Michikawa, Y., Mazzucchelli, F., Bresolin, N., Scarlato, G., Attardi, G., 1999. Agingdependent large accumulation of point mutations in the human mtDNA control region for replication. Science 286, 774–779. Mindell, D.P., Sorenson, M.D., Dimcheff, D.E., Hasegawa, M., Ast, J.C., Yuri, T., 1999. Interordinal relationships of birds and other reptiles based on whole mitochondrial genomes. Syst. Biol. 48, 138–152. Mishmar, D., et al., 2003. Natural selection shaped regional mtDNA variation in humans. Proc. Natl. Acad. Sci. USA 100, 171–176.

Mitani, H., et al., 2006. The Medaka genome: why we need the multiple fish models in vertebrate functional genomics. In: Volff, J.N. (Ed.), Structure and Evolution of Vertebrate Genomes. Genome Dynamics, vol. 2. Karger Publishers, Basel, pp. 165–182. Miya, M., Nishida, M., 2000. Use of mitogenomic information in teleostean molecular phylogenetics: a tree-based exploration under the maximum parsimony optimality criterion. Mol. Phylogenet. Evol. 17, 437–455. Miya, M., et al., 2003. Major patterns of higher teleostean phylogenies: a new perspective based on 100 complete mitochondrial DNA sequences. Mol. Phylogenet. Evol. 26, 121–138. Moilanen, J.S., Majamaa, K., 2003. Phylogenetic network and physiochemical properties of nonsynonymous mutations in the protein-coding genes of human mitochondrial DNA. Mol. Biol. Evol. 20, 1195–1210. Moritz, C., Dowling, T.E., Brown, W.M., 1987. Evolution of animal mitochondrial DNA: relevance for population biology and systematics. Ann. Rev. Ecol. Syst. 18, 269–292. Nesbo, C.L., Arab, M.O., Jakobsen, K.S., 1998. Heteroplasmy, length and sequence variation in the mtDNA control regions of three percid fish species (Perca fluviatilis, Acerina cernua, Stizostedion lucioperca). Genetics 148, 1907–1919. Oleksiak, M.F., Roach, J.L., Crawford, D.L., 2005. Natural variation in cardiac metabolism and gene expression in Fundulus heteroclitus. Nat. Genet. 37, 67–72. Rand, D.M., 1993. Endotherms, ectotherms, and mitochondrial genome-size variation. J. Mol. Evol. 37, 281–295. Roberti, M., Musico, C., Polosa, P.L., Milella, F., Gadaleta, M.N., Cantatore, P., 1998. Multiple protein-binding sites in the TAS-region of human and rat mitochondrial DNA. Biochem. Biophys. Res. Commun. 243, 36–40. Ross, O.A., et al., 2001. Mitochondrial DNA polymorphism: its role in longevity of the Irish population. Exp. Gerontol. 36, 1161–1178. Ruiz-Pesini, E., Mishmer, D., Brandon, M., Procaccio, V., Wallace, D.C., 2004. Effects of purifying and adaptive selection on regional variation in human mtDNA. Science 303, 223–226. Sakaizumi, M., Moriwaki, K., Egami, N., 1983. Allozymic variation and regional differentiation in wild population of the fish Oryzias latipes. Copeia 1983, 311–318. Setiamarga, D.H., Miya, M., Yamanoue, Y., Azuma, Y., Inoue, J.G., Ishiguro, N.B., Mabuchi, K., Nishida, M., 2009. Divergence time of the two regional medaka populations in Japan as a new time scale for comparative genomics of vertebrates. Biol. Lett. 5, 812–816. Shadel, G.S., Clayton, D.A., 1997. Mitochondrial DNA maintenance in vertebrates. Annu. Rev. Biochem. 66, 409–435. Shima, A., Shimada, A., Sakaizumi, M., Egami, N., 1985. First listing of wild stocks of the medaka Oryzias latipes currently kept by the Zoological Institute, Faculty of Science, University of Tokyo. J. Fac. Sci. Univ. Tokyo Sec. IV 16, 27–35. Skrypina, N.A., Timofeeva, A.V., Khaspekov, G.L., Savochkina, L.P., Beabealashvilli, R.Sh., 2003. Total RNA suitable for molecular biology analysis. J. Biotechnol. 105, 1–9. Somero, G.N., 1995. Proteins and temperature. Annu. Rev. Physiol. 57, 43–68. Somero, G.N., 2003. Adaptation of enzymes to temperature: searching for basic “strategies”. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 39, 321–333. Stamatakis, A., 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 22, 2688–2690. Swofford, D.L., 2002. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods) version 4.0b10. Sinauer Associates, Sunderland. Takehana, Y., Nagai, N., Matsuda, M., Tsuchiya, K., Sakaizumi, M., 2003. Geographic variation and diversity of the cytochrome b gane in Japanese wild populations of medaka, Oryzias latipes. Zool. Sci. 20, 1279–1291. Takehana, Y., Uchiyama, S., Matsuda, M., Jeon, S.R., Sakaizumi, M., 2004. Geographic variation and diversity of the cytochrome b gene in wild populations of medaka (Oryzias latipes) from Korea and China. Zool. Sci. 21, 483–491. Takehana, Y., Naruse, K., Sakaizumi, M., 2005. Molecular phylogeny of the medaka fishes genus Oryzias (Beloniformes: Adrianichthyidae) based on nuclear and mitochondrial DNA sequences. Mol. Phylogenet. Evol. 36, 417–428. Tanaka, M., Gong, J.S., Zhang, J., Yoneda, M., Yagi, K., 1998. Mitochondrial genotype associated with longevity. Lancet 35, 185–186. William, J., Ballard, O., Kreitman, M., 1995. Is mitochondrial DNA a strictly neutral marker? Trends Ecol. Evol. 10, 485–488. Yamahira, K., Takeshi, K., 2008. Variation in juvenile growth rates among and within latitudinal populations of the medaka. Popul. Ecol. 50, 3–8. Yamahira, K., Kawajiri, M., Takeshi, K., Irie, T., 2007. Inter- and intrapopulation variation in thermal reaction norms for growth rate: evolution of latitudinal compensation in ectotherms with a genetic constraint. Evolution 61, 1577–1589. Yang, Z., 1997. PAML: a program package for phylogenetic analysis by maximum likelihood. CABIOS 13, 555–556.