Mapping and Meta-Analysis of Height QTLs in Soybean  

Ya'nan Sun1 , Huaihai Luan2 , Zhaoming Qi1 , Dapeng Shan4 , Chunyan Liu1,2 , Guohua Hu2,3 , Qingshan Chen1
1. Northeast Agricultural University, Harbin, 150030, P.R. China
2. The Crop Research and Breeding Center of Land-Reclamation, Harbin, 150090, P.R. China
3. The National Research Center of Soybean Engineering and Technology, Harbin, 150030, P.R. China
4. HeilongjiangAcademy of Agricultural Sciences, Suihua Institute, Suihua, 152052, P.R. China
Author    Correspondence author
Legume Genomics and Genetics, 2012, Vol. 3, No. 1   doi: 10.5376/lgg.2012.03.0001
Received: 19 Dec., 2011    Accepted: 13 Jan., 2012    Published: 16 Jan., 2012
© 2012 BioPublisher Publishing Platform
This article was first published in Molecular Plant Breeding in Chinese, and here was authorized to translate and publish the paper in English under the terms of Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Preferred citation for this article:

Sun et al., 2012, Mapping and Meta-Analysis of Height QTLs in Soybean, Legume Genomics and Genetics, Vol.3, No.1 1-7 (doi: 10.5376/lgg.2012.03.0001)

Abstract

Plant height is one of very complicated quantitative traits of soybean yield and the phenotype data is variable to environment. Since molecular marker technology appeared, many studies about soybean height QTL have been reported. In this study, we analyzed QTLs data of soybean height between 2006 and 2008 with a recombination inbred lines (RIL) population derived from a cross between Charleston and Dongnong 594 by mixed linear model approach. The results showed that 15 QTLs for plant height were mapped in the linkage group (LG) B1, LG D1a, LG G by software Windows QTL Cartographer Ver. 2.5. Furthermore, another 78 QTLs of plant height mapped in different populations and environments were collected, combined with the QTLs above, projected and integrated the reference map with the software BioMercator2.1. Finally, a total of 12 consensus QTLs and their corresponding markers were obtained respectively. The minimum CI was shrunk to 0.24 cM. These results would lay the foundation for mapping plant height QTL genes precisely in soybean in future.

Keywords
Soybean [Glycine max (L.) Merr.]; Plant Height; QTL mapping; Meta-analysis

Soybean (Glycine max (L.) Merr.), grown worldwide as the major crop, enrich seed protein and oil. Nowadays, obtaining high and steady yield is the main goal to soybean researchers. Molecular markers have provided a powerful new tool for breeders to search for new sources of variation and to investigate genetic factors controlling quantitatively inherited traits. Plant height is an important factor of yield, Wilcox et al (1981) found plant height have greater effect on yield. Luo (2001) confirmed the lower plants were stronger resistance to lodging and high-yielding in the research of soybean density cultivation. Wang et al (2004) detected 5 soybean height QTLs from a BC2F4 lines derived from the cross of IA2008×PI468916, 4 height QTLs of position matched to yield, further verified that plant height related to yield. Currently, mapping QTLs tend to use the data of field trials with multiple years and sites, but the stability and validity of the QTL mapping still stay in the traditional quantitative genetics research stage. Therefore, research based on stability analysis of traits and quantitative genetics is necessary (Wang et al., 2007).

With the development of molecular genetics and molecular marker technology, some high-density integrated genetic maps, such as the map by Gong (2006), have been constructed, especially soymap2 with high density (Song et al., 2004), lay the foundation for soybean height QTL from fine-mapping to gene mapping. The information of 78 soybean plant height QTLs by using CIM and MIM methods were listed in the public database Soybase (http://soybase.org/) and published in recent published papers. Sun (2006) investigated the developmental behavior of plant height QTL by unconditional and conditional QTL mapping methods (Zhu, 1995). But the reports on matured soybean plant height QTL under different environment were rare.

Mining the major QTL has become a hot spot recently years. Psychologist Glass (1976) firstly widely used Meta-analysis in physic, sociology, and behavioral science, integrated and obtained the better result conveniently. Goffinet and Gerber (2000) used meta-analysis to increase their precision and validity by mathematical model to refine integrated QTLs. In soybean, Guo (2006) integrated 62 QTL associations for resistance to soybean cyst nematode in soybean and obtained 10 consensus QTLs and the corresponding markers. Qi (2009) integrated and analyzed 65 QTLs for soybean 100-SW, and 12 consensus QTLs and the corresponding markers have been obtained. These results laid the foundation for marker-assisted selection (MAS) and gene cloning in soybean.

In the present study, we aimed to detect real QTLs under different environments. Using genetic statistics software Windows QTL Cartographer Ver. 2.5 (Wang et al., 2007), we constructed a set of soybean height from the RILs in 2006~2008, which was used to study the consistency of QTLs across places and years as well as to detect real QTLs. We collected the information of mapping QTLs from the present study from many different populations and environments. Finally, we integrated these data to the reference map, 12 consensus QTLs of soybean height were obtained, which could be used for MAS of soybean height.

1 Results and Analysis
1.1 The phenotypic information in RIL
The phenotypic results of RIL population between 2006 and 2008 were shown in table 1. The difference of maternal parent in height was significant. All plant height data displayed a typical quantitative genetic model-approximate normal distribution, and it was suitable for QTL mapping (Figure 1).

 

Table 1 The height of the parents of RIL population between 2006 and 2008


 

Figure 1 The frequency distribution of soybean height of RIL population in three years


1.2 QTL mapping of height in RILs
Based on the mixed linear model approach, we analyzed the data of soybean height in 2006~2008, including CIM (Composite Interval Mapping) and MIM (Multiple Interval Mapping) methods. In this population, 4 QTLs were detected by CIM, and another 11 QTLs were detected by MIM (Table2, Table 3). Compared with different mapping methods, CIM and MIM detected twice or more times QTLs on LG B1, LG D1a, LG G, the CI deviation was small, so the real QTLs could exist on these LGs. The explained variation of all QTLs ranged from 4.00% to 59.70% with a LOD-value between 3.04 and 10.1. On LG B1, QphB1-1 was mapped between Satt509 and Satt229, its CI was 192.9~254.8 cM, explaining about 39.00% of the variation with an additive effect of 9.69; QphB1-2 and QphB1-3 were mapped between Satt251 and Satt229, CIs of them were 192.6~261.5 cM and 261.5~329 cM, explaining about 59.70% of the variation with an additive effect of 12.20; QphB1-4, QphB1-5, and QphB1-6 were mapped between Satt509 and Satt197, Satt197 and Satt251, Sat_099 and Sat_113, CIs of them were 83~121.6 cM, 135~185.9 cM, and 356.5~422.2 cM, with the same additive effect of 12.20. On LG D1a, QphD1a-1 was mapped between Sat_062 and Sat_106, CI was 198.9~206.3 cM, explaining about 10.0% of the variation with an additive effect of -10.53; QphD1a-2 was mapped between Sat_062 and Sat_106, its CI was 198.4~207.7 cM, explaining about 9.10% of the variation with an additive effect of -13.18; QphD1a-3 was mapped between Satt220 and Sat_062, its CI was 140.9~170.0 cM, explaining about 4.00% of the variation with an additive effect of -4.53, three QTLs were detected overlap of the confidence interval and QTL of soybean height maybe exist on LG D1a. On LG G, QphG-1 was mapped between Satt199 and Sat_094, its CI was 0.4~6.6 cM, explaining about 11.00% of the variation with an additive effect of -7.57; QphG-2 was mapped between Satt505 and Satt288, its CI was 0.2~17.8 cM, explaining about 13.00% of the variation with an additive effect of -4.98; QphG-3 was mapped between Satt199 and Sat_094, its CI was 0~15.7 cM, explaining about 13.80% of the variation with an additive effect of -7.16. The major QTLs of height maybe exist on LG G, the additive effect was between -4.98 and 7.57, explaining from 11.0% to 13.80% of the variation.

 

Table 2 The results of plant height QTL by CIM


 

Table 3 The results of plant height QTL by MIM


Based on the data of height in three years, for minimizing the error, we selected the QTLs which were identified twice or more times, the QTL on LG B1, LG D1a and LG G were the major QTL of soybean height which were less affected by environment.

1.3 Collection of soybean height QTLs and Meta-analysis
We collected the information of 78 QTLs, including parents, population size, analysis method and population type from previous studies (Table 4), integrated including 36 marker intervals and 20 marker loci.

 

Table 4 Information of soybean height QTLs


With the method of meta-analysis, 12 consensus QTLs of soybean height were obtained and located on LG B1, LG C2, LG D1a, LG F, LG G, LG K, and LG M (Table 5). As Figure 2 showed, QTLs were integrated on LG B1, three consensus QTLs were analyzed from three, two, three original QTLs, respectively, the minimum CI could reach 0.24 cM.

 

Table 5 Meta-analysis of soybean height QTLs


 

Figure 2 Meta-analysis of height QTLs on LG B1


2 Discussion
In recent years, many studies on QTL mapping of soybean height used one-year phenotypic data, but little based on the multiple years and sites with different methods. From a statistical point of view, the accuracy of QTL position and affect could be detected if the database were analyzed in different locations and years (Jansen et al., 1995). In this research, the phenotype of Charleston has a distant genetic relationship and obvious differences result from different environment and the selected randomness of phenotype measuring. The phenotypic data of parents didnt change the result of QTL mapping, which based on the linkage map drawn by Zhang (2004) and phenotypic data were measured from 149 individuals. It is better to try and replenish different mapping method, different environments of the same material QTL mapping, QTL×Environment (Q×E) interaction, to obtain the real QTL with genetic stability and explaining high phenotypic variation, to improve the traits and accelerate the development of breeding (Guo et al., 2003). Epistatic effects of QTLs and QE interaction effects were considered on protein and oil content in soybean by mixed linear approach (Shan et al., 2008; 2009). Based on different model, CIM and MIM were different mapping methods (Zeng, 1994; Kao, 1999). In the present study, 15 QTLs of soybean height were detected by CIM and MIM in different environments. The objective of planting under different environments is to minimize the error, QTLs were detected twice or more times would be the major QTL.

Wang and Paigen (2002) found that 18 of the 22 human HDL-C QTLs were within the murine HDL-C QTL for guiding future research on the genetic regulation of HDL concentrations and for finding gene targets for up regulating HDL levels in mice and humans. 127 plant height QTLs of maize were refined by means of "overview" analysis. At last, 40 "real" QTLs were identified (Wang et al., 2006). CI of QTLs was greatly narrowed down by the method of Meta-analysis. The “real” QTLs with precise locations were then compared with genes affecting plant height in maize and rice, and a number of candidate genes for plant height QTG were found. Based on data of soybean height from different populations and mapping methods, the method of meta-analysis was used to integrate and obtain the major QTL, shrink the confidence interval, and improve the accuracy and validity of QTL position. In present study, 78 QTLs of height that have been reported were collected, and 56 QTLs were integrated. In the end, 12 consensus QTLs and their corresponding markers were obtained respectively. 3 QTLs were detected twice or more times by CIM and MIM on LG B1, LG D1a and LG G and furthermore used in Meta-analysis. The minimum CI was shrunk to 0.24 cM and false positive decreased greatly. According to this result, it is better to select markers covering the target segment for fine-mapping and gene mining, the real QTL is important for researching on improving soybean yield.

3 Materials and Methods
3.1 Materials and methods for QTL mapping
3.1.1 Plant materials and experimental design
A total of 149 F2:14-F2:16 population derived recombinant inbred lines (RILs) were advanced by single-seed-descent from crosses between Charleston (provided by Dr. Lijuan Qiu, Chinese Agricultural Academy of Science, Beijing, China) and Dongnong594 (developed by Northeast Agricultural University, Harbin, Heilongjiang, China). RILs were planted in experimental land of the Crop Research and Breeding Center of Land-Reclamation in Harbin between 2006 and 2008. The plants were cultured in a randomized experimental single-row plot (1 m in length and 0.5 m in width) for two or three repeats. The timing and frequency of cultural management procedures for the trials were consistent with soybean production practices under their respective environments. In our study, the RIL population was performed in duplicate in 2006 and 2007, whereas it was performed in triplicate during 2008.

3.1.2 Measurement of soybean height
Five plants of each line were randomly harvested and used for the phonotype values calculation. According to the "Bean Germplasm Description Specifications and Data Standards", measurement of plant height was performed by distance from cotyledon section to the top of the main stem each plant in soybean maturity, and an average of five replications was used.

3.1.3 Statistical analysis
The linkage map of RIL population was constructed by Zhang (2004). QTLs under different environments were separately mapped by software Windows QTL Cartographer Ver. 2.5. QTL position was determined with P<0.005 and LOD>3.0 as the threshold.

3.2 Materials and methods of Meta-analysis
Mapping information of soybean height QTLs, including names, traits, chromosomes, markers LGs, marker and population was collected from the public database Soybase (www.soybase.org) and recent reports. QTLs of one population for a given trait under a given environment were integrated as one experiment, and QTLs analyzed by the mean of several locations, were also integrated as one experiment. Two important parameters of QTL were the map position (most likely position and CI) and explained proportion of phenotypic variance, respectively. When the CI for QTL position was not available in literature, the formula proposed by Darvasi (1993; 1997) was used to calculate 95% CI.

Although the QTLs from different populations had different backgrounds and obtained by different methods, all QTLs were projected from original maps to the reference map, soymap2 (Song et al., 2004), by most likely position and CI with homothetic function. Calculations referred to reference (Chardon, 2004). Most projected QTLs were in a cluster in the integrated map. Based on the analysis principle of the software BioMercator 2.1, the QTL cluster was that the CI of original QTL covered half with each other, or the CI of original QTL covered another, reference to Qi (2009).

Meta-analysis was used to estimate the existence of the real QTL and CI. The basic process of meta-analysis was as follows. QTLs from many independent experiments associated on the same LG at neighboring interval were used to calculate a real QTL. This QTL could give five models, and the minimum Akaike Information Criteria (AIC) value was used to choose the best model QTL, called real QTL. The formula was estimated and described by Goffinet and Gerber (2000).

The height QTL was analyzed by software BioMercator ver. 2.1-tools-Meta analysis. The consensus QTLs were determined by the smallest AIC value.

Authors' contributions
Ya'nan Sun is first writer who perform this experiment, data analysis and write manuscript; Huaihai Luan and Chunyan Liu participated the experiment design and conducted this experiment; Zhaoming Qi and Dapeng Shan participated the field experiment; Guohua Hu and Qingshan Chen are the persons who take charge of this project, including experiment design, data analysis, writing and modifying of the manuscript, and corresponding authors of this manuscript. All authors have read and approved the final manuscript.

Acknowledgments
This study was partially supported by the Transgenic Specific Technology program (2009ZX08009-013B), the Chinese Introducing International Super Agricultural Science and Technology program (2006-G1(A)) and the Public Agricultural research special funds projects (200903003).

Reference
Chapman A., Pantalone V.R., Ustun A., Allen F.L., D. Landau-Ellis, Trigiano R.N., and Gresshoff P.M., 2003, Quantitative trait loci for agronomic and seed quality traits in an F2 and F4:6 soybean population, Euphytica, 129(3): 387-393 
http://dx.doi.org/10.1023/A:1022282726117

Chardon F., Virlon B., Moreau L., Falque M., Joets J., Decousset L., Murigneux A., and Charcosset A., 2004, Genetic architecture of flowering time in maize as inferred from quantitative trait loci meta-analysis and synteny conservation with the rice genome, Genetics, 168(4): 2169-2185
http://dx.doi.org/10.1534/genetics.104.032375 PMid:15611184 PMCid:1448716

Chen Q.S., Zhang Z.C., Liu C.Y., Xin D.W., Shan D.P., Qiu H.M., and Shan C.Y., 2007, QTL analysis of major agronomic traits in soybean, Scientia Agricultura Sinica, 40(1): 41-47

Darvasi A., and Soller M., 1997, A simple method to calculate resolving power and confidence interval of QTL map location, Behavior Genetics, 27(2): 125-132 
http://dx.doi.org/10.1023/A:1025685324830 PMid:9145551

Darvasi A., Weinreb A., Minke V., Weller J.I., and Soller M., 1993, Detecting marker-QTL linkage and estimating QTL gene effect and map location using a saturated genetic map, Genetics, 134(3): 943-951 
PMid:8349116 PMCid:1205528

Glass G.V., 1976, Primary, secondary, and meta-analysis of research, Educational Researcher, 5(10): 3-8 
http://dx.doi.org/10.2307/1174772 http://dx.doi.org/10.3102/0013189X005010003

Goffinet B., and Gerber S., 2000, Quantitative trait loci: a meta-analysis, Genetics, 155(1): 463-473 PMid:10790417 PMCid:1461053

Gong P.T., Mu J.G., Zhao J.R., Wang X.L., Bai Y.N., and Fang X.J., 2006, An Integrated Soybean Genetic Linkage Map Comprising 315 SSRs and 40 AFLPs, Molecular Plant Breeding, 4(3): 309-316

Guan R.X., 2004, QTL mapping of soybean agronomic characters and genetic diversity analysis of soybean cultivars from China and Japan, Thesis for Post-D., Beijing: Chinese Academy of agricultural sciences, Supervisors: Chang R.Z., Qiu L.J., pp.27-29

Guo B., Sleper D.A., Lu P., Shannon J.G., Nguyen H.T., and Arelli P.R., 2006, QTLs associated with resistance to soybean cyst nematode in soybean: Meta-analysis of QTL location, Crop Sci., 46: 595-602 
http://dx.doi.org/10.2135/cropsci2005.04-0036http://dx.doi.org/10.2135/cropsci2005.04-0036-2

Guo L.B., Luo L.J., Xing Y.Z., Xu C.G., Mei H.W., Wang Y.P., Zhong D.B., Qian Q., Ying C.S., and Shi C.H., 2003, Dissection of QTLs in two years for important agronomic traits in rice (Oryza sativa L.), Chinese J. Rice Sci., 17(3): 211-218

Jansen R.C., Van Ooijien J.M., Stam P., Lister C., and Dean C., 1995, Genotype-by-environment interaction in genetic mapping of multiple quantitative trait loci, Theor. Appl. Genet, 91(1): 33-37 

http://dx.doi.org/10.1007/BF00220855


Kabelka E.A., Diers B.W., Fehr W.R., LeRoy A.R., Baianu I.C., You T., Neece D.J., and Nelson R.L., 2004, Putative alleles for increased yield from soybean plant introductions, Crop Sci., 44: 784-791 
http://dx.doi.org/10.2135/cropsci2004.0784

Kao C.H., Zeng Z.B., and Robert D., Teasdale R.D., 1999, Multiple interval mapping for quantitative trait loci, Genetics, 152(3): 1203-1216 
PMid:10388834 PMCid:1460657

Li D.D., Pfeiffer T.W., and Cornelius P.L., 2008, Soybean QTL for yield and yield components associated with Glycine soja alleles, Crop Sci., 48(2): 571-581 
http://dx.doi.org/10.2135/cropsci2007.06.0361

Luo G.Y., Zhan Y., Liu S.L., Kong X., Wang S.M., Sun D.M., and Gai J.Y., 2001, New record of high yield in Xindadou No.1 and Shidadou No.1, Soybean Science, 20(4): 270-273

Mansur L.M., Lark K.G., Kross H., and Oliveira A., 1993, Interval mapping of quantitative trait loci for reproductive, morphological, and seed traits of soybean, Theor. Appl. Genet., 86(8): 907-913
http://dx.doi.org/10.1007/BF00211040

Qi Z.M., Sun Y.N., Chen L.J., Guo Q., Liu C.Y., Hu G.H., and Chen Q.S., 2009, Meta-analysis of 100-seed weight QTLs in soybean, Scientia Agricultura Sinica, 42(11): 3795-3803

Shan D.P., Qi Z.M., Qiu H.M., Shan C.Y., Liu C.Y., Hu G.H., and Chen Q.S., 2008, Epistatic effects of QTLs and QE interaction effects on oil content in soybean, Acta Agron. Sin., 34(6): 952-957
http://dx.doi.org/10.3724/SP.J.1006.2008.00952

Shan D.P., Zhu R.S., Chen L.J., Qi Z.M., Liu C.Y., Hu G.H., and Chen Q.S., 2009, Epistatic effects and QE interaction effects of QTLs for protein content in soybean, Acta Agron. Sin., 35(1): 41-47 
http://dx.doi.org/10.3724/SP.J.1006.2009.00041

Song Q.J., Marek L.F., Shoemaker R.C., Lark K.G., Concibido V.C., Delannay X., Specht J.E., and Cregan P.B., 2004, A new integrated genetic linkage map of the soybean, Theor. Appl. Genet., 109(1): 122-128 
http://dx.doi.org/10.1007/s00122-004-1602-3 PMid:14991109

Sun D.S., Li W.B., Zhang Z.C., Chen Q.S., and Yang Q.K., 2006, Analysis of QTL for plant height at different developmental stages in soybean, Acta Agron. Sin., 32(4): 509-514

Wang D., Graef G.L., Procopiuk A.M., and Diers B.W., 2004, Identification of putative QTL that underlie yield in interspecific soybean backcross populations, Theor. Appl. Genet., 108(3): 458-467
http://dx.doi.org/10.1007/s00122-003-1449-z PMid:14504749

Wang S.C., Christopher J. Basten, and Zeng Z.B., 2007, Windows QTL cartographer V2.5, Program in Statistical Genetics, North Carolina State University, Raleigh, NC. (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm)

Wang X., and Paigen B., 2002, Quantitative trait loci and candidate genes regulating hdl cholesterol: a murine chromosome map, Arteriosclerosis, Thrombosis, and Vascular Biology, 22(9): 1390-1401
http://dx.doi.org/10.1161/01.ATV.0000030201.29121.A3 PMid:12231556

Wang X.Z., 2008, Inheritance, stability analysis and QTL mapping of yield related traits in soybean, Thesis for Ph.D., Beijing: Chinese Academy of agricultural sciences, Supervisors: Zhou X.A., pp.48-51

Wang X.Z., Zhang X.J., Zhou R., Sha A.H., Wu X.J., Cai S.P., Qiu D.Z., and Zhou X.A., 2007, QTL analysis of seed and pod traits in soybean ril population, Acta Agron. Sin., 33(3): 441-448

Wang Y., Yao J., Zhang Zh.F., and Zheng Y.L., 2006, The comparative analysis based on maize integrated QTL map and meta-analysis of plant height QTLs, Chinese Science Bulletin, 51(18): 2219-2230
http://dx.doi.org/10.1007/s11434-006-2119-8

Wang Z., 2004, Construction of soybean SSR based map and QTL analysis important agronomic traits, Thesis for M.S., Guangxi: Guangxi University, Supervisors: Fang X.J., pp.66-72

Wilcox J.R., and Tuneo S., 1981, Interrelationships among height, lodging and yield in determinate and indeterminate soybeans, Euphytica, 30(2): 323-326
http://dx.doi.org/10.1007/BF00033993

Wu X.L., Wang Y.J., He C.Y., Chen S.Y., Gai J.Y., and Wang X.C., 2001, QTL analysis of major agronomic traits in soybean, Genetics, 28(10): 947-955

Zeng Z.B., 1994, Precision mapping of quantitative trait loci, Genetics, 136(4): 1457-1468 PMid:8013918 PMCid:1205924

Zhang W.K., Wang Y.J., Luo G.Z., Zhang J.S., He C.Y., Wu X.L., Gai J.Y., and Chen S.Y., 2004, QTL mapping of ten agronomic traits on the soybean (Glycine max L. Merr.) genetic map and their association with EST markers, Theor. Appl. Genet, 108(6): 1131-1139 
http://dx.doi.org/10.1007/s00122-003-1527-2 PMid:15067400

Zhang Z.C., Zhan X.L., Chen Q.S., Teng W.L., Yang Q.K., and Li W.B., 2004, QTL mapping of seed oil and protein content of soybean, Soybean Science, 23(2): 81-85

Zhu J., 1995, Analysis of conditional genetic effects and variance components in developmental genetics, Genetics, 141(4): 1633-1639 PMid: 8601500 PMCid:1206893

Zhu X.L., 2006, Constructing of genetic linkage map and QTL mapping of important agronomic traits in two soybean populations, Thesis for M.S., Harbin: Northeast Agricultural University, Supervisors: Hu G.H., pp.23-36

Legume Genomics and Genetics
• Volume 3
View Options
. PDF(834KB)
. FPDF
. HTML
. Online fPDF
Associated material
. Readers' comments
Other articles by authors
. Ya'nan Sun
. Huaihai Luan
. Zhaoming Qi
. Dapeng Shan
. Chunyan Liu
. Guohua Hu
. Qingshan Chen
Related articles
. Soybean [ Glycine max (L.) Merr. ]
. Plant Height
. QTL mapping
. Meta-analysis
Tools
. Email to a friend
. Post a comment

ass fuck sexo anal Gruppen Pornos Blondine Pornos inzest porn hd
ankara escort
casumo
escort bodrum
erotik filmler
beykent spor salonu
free brazzers
istanbul escort