Genetic signature of strong recent positive selection at interleukin-32 gene in goat.

OBJECTIVE: Identification of the candidate genes that play key roles in phenotypic variations can provide new information about evolution and positive selection. Interleukin (IL)-32 is involved in many biological processes, however, its role for the immune response against various diseases in mammals is poorly understood. Therefore, the current investigation was performed for the better understanding of the molecular evolution and the positive selection of single nucleotide polymorphisms in IL-32 gene.
METHODS: By using fixation index (FST ) based method, IL-32 (9375) gene was found to be outlier and under significant positive selection with the provisional combined allocation of mean heterozygosity and FST . Using nucleotide sequences of 11 mammalian species from National Center for Biotechnology Information database, the evolutionary selection of IL-32 gene was determined using Maximum likelihood model method, through four models (M1a, M2a, M7, and M8) in Codeml program of phylogenetic analysis by maximum liklihood.
RESULTS: IL-32 is detected under positive selection using the FST simulations method. The phylogenetic tree revealed that goat IL-32 was in close resemblance with sheep IL-32. The coding nucleotide sequences were compared among 11 species and it was found that the goat IL-32 gene shared identity with sheep (96.54%), bison (91.97%), camel (58.39%), cat (56.59%), buffalo (56.50%), human (56.13%), dog (50.97%), horse (54.04%), and rabbit (53.41%) respectively.
CONCLUSION: This study provides evidence for IL-32 gene as under significant positive selection in goat.

INTRODUCTION

The goat (Capra hircus) was first domesticated about 10,000 years ago [1,2]. China has a prolonged history breeding, production and domestication of various breeds of goats [3]. These native goat breeds are spread from the harsh, cold and dry Qinghai-Tibet Plateau to the humid and warm Southern regions. Chinese goats are classified into five groups on the basis of geographical and ecological distribution and 25 breeds of goats that are reared and domesticated in Southeast China [2]. However, until now their phylogeny and origin is poorly described.

Cytokines help to communicate between different immune cells during both acquired and innate immune reactions [4]. Cytokines modulate host immune responses to various kinds of infections, inflammations, injuries and diseases of unknown etiology [5]. The genetic variants at many loci of cytokines in humans, as well as in many livestock species, have been reported to be susceptible towards some particular pathogens, in addition to immunopathology of the host [6]. Therefore, the genetic variation at some specific loci of cytokines might cause the significant differences among different individuals in their immune reactions or resistance to certain pathogens.

Interleukins are immunomodulatory proteins that assist to elicit various responses in tissues and cells and modulate activation, growth and differentiation during different immunological reactions [7,8]. Interleukin (IL)-32 is recognized as pro-inflammatory cytokine that plays significant role in many biological processes. It is reported that NK4 (natural killer cell transcript-4) stimulates numerous pro-inflammatory cytokines, namely tumor necrosis factor-α as well as IL-8, therefore the name of NK4 was changed as IL-32 [9,10]. It is worth mentioning here that IL-32 has no similarities with other identified ILs, however it possess characteristic pro-inflammatory properties [11]. Interestingly, IL-32 gene is present in most of the mammalian species excluding rodents. The absence of IL-32 gene in rodents restricts in vivo animal studies and further progress on IL-32 research for medical applications [10,12]. Therefore, this limitation can be solved by evolutionary and positive selection studies of IL-32 gene in animals other than rodents.

Different mutations are significant for evolution to arise through natural selection of different traits resulting in population changes due to environmental and genetic influences. Conversely the changes in single trait are frequently channeled by its co variation with different traits. These variations may result in different directions in a phenotypic space with a small to no conflict even when individual traits do differ [13,14]. Therefore, consideration of the inherent factors that structure the variation within a population, and how this variation evolves, is critical for understanding the probable evolution of that population.

Several researchers have reported the significance of molecular markers for investigating genetic studies among various livestock species [15]. Molecular markers have been identified as a leading tool for calculating genetic variation and for determining the genetic gap among strongly related populations of ruminants due to their random division in the genome, higher polymorphism and availability of automated genotype scoring [16,17].

Efforts to identify the signature of selection can be helpful to determine which artificial selection leads to positive genetic gain. Since genetic selection can enhance the advantageous allele’s frequency as well as sequence variations to the allele, positive selection signals can be found by identifying the sequence variations near suitable mutations [18,19].

New approaches are needed to increase knowledge of the limited genetic factors from the association studies. One attempt might be to identify selection signatures in which positive selection has occurred by artificial selection on the meat quality. Because the selection could increase beneficial allele frequency and further augment allele frequencies of sequence variants in the vicinity, signals for the positive selection can be identified by harboring distribution of the sequence variants around a favorable mutation [18,19].

The evolutionary positive selection studies on goat IL-32 in scientific literature are still scarce. Therefore, the present investigation was performed to investigate the genetic and nucleotide variations to better understand the molecular evolution and the positive selection of IL-32 gene in Chinese local goat breeds. This study aimed to provide novel information about the nucleotide variations and recent positive selection in IL-32 gene in goat as compared to other vertebrate species using nucleotide sequences from National Center for Biotechnology Information (NCBI) database.

MATERIALS AND METHODS

Ethics statement

All the experimental procedures of this study were approved by the Law of Animal Husbandry in People’s Republic of China (Dec 29, 2005). The experimental protocols for collection of ear tissue samples were reviewed and permitted by the Biological Studies Animal Care and Use Committee of National Animal Husbandry Service, Hubei, PR China. Measures were taken to reduce any distress to experimental animals during collection of samples.

Animal selection and extraction of DNA

For this study, 200 goats of four local breeds (Enshi black, Hybrid white yellow, Nanjing yellow and Yichang white) from southern districts of China were selected. The genomic DNA extraction from ear tissues of these goats was extracted using Genomic DNA kit TIANamp (TianGen, Beijing, China).

Genotyping and sequencing

The genomic DNA for the extension of eight anticipated single nucleotide polymorphisms (SNPs) loci was used for screening and characterization of these SNPs. Using the SeqMan program the sequences of these identified SNPs were aligned and genotyped in 200 goats using Matrix-assisted laser desorption/ionization time of flight mass spectrometry (MALDI-TOF) assay (SquenomMassARRAY: emoji:, BGI company, Shenzhen, China).

Frequency distribution analysis

Selection effects were assessed and analyzed by Lositan software [20]. The fixation index F and p values for each locus were calculated by allele frequencies based on heterozygosity. Simulations were comprised of four populations, 200 individuals, 8 loci, 17 SNPs and an expected 0.102 F value. This technique enhanced the divergent selection by identifying the outliers with F values more than estimated, overcoming for heterozygosity [21]. On real data, 100,000 simulations were employed to build the datasets of population. Quantiles were supposed for provisional F joint distribution against mean heterozygosity with confidence limits 95%. Loci that expressed outside the simulated neutral distribution with a characteristic differentiation behavior (i.e. F) were finally regarded as outliers.

Evolution analysis

Coding sequences of 11 mammalian species were retrieved from Gene Bank and concluded sequences of proteins were aligned by MEGA6.0 program [22] using parameters of default alignment, followed by manual adjustment. Maximum likelihood technique was employed for the determination of positive selection and amino acid sites under selection. For this study, four models (M1a, M2a, M7, and M8) were used in Codeml program of phylogenetic analysis by maximum liklihood (PAML) package in two steps. First, likelihood ratio test (LRT) was applied to measure the positive selection for the presence of sites with ω>1. Then, comparison was done with a null model to obtain this, that did not permit sites with ω>1 and a more generalized model that did. In the second step, LRT comparison was made with a null model M7, assuming a beta distribution B (p, q) for x ω (in the interval 0<ω<1, where 0 indicated complete constraint and 1 was the expectation under no selection pressure), and another M8 model using an additional class of sites with ω estimated was included. Consequently, a proportion of sites (p0) came from the beta distribution B (p, q) and the remaining sites (p1 = 1–p0) had a ω ratio estimated from data that was >1 in M8 model. It was assumed that M8 model fit the sequences better than M7, if the LRT was statistically significant (p<0.05), as a result the positive selection was accepted in M8 model.

Phylogenetic analysis

The sequences of the nucleotides of goat IL-32 gene were compared with that of sheep, dog, human, rabbit, cattle, buffalo, bison, horse, camel, and cat for evolutionary study. After retrieving the nucleotide sequences from NCBI, MEGA6 software package the Neighbor Joining method was used for constructing phylogenetic tree. The genetic sequence analysis of similarity among eleven species with nucleotide sequences of IL-32 gene was done using ClustalW software [23].

Codon-based nonsynonymous- synonymous and neutrality test

Maximum Likelihood computations of dN and dS were conducted using HyPhy software package (MEGA6, Hachioji, Tokyo, Japan) [24]. There were 273 positions in the final dataset. The nonsynonymous- synonymous (dN-dS) ratio is useful for determining the codons that had undergone positive selection, where dS denotes a synonymous substitutions number in each site, and dN a nonsynonymous substitutions number in each site, where a higher number of nonsynonymous substitutions is an indication of positive values. For neutrality test, the variance of the difference was computed using the analytical method. Analyses were conducted using the Nei-Gojobori method (Nei and Gojobori, 1986). There were of 91 positions in the final dataset. Both the analyses involved 11 nucleotide sequences. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA6 [22].

RESULTS

The sequence description among the four goat breeds unveiled 17 SNPs from eight loci and these SNPs were genotyped for further analysis. F based technique recognized IL-32 gene loci as in selection sweep in the studied breeds. Positive selection of IL-32 gene by frequency distribution (FDIST) method, via the outliers approach, IL-32 gene was present outside the 95% confidence interval of the conditional combined allocation of mean heterozygosity and F by Lositan FDIST examination (Figure 1). IL-32 gene was found to be significantly in the region resulting from (p<0.05) positive selection (Table 1).

Figure 1

Candidate gene interleukin-32 under positive selection keeping the 95% confidence interval. FST, fixation index; He, heterozygosity.

Table 1

Locus, heterozygosity (He) and fixation index (FST) for each of 17 genotyped SNPs

Locus	Het	FST	p value
IL-3	0.292	0.0087	0.526
IL1-α	0.487	−0.005	0.251
IL1-α	0.362	−0.001	0.354
IL1-β	0.479	−0.006	0.236
IL1-β	0.481	−0.006	0.244
IL1-β	0.412	−0.001	0.334
IL-4	0.429	0.04	0.884
IL-8	0.517	0.057	0.92
IL-8	0.534	0.013	0.589
IL-31	0.497	−0.005	0.251
IL-33	0.492	−0.009	0.168
IL-33	0.201	0.06	0.951
IL-33	0.408	0.013	0.606
IL-32(1158)	0.438	0.043	0.893
IL-32(2017)	0.465	−0.008	0.154
IL-32(7638)	0.485	0.022	0.735
IL-32(9375)	0.487	0.0861	0.985*

SNP, single nucleotide polymorphisms; IL32, interleukin. p*, simulated FST

Positive selection of IL-32 gene by FDIST analysis

By keeping 95% confidence interval in Lositan FDIST analysis, IL-32 (9375) gene was present in outlier area with the provisional combined allocation of F and mean heterozygosity (Figure 1). Furthermore, positively selected IL-32 gene was significantly (p<0.05) present in corresponding positive selection area (Table 1).

Evolutionary analysis of positive selection

Excess of non-synonymous over synonymous substitutions indicates molecular evidence for positive selection. Two paired models (M1–M2; M7–M8) were used for the database analysis and ωs and log-likelihood values were obtained (Table 2). Likelihood logs were used for the LRT determination, with the intention to investigate whether the substitute models (M2:M8) were better than the null models (M1:M7) in fitting sequences. In the first LRT, M1 model provided a log-likelihood of LnL = −2,445.220, with the estimated value ω0 = 0.317 and M2 model presented LnL = −2,444.825, with the estimated value ω3 = 1.749. The statistical value for this test was 2ΔLnL = 0.791, with df = 2 and p<0.001. We neglected the M1 and accepted the M2 model. According to M2 model, 64% sites were found under purifying selection, 28% sites were under neutral evolution and only 6% sites were protected by positive selection. The values of ωs were found as 0.336, 1.000, and 1.749, respectively. In the second LRT, the comparison was made between M8 and M7 models (more parameter-rich models), and we found 4.318 as the value of 2ΔLnL. After its comparison with x2 distribution value (df = 2), we discarded the M7 model, but accepted the M8 model (p< 0.001). Hence, M8 proved a better fit than the M7 model for the sequences. According to M8 model, the ωs value of 87% sites were best fit in the beta distribution B (1.798, 0.124), while only 12% sites showed ω1 = 1.628. Hence, the evolutionary study of the IL-32 gene in 11 mammalian species indicated the gene was under purifying selection where previously no positive selection had been reported.

Table 2

Results of positive selection tests for interleukin-32 gene

Model	Parameter estimates1)	lnL	LRTs	Positive selection sites
Model 1	p0 = 0.638, p1 = 0.361ω0 = 0.317, ω1 = 1.000	−2,445.220	0.791	NEB: Not allowed
Model 2	p1 = 0.649, p2 = 0.287, p3 = 0.062ω1 = 0.336, ω2= 1.000, ω3= 1.749	−2,444.825		BEB:
Model 7	p = 1.257, q = 1.077	−2,443.730	4.318	NEB: Not allowed
Model 8	p0 = 0.875, p = 1.797, q = 2.009 (p1 = 0.124) ω = 1.628	−2,441.571		BEB: 2C, 6A, 10L, 11K, 23Q, 30S, 49S, 55D 75E, 89S

LnL, log likelihood; LRT, likelihood ratio; NEB, naive empirical bayes; BEB, bayes empirical bayes.

ω, the ratio of nonsynonymous-to-synonymous substitutions; p0, purifying selection; p2, neutral selection, p3, positive selection substitution ratio for all sites; p and q, the β distribution parameters.

p>95%;

p>99%.

Phylogenetic relationship of IL-32 gene among species

We collected complete coding sequences of IL-32 gene from 11 different mammalian species while keeping the goat as reference species from the public database, Gene Bank. MEGA6.0 software was used for the alignment of all the obtained sequences along with their respective amino acids (23). Before employing this software for phylogenetic analysis, manual editing of all the obtained sequences was performed. The caprine IL-32 gene was found closely associated with that of sheep, camel, bison, cattle, buffalo, cat, horse, human, rabbit, and dog (Figure 2). The genetic similarity of nucleotide sequences of IL-32 gene of goat was more closer to sheep (96%) than camel (58%) bison (91%) and cat, buffalo, and human (56%), dog (50%), horse (54%), rabbit (53%), respectively Table 3.

Figure 2

Phylogenetic relationship of IL-32 gene in different species.

Table 3

Coding sequences similarity analysis and values of pair wise comparison of interleukin-32 gene among different species

	Goat	Sheep	Bison	Camel	Buffalo	Human	Cat	Cow	Rabbit	Dog	Horse
Goat	-	96.54	91.97	58.39	56.50	56.13	56.59	55.95	53.41	50.97	54.04
Sheep	96.54	-	91.97	50.87	52.81	43.94	55.41	51.43	49.35	49.13	55.19
Bison	91.97	91.97	-	48.91	52.55	45.99	54.50	50.12	50.36	49.88	54.99
Camel	58.39	50.87	48.91	-	60.87	53.76	57.56	56.90	51.14	52.38	53.62
Buffalo	56.50	52.81	52.55	60.87		55.91	60.98	81.43	53.66	55.49	54.27
Human	56.13	43.94	45.99	53.76	55.91	-	58.28	54.29	52.47	53.76	57.42
Cat	56.59	55.41	54.50	57.56	60.98	58.28	-	59.29	55.43	63.18	57.79
Cow	55.95	51.43	50.12	56.90	81.43	54.29	59.29	-	48.81	56.90	55.00
Rabbit	53.41	49.35	50.36	51.14	53.66	52.47	55.43	48.81	-	52.71	49.90
Dog	50.97	49.13	49.88	52.38	55.49	53.76	63.18	56.90	52.71	-	55.42
Horse	54.04	55.19	54.99	53.62	54.27	57.42	57.79	55.00	49.90	55.42	-

Similar values in the under triangle (%).

Determination dN-dS ratio and neutrality test

Using coding sequences, non-synonymous/synonymous rate ratios (dN/dS) were calculated and positively selected sites (having values >1) are shown in Table 4 Each codon, and estimates of the number of inferred synonymous (s) and nonsynonymous (n) substitutions are also presented along with the numbers of sites that are estimated to be synonymous (S) and nonsynonymous (N). Using the coding sequences of species, the neutrality analysis was performed and the results indicated that the <0.05 values were considered as significant Table 5.

Table 4

Maximum likelihood analysis of natural selection codon-by-codon

Codon#	Codon Start	Triplet	Syn (s)	Nonsyn (n)	Syn sites (S)	Nonsyn sites (N)	dS	dN	dN-dS
18	871	ACA	1	6	0.851	2.128	1.174	2.818	1.643
20	892	CGC	1	6	0.715	2.176	1.397	2.756	1.358
22	898	GTC	1	5	0.953	2.046	1.048	2.443	1.395
33	967	GAC	1	5	0.914	2.085	1.093	2.397	1.304
35	973	TCC	0.5	7.5	0.386	2.362	1.292	3.174	1.882
41	994	GAC	0	4	0.418	2.553	0	1.566	1.566
46	1021	TTC	0	6	0.227	2.766	0	2.169	2.169
51	1036	CAG	0	5	0.325	2.394	0	2.088	2.087
67	1102	CGC	0	4	0.538	2.280	0	1.754	1.754
71	1114	GAG	0	3	0.298	2.354	0	1.273	1.273
85	1171	CAG	0	4	0.335	2.363	0	1.692	1.692
88	1183	AGC	0.7	7.3	0.468	2.473	1.423	2.964	1.540
89	1186	CTC	2	9	0.873	2.115	2.289	4.254	1.965
91	1192	CAG	0.7	7.3	0.422	2.338	1.579	3.135	1.556

dN, nonsynonymous; dS, synonymous. dN-dS values >1 indicating significance.

Table 5

Codon-based test of neutrality for analysis among species

	Goat	Sheep	Bison	Camel	Buffalo	Human	Cat	Cow	Rabbit	Dog	Horse
Goat	-	−0.125	−0.466	0.371	−0.347	−0.363	0.137	0.114	−1.688	−0.365	0.144
Sheep	0.901	-	−0.428	−0.042	−0.039	−0.967	−0.76	−0.35	−2.249	−1.218	−0.583
Bison	0.642	0.669	-	−0.124	−1.355	−0.039	−1.008	−0.669	−3.076	−1.22	−0.459
Camel	0.711	0.967	0.901	-	−3.17	−2.028	−2.918	−3.214	−3.32	0.291	−3.089
Bifflao	0.73	0.969	0.178	0.002*	-	−2.114	−2.555	−0.668	−2.19	−0.728	−4.236
Human	0.717	0.335	0.969	0.045*	0.037*	-	−2.143	−1.681	−2.834	−0.939	−3.076
Cat	0.892	0.449	0.316	0.004*	0.012*	0.034*	-	−1.127	−1.855	−2.257	−3.016
Cow	0.909	0.727	0.505	0.002*	0.505	0.095*	0.262	-	−2.094	−0.996	−2.642
Rabbit	0.094	0.026*	0.003*	0.001*	0.03*	0.005*	0.066*	0.038*	-	−0.947	−3.419
Dog	0.716	0.226	0.225	0.772	0.468	0.350	0.026*	0.321	0.345	-	−0.594
Horse	0.886	0.561	0.647	0.002*	0.00	0.003*	0.003*	0.009*	0.001*	0.553	-

Represent significance respectively.

DISCUSSION

The ongoing development in the significant directories of genetic variation has made novel suggestions for the discovery of positive selection targets, which ultimately would be helpful to explain the roles of drift and selection in evolutionary development. Furthermore, signatures of positive selection hamper the genome regions that functionally play key roles. Therefore, discovering such genomic regions would assist in the identification of genetic deviations that interrupt these functional genomics regions. Consequently, the history of a population can be created from the chronological events and information of the targets for positive selection.

The genetic basis of various traits in many species has been studied by candidate gene approach. The identification of the genes that play a key role in phenotypic variation in natural livestock populations can provide new information about the evolutionary process and enables direct selection [25]. In the current investigation, we identified IL-32 gene by operating the mean heterozygosity and F because this approach can identify genetic discrimination for each locus.

Previously, many genes related to positive selection were identified by LRT. The genes inferred by LRT are considered true cases of adaptation that might have been involved in positive selection [26]. In the current study, we used likelihood logs for the LRT determination. Both the alternative models M2 and M8 were implemented. According to M2 model, 65% sites were found under purifying selection, 28% under neutral evolution and only 6% were protected by positive selection. According to M8 model, the ωs value of 87% sites best fit in the beta distribution B (1.797, 0.12486), while only 12% sites showed a ω1 = 1.628 value and were concerned with positive selection. Hence, the evolutionary study of the IL-32 gene in 11 species of vertebrates found the gene to be under purifying selection and previously no positive selection has been reported.

Due to the co-evolution of pathogens with host, the genes related to the host immune system are under constant selection pressure for adaptation. In this way, valuable mutations are greatly favored and harmful mutations are promptly eliminated [7,27, 28]. In this study, IL-32 was found to be a gene of the immune system under positive selection in many mammalian species. Keeping this outcome in mind, with the help of maximum-likelihood approaches, we performed this investigation for signatures of positive selection at IL-32 gene in different mammals.

Directional selection increases allele frequency of suitable mutations and their immediately adjacent sequence variations, which stay in the genome as signatures of selection [29-31]. Previous studies support our findings in which a few quantitative trait loci linked to economically beneficial traits corresponded to strong selection regions in genome as determined by Linkage disequilibrium [32] or allele frequency [19]. In this study, we used log likelihood method of the allele frequencies of sequence variations and Lositan software against identified SNPs for detection of selection signatures. These signatures of selection can improve the penetration of the gene and can also play a role in many other physiological functions as well. Focus on such studies with a large sampling size can enhance the statistical effectuality for identifying signatures of selection.

The phylogeny of goat IL-32 is poorly described in scientific reports. Therefore, complete coding sequences were compared in different mammalian species and it was found that the IL-32 gene of goat shared 96% identity with sheep. After construction of a phylogenetic tree, it was found that goat IL-32 gene had a close resemblance with sheep IL-32 gene and therefore they were placed in the same phylogenetic group. Hence, our study provided novel information about the nucleotide variations in IL-32 gene for the genetic control of immune responses in goat breeds as well as in other analyzed mammalian species.

Bustamante observed in (2005) [33] that genes having lower non synonymous mutations were significantly associated with Mendelian diseases, moreover, several genes known to be possible targets of positive selection were implicated in cancer causing processes [34]. Genes that are related to complex diseases, like cytochrome P450 (CYP3A) [35], asthma (IL-13 [36], IL4 [37,38] and IL1A [39]), cardiovascular disease (matrix metallopeptidase 3 [MMP3] [40] and angiotensinogen [AGT] gene [41] type 2 diabetes calpain-10 [CAPN10] [36]) have also been reported to have signatures of selection. So signatures of selection might play a role in disease mapping via genome wide analysis.

Positive selection might play an important role in sustaining the immune response in the host against various infectious agents and environmental toxicants [28]. In the current investigation, we determined the positive selection for IL-32 along with few of the identified codons related to dN-dS. Incomplete information about the functions of specific codons in IL-32 of the mammalian species limits the absolute understanding of our current interpretation. However, for complete comprehension of the function of the observed differences in mammalian IL-32, further structural and functional studies involving mutagenesis and crystallographic approaches are needed.

CONCLUSION

As a whole, all the approaches employed to investigate positive genetic selection and nucleotide variations in this study lead to better understanding of the molecular evolution and positive selection of SNPs in IL-32 gene. However, the distribution of genetic variation as a result of selection could be supported by demographic records of the population.