Phylogenetic and synteny analysis
To identify the visual opsin gene repertoire of turbot, zebrafish opsin sequences were used as BLASTp query sequences with e-value < 10–10, and they were downloaded from the chromosome of the reference genome. Opsin sequences of four flatfish (spotted halibut, Verasper variegatus; Japanese flounder, Paralichthys olivaceus; barfin flounder, Verasper moseri, and Atlantic halibut, Hippoglossus hippoglossus) and four freshwater species living in shallow water (zebrafish, Danio rerio; medaka, Oryzias latipes; guppy; and zebra mbuna, Maylandia zebra) were obtained from GenBank. For accession numbers, see previous research [11, 28, 44, 45]. Phylogenetic relationships among the visual opsin nucleotide sequences were inferred using MEGA 7 software [46] by applying the neighbor-joining method [47] and maximum composite likelihood model algorithm [48]. The reliability of tree topology was evaluated by bootstrap analysis with 1000 replications and uniform rates among sites. Genomicus synteny [49] and Ensembl genome browsers were used to assess the syntenic regions between turbot and five other teleost genomes (tongue sole, Cynoglossus semilaevis; Japanese medaka HNI; zebra mbuna; zebrafish, and guppy).
Branch and branch-site test of selection
To estimate the differences in selection pressure between five flatfish and four freshwater fish, branch-specific models and branch-site models of maximum likelihood were implemented in the CODEML program of PAML4.9i [50]. Tree topologies of each gene were the same as those obtained in section 2.1. First, we compared the one-ratio model and free-ratios model in the branch test of selection [24]. We then compared null Model A against the alternative Model A in the branch-site test to determine whether positive selection was present; different sets of foreground lineages were marked with letters in each gene (Fig. 1) [51]. The likelihood ratio test (LRT) was used to compare all alternative models and their corresponding null models. The Bayes Empirical Bayes method was used to obtain the posterior probability of sites under positive selection [52].
Analysis of spectral tuning sites
Representative spectral tuning sites of amino acids were compared among the opsins of five flatfish and four freshwater species. The sites were referenced from a previous study [10]. The amino acid sequence alignments were accomplished using ClustalX [53]. The number of amino acid sites was standardized to bovine rhodopsin, except for blue opsin, which was standardized to barfin flounder SWS2A.
Estimation of turbot RH2 divergence times
The divergence time of turbot RH2 genes was estimated by MCMCTREE within PAML4.9i [50]. The neighbor-joining tree of six species (turbot; Atlantic salmon, Salmo salar; fugu, Takifugu rubripes; common carp, Cyprinus carpio; lungfish, Neoceratodus forsteri, and rainbow trout) was acquired by MEGA 7 using the same method as described above. The fossil calibrations were adopted from TimeTree [54]. The tree topology and fossil calibrations were set as (((((((S.maximus_RH2B1, S.maximus_RH2B2), S.maximus_RH2C), T.rubripes_RH2), (S.maximus_RH2A1, S.maximus_RH2A2)), (S.salar_RH2, O.mykiss_RH2)) '>1.86<2.27', (C.carpio_RH2-1, C.carpio_RH2-2)) '>2.05<2.55', N.forsteri_RH2). The RootAge was set as '<6.0'.
Quantification of turbot opsin RNA expression
Laval and juvenile turbot were bred at Shenghang Aquatic Science and Technology Company (Weihai, Shandong Province, China). The ages of the individuals used for analysis of heterochronic shifts in opsin gene expression were: 0.5 month (15 days post hatch), 1 month, 4 months, 9 months, and 18 months. All individuals for mRNA expression analysis were euthanized with 300 mg/L of MS-222 (Sigma, Shanghai, China) before being decapitated. The eyes were removed and immediately stored in liquid nitrogen until analyzed.
Total RNA was extracted using the RNA Isolation Kit (Vazyme Biotech Co, Nanjing, China). RNA purity and concentration were examined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Shanghai, China), and RNA integrity was verified by gel electrophoresis. The PrimeScript RT reagent kit with gDNA Eraser was used to synthesize first-strand cDNA from 0.5 μg of total RNA (TaKaRa, Dalian, China). Using the Bio-Rad CFX Connect™ Real-Time PCR System (Bio-Rad, Hercules, CA, USA), quantitative real-time PCR (qPCR) was conducted with TB Green Premix Ex Taq (TaKaRa) following the manufacturer's protocol. Melting curves were plotted to confirm amplification specificity. Using the 2−ΔΔCt method, the relative expression level of each opsin gene was normalized to β-actin, which was selected as the internal reference gene after evaluating the expression pattern of eight commonly used housekeeping genes [55]. Table S1 shows the specific primers for qPCR. The PCR mixture (20 µl) contained 10 μl of TB Green Premix Ex Taq, 7.6 μl of PCR-grade water, 1.6 μl of cDNA, and 0.4 μl of each of the primers. The qPCR reaction was performed in triplicate with the program of 95 ºC for 3 min, followed by 40 cycles at 95 ºC for 10 s, 57/58 ºC for 30 s, and 72 ºC for 30 s. Proportional opsin expression was determined as a fraction by calculating the proportion of each cone opsin (Ti) relative to the total cone opsin expression (Tall) as follows:
where Ei represents the PCR efficiency for each pair of primers and Cti is the critical cycle number for each gene [56, 57].
To assess significance of the change in opsin expression between different stages, the least significant difference (LSD) post hoc test with 95% confidence level was used. Data were analyzed using SPSS 23.0 software. Relative expression data are shown as the mean ± standard deviation, and proportional expression data are shown as averages.