Association between the Polymorphisms of fads2a and fads2b and Poly-Unsaturated Fatty Acids in Common Carp (Cyprinus carpio)

Simple Summary Fishes are the major dietary source of polyunsaturated fatty acids (PUFAs) for humans. The limited availability of PUFAs derived from fish represents a critical bottleneck in food production systems, one that numerous research institutions and aqua-feed companies in this field are trying to overcome. This problem could be minimized by select-bred fish to be capable of more effectively producing endogenous PUFAs. Fatty acid desaturase 2 (fads2) is one of the rate-limiting enzymes in the synthesis of PUFAs. The common carp, one of the most important food and ornamental allo-tetraploid fish, encodes two fads2 genes (fads2a and fads2b). The PUFA contents among common carp individuals were numerous, suggesting that there might exist polymorphisms in fads2a and fads2b. However, selective breeding of common carp with high PUFA contents was hindered due to a lack of effective molecular markers. This study investigated the contents of PUFAs in common carp and identified polymorphisms in the CDS regions of fads2a and fads2b. The association study identified three cSNPs associated with the PUFA contents and suggested that fads2b might be the major gene responding for common carp PUFA contents. These cSNPs would be potential markers for future selection to improve the PUFA contents in common carp. Abstract Fatty acid desaturase 2 (fads2) is one of the rate-limiting enzymes in PUFAs biosynthesis. Compared with the diploid fish encoding one fads2, the allo-tetraploid common carp, one most important food fish, encodes two fads2 genes (fads2a and fads2b). The associations between the contents of different PUFAs and the polymorphisms of fads2a and fads2b have not been studied. The contents of 12 PUFAs in common carp individuals were measured, and the polymorphisms in the coding sequences of fads2a and fads2b were screened. We identified five coding single nucleotide polymorphisms (cSNPs) in fads2a and eleven cSNPs in fads2b. Using the mixed linear model and analysis of variance, a synonymous fads2a cSNP was significantly associated with the content of C20:3n-6. One non-synonymous fads2b cSNP (fads2b.751) and one synonymous fads2b cSNP (fads2b.1197) were associated with the contents of seven PUFAs and the contents of six PUFAs, respectively. The heterozygous genotypes in both loci were associated with higher contents than the homozygous genotypes. The fads2b.751 genotype explained more phenotype variation than the fads2b.1197 genotype. These two SNPs were distributed in one haplotype block and associated with the contents of five common PUFAs. These results suggested that fads2b might be the major gene responding to common carp PUFA contents and that fads.751 might be the main effect SNP. These cSNPs would be potential markers for future selection to improve the PUFA contents in common carp.


Introduction
Fish provide polyunsaturated fatty acids (PUFAs) for human. PUFAs are necessary for important biological functions in human, such as regulating lipid metabolism, stimulating growth and development, exerting anticancer and antiaging effects, enforcing immuneregulation, promoting cardiovascular health, and aiding weight loss [1][2][3]. PUFAs in the diet also help to elevate serum peroxides and antioxidant reserves and to lower plasma triglycerides, thus helping to ameliorate cardiovascular diseases [4].
The common carp (Cyprinus carpio) is cultured commercially worldwide, generating numerous geographical populations. Domestication has also accelerated the phenotypic diversity of the common carp [21]. This fish is able to convert dietary C18 FAs to PUFAs [6,11]. Due to its allo-tetraploiy [22], the common carp encodes one more fads2 gene than most diploid freshwater fish. Two fad cDNAs, fads2a and fads2b, were previously identified [23]. The PUFA contents among common carp individuals were numerous [24], suggesting that there might exist polymorphisms in fads2a and fads2b. Indeed, in common carp var. Jian, one intronic SNP in fads2a and one intronic SNP in fads2b were identified [23]. However, whether these two intronic SNPs were associated with the PUFA contents has not been studied. The question of whether there are polymorphisms in these two genes in other varieties and the association between polymorphisms and the PUFA contents also require studying. Furthermore, which of these two genes is the major gene responsible for PUFA biosynthesis is still unknown.
Genomic markers have been demonstrated to be efficient in the representation of the associated phenotypes and aquaculture breeding [25][26][27]. Selective breeding of common carp with high PUFA contents is hindered due to a lack of effective molecular markers. To answer the above questions and to develop SNP markers for selective breeding of PUFAs, in this study, we measure the PUFA contents and identified more polymorphisms in common carp fads2a and fads2b than before. The association between certain genotypes and the contents of multiple PUFAs in the common carp population are performed. Three significantly associated SNPs are useful in the molecular-assisted selection of high PUFA contents in common carp.

Ethics Statement
This study was conducted following the recommendations of the Care and Use of Animals for Scientific Purposes established by the Animal Care and Use Committee of the Chinese Academy of Fishery Sciences (ACUC-CAFS). Before collecting the tissues, all fishes were euthanized in MS222 solution.

Sampling
To increase the diversities of HUFA contents, in 2018, we collected the juvenile common carp including the var. Huanghe (HHC, at Donge, Shandong province, China), the var. Furui (FRC), and the var. Jian (JC, at Wuxi, Jiangsu province, China) (Supplementary Figure S1). The samples were maintained at the experimental fish farm of Chinese Academy of Fishery Sciences (Fangshan, Beijing, China) and fed with the same commercial diet (Tongwei, China) in the culture ponds. This commercial diet contained fish meal, fish oil, colza oil, and soybean meal. In 2019, a total of 269 one-year-old common carp were randomly selected, including 124 FRC individuals, 47 HHC individuals, and 98 JC individuals. The liver and muscle of each individual were dissected, immediately frozen in liquid nitrogen, and stored at −80 • C.
2.3. RNA Extraction, cDNA Synthesis, Sequencing, and Genotyping RNA was extracted from the liver of each individual using TRIzol®(100 mg sample /1 mL TRIzol, Invitrogen, Waltham, MA, USA). The concentration of RNA was determined using the NanoDrop 2000 spectrophotometer (Thermo Scientific™, Wilmington, DE, USA), and its quality was determined by 1% agarose gel electrophoresis. Next, 1 µg RNA was reverse transcribed to cDNA using the SuperScript III RNase H-Reverse Transcriptase Kit (Invitrogen, Carlsbad, CA, USA) with oligo-dT. Based on the reference full-length sequence of fads2a and fads2b in common carp (GenBank accessions: MK852165 and MK852166), genespecific primers were designed to amplify the complete coding sequence (CDS) regions of these two genes (Table 1).   (Table 1), 1 min at 72 • C, and a final extension of 10 min at 72 • C. Finally, the PCR products were sequenced with the Sanger method. Since the CDS of each gene was less than 1350 bp, two ends of one round of Sanger sequencing covers the complete CDS.
Since fads2a and fad2b have high nucleotide sequence similarities at 89.86%, to avoid misalignments and putative artificial coding SNPs (cSNPs), all of the amplified sequences were aligned to fads2a and fads2b using blastn. For one sequence expected to be from fads2a, if it had a higher alignment score for fads2b than fads2a, then this sequence should be aligned to fads2b to call for SNPs. The same operation was performed on all sequences expected from fads2b. To detect SNPs and to perform genotyping in each sample, the sequences were aligned to their corresponding reference genes in the novoSNP software [28] with the position ranging from 100 bp to 700 bp. If one locus had only one peak, this locus was considered a homozygote. If the base in this locus was different from the reference, then this locus was a cSNP. For one locus with two peaks, the genotype was auto-detected with the software and considered a heterozygous site. If one site had an F-score over 30 measured by novoSNP, it was of high quality.

Genetic Diversities of Common Carp fads2a and fads2b
We investigated the genetic distances and structures among three varieties with all retained genotypes in fads2a and fads2b. The genetic distances were estimated using Principle Component analysis (PCA) in the Tassel 5.0 software [29]. The first two eigenvectors were plotted. The population structure among three varieties we analyzed using Admixture 1.3.0 [30]. K was set from 2 to 6. The population structures and the genetic compositions of each individual were visualized by pophelper v2.3.1 [31].
As three varieties were grouped together in the PCA analysis and had similar structures, we combined these varieties into one population. The observed heterozygosity (Ho), expected heterozygosity (He), genotype frequency, and allele frequency in the population were calculated using the Genepop software [32]. If the frequency of one SNP in the population was higher than 0.04 (at least 5 out of 122 individuals with fads2a genotype information), it was retained in the downstream analysis. The polymorphism information content (PIC) for each retained SNP locus was estimated using PICcalc 0.6 [33]. The effects of SNPs to the coding sequences were classified into synonymous substitution, non-synonymous substitution, stop gain, and stop loss.

PUFA Content Comparison
Before the following association study, we first examine whether there existed significant differences in the PUFA types and contents among three varieties. With the contents of 12 PUFAs, the distances measured by the PCA analysis were estimated with Tassel 5.0.
The Spearman correlation coefficient of the contents of any two PUFAs in all individuals was calculated using the R 'cor.test' function. Using the TBtools software [35], we clustered 12 PUFAs on the basis of the correlations using the 'Median method' and visualized all correlation coefficients using the R function 'heatmap'.

Associations of Polymorphisms in fads2a and fads2b with the Contents of 12 PUFAs
Associations between the contents of each PUFA and the genotypes were studied using the mixed linear model (MLM) and the analysis of variance (ANOVA), respectively. The MLM model, including the PCA and the kinship matrix, which is MLM (PCA + K), was applied into the association study with Tassel 5. This method was widely used in the association study between the polymorphisms in candidate genes and phenotypes [36][37][38]. The percentage of phenotypic variation explained by each marker was calculated with Tassel 5.
For each SNP locus, we grouped the individuals on the basis of their genotypes. Then, we performed pairwise group comparison with ANOVA to estimate whether there existed significant content difference in each PUFA between the two groups. All p values were corrected using the false discovery rate (FDR) method for multiple hypothesis testing.
As for each PUFA, one SNP locus was considered significantly associated with the content if it had a p value < 0.05 in the MLM method and a FDR-corrected p value < 0.05 in the ANOVA. In each gene, the pairwise linkage disequilibrium (LD) between any two SNPs was measured by D' calculated with PLINK v1.90 [39] and visualized with the R package LDheatmap [40]. The haplotype block was defined if the average D' of any two compared SNPs in this region was greater than 0.6.

Genotyping Common Carp fads2a and fads2b
The CDS sequences of fads2a and fads2b were highly similar, with a nucleotide identity of 89.86% and a protein identity of 92.06%. The fads2a CDS spanned 12 exons, and the fads2b CDS covered 12 exons. For fads2a, base contents of A, G, C, and T in the sequence were 25.0%, 26.2%, 24.6%, and 24.2%, respectively; the A + T content (49.2%) was lower than the G + C content (50.8%). For fads2b, the average A, G, C, and T base contents in amplified sequences were 25.4%, 26.1%, 24.5%, and 24%, respectively; the A + T content (49.4%) was lower than the G + C content (50.6%). All of these data suggest high similarity between the CDS of these two genes. After amplification and validating the real source genes, we successfully sequenced the entire CDS of fads2a in 122 individuals and the entire CDS of fads2b in 237 individuals.
We obtained ten fads2a cSNPs and 27 fads2b cSNPs (Supplementary Table S1). After filtering low-frequency cSNPs with a MAF threshold of 0.04, five cSNPs, including one nonsynonymous SNPs and four synonymous SNPs, were identified in fads2a from 122 samples, corresponding to 14 genotypes (Table 2). These cSNPs were distributed in four exons. Likewise, eleven cSNPs, including six non-synonymous SNPs and five synonymous SNPs were identified in fads2b from 237 samples ( Table 2), being distributed over exons 1 to 11. In these loci, 30 genotypes were detected. Taking the CDS lengths into consideration, the SNP loci only accounted for 0.37% in fads2a and 0.82% in fads2b, suggesting low polymorphism in these two genes. However, fads2b had higher polymorphisms than fads2a.

Genetic Diversities of fads2a and fads2b
PCA based on the fads2a genotypes clustered three common carp varieties together ( Figure 1A). The first two principal components explained 30% and 24% of the total genetic variances. A similar phenomenon of clustering together was observed with the fads2b genotypes ( Figure 1B). The first two principal components explained 60.1% and 11.14% of the total genetic variances, respectively. These data proved that these three varieties had similarity genetic background in both fads2a and fads2b. These findings were further supported by the population structure analysis with fads2a cSNPs and fads2b cSNPs, respectively (Supplementary Figure S2). In the admixture plots with different K values, three varieties had similar genetic components. There was no variety-specific SNP in these two genes, as evidenced by both PCA and population structure analysis. Therefore, we combined these three varieties into one population in the following analysis. The minor allele frequencies (MAF) of fads2a cSNPs ranged from 0.05 to 0.28. The MAFs of fads2b cSNPs ranged from 0.04 to 0.47. No significant MAF distribution differences between these two genes were observed (Mann-Whitney U test, p value = 0.377). Ho in fads2a ranged from 0.1 to 0.51 and He ranged from 0.09 to 0.4. PIC analysis indicated that the fads2a.288 locus displayed a low polymorphism level (PIC < 0.1) compared with the other four loci. Among the fads2b cSNPs, Ho ranged from 0.02 to 0.38 and He ranged from 0.07 to 0.49. Two loci displayed low polymorphism levels (PIC < 0.1) compared with the remaining nine loci. All fads2a cSNPs and fads2b cSNPs significantly deviated from the Hardy-Weinberg equilibrium ( Table 2, chi-square p value < 0.05).

Diversities of PUFA Contents in Common Carp
We obtained the contents of 12 PUFAs in 269 common carp samples (Supplementary Table S2). PCA analysis with the contents of 12 PUFAs demonstrated that there were no variety-enriched PUFAs (Supplementary Figure S3), similar to no variety-specific SNPs in these two genes. Hence, in the following phenotype comparison and association study, we also grouped three varieties into one population. The n-6 PUFA contents (44.68 ± 16.81%) were higher than the n-3 PUFA contents (6.76 ± 3.58%). The most abundant PUFA was C18:2n-6 (24.33 ± 11.03%), followed by C20:4n-6 (8.17 ± 4.96%) and C18:3n-6 (6.36 ± 6.75%) ( Table 3). The content distributions of 12 PUFAs were not subject to the normal distribution (Kolmogorov-Smirnov test p value < 0.05, Supplementary Figure S4). A correlation matrix of pairwise content comparisons between two PUFAs ( Figure 2) exhibited two major clades. The first clade including PUFAs with 18 carbons. The components in the second clade mainly included at least 20 carbons (except C18:3n-6), which were the derivates of the first clade. Hence, it is reasonable that there exist negative correlations between these two clades. Based on the correlations among the components, the second clade was further classified into two groups: n-6 PUFAs (except C20:3n-3) and n-3 PUFAs. The second-clade classification was consistent with the classification based on the location of the last double bond with reference to the terminal methyl end of the molecule (n-3 and n-6) [11]. The clustering suggested that correlations among the PUFA products were higher than those between the products and substrates. 3 PUFAs. The second-clade classification was consistent with the classification based on the location of the last double bond with reference to the terminal methyl end of the molecule (n-3 and n-6) [11]. The clustering suggested that correlations among the PUFA products were higher than those between the products and substrates.

Three cSNPs Associated with the PUFA Contents
Out of five fads2a cSNPs, only one synonymous cSNP (fads2a.624) was found to be significantly associated with the content of C20:3n-6 by both MLM (p value of 0.025) and ANOVA (corrected p value of 0.0037, Table 4, and Figure 3A). The homozygous genotype (GG) of this locus was associated with higher C20:3n-6 content than the heterozygous genotype (AG), with a fold change of 1.50.

Comm per sh
Comm nus si

Three cSNPs Associated with the PUFA Contents
Out of five fads2a cSNPs, only one synonymous cSNP (fads2a.624) was found to be significantly associated with the content of C20:3n-6 by both MLM (p value of 0.025) and ANOVA (corrected p value of 0.0037, Table 4, and Figure 3A). The homozygous genotype (GG) of this locus was associated with higher C20:3n-6 content than the heterozygous genotype (AG), with a fold change of 1.50.   Among 11 fads2b cSNPs, two cSNPs were significantly associated with the contents of multiple PUFAs by both MLM and ANOVA. One non-synonymous cSNP (fads2b.751), resulting in the conversion of valine to isoleucine, was associated with the contents of four n-3 PUFAs (C18:3n-3, C20:3n-3, C20:4n-3, and C22:5n-3) and three n-6 PUFAs (C20:3n-6, C22:4n-6, and C22:5n-6). In this locus, two genotypes including the homozygous genotype (GG) and the heterozygous genotype (AG) were observed while the homozygous genotype of AA was not identified. Interestingly, in all observed associated PUFAs, the individuals with the heterozygous genotype had significantly higher contents than those with the homozygous genotype (Table 4, Figure 3B-H). The fold changes of mean PUFA contents between two groups ranged from 1.4 to 2.1.
The other cSNP (fads2b.1197) was synonymous and associated with the contents of three n-3 PUFAs (C20:3n-3, C20:4n-3, and C22:5n-3) and three n-6 PUFAs (C18:2n-6, C20:3n-6, and C22:5n-6). In this locus, two genotypes including the homozygous genotype (CC) and the heterozygous genotype (CT) were observed while the homozygous genotype of TT was not identified. Likewise, in all observed associated PUFAs, the individuals with the heterozygous genotype had significant higher contents than those with the homozygous genotype (Table 4, Figure 3I-N). The fold changes of mean PUFA contents between two groups ranged from 1.5 to 2.0.

Discussion
The limited availability of PUFAs derived from fish represents the critical bottleneck in food production systems, one that numerous research institutions and aqua-feed companies in this field are trying to overcome. Increasing PUFA in fish may have important impacts on human heath, especially in a context where the traditional sources of PUFA, wild fish, and fish oils, are limited by the stagnation of fisheries. Although producing aquaculture fish with more PUFA could be achieved by feeding different fish ingredients [41,42], this goal could be minimized by select-bred fish to be capable of producing endogenous PUFAs more effectively.
The fads2 gene is a rate-limiting enzyme catalyzing the desaturation of PUCAs, and the overexpression of salmon fads2 gene in zebrafish could improve PUFA contents [43]. Two fads2 genes were identified in the common carp and had differential expression patterns during the development from embryo to larval [23]. Our study intends to find SNPs linked to the PUFA contents in common carp, using a candidate gene approach, by searching associations between polymorphisms in the coding region of fads2a/fads2b and the PUFA contents. The PUFA content diversities in common carp suggest the possible existence of SNPs in fads2a/fads2b. Since the genic regions of fads2a and fads2b were longer than 11,500 bp [23], it is laborious to amplify the whole genic regions. Therefore, we investigated the presence of SNPs in the CDS regions of fads2a and fads2b. The lower amplification efficiency of fads2a than fads2b was possibly due to the low gene expression of fads2a in liver.
Although under long-term selection (both natural and artificial) common carp acquired genetic diversity and phenotypic diversity, the three varieties studied had no significant differences in the genetic diversities of these two genes and the PUFA content diversities. Three varieties of common carp were grouped into one population to increase the sample size, the genetic diversities of fads2a and fads2b, and the PUFA content diversities. It is still unknown which gene makes more contributions to the PUFA biosynthesis. More associated PUFAs and higher explained percentage of phenotypic variations in fads2b than fads2a suggest that fads2b might be the major effect gene for PUFA biosynthesis. More functional studies on these two genes are required for further validation [12,13,15,16].
The non-synonymous fads2b cSNP might improve desaturase efficiency on the n-3 and n-6 substrates, leading to higher PUFA contents. We also observed one synonymous cSNP associated with the PUFA contents. One possible reason is the tight linkage between the non-synonymous cSNP and the synonymous cSNP, and the synonymous cSNP is a possible biomarker representing the non-synonymous cSNP. The other possible reason is that the synonymous cSNP might affect the DNA methylation level of fads2b. In humans, the DNA methylation level in the fads1/2/3 gene cluster was linked with genetic variants and desaturase activities and thus affected fads2 gene expression [44,45]. We predicted two potential CpG islands (from 157 to 731 bp, and from 1015 to 1237 bp) on the fads2b CDS region using the method described by Gardiner-Garden and Frommer (1987) with the parameters of GC content over 50% and of the observed/expected ratio over 60% [46]. The synonymous cSNP of fads2b.1197 is located on the second CpG island. Whether this cSNP affects the fads2b DNA methylation level and then its expression requires future study.
We observed the eterozygote advantage with the PUFA contents. The heterozygote advantage was widely observed in the economic traits including the growth rates [47]. In two cSNPs, fads2b.751 and fads2b.1197, the heterozygous genotypes were associated with higher PUFA contents than the homozygous genotypes, suggesting that the heterozygosity of fads2b might increase the desaturase activity. Steer et al. [48] reported heterozygotes of fads2 in humans to be increasingly capable of synthesizing docosahexaenoic acid from alphalinolenic acid, thereby suggesting that heterozygotes can achieve nutritional adequacy.
We identified cSNPs in fads2a and fads2b associated with the PUFA contents. These cSNPs will be applied into future marker-assisted selective breeding of common carp to improve PUFA contents. The question of whether there exist polymorphisms in the promoters, miRNA binding sites, and other functional elements of fads2a and fads2b will be studied in the future.

Conclusions
We detected many polymorphisms in the coding regions of fads2a and fads2b in common carp and found higher polymorphisms in fads2b than fads2a. We identified the contents of 12 PUFAs in common carp muscle, which exhibited diversities among individuals. The association studies between the polymorphisms of these two genes and the PUFA contents in common carp identified three cSNPs linked to the PUFA contents. The gene fads2b might be the major gene responding for common carp PUFA contents, and fads.751 might be the main effect SNP. Taken together, our findings highlight the importance of polymorphisms in fads2a and fads2b, critical to the biosynthesis of n-6 and n-3 PUFA levels in common carp. These cSNPs would be potential markers for future selection to improve PUFA contents in common carp.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/ani11061780/s1, Table S1: The genotypes of fads2a and fads2b, Table S2: The contents of 12 PUFAs in 269 samples, Figure S1: The images of three varieties of common carp used in this study, Figure S2: Population structures using the genotypes of fads2a and fads2b, Figure S3: PCA plot clustering three common carp varieties with the contents of 12 PUFAs, Figure S4: Comparing the content distributions of 12 PUFAs with the normal distribution.

Data Availability Statement:
The data presented in this study are available in supplementary material.

Conflicts of Interest:
The authors declare no conflict of interest.