Genome-Wide Association Study for Squalene Contents and Functional Haplotype Analysis in Rice.

Squalene is an isoprenoid compound that acts as the intermediate metabolite in cholesterol synthesis. Squalene is not very susceptible to peroxidation, and it quenches singlet oxygen in the skin, which is caused by UV exposure and other ionizing radiation sources. Squalene is a precursor to phytosterol synthesis, and it has been widely studied for its ability to reduce oxidation, cancer activity, and cholesterol levels. We performed a genome-wide association study for squalene in rice using 1.6 million high-quality SNPs extracted from 295 accessions' resequencing data. The candidate gene locus Os09g0319800-an orthologue of terpene synthase in Arabidopsis-showed up as the most likely candidate gene amongst the identified loci. Nucleotide variations in the promoter were associated with squalene content variations within the japonica group. The results of this study can provide clues for understanding the mechanisms of squalene biosynthesis in rice.

Introduction

Squalene, which is beneficial to human health, is contained in rice that is consumed by a large population. Genetic studies related to the squalene content of rice are valuable for high squalene-containing rice breeding. Squalene is a natural triterpene hydrocarbon known for its key role as a precursor in phytosterol or cholesterol synthesis, as well as all steroids in plants, animals, and humans.[1] Squalene can be obtained from many plants, including olive oil, palm oil, wheat-germ oil, soybean oil, rice, Amaranthus grains, and even the livers of sharks.[2−5]

Squalene has received much scientific attention due to the beneficial effects of its natural products in the context of human health. It has been used in cosmetics, medicine, moisturizing ingredients, functional food applications, as a immunologic adjuvant, and for protection against carcinogenesis, cardiovascular diseases, and cancer.[6−8] Squalene has been identified in vitro and in animal studies to have anticancer, antioxidant, drug carrier, detoxification, and skin hydration effects.[9] Devarenne et al.[10] reported that the highest level of squalene synthase is found in the apical meristem of tobacco. Compounds that perform important biological functions structurally related to squalene also exist in nature.[11] For example, animals use quinyl groups to form ubiquinone side chains, in which ubiquinone coenzyme Q10 (CoQ10) has an antioxidant effect in the human body.[12] Similar to squalene, substances that require a prenyl group include carotene; vitamins A, K, D, and E; tocopherols; tocotrienols; cyclic terpenoid compounds; and dolichol.[13] The formula for squalene structures is C30H50, and it is enzymatically collapsed into three interconnected, closed, six-carbon rings attached to a five-carbon ring, which includes a prenyl side-chain. This squalene synthesis starts with the production of 3-hydroxy-3-methylglutaryl coenzyme A (HMG CoA). From the onset, the reduction of HMG CoA leads to the formation of mevalonate.[13] Squalene also plays an important role as a precursor to phytosterol, which has a positive human physiological effect. The three phytosterols that are commonly found in plants, are β-sitosterol, campesterol, and stigmasterol.[14] They are essential natural plant components that are rich in vegetable oils, nuts, seeds, and grains and exist in plant sterol-enriched margarine.[15,16] Phytosterols have also been suggested to have anti-oxidative, anti-inflammatory, and anticarcinogenic benefits, and their cholesterol-lowering abilities are most widely researched.[17] Numerous studies indicate that plant sterols hinder cholesterol absorption in the intestines and thereby lower the levels of total low-density lipoprotein (LDL) and plasma cholesterol.[17] The intestines absorb plant sterols poorly (0.4–3.5%) and plant stanols even more so (0.02–0.3%).[18] Campesterol in brassinosteroids has been identified as a precursor for biosynthesis in developing seeds, a phytohormone class that manages integral embryonic steps.[19] In addition, phytosterols from rice bran appeared at over 40%, impacting their use in the food industry and for human health purposes.[20] Apart from their effects on lipid metabolism, plant phytosterols can affect other metabolic processes. While concentrations of campesterol are approximately three times higher than those of sitosterol, the focus of most studies has remained on the metabolic effects of sitosterol. Phytosterols are a group of triterpenes from lipids (C30), present in all plant tissues.[21,22] Phytosterols in cultivated rice were found to possess the highest concentration of sitosterol, whereas stigmasterol and campesterol were detected at lower levels.[23] According to Kumar et al.,[24] β-sitosterol was found in N22 (drought tolerant) rice seedlings.

The fact that squalene is included in rice is important because more than half of the global population relies on rice, which is a primary staple in several parts of the world, particularly in Asia.[25,26] Rice’s nutritional qualities represent a significant factor in the consumer choice, especially in developing countries, where rice is the primary source of some essential nutrients.[27] On a global scale, rice provides energy (21%) and protein (15%) to humans. Improving its valuable nutrient contents, which would have important health benefits for humans, has long been a target of rice breeding programs.

However, only a limited number of genetic studies have addressed squalene content in rice. Most research has been conducted using the plant model, Arabidopsis. Fortunately, today, many techniques have been developed to search for squalene-releated genes by correlating phenotypical traits with single nucleotide polymorphisms (SNPs) in the genome. The genome-wide association study (GWAS) is one of the techniques that is utilized to efficiently produce analysis of the genetic diversity of complex traits’ architecture and provide valuable initial insights into candidate genes or loci for subsequent validation.[28] Han et al.[29] also found that the GWAS represents an effective means for the exploration of the genetic architecture of phenolic compounds in both cultivated and Tibetan wild barley. In addition, Volante et al.[30] utilized GWAS to study japonica rice fields under both permanent flooding (PF) and limited water (LW) conditions. The GWAS was presented to map genes/loci resistant to rice blast in open-access rice diversity panel 1 (RDP1).[31] Hence, the dissection of plant biosynthetic pathways has experienced rapid advancement due to a wide range of valuable tools: whole-genome sequencing; high throughput physical mapping; global gene expression analysis; and metabolite profiling in a variety of organisms. Combining these approaches allows for fast access to and the development of essential biochemical information in model systems, as well as the efficient bridging of the information/research gap between agricultural crops and engineering pathways, i.e., those for micronutrients, in global staple crops. Wang et al.[32] performed the GWAS of vitamin E for methyltransferase (OsγTMT) and α-tocopherol, which were found in 13 rice candidate genes. Furthermore, Li et al.[33] used the GWAS for SNP identification and insertion/deletions (indels) linked to α-tocopherol content in kernels of maize.

We investigated the squalene content of 295 rice core sets, performed the GWAS using SNPs from resequencing, and found candidate genes that were associated with squalene content. Of these, Os09g0319800, an orthologue of terpene synthase in Arabidopsis, was found to be the most possible candidate gene amongst the identified loci. Our GWAS on squalene will contribute to understanding the mechanism of nutrient biosynthesis and providing a foundation to future studies on breeding rice with high squalene content.

Results and Discussion

Phenotypic Evaluation and Correlation among Individuals

We observed significant phenotype variation between the japonica and indica varieties in terms of the tochromanol content. We considered whether there was a similar pattern of correlation between the ecotype and the nutrient content in terms of squalene and its derivatives, the phytosterols. It was surprising to discover that there were distinctive differences between the japonica and indica groups with respect to their squalene and phytosterol contents, except for campesterol, with p < 0.01. The other three components all had p values lower than 0.001 (Figure ).

Figure 1

Statistical analysis of squalene and phytosterol contents between japonica and indica rice from 295 accessions. **: significant difference at p < 0.01, and ***: significant difference at p < 0.001.

Population Structure

The 295 accessions and population structure were identical to those identified in ref (32). For the 295 accessions, four subgroups, temperate japonica, tropical japonica, aus, and indica, were separated with some admixed accessions, where K = 4 (Figure S1). At K = 2, the 295 accessions were separated into two groups represented by indica and japonica. The japonica group was divided into temperate japonica and tropical japonica at K = 3. Temperate japonica was divided into two parts, one composed of bred varieties at K = 4. The aus type was separated from indica at K = 5, and the four subgroups were finally well separated from each other. For the 295 accessions together, the same patterns were observed as in the 137 accessions. It is worth mentioning that our temperate japonica varieties were divided into two parts at K = 4, which indicated a relatively complicated genetic background within our japonica varieties.

Principal Component Analysis (PCA)

Based on phylogenetic tree and structure results, the four major ecotype groups separated well from each other.[32] Utilizing 1.6 million high-quality SNPs, we performed a PCA. According to the parallel analysis, the scree plot suggested a component number of two. Based on the first two PCs, we observed clear subpopulation structures. The two subpopulations, indica and japonica, were well differentiated, with clear clusters based on the first two PCs and with the admixture accessions situated between the two groups.[32] Furthermore, temperate and tropical japonica clustered separately, as did aus and indica. For the phenotype, we examined the correlation between the different subgroups using squalene and phytosterol contents. We applied the R package “psych” for a PCA of squalene and phytosterol contents in 295 accessions.[34] Parallel analysis suggested that the number of factors and components were three and two, respectively (Figure S2). We therefore performed the PCA with two components. Based on squalene and phytosterol contents, all rice accessions in our study clustered into two groups, composed of indica and japonica subspecies (Figure S3). These genetic and phenotypic variations suggested that genetic factors underlay squalene and phytosterol content-level differences between the indica and japonica groups.

GWAS of Squalene Content in indica, japonica, and the Entire Population using 295 Accessions

We studied genome-wide association with squalene content, an important precursor of phytosterols. The GWAS was performed using the 137 accessions and 158 Korean-bred varieties together. We used the same 1.6 million high-quality SNPs without missing data across all this study’s accessions. Similarly, to determine the existence of different allelic networks associated with trait variations within the differing subpopulations, we conducted the GWAS to analyze association signals in indica (72 total accessions), japonica (217 total accessions), and entire rice populations. Since most accessions from the 158 Korean-bred varieties (temperate japonica: n = 133, tropical japonica: n = 3, indica:n = 21, and admixture: n = 1) belonged to the japonica type, more than 130 japonica were added into the panel, while only 21 indica were added. We assumed that the 295 accessions’ association panel was closer to a japonica panel and therefore japonica varieties’ mapping power was increased. For the entire population, the genome-wide significance threshold was set as p < 1 × 10–6 considering the sample size. There were three significantly associated signals—located on chromosomes 6, 9, and 10—in the entire population panel (Figure ). In the japonica population, the genome-wide significance threshold was set as p < 1 × 10–5, and two associated signals—which had identical genomes located to signals identified from the entire population on chromosomes 9 and 10—were found. We did not find significant association signals from the indica panel (Figure ). Only squalene GWAS results showed significant association signals in this study. The candidate gene locus Os09g0319800—an orthologue of terpene synthase in Arabidopsis—was the most possible factor among the identified loci. Squalene is an isoprenoid compound, which acts as an intermediate metabolite for cholesterol and other steroid synthesis.[9] Squalene is not very vulnerable to peroxidation and quenches singlet oxygen species in human skin, thereby protecting it from lipid peroxidation due to exposure to UV and other ionizing radiation sources.

Figure 2

Genome-wide association study of squalene content in the entire (295 accessions), japonica, and indica populations using a compressed mixed linear model. The −log10(p) values from a genome-wide scan are plotted against the position on each of the 12 chromosomes. (a–c) Manhattan plots for squalene content detected in the entire, japonica, and indica populations, respectively.

We investigated significantly associated signals on chromosomes 6, 9, and 10. Signals located on chromosome 6 were only identified in the entire population and appeared to be responsible for phenotype variations amongst the indica and japonica groups. However, the false discovery rate (FDR)-adjusted p value was too high to be treated as an authentic candidate gene locus (Table ). Therefore, we excluded association signals on chromosome 6 for further analysis. However, signals on chromosomes 9 and 10 were authentic according to the results. We found certain gene loci, which could be used for further analysis.

Table 1

Genome-wide Significant Association Signals of Squalene Content in the Entire Population using Compressed Mixed Linear Model

chromosome	positiona	Pvalue	MAFb	FDRc	gene	annotation
6	25 766 225	7.42 × 10–4	0.403846	0.99458	Os06g0635300	GTLd
6	27 534 120	5.24 × 10–5	0.110577	0.99458	Os06g0665900	TDPe
9	9 083 379	9.55 × 10–7	0.387500	0.021728	Os09g0319701	CHPf
9	9 097 404	7.63 × 10–7	0.383333	0.020960	Os09g0319800	TSDPg
10	10 429 128	1.30 × 10–6	0.316667	0.039441	Os10g0347000	X8-DPh
10	10 397 716	2.15 × 10–6	0.254167	0.039441	Os10g0346300	CDPi
10	10 417 578	2.88 × 10–6	0.254167	0.039748	Os10g0346600	SVRPj

Position: the position of SNP.

MAF: minor allele frequency.

FDR: false discovery rate.

GTL: gastric triacylglycerol lipase.

TDP: thioredoxin domain-containing protein.

CHP: conserved hypothetical protein.

TSDP: terpene synthase-like domain-containing protein.

X8-DP: X8 domain-containing protein.

CDP: cupredoxin domain-containing protein.

SVRP: similar to vacuolar sorting receptor 4 precursor.

Identification of Putative Nucleotide Polymorphisms for the Squalene Content in Rice

Association signals derived from the results for the 295 accessions’ panel, in which two of the most significant signals were on chromosomes 9 and 10 (Figure ). The peak SNP on chromosome 9 was chr09:9097404, and chr10:10429128 on chromosome 10. We screened ∼50 kb around the peak SNP region for candidates, using the International Rice Genome Annotation database. Two gene loci around the peak SNP on chromosome 9: Os09g0319701 and Os09g0319800. No detailed information on Os09g0319701 has been reported, and it is annotated as a conserved hypothetical protein. Os09g0319800 is a terpene synthase-like domain-containing protein, involved in squalene’s biosynthesis pathway.[35] The chr09:9097404 peak is located in Os09g0319800 and is associated with squalene biosynthesis. Therefore, we chose Os09g0319800 as our candidate gene. We identified three gene loci around the SNP on chr10:10429128: Os10g0346600, Os10g0347000, and Os10g0347800. Os10g0346600 is annotated as a vacuolar sorting receptor 4 precursor, which is related to calcium ion binding in Arabidopsis.[36] Os10g0347000 is annotated as an X8 domain-containing protein, which displayed functions as glycosyltransferases.[37] Os10g0347800 is annotated as a forminlike protein, associated with actin filament formation and regulation in Arabidopsis.[38] Unfortunately, none of the gene loci on chromosome 10 showed a direct connection with squalene biosynthesis, and thus, only Os09g0319800 was chosen for further analysis.

Figure 3

Association mapping results and genomic locations of significant polymorphisms for squalene in the 295 accessions. The dashed line indicates significantly associated regions on chromosomes 9 and 10. Candidate gene names are highlighted in bold.

Comparative Sequence Analysis for Os09g0319800

We sought to identify if the homologous genes of Os09g0319800, from the near range, exists in indica, which is another subspecies of japonica, and wild rice (Oryza genus). We also extended the search to the monocotyledons (class Liliopsida) and dicotyledons (class Magnoliopsida)

First, we identified whether Os09g0319800 exists in indica through BLAST and a public homologs database.[39] The nucleotide BLAST querying the genomic DNA sequence of Os09g0319800 to the nucleotide collection (nr/nt) database (http://blast.ncbi.nlm.nih.gov) showed that the best hit in indica was of E value = 4 × 10–77 and identity = 81.68%, but query coverage = 14%. This result meant that the gene corresponding to Os09g0319800 was not found in the indica genome in genomic DNA. The amino acid BLAST query of the non-redundant protein sequence (nr) database with the amino acid sequence of Os09g0319800 revealed hypothetical protein OsI_15038 with E value = 7 × 10–97, identity = 64.44%, and query coverage = 73%.[40] However, unlike Os09g0319800, which is located on chromosome 9, the gene was located on chromosome 4. In the amino acid BLAST query of Os09g0319800, the japonica genes were also hit, including Os04g0178300 in chromosome 4 (E value = 4 × 10–98, identity = 60.47%, and query coverage = 83%). Thus, OsI_15038 is considered to be a paralogous gene of Os09g0319800. In another public database, Plants Ensembl (http://http://plants.ensembl.org), the reconciliated phylogenetic trees of the homologous genes of Os09g0319800 showed no homologs in indica (Figure S4).[41] Therefore, Os09g0319800 is assumed to be a gene that does not exist in indica. It is interesting to note that the orthologs of Os09g0319800 exist in many wild rice plants (Oryza rufipogon, Oryza nivara, Oryza barthii, and Oryza glumipatula) (Figure S4).

Next, we identified the othologous genes of Os09g0319800 among the plants in the class Liliopsida and Magnoliopsida using Plants Ensembl. The result showed that Os09g0319800 also had a total of nine high-confidence orthologous genes among four species in Liliopsida and for the species in Magnoliopsida, including BnaAnng04000D from Brassica napus, which was present[42] (Table S2). In Liliopsida, all of the species were wild rice (Oryza genus) except Eragrostis tef, a cereal crop in Africa, that had recently completed the whole-genome sequencing and phylogenetic tree reconstruction and been identified as a plant belonging to Chloridoideae.[43]

Functional Haplotype Analysis for Os09g0319800

Using our resequencing data, we generated the promoter region and coding region haplotypes for Os09g0319800. We identified a total of six haplotypes combining the promoter and coding region sequences (Figure ). The indica type was absent from the haplotype analysis, since all indica accession sequences were missing in this genomic region. No SNP was present in the coding region of japonica varieties compared to the reference sequence. Two types of coding regions, constituted by a total of nine variant positions—seven of which lead to amino acid substitution—were identified in aus varieties and some admixtures. The squalene content varied between different haplotypes, types 1–3 composed of japonica type varieties had much higher squalene content than accessions in types 4 and 5, which represented the aus type and some admixtures. There were noticeable squalene content variations between types 1 and 2, and only one SNP was found in the promoter part between the two types. Type 3 carried five SNPs in the promoter region and showed significant squalene content variations compared to other japonica varieties. Therefore, we think that promoter region SNPs in types 2 and 3 could be causing phenotypic variations. Nucleotide variations in the promoter and coding regions of aus type, and some of the admixtures, could be responsible for squalene content differences compared to the japonica group.

Figure 4

Nucleotide polymorphism of Os09g0319800 promoter and coding regions, and their relationship with squalene content. (a) Haplotypes of Os09g0319800 promoter and coding regions. (b) Squalene content among haplotypes. (c) Squalene content among ecotypes. The error bars represent one standard deviation. TEJ: temperate japonica, TRJ: tropical japonica, AUS: aus, and Admix: admixture.

Asia, especially East and Southeast Asia, are the major rice cultivation areas in the world. Rice is the major energy source for most people in the developing countries of the Southeast Asian Nations (ASEAN). Nutrient intake from rice represents the total dietary intake amount for people in these countries. It was surprising to find significant differences in squalene, phytosterols, and tocochromanol contents between the different ecotypes. According to our results, compared to the indica varieties, japonica varieties contained a much higher phytonutrients’ content, which are indispensable for human health. Genetic factors must be underlying this phenotypic variation, since there is a substantial difference between japonica and indica varieties regarding their genetic background. Furthermore, noticeable variations among subgroups (temperate japonica, tropical japonica, indica, and aus) were observed. Dissection of variations in genes related to nutrient biosynthesis in rice can facilitate the process of improving valuable compounds through breeding activities, which could yield incredible benefits for human health.

Experimental Section

Plant Materials

A total of 4406 varieties of rice exist worldwide, collected by the National Genebank of the Rural Development Administration (RDA-Genebank, Republic of Korea) and developed into a core set of 166 rice accessions using the program PowerCore.[44] Of these accessions, 137 were selected for the resequencing project (Table S1), in addition to 158 Korean-bred varieties. The experiment was conducted during the rice-growing seasons of 2016 and 2017, at Kongju National University. Each accession was transplanted in two rows, with 15 cm separating each plant, and 30 cm separating each row. Field management followed the standard agricultural practice, with each field replicated twice in a completely randomized design. Young plant leaves were sampled and then stored at −80 °C, until their use for genomic DNA extraction using the CTAB method.[45]

Analysis of Squalene and Phytosterol Content

Squalene and phytosterol samples were prepared for analyses using the procedures described by Bhandari et al.,[46] albeit with some modifications. Fresh rice samples (5 g) were mixed with 20 mL of ethanol and 0.5 g of ascorbic acid as an antioxidant, followed by shaking in a water bath at 80 °C for 10 min. Then, 600 μL of 44% KOH were added in the saponification step, followed by shaking at 80 °C for 18 min. The tubes were incubated on ice. Then, 10 mL of n-hexane and 10 mL of distilled water were added and mixed, and the sample was centrifuged for 5 min at 3,000 rpm, whereby we collected the upper hexane layer. This process was repeated three times, and the collected hexane layers were pooled and washed with 10 mL of distilled water a total of three times. The samples were then passed through anhydrous Na2SO4 to remove water, concentrated using a rotary evaporator, dissolved in 1 mL of isooctane, and injected into a GC instrument (GC 400 Varian) with a flame ionization detector and a capillary column (CP-SIL 8CB, 30 m × 0.25 mm, 0.4 μm film thickness). The injector and detector (FID) temperatures were set at 290 °C, and the injection volume was 1 μL, with a split ratio of 1:20 under a constant column flow (1.0 mL min–1) of helium gas. The oven temperature was initially set at 220 °C for 2 min; increased to 290 °C with 5 min–1 and held for 14 min; and then increased to 310 °C with 10 °C min–1 and held for 18 min.

Statistical Analysis

Statistical analyses were performed utilizing Microsoft Excel 2010 and SPSS (version 20.0). The Pearson correlation coefficients between squalene and phytosterol contents, as well as the correlations between the 2-year phenotype data, were analyzed using SPSS (version 20.0). Statistical differences in squalene and phytosterol contents between the japonica and indica groups were determined using Sigmaplot (version 12).

Whole-Genome Resequencing and Single Nucleotide Polymorphism (SNP) Calling

Each accession used genomic DNA from a single plant. Total genomic DNA was extracted from young leaf tissue with the DNeasy Plant Mini Kit (Qiagen). Rice varieties were sequenced to obtain genomes with ∼8× coverage to detect SNPs. We generated a total of 3900 million raw reads, which were aligned against the rice reference genome (IRGSP 1.0) for SNP identification and genotype calling. In total, 12.23 million SNPs and indels were obtained through SNP calling. The genotype data set was generated for GWAS. Only SNPs with minor allele frequency >5% that contained genotype calls for all 137 accessions were used. Finally, GWAS utilized ∼1.6 million high-quality SNPs.

Population Structure Analysis

We conducted principal component analysis (PCA), constructing a phylogenetic tree based on the 1.6 million high-quality SNPs.[47,48] The population phylogenetic tree was constructed using PHYLIP (version 3.68), and MEGA7 (version 7.0)[49] was employed to display the tree.[50,51] Furthermore, we used the program ADMIXTURE (version 1.3.0), which estimates individual ancestry and admixture proportions, assuming that K populations exist based on the maximum-likelihood method.[52] We also applied PCA to the phenotype data from 2013 using the R package psych.[34]

Genome-Wide Association Study (GWAS)

We performed the GWAS utilizing the software genome association and prediction integrated tool (GAPIT) (version 3.0).[48] We then selected the compressed mixed linear model (CMLM), which uses a kinship (K) matrix as the variance–covariance matrix between individuals. An optimized number of principal components (PCs) was obtained from the PCA results. We performed GWAS for 295 accessions from the 137 core collection and 158 Korean-bred varieties.

Functional Haplotype Analysis for Os09g0319800

Whole-genome resequencing data for 295 rice accessions were imported into the TASSEL (version 5.0)[53] software to extract nucleotide polymorphism. IRGSP 1.0 was employed as a reference genome for functional haplotype analysis of the Os09g0319800 gene on chromosome 9—positioned at 9 096 555–9 101 312 (+strand). Sequences were aligned using MAGA7 (version 7.0).[49] Haplotype diversity was analyzed utilizing DnaSP (version 6) software.