Phylogenetic analysis and ontogenetic changes in the cone opsins of the western mosquitofish (Gambusia affinis).

To convert external light into internal neural signal, vertebrates rely on a special group of proteins, the visual opsins. Four of the five types of visual opsins-short-wavelength sensitive 1 (Sws1), short-wavelength sensitive 2 (Sws2), medium-wavelength sensitive (Rh2), and long-wavelength sensitive (Lws)-are expressed in cone cells for scotopic vision, with the fifth, rhodopsin (Rh1), being expressed in rod cells for photopic vision. Fish often display differing ontogenetic cone opsin expression profiles, which may be related to dietary and/or habitat ontogenetic shift. The western mosquitofish (Gambusia affinis) is an aggressive invader that has successfully colonized every continent except Antarctica. The strong invasiveness of this species may be linked to its visual acuity since it can inhabit turbid waters better than other fishes. By genome screening and transcriptome analysis, we identify seven cone opsin genes in the western mosquitofish, including one sws1, two sws2, one rh2, and three lws. The predicted maximal absorbance wavelength (λmax) values of the respective proteins are 353 nm for Sws1, 449 nm for Sws2a, 408 nm for Sws2b, 516 nm for Rh2-1, 571 nm for Lws-1, and 519 nm for Lws-3. Retention of an intron in the lws-r transcript likely renders this visual opsin gene non-functional. Our real-time quantitative PCR demonstrates that adult male and female western mosquitofish do not differ in their cone opsin expression profiles, but we do reveal an ontogenetic shift in cone opsin expression. Compared to adults, larvae express proportionally more sws1 and less lws-1, suggesting that the western mosquitofish is more sensitive to shorter wavelengths in the larval stage, but becomes more sensitive to longer wavelengths in adulthood.

Introduction

Vertebrates rely on a group of specialized photoreceptor cells in the retina of eyes to convert external light into neural signals. The visual pigment located in the photoreceptor outer segment is responsible for deciphering light information and it comprises a light-absorbing protein (visual opsin) and a vitamin A-derived chromophore. Visual opsins are membrane-bound G-protein coupled receptors (GPCRs), and vitamin A-derived chromophores include vitamin A1 (retinal) or A2 (3,4-dehydroretinal). Differences in phylogenetic classes of visual opsin genes or in the types of chromophore can alter the maximal absorbance wavelength (λmax) of visual pigments [1]. Photoreceptor cells can be categorized as rod cells (scotopic or low-light vision) or cone cells (photopic or well-lit vision), with cone cells being further divided into single-, double-, twin-, triple-, or quadruple-cone cells [2]. Based on their expression of visual opsin genes, photoreceptor cells can also be categorized into five main classes: Rh1, Sws1, Sws2, Rh2, and Lws [3]. Most vertebrates only express one kind of rhodopsin rh1 gene in rods [4], and the other four classes of visual opsins are all expressed in cones with different spectral sensitivities: Sws1 is sensitive to UV and violet light, Sws2 to blue, Rh2 to green, and Lws to red [1].

Habitat has a profound effect on the visual system of the animals living in it [5, 6]. For example, the cichlids in Lake Victoria have distinct lws alleles with different λmax values that correspond to varied lighting environments with depth [7], and the opsin expression profiles of bluefin killifish (Lucania goodei) differ between clear springs and tannin-stained swamps [8, 9]. Apart from their lighting environments, habitats can also exert a selection pressure on visual systems in other ways. For example, predation risk variation in guppy has been demonstrated to alter lws allele frequencies and opsin expression profiles among populations [10].

Aquatic habitats provide more diverse lighting environments than terrestrial habitats since water itself selectively absorbs light and suspended or dissolved substances can absorb and/or scatter light [11, 12]. Multiple independent opsin gene duplication events have occurred over the evolutionary history of teleosts [13]. For example, the haplochromine cichlids have up to six cone opsins with distinct λmax values, including one sws1, two sws2, two rh2, and one lws [14]. Notably even for animals with multiple cone opsins, these opsins are not expressed equally or constantly [15-17]. Lighting environments, diurnal periodicity, and diet all contribute to changes in cone opsin expression profiles [6, 15, 16]. Many fishes exhibit the ability to adjust their cone opsin expression profiles, termed “phenotypic plasticity”. One type of phenotypic plasticity is changes in expression profile as an individual grows. Such ontogenetic changes in opsin profiles have been reported for giant mottled eel (Anguilla marmorata), milkfish (Chanos chanos), barfin flounder (Verasper moseri), and spotted unicornfish (Naso brevirostris) [17-20]. Two other types of phenotypic plasticity are developmental and adult plasticity; with respect to opsins, the former pertains to how different rearing conditions can result in distinct cone opsin expression profiles [21, 22], and the latter explains how adult fish can alter their expression profiles to adapt to local lighting environments [23-25]. Both developmental and adult plasticity may enable fishes to flourish in turbid waters [22, 25].

Currently, the genomic arrangement of cone opsin genes has been fully investigated for two poeciliid fishes, i.e., guppy and green swordtail (Xiphophorous helleri) [26, 27]. Both of these poeciliids possess nine cone opsin genes, including one sws1, two sws2, two rh2, and four lws. The duplicated opsin paralogs are tandemly arranged, except for lws-r, with the two sws2, three lws, and two rh2 opsins forming a linkage group [28]. The evolutionary trajectories of lws genes in poeciliid fishes reflect interactions between genomic organization and sexual selection [29].

Not all poeciliid fishes have had their cone opsin genes investigated. Although the western mosquitofish (Gambusia affinis) (Poeciliidae, Cyprinodotiformes) has been introduced in many countries to control mosquito populations since the early 1900s [30], becoming one of the worst invasive species in the process [31], our understanding of its cone opsins remains limited [32]. To date, research has only demonstrated that the western mosquitofish possesses all four phylogenetic classes of cone opsin genes, including two sws2 genes (sws2a and sws2b), but exactly how many cone opsins it has remains uncertain [32]. Moreover, previous studies indicate that the western mosquitofish may display ontogenetic changes in cone opsin expression profiles because it exhibits an ontogenetic dietary shift, colors are critical to it in the mating choice of both sexes, and sex steroid estrogen influences its cone opsin expression profiles [32-36].

Sensory-based approaches can be applied to protecting endangered species and controlling alien ones [37]. A comprehensive understanding of the visual spectral characteristics of western mosquitofish could reveal how this species is so aggressively invasive and provide insights into the sensory biology of other invasive species. Here, we sought to determine how many cone opsin genes are present in the western mosquitofish and to predict the λmax values of the respective opsins. Moreover, since longer-wavelength signals, i.e., yellow, are more attractive to adult individuals [34, 35], we also hypothesized that the western mosquitofish would display an ontogenetic change in cone opsin expression profiles, with the visual spectrum shifting to red during development.

Materials and methods

Overall, a total of 84 western mosquitofish were used in this study. One adult specimen was used for transcriptome analysis, ten adult specimens were used for cloning cone opsin genes, 30 adult specimens were used for reproductive experiments, and the resulting 43 laboratory-raised specimens (comprising 30 larva, seven adult males, and six adult females) were used to quantify ontogenetic cone expression.

Fish collection and experimental setup

Western mosquitofish samples were caught using a hand-net at two different localities, Tunghai Lake at Tunghai University (Taichung City, Taiwan) and Linyuan Ocean Wetland Park (Kaohsiung City, Taiwan). Adult fish were caught from 9:00 to 11:00 AM and then euthanized using 0.025% MS-222 (Ethyl 3-aminobenzoate, methanesulphonic acid salt). After each specimen was rendered comatose, its two eyes were removed and stored in RNAlater solution (AM7020, Invitrogen), with the rest of the body being preserved in 100% ethanol at -20˚C until DNA extraction. Eye samples for cone opsin gene cloning and RNA-seq were first kept at 4˚C for 24 hours and then stored at -80 ˚C until RNA extraction. Overall, three males and three females from Tunghai University and two males and two females from Linyuan Ocean Wetland Park were prepared for DNA and RNA extraction.

A further 30 live adult specimens from Tunghai University were transported to the Marine Research Station of Academia Sinica (Yilan, Taiwan). These specimens were housed indoors in a glass aquarium (90 x 40 x 40 cm length:width:height) with continuous dripping influx of dechlorinated tap water to maintain water conditions as clean as possible. A full-spectrum LED grow light (Hydroponics Lin Agritech) was used as light source and placed approximately 28 cm above the water’s surface. The photoperiod of the climate-controlled room was set to 14:10 h light: dark, and the spectral irradiance of the light source was measured using a spectrometer (MK350N, Gamma Scientific) (S1 Fig). Fish were fed ad libitum with artificial fish feed (Otohime B2, Marubeni Nisshin Feed) twice a day. After one month, we prepared three floating fish-breeding cages in the aquarium to host pregnant females and to collect F1 offspring. When we had collected more than 50 F1 individuals, all parental individuals were removed from the aquarium and all F1 offspring were released from the cages into the main aquarium. These F1 offspring then mated to produce F2 individuals under the same lighting conditions as their parents. The floating fish-breeding cages were again employed for pregnant F1 females. The F1 adults and newborn F2 larvae (less than 10 days old) were used for ontogenetic study. In order to limit the effect of diurnal variation in opsin expression, fish were euthanized (0.025% MS-222) between 10:00 and 11:00 AM. After each specimen was rendered comatose, larval specimens were directly preserved in RNAlater solution (AM7020, Invitrogen) whereas, for adult specimens, body weight and standard length were first determined before removing both eyes and storing them in RNAlater solution (AM7020, Invitrogen). These samples were also kept at 4˚C for 24 hours before being transferred to -80˚C until RNA was extracted. All experiments were performed and specimens handled with approval (107–28) from the Institutional Animal Care and Use Committee (IACUC) of Tunghai University.

DNA and RNA extraction, and cDNA preparation

DNA samples were extracted from fin tissues using a DNA extraction kit (Cat No./ID: GS100, Geneaid). An RNeasy Plus Universal Kit (Cat No./ID: 73404, QIAGEN) was used to isolate total RNA according to the manufacturer’s protocol. One eyeball from each adult specimen or five larval heads as a pooled sample was placed in a 2 ml microcentrifuge tube with stainless steel beads for tissue homogenization. Cone opsin genes are exclusively expressed in the retina, so RNA extraction from whole larval heads does not bias quantification of cone opsins [38]. The tissue was homogenized using a Roche MagNA Lyser system (Roche). Total RNA content and quality were measured using a NanoDrop 1000 (Thermo Scientific). Two micrograms (μg) of total RNA were reverse-transcribed using a Verso cDNA Synthesis Kit (Cat No. 00764129, ThermoFisher Scientific) with Random Hexamer and Anchored Oligo-dT 3:1 (v/v) in a final volume of 20 μl. The remaining RNA was preserved at -80˚C after adding RNase inhibitor (RG90925, Lucigen).

Transcriptomics

A female specimen collected from Tunghai University was subjected to RNA sequencing (RNA-seq) by Genomics Tech. (Taipei, Taiwan). The mRNA was purified from total RNA using poly-T oligo-attached magnetic beads and then fragmented by heating. The RNA library was prepared with the TruSeq® Stranded mRNA Library Prep system (Illumina, San Diego, CA, USA). An Agilent Bioanalyzer 2100 system and a DNA High Sensitivity Chip was employed for quality control. Then, the RNA library was sequenced on an Illumina NovaSeq 6000 platform using a 150-bp paired-end strategy with an output of 6 Gb.

The raw reads were trimmed using Trimmomatic v0.36 [39] and then de-novo assembled in Trinity v2.8.4 [40]. Identical transcripts were eliminated using CD‐HIT‐EST v4.8 [41, 42], with a 95% clustering cut-off. Then, the trimmed reads were aligned onto the transcripts to estimate the sequencing depth of each transcript using Bowtie2 v2.3.4.3 [43] and RSEM v1.2.28 [44]. The average Fragments Per Kilobase Million (FPKM) was 52.67. Open reading frames (ORFs) were predicted in Transdecoder v5.3.0 (https://github.com/TransDecoder/TransDecoder) and annotated in BLASTX v2.5.0 [45]. Finally, the transcripts were translated into amino acid sequences and annotated using BLASTP v2.5.0 [45].

Identification of cone opsin genes

Known guppy cone opsin genes––sws1, sws2a, sws2b, rh2, and lws (Table 1)––were used as references to query and annotate the cone opsin genes from the western mosquitofish draft genome [46] and our draft transcriptome using BLASTN 2.8.1+ [45]. We designed specific pairs of primers for each of the identified cone opsin genes. Furthermore, we aligned the complete β-actin gene sequence of guppy (EU143771) and the partial β-actin gene sequence (KP284099) of mosquitofish to design a pair of primers spanning the β-actin exon and intron regions. All primer sequences are presented in Table 2.

Table 1

List of cone opsin genes from cyprinodontiform fishes, medaka (Oryzias latipes), and zebrafish (Danio rerio).

Scientific name	tmt-opsin	val-opsin	sws1	sws2	rh2	lws
Cypriniformes
Danio rerio	KT008411	NM_131586	BC060894	NM_131192	NM_131253	NM_001002443
			(sws1)	(sws2)	(rh2-1)	(lws-1)
					NM_182891	NM_001313715
					(rh2-2)	(lws-2)
					NM_182892
					(rh2-3)
					NM_131254
					(rh2-4)
Beloniformes
Oryzias latipes			AB223058	AB223056	AB223053	AB223051
			(sws1)	(sws2a)	(rh2-a)	(lws-a)
				AB223057	AB223054	AB223052
				(sws2b)	(rh2-b)	(lws-b)
					AB223055
					(rh2-c)
Cyprinodontiformes
Poeciliidae
Poecilia reticulata			HQ260685	HQ260684	HQ391990	HQ260679
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				HQ391991	HQ260683	HQ260680
				(sws2b)	(rh2-2)	(lws-2)
						HQ260681
						(lws-3)
						HQ260682
						(lws-r)
Poecilia caymanensis			KX768650	KX768635	KX768663	KX768609
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				KX768622		KX768598
				(sws2b)		(lws-2)
						KX768587
						(lws-3)
						KX768576
						(lws-r)
Poecilia velifera			KX768639	KX768624	KX768654	KX768600
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				KX768611		KX768589
				(sws2b)		(lws-2)
						KX768578
						(lws-3)
						KX768567
						(lws-r)
Limia nigrofasciata			KX768649	KX768634	KX768662	KX768608
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				KX768621		KX768597
				(sws2b)		(lws-2)
						KX768586
						(lws-3)
						KX768575
						(lws-r)
Heterandria formosa			KX768651	KX768636	KX768664	KX768610
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				KX768623		KX768599
				(sws2b)		(lws-2)
						KX768588
						(lws-3)
						KX768577
						(lws-r)
Micropoecilia bifurca			KX768646	KX768631	KX768659	KX768605
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				KX768618		KX768594
				(sws2b)		(lws-2)
						KX768583
						(lws-3)
						KX768572
						(lws-r)
Pamphorichthys minor			KX768647	KX768632	KX768660	KX768606
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				KX768619		KX768595
				(sws2b)		(lws-2)
						KX768584
						(lws-3)
						KX768573
						(lws-r)
Xiphophorus hellerii			GU454734		GU454732	DQ075246
			(sws1)		(rh2-1)	(lws-1)
					GU454733
					(rh2-2)
Gambusia affinis			MN817658	MN817660	MN817659	MN817662
			(sws1)	(sws2a)	(rh2-1)	(lws-1)
				MN817661		MN817663 (lws-3)
				(sws2b)
						Hoffberg, Troendle [46]
						(lws-r)

Two non-visual opsins—tmt-opsin and val-opsin—from zebrafish are used as outgroups for cone opsins.

Table 2

Sequences of the primers used for amplifying β-actin and cone opsin genes and for quantifying expression levels of cone opsin genes in the western mosquitofish.

Gene	Primer sequence (5'-3')	Length of Complete cDNA (bp)	Length of target cDNA segment (bp)	Reference
Gene amplification
sws1	For: TGAGGTCGGAGGAACTCA	1008	1068	this study
	Rev: TGGATCCTCTGAACTCAGATC
sws2a	For: AATTCAACCGAGTCGTAGAC	1056	1061	this study
	Rev: GTTTGGACATAGCATTTACGA
sws2b	For: TATACCGCATCCAGCTCGTA	1065	1143	this study
	Rev: ACACGTCAATGTTGGTCATTC
rh2-1	For: CAGTGAAATCTCCCACTCAC	1059	1081	this study
	Rev: CAGACTTGAAGACGCGATTA
lws-1	For: GAGGAAGGTCTGAGAGCTGA	1071	1040	this study
	Rev: TGATGTGGACACTTCAGAGC
lws-3	For: GCCGAGAGAAGATCTCCAGA	1071	1040	this study
	Rev: TGATGTGGACACTTCAGAGC
lws-r	For: ATCAAGCAGCTCAGATCGTC	*	1310	this study
	Rev: TGATGTGGACACTTCAGAGC
β-actin	For: CCTGTACGCTTCTGGTCGTA	1128	566 (776a)	this study
	Rev: CCTCCAATCCAGACAGAGTA
qPCR
sws1 (Ei = 96.65%)	For: CCGACTGGTACACCCACAAT		56	Friesen et al. 2017 [32]
	Rev: TGGACAGAGGCATGCAGAAG
sws2a (Ei = 100.05%)	For: TCTGCTTCTGCTTTGCCGTA		92	Friesen et al. 2017 [32]
	Rev: TGGTCACCTCCTTTTCAGCC
sws2b (Ei = 100.04%)	For: TTATCTGCAAGCCACTCGGG		104	Friesen et al. 2017 [32]
	Rev: AACACTGCATTCCCTCAGGG
rh2-1 (Ei = 96.45%)	For: TCAGAACAAAAAGCTGCGCC		46	Friesen et al. 2017 [32]
	Rev: TGAAACCAAAGGCGCACATG
lws-1 (Ei = 97.95%)	For: CTGAGAGCTGAGAAACCTTC		54	this study
	Rev: CTCTTGTTGTATCTTCGTGC
lws-3 (Ei = 97.10%)	For: AGAGAAGTGTGAGGTGAAGC		112	this study
	Rev: CTCTTGTTGTATCTTCGTGC

Ei = efficiency of qPCR primers.

*lws-r is likely a pseudogene.

aLength of target DNA segment.

PCR amplifications of opsin genes were performed in a final reaction volume of 25 μL, containing 2 ng cDNA, 1 μl each of the forward and reverse primers (10 μM), 12.5 μl of Fast-RunTM Advanced Taq Master Mix (ProTech, Taipei, Taiwan), and distilled water. The thermal cycling protocol was as follows: one cycle at 94°C for 4 min; 35 cycles of denaturation at 94°C for 30 sec, 55–60°C for 30 sec, and 72°C for 1 min; and a final single extension step at 72°C for 5 min. Sequencing was performed by Mission Biotech Inc., Taipei, Taiwan. The contig sequences were constructed using the program CodonCode Aligner 6.0.2, and the results were identified by BLAST analysis against the database of the National Center for Biotechnology Information (NCBI).

Phylogenetic analyses

The cone opsin gene sequences of other poeciliid fishes, medaka (Oryzias latipes), and zebrafish (Danio rerio) were downloaded from NCBI, and two non-visual opsin genes of the zebrafish—val-opsin and tmt-opsin—were selected as outgroups [38]. All opsin gene sequences included in our phylogenetic analysis are shown in Table 1. Gene sequences were aligned using MACSE v2, which is designed to align protein-coding nucleotide sequences based on their corresponding amino acid translations [47]. Phylogenetic analyses were performed using a partitioned Maximum Likelihood (ML) approach. RAxML 8 [48] was used for ML analyses (MLA). Partitions were set with respect to codon position, and the GTR+G+I model (with four discrete rate categories) was adopted for each partition [49]. The ML tree was obtained by performing 100 different runs using the default algorithm. The best ML tree was chosen from likelihood scores among suboptimal trees from each run. Nodal support for MLA was determined from bootstrap analysis using RAxML [50], employing 10,000 non-parametric bootstrap replications with the ML criterion.

Apart from targeting the coding region of opsins, we also conducted phylogenetic analysis on the 5’ untranslated region (UTR) of lws from many poeciliid fishes given that Watson et al. (2011) [26] and Sandkam et al. (2017) [29] have demonstrated that gene conversion between poeciliid lws paralogs is more frequent in coding regions than in non-coding regions. Upon adding our 5’ UTR data from western mosquitofish to the poeciliid data of Sandkam et al. (2017) [29], the sequences were aligned using MUSCEL in Seqotron [51]. MLA was performed with the GTR+G+I model (with four discrete rate categories) [49] in RAxML 8 [48]. The ML tree was obtained by performing 100 different runs using the default algorithm. The best ML tree was chosen from likelihood scores among suboptimal trees from each run. Nodal support for MLA was determined from bootstrap analysis with RAxML [50], employing 10,000 non-parametric bootstrap replications with the ML criterion.

Prediction of the λmax values of the cone opsins

To predict the λmax values of western mosquitofish cone opsins, we relied on inferences based on its closest evolutionary relative, the guppy, for which λmax values were determined previously from in vitro expression data [52]. Moreover, given that some amino acid residues of a visual opsin have been demonstrated to have more pronounced effects on the λmax value than other residues, for this study, not only did we target the five residues (the 180th, 197th, 277th, 285th, and 308th amino acid residues corresponding to the human M/LWS pigment) based on the “5-sites rule” [1, 53], but also an additional five residues (the 110th, 112th, 138th, 223th, and 227th amino acid residues corresponding to the human M/LWS pigment), which also influence the λmax value of Sws2a in barfin flounder (Verasper moseri) and Rh-2 opsins in Oryzia [20, 54]. We assumed that the λmax values of two opsins were identical when a western mosquitofish cone opsin presented the same ‘10 key-sites’ composition as the orthologous opsin in guppy.

Real-time quantitative PCR (qPCR)

The specific primers designed for qPCR were based on the annotated sequences from our whole-genome BLAST and our draft transcriptome. The amplification efficiency and melting curve of each qPCR primer pair was tested by 5-fold serial dilutions of the templates, with three replicates for each gene and sample. Each qPCR primer pair was adopted only when its amplification efficiency fell between 90% and 110% and the melting curve analysis revealed each pair of specific primers only generated a single product. Expression of opsin genes was determined by qPCR in a Roche LightCycler480 system (Roche). Each reaction contained 10 μl of Roche LightCycler480 SYBR Green I Master (Roche Applied Science, Indianapolis, IN), 50 ng of cDNA, and 1 μl of each primer (10 μM) (Table 2) in a final volume of 20 μl. The qPCR reactions were performed in a LightCycler 480 Multiwell Plate system (Roche 04729692001) with Optical Adhesive Film (Applied Biosystems Ref 4360954). The following thermal cycles were performed: one cycle of 50°C for 2 min and 95°C for 10 min; followed by 45 cycles of 95°C for 10 sec, 60°C for 10 sec, and 72°C for 10 sec; and then one cycle of 95°C for 5 sec and 65°C for 1 min. qPCR products were qualified according to a melting-curve analysis. Additionally, representative samples were electrophoresed to verify that only a single product (band) was present. RNA-free water was used as a template in the control reactions to determine non-specific primer amplification background levels. Three replicates were performed for each cone opsin gene for each specimen.

Proportional expression of each cone opsin was calculated using the following equation: where T/T is the proportional expression of a given opsin gene i, and E is the amplification efficiency for each pair of opsin primers.

Expression of each cone opsin was first evaluated using Levene’s test and then compared using one-way ANOVA with a post-hoc Tukey HSD when Levene’s test was not significant or using Welch’s ANOVA with a Games-Howell post-hoc test when Levene’s test was significant. Statistical tests were performed in R version 3.6.0 (R Foundation for Statistical Computing, Vienna, Austria).

Results

Identification of cone opsins

The naming of poeciliid cone opsins has varied through time, but was recently reviewed by Sandkam et al. (2018) [28]. Based on that review, we have labelled the cone opsins we identified in western mosquitofish according to the orthologous opsins in guppy. We identified five cone opsin genes from the published draft genome of western mosquitofish [46]: sws1, sws2a, sws2b, rh2-1, and lws-r. We annotated six cone opsins from our draft transcriptome, including sws1, sws2a, sws2b, rh2-1, lws-1, and lws-3. These newly identified cone opsin genes (apart from lws-r, see below) were amplified from cDNA using specific primers (Table 2), with accession numbers MN817658 to MN817663 (Table 1).

The size of β-actin fragments amplified using DNA or cDNA as template differed (Fig 1B), with DNA templates producing larger fragments. However, amplified lws-r fragments were all identical in size no matter which template was employed (Fig 1A). Compared to the cDNA sequence of lws-r in guppy (KX768568), the western mosquitofish lws-r transcript contains an intron and also has a 46-basepair deletion in the coding region (equivalent to positions 1166 to 1211 of the guppy complete lws-r sequence). This intronic retention in the western mosquitofish lws-r transcript would result in an early stop codon during translation (S2 Fig), so it was excluded from further analyses apart from our phylogenetic assessment of the 5’ UTR region of lws.

Fig 1

PCR-amplified fragment sizes of lws-r and β-actin.

(A) lws-r opsin, and (B) β-actin. The initial character of each lane label represents DNA (D) or cDNA (R), and the number represents a specific specimen: 1 to 3 are female specimens collected at Tunghai University; 4 to 6 are male specimens collected from Tunghai University; 7 and 8 are female specimens collected at Linyuan Ocean Wetland Park; and 9 and 10 are male specimens collected at Linyuan Ocean Wetland Park.

Phylogenies of western mosquitofish cone opsins

The cone opsin gene dataset comprised a total of 85 taxa and 1344 aligned nucleotide positions, with 983 variable sites and 848 parsimony-informative sites (S1 Appendix). In our cone opsin ML tree, the sws1, sws2, rh2, and lws families are all respectively clustered with high statistical support (bootstrap values ≥ 70), with the lws family lying basal to these other three cone opsin families (Fig 2). This cone opsin ML tree also demonstrates that all six cone opsins from the western mosquitofish cluster with their corresponding cone opsin genes from other poeciliid fishes with high statistical support. The two western mosquitofish lws genes (lws-1 and lws-3) are clustered together and are paralogous to the corresponding genes of other poeciliid fishes (Fig 2E), and the other four cone opsin genes of the western mosquitofish (sws1, sws2a, sws2b, and rh2-1) are all orthologous to their corresponding genes in other poeciliids (Fig 2B to 2D).

Fig 2

Phylogenetic tree of the coding regions from cone opsin genes based on partitioned maximum likelihood analysis (three partitions).

Branch colors represent different cone opsin families. Phylogenies for each opsin family are shown in panels B to E. Solid circles on branch nodes indicate statistically robust nodes with bootstrap values ≥ 70.

The 5’ UTR dataset of lws genes comprised a total of 64 taxa and 140 aligned nucleotide positions (S2 Appendix), with 111 variable sites and 84 parsimony-informative sites. From our lws 5’ UTR ML tree (Fig 3), it is clear that lws genes from poeciliid fishes can be divided into four strongly supported groups: lws-1 group, lws-2 group, lws-3 group, and lws-r group. The three lws genes possessed by western mosquitofish (lws-1, lws-3, and lws-r) are all orthologous to the corresponding lws genes of other poeciliid fishes. Notably, our genomic (BLAST) and transcriptomic analyses reveal lws-2 to be absent from western mosquitofish.

Fig 3

Phylogenetic tree of the 5’ untranslated region (UTR) of poeciliid lws opsin genes based on maximum likelihood analysis.

Oreochromis niloticus was used as outgroup. Solid circles on branch nodes indicate statistically robust nodes with bootstrap values ≥ 70.

Predicted λmax values of western mosquitofish cone opsins

The λmax values of guppy cone opsins are 353 nm (Sws1), 438 nm (Sws2a), 408 nm (Sws2b), 516 nm (Rh2-1), 571 nm (Lws-1), and 519 nm (Lws-3) [52]. Apart from Sws2a, all western mosquitofish cone opsins have the same ‘10-sites’ composition as their orthologs in guppy, so we infer the λmax values of western mosquitofish cone opsins to be as follows: Sws1 = 353 nm, Sws2b = 408 nm, Rh2-1 = 516 nm, Lws-1 = 571 nm, and Lws-3 = 519 nm. Compared to western mosquitofish Sws2a, guppy Sws2a has an Ala-to-Thr mutation at residue 110 (Ala110Thr) and a Thr-to-Ala mutation at residue 285 (Thr285Ala). Comparing the two orthologous Sws2a opsins from guppy and medaka, we assume that the effect of Ala110Thr on the λmax value is -1 nm. Moreover, Cowing et al. (2002) [55] inserted the Thr285Ala mutation into the Sws2 opsin of Cottus gobio and observed an approximatley -10 nm shift in the λmax value. Accordingly, we infer a λmax value for western mosquitofish Sws2a of 449 (= 438+1+10) nm (Table 3).

Table 3

Ten-site amino acid haplotypes for cone opsin genes in guppy and western mosquitofish.

Scientific name	Opsin	Key site position										Measured λmax (nm)	Predicted λmax (nm)
		110	112	138	180	197	223	227	277	285	308
Oryzias latipes	Lws-A	S	C	I	S	H	L	C	Y	T	A	561a
	Lws-B	S	C	I	S	H	L	C	Y	T	A	562a
Poecilia reticulata	Lws-1	S	C	I	A/ S	H	L	C	Y	T	A	A: 562b/ S: 571b
	Lws-2	S	C	I	P	H	L	C	F	A	A	516b
	Lws-3	S	C	I	S	H	L	C	Y	T	A	519b
	Lws-r	S	C	I	S	H	L	C	Y	T	A
Gambusia affinis	Lws-1	S	C	I	S	H	L	C	Y	T	A		571
	Lws-3	S	C	I	S	H	L	C	Y	T	A		519
Oryzias latipes	Rh2-A	C	Y	Q	A	E	M	H	F	A	S	452a
	Rh2-B	T	T	E	A	E	M	H	F	A	A	516a
	Rh2-C	T	N	Q	A	E	M	H	F	A	A	492a
Poecilia reticulata	Rh2-1	T	T	E	A	E	M	H	F	A	A	516b
	Rh2-2	T	Y	Q	A	E	L	H	F	A	A	476b
Gambusia affinis	Rh2-1	T	T	E	A	E	M	H	F	A	A		516
Oryzias latipes	Sws2a	A	C	M	A	E	L	C	F	A	S	439a
Poecilia reticulata	Sws2a	T	F	M	A	E	L	C	F	A	S	438b
Gambusia affinis	Sws2a	A	C	M	A	E	L	C	F	T	S		449c
Oryzias latipes	Sws2b	C	V	M	G	E	L	C	F	A	S	405a
Poecilia reticulata	Sws2b	C	Y	M	A	E	L	C	F	A	S	408b
Gambusia affinis	Sws2b	C	Y	M	A	E	L	C	F	A	S		408
Oryzias latipes	Sws1	V	L	L	G	E	L	C	F	A	A	356a
Poecilia reticulata	Sws1	V	V	L	G	E	L	C	F	A	A	353b
Gambusia affinis	Sws1	V	V	L	G	E	L	C	F	A	A		353

Measured λmax of guppy refers to reconstituted opsins. Predicted λmax of western mosquitofish refers to predictions based on the orthologous opsin in guppy.

aλmax value from Matsumoto et al. (2006) [56]

bλmax value from Kawamura et al. (2016) [52].

cInferred from the λmax value of the orthologous opsin in guppy and the putative -10 nm effect of the Thr-to-Ala mutation at residue 285 [55] and the putative -1 nm effect of the Ala-to-Thr mutation at residue 110.

Cone opsin gene expression

The rh2-1 cone opsin gene exhibited the greatest expression in all developmental stages we tested and for both sexes, whereas proportional expression levels of both the sws2a and lws-3 opsin genes were always lower than 1% (Fig 4). Apart from lws-3 (ANOVA F2, 16 = 1.3571, p > 0.05), proportional expression of all of the western mosquitofish cone opsin genes was significantly different among larval, adult male, and adult female specimens (Fig 4). One-way ANOVA and post-hoc Tukey HSD revealed that larvae expressed more Sws1 and Sws2a opsins than adults of both sexes (Sws1 F2, 16 = 4.9853, p < 0.05; Sws2a F2, 16 = 13.529, p < 0.05). A Welch’s ANOVA test also supported that expression levels of the sws2b (F2, 9.108 = 6.0534, p < 0.05), rh2-1 (F2, 10.672 = 13.231, p < 0.05), and lws-1 (F2, 9.81 = 16.72, p < 0.05) genes were all significantly different among larval, adult male, and adult female specimens. A post-hoc Games-Howell test revealed that larvae expressed less Sws2b opsin than adult males, less Lws-1 than adults of both sexes, but more Rh2-1 opsin than adults of both sexes.

Fig 4

Box-plot showing the median (—), 25th and 75th percentiles (box), 95% range (|), and outliers (•) of cone opsin proportional expression values for adult and larval specimens.

An asterisk (*) indicates p < 0.05.

Discussion

We identified a total of seven cone opsin genes, including one sws1, two sws2, one rh2, and potentially three lws, in this study. In teleosts, a visual opsin genic tandem array is commonly observed for sws2, rh2, and lws paralogs [13], and the sws2-lws synteny is a synapomorphy from the common ancestor of mammals, birds, reptiles, and fish [57]. However, the western mosquitofish draft genome [46] only revealed that two sws2 genes are tandemly arrayed, indicating that the draft genome may not be complete. Lin et al. (2017) [13] considered one of the lws genes, lws-r, to be apomorphic in the Cyprinodontoidei since it lacks intron II-V but has a highly conserved exon/intron I. The lws-r gene is also thought to have arisen by retrotransposition and is located in an intron of gephyrin (gphn) [26, 58]. Here, we confirm by whole-genome BLAST analysis that the western mosquitofish has the lws-r gene and, as for guppy and green swordtail, its lws-r also contains only one intron and is located in a gphn gene intron.

We identified a further two western mosquitofish lws opsin transcripts (lws-1 and lws-3) by RNA-seq. In our phylogenetic tree of lws 5’ UTRs, these two lws genes cluster separately and with high statistical support with other poeciliid lws-1 and lws-3 genes, respectively (Fig 3). However, the phylogenetic analysis based on the coding regions revealed a very different topology in that these two lws genes grouped together with high statistical support rather than separately clustering with their own orthologous lws genes (Fig 2E). Thus, the lws-1 and lws-3 genes appear to group by species as if independent gene duplication events have occurred in poeciliid fishes. Sandkam et al. (2017) [29] stated that the genomic environment in poeciliid fishes tends to homogenize their lws-1 and lws-3 genes through gene conversion because of their genomic proximity, gene structure, and orientation, but that this gene conversion may be curbed in some lineages where bright red/orange male coloration is critical to sexual selection. Moreover, the distinct topologies of the coding region-based tree and the 5’ UTR-based tree is consistent with gene conversion being more frequent in coding regions than in non-coding ones [62]. Although the draft western mosquitofish genome does not reveal if lws-1 and lws-3 are in close proximity, have a similar gene structure, or are in inverted orientations, our coding region-based phylogenetic tree may indicate that gene conversion has homogenized the western mosquitofish lws-1 and lws-3 genes.

Compared to the guppy (nine cone opsins), western mosquitofish appear to lack two cone opsin genes, i.e., rh2-2 and lws-2. It is possible that these two opsins have been lost from the western mosquitofish genome, which is not unusual among the Cyprinodontoidei [13]. Alternatively, these two opsins do exist in the western mosquitofish genome but they are not (or barely) expressed in the larval/adult stages we sampled.

That the lws-r gene is located in the gphn intronic region raises two questions. First, does lws-r interfere with the expression of gphn and, second, is a hybrid transcript generated [27]? Based on our draft transcriptome, we can confirm that the western mosquitofish expresses gphn as well as the green swordtail, and can also rule out gene fusion between these two genes. However, although western mosquitofish lws-r is present in the draft genome and we successfully amplified it from cDNA, it is absent from our transcriptome assembly. Interestingly, the lws-r sequences from DNA and cDNA template (Fig 1A) were of the same length and both contained the intron region (S2 Fig). This outcome is not likely due to DNA contamination because β-actin fragments amplified from cDNA samples exhibited a single size and were shorter than those of DNA samples (Fig 1B).

In the process of library preparation for RNA-seq, we used poly-T oligo-attached magnetic beads to purify mRNA. Thus, the absence of lws-r in our RNA-seq dataset may be attributable to lws-r transcript not undergoing full post-transcriptional modification so that the poly(A) tail had not been added. Alternatively, transcription of the gphn gene may enable lws-r in guppy and green swordtail to manifest [27, 58]. This “hijack” hypothesis may explain why the premature lws-r, containing the intronic sequence and lacking the poly(A) tail, could be amplified from cDNA template. In our study, cDNA samples were reverse-transcribed using both Hexamer and Anchored Oligo-dT, so the spliced intronic fragments from post-transcriptional gphn modification could be converted into cDNA. It should be noted that we did not detect any “mature” mRNA of lws-r from either larval or adult samples, so it is doubtful that Lws-r exerts a light-sensing function in western mosquitofish. Yokoyama et al. (2014) [59] reported that the sws1 opsin gene of the pearleye (Benthalbella dentata) became a pseudogene about 130 million years ago, but it still produces transcripts and has a low nucleotide substitution rate, perhaps indicating that it retains some unknown functions. It would be interesting to determine if lws-r in the western mosquitofish is a functional visual opsin, has other unknown functions, or its transcript is just a side product of gphn expression.

Apart from Sws2a, we found that the other five identified cone opsins in western mosquitofish possess the same “10-sites” composition as their guppy orthologues. Accordingly, we inferred that these five cone opsins have the same λmax values as defined for guppy. A shift in the λmax value of a visual opsin can alter the evolutionary trajectory of another one. For instance, in many species of flounder, a shift to green of Sws2a has been suggested as having compensated for the functional loss of one of their Rh2 paralogs [20, 60]. This green-shift phenomenon might also occur in western mosquitofish. Compared to guppy Sws2a, the Ala110Thr and Thr285Ala mutations in western mosquitofish Sws2a may induce an 11-nm shift in long-wavelength [55, 56]. By incorporating more closely-related poeciliids in further genomic studies, it should be possible to determine if the greenward shift of Sws2a is a response to loss of the rh2-2 gene from the western mosquitofish genome, as has been reported for Verasper fishes [60].

Our qPCR results reveal that apart from lws-3 with less than 0.05% proportional expression for all developmental stages and sexes, the remaining five cone opsins all display significant differences in proportional expression between larval and adult stages. The significant ontogenetic decrease in sws1 expression and increase in lws-1 suggest that the visual spectrum of adult western mosquitofish is more sensitive to longer wavelengths than that of larvae. UV sensitivity is critical for detecting zooplankton [61-66]. Many fishes display an ontogenetic long-wavelength shift in spectral sensitivity when they no longer feed on zooplankton [67-70]. Since Gambusia fishes exhibit an ontogenetic dietary shift from micro to large invertebrates [33, 71], the decline in ultraviolet sensitivity of western mosquitofish is consistent with its ontogenetic dietary shift.

The lws-1 allele of western mosquitofish is 180S, which is a more red-sensitive allele (180S, λmax = 571 nm, compared to 180A λmax = 562 nm; Table 3). Sandkam et al. (2015) [10] demonstrated that female guppy carrying the 180S allele prefer males with more abundant red/orange coloration than those carrying the 180A allele. Therefore, allele type and expression levels of lws-1 could not only explain why adult male western mosquitofish are strongly attracted to the yellow spot around the urogenital opening of mature females [34], but also why female western mosquitofish are more inclined to breed with male guppies than conspecific males [72].

In conclusion, we have verified that the western mosquitofish has seven cone opsin genes, including one sws1, two sws2, one rh2, and three lws, but whether the lws-r gene is functional or not remains unclear. The predicted λmax value for Sws1 is 353 nm, it is 449 nm for Sws2a, 408 nm for Sws2b, 526 nm for Rh2-1, 571 nm for Lws-1, and 519 nm for Lws-3. The mutations at residues 110 and 285 of Sws2a have likely induced a green-shift in the λmax value of that opsin, which may be associated with loss of the rh2-2 gene. We have also revealed an ontogenetic red shift in the visual spectrum of western mosquitofish. The ontogenetic decline in UV-sensitivity could be correlated with a dietary shift. Furthermore, lws-1 allele type and expression level appears to reflect the asymmetric negative reproductive interference reported between sympatric guppy and western mosquitofish. The comprehensive understanding of cone opsins in western mosquitofish from this study will facilitate further explorations of the mechanisms, such as developmental and/or adult plasticity of opsin expression, that enable western mosquitofish to flourish better in turbid waters compared to other fishes [73].

The lighting spectrum of the full-spectrum LED light deployed above the aquarium in this study.

(JPG)

Click here for additional data file.

Nucleotide sequence alignment of the western mosquitofish lws-r fragment from the published draft genome [51], the sequences we generated from DNA and cDNA templates in this study, and the published lws-r complete gene and cDNA sequences of guppy (KX768568).

(PDF)

Click here for additional data file.

(DOCX)

Click here for additional data file.

(DOCX)

Click here for additional data file.

(JPG)

Click here for additional data file.

(JPG)

Click here for additional data file.