BACKGROUND
cotton is the most important natural textile fibre and is a significant oilseed crop, providing approximately 35% of the total fibre used in the world  <cit> , profoundly affecting the world economy and changing human daily life.
china is the world’s largest raw cotton producer and consumer, providing two-thirds of the world’s cotton together with india and the united states  <cit> .
the cotton genus consists of over  <dig> species distributed in arid to sub-arid regions of the tropics and sub-tropics  <cit> .
among the cotton family, there are  <dig> traditional ad tetraploid species that are inter-crossable to various degrees, and  <dig> diploid species that are divided into  <dig> diploid genome groups   <cit> .
there are four species have been domesticated for spinnable-fibre and collectively provide the world’s cotton production, including two diploids of g. herbaceum  and g. arboreum , which are old world cottons from asia-africa, and two new world tetraploid cottons from the americas, which are g. hirsutum or upland cotton  and g. barbadense or extra-long staple fibre cotton   <cit> .
tetraploid cottons are presumably derived from a single polyploidization event that occurred 1– <dig> million years ago  between ancestors of g. arboreum  and g. raimondii   <cit> , accounting for 98% of the world’s cotton production.
upland cotton is the primary source for cotton production, representing over 95% of the cotton fibre produced in the world  <cit> ; however, its genetic diversity is narrow.
g. barbadense has been used to improve fibre quality in upland cotton, as the second domesticated tetraploid cotton.
however, it has some limitations, such as low yield and limited adaption.
on the contrary, wild cotton germplasms harbour extensive genetic diversity and potential practicability, and they have rich sources of novel traits and are currently being mined to improve many beneficial agronomic traits.
for example, g. tomentosum contains many unique agronomic traits, including insect-pest resistance, salt tolerance, heat tolerance, drought tolerance, nectarilessness and lint colour.
g. mustelinum has fine fibre quality and verticillium wilt resistantce.
g. darwinii has many excellent traits, such as drought tolerance, finer fibre fineness, fusarium wilt and verticillium wilt resistance  <cit> .
in addition, with the rapid development of next-generation sequencing  technologies, the recent publications on the genome sequence of the tetraploid ad1-genome  <cit>  and ad2-genome  <cit>  and the diploid a-genome  <cit>  and d-genome  <cit>  have improved the development of new analyses and comparative approaches for the genomics of both diploid and polyploid cottons.
single nucleotide polymorphisms  are known to be of considerable importance because they have a much higher abundance in the genome and are used to determine the population structure and for linkage disequilibrium  analysis  <cit> .
similarly, insertions/deletions  have been used in rice, cotton and chickpeas for fine mapping and marker-assisted selection .
furthermore, gene expression and function could be affected by the position of snps and indels within a genome  <cit> .
for example, variations present in coding regions and regulatory sequences may change protein functions and induce/repress gene expression.
therefore, the discovery of polymorphisms is very important in the study of genomic variation in crop species.
to date, snp discoveries have been applied in other crops beyond cotton  <cit> .
snp discoveries in cotton have progressed using different methods, such as bac-end sequences  <cit> , single copy sequences  <cit> , transcriptome sequencing  <cit> , reduced representation libraries  techniques  <cit> .
nevertheless, based on these techniques, access to information on snps is still very limited.
in addition, the resources of wild cotton have not yet been excavated due to the large and complex genome, which greatly hinders the genetic research and cotton breeding.
fortunately, with the rapid developments and applications of ngs technologies, many complexity reduction approaches have been developed based on ngs platforms, such as genotype-by-sequencing , 2b-rad and rrls.
therefore, with the successful application of these methods, a large number of sequence polymorphisms have been detected, including snps and indels  <cit> .
recently, based on double barcode genotyping systems and deep sequencing, slaf-seq was developed, which is an accurate and cost-effective high-throughput sequence-based technology  <cit> .
more importantly, the slaf-seq does not depend on the sequence of the reference genome, and it reduces the complexity of the reference genome, which is particularly important for species with an unknown genome  <cit> .
here, the slaf-seq method was used to identify dna polymorphisms between the reference   <cit>  and  <dig> tetraploid cotton species, i.e., g. hirsutum cv.
emian <dig> , g. barbadense acc.
3– <dig> , g. tomentosum , g. mustelinum  and g. darwinii  genomes.
using the slaf-seq data, an integrated analysis was carried out.
therefore, there were three main objectives in this study.
first, interspecific variations were explored and characterized.
second, the phylogenetic trees of tetraploid cotton species were reconstructed.
third, the interspecific snp and indel markers were developed for future introgression breeding.
in summary, the new dna polymorphisms presented in this study will dramatically increase the efficiency for future research into high-density interspecific mapping, introgression breeding, genetic dissection and gene utilization in cotton.
RESULTS
build variation resources of tetraploid cottons
in the present study, a new variation resource was generated that had direct relevance to the current scenario of cotton breeding strategies.
using an illumina high-throughput sequencing platform, a total of  <dig> , <dig> 80-bp long paired-end reads were generated from the five tetraploid cotton species, including two cultivars and three wild species .
after quality filtering,  <dig> , <dig> high-quality reads were obtained, which varied from  <dig>  to  <dig>  million for the different materials.
nearly 80% of the high-quality reads were mapped to the cotton reference genome, which covered approximately  <dig> % of the total genome for each material .
we also developed  <dig>  slaf tags among the five cotton species.
slaf tags were compared, polymorphic slaf tags were identified , and it was observed that the number of polymorphic slaf tags were unevenly distributed between the a subgenome  and d subgenome .
the average depth for g. hirsutum cv.
emian <dig> , g. barbadense acc.
3– <dig> , g. tomentosum , g. mustelinum  and g. darwinii  were  <dig> ,  <dig> ,  <dig> ,  <dig> , and  <dig> , respectively.table  <dig> summary of sequence data and mapping statistics on the tm- <dig> genome
mapq <dig> mapping quality of 30
identification and characteristics of snps and indels
totally,  <dig> , <dig> snps and  <dig>  indels were identified.
additionally, after filtering out the low confidence, a total of  <dig>  polymorphisms, including  <dig>  snps and  <dig>  indels, were detected between the reference tm- <dig> genome and the five tetraploid cotton species.
eventually, the total number of homozygous dna polymorphisms was  <dig>  for subsequent analysis, which included  <dig>  snps and  <dig>  indels .
further, the snps and indels between each of the five tetraploid cotton species, and the reference tm- <dig> genome were also detected and showed that the gb_3- <dig> had the highest dna polymorphisms than other species, followed by gd; the dna polymorphisms of gh_e <dig> were the least .
notably, the overlap of snps and indels revealed a moderate degree of overlap among the five tetraploid cotton species .fig.
<dig> distribution and number of snps and indels detected in each of cotton species.
a, b the number of snps  and indels  detected in each of cotton species are illustrated in the bar graphs by different colours as indicated.
total number in each of cotton species is indicated by red bar graphs.
the homozygous number in each cotton species is indicated by green bar graphs.
c, d overlap of snps  and indels  across five species.
e, f the distribution of snps  and indels  detected in five cotton species .
i: the chromosomes; ii-vi: the total number of homozygous snps and indels of gd, gm, gt, gb_3– <dig> and gh_e <dig>  respectively
the snps and indels that were detected in each of the five cotton species were further analysed.
the snps were classified as transitions  and transversions  based on nucleotide substitutions .
among the transitions, the proportions of a/g were slightly higher than c/t transitions in gh_e <dig>  gt, and gm, while the proportions of a/g was slightly lower than c/t transitions in gb_3- <dig> and gd.
for transversions, the percentage of t/a transversions was relatively higher than others, namely, a/c, g/t and c/g.
the ratio between transitions and transversions  for gh_e <dig>  gb_3- <dig>  gt, gm and gd were  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig> , respectively.
for indels, the length of the insertions ranged from 1 bp to 20 bp; however, the length distribution of deletions was 1 bp to 35 bp .
the majority of indels were single nucleotides, accounting for  <dig> %, di- to tetra-nucleotides were  <dig> %, and the remaining  <dig> % were ≥ 5 bp.fig.
<dig> percentage of snps substitution types, the number of at-to-gc conversion in transposons, and length distribution of indels among five cotton species.
a percentage of different snps substitution types in five cotton species.
b the number of at-to-gc conversion in transposons in five cotton species.
tes are indicated by blue bar graphs; non-tes are indicated by red bar graphs.
c length distribution of indels among five cotton species.
the number of insertions and deletions were shown on the y-axis in the bar graph.
the numbers of various lengths  are shown on the x-axis in the bar graph.
gh_e22/ref_tm- <dig>  gb/ref_tm- <dig>  gt/ref_tm- <dig>  gm/ref_tm- <dig> and gd/ref_tm- <dig>  respectively indicates dna polymorphisms identified between gh_e <dig>  gb_3- <dig>  gt, gm, gd and the reference genome g. hirsutum acc.
tm-1
the dna polymorphisms in transposable elements  were also explored in the five cotton species.
the percentages of snps/indels in tes for gh_e <dig>  gb_3- <dig>  gt, gm and gd were  <dig> %/ <dig> %,  <dig> %/ <dig> %,  <dig> %/ <dig> %,  <dig> %/ <dig> %, and  <dig> %/ <dig> %, respectively.
the number of snps and indels was calculated in te and non-te regions .
further, through the analysis of snps and indels in tes, it was found that the dna polymorphisms mainly existed in long terminal repeat  retrotransposons, particularly gypsy-type .
additionally, at-to-gc conversion was larger in the te than non-te regions, which was significantly outnumbered in gb_3- <dig> than in other cotton species .
genomic distribution of snps and indels across the reference tm- <dig> genome
to survey the landscape of dna polymorphisms for each of the five cotton species, the distributions of dna polymorphisms were obtained.
the number of homozygous snps and indels was detected for subsequent analysis, and their frequency  varied across all  <dig> chromosomes.
herein, snps were observed to be unevenly distributed among the chromosomes and between the at and dt in each of cotton species, with more snps in the at than the dt .
likewise, the distribution of indels among the chromosomes and between the at and dt was also not uniform .
in addition, the total number of dna polymorphisms and their density  varied across the  <dig> cotton chromosomes in the five cotton species, and the total number of dna polymorphisms was found to be positively correlated with the chromosome length .
among them, the maximum number of homozygous snps was found on chromosome a <dig> , and the minimum number was found on chromosome d <dig> .
however, for indels, the maximum number was  <dig> on chromosome a <dig>  and the minimum number was  <dig> on chromosome d <dig> .
the distribution of gene densities across chromosomes from the tm- <dig> reference genome was also calculated, which showed the lowest gene density was in the middle of the chromosome .
in addition, the genomic distribution of dna polymorphisms was determined by calculating the frequency of occurrence .
the highest frequency of homozygous snps was on chromosome a <dig> , and the lowest frequency of snps was on chromosome d <dig> .
similarly, the homozygous indel frequency was the highest on chromosome a <dig> , and the lowest on chromosome d <dig> .
similarly, the total numbers of snps and indels were unevenly distributed amomg the chromosomes and between at and dt, with fewer dna polymorphisms in the d-subgenome than the a-subgenome .
furthermore, a total of  <dig> pairs of homeologous genes containing snps were identified between at and dt in the reference tm- <dig> genome, including  <dig> pairs of homeologous genes with one containing snps in at, and the other gene did not contain snps in dt;  <dig> pairs of homeologous genes with one containing snps in dt, and the other gene did not contain snps in at; and  <dig> pairs of homeologous genes with genes containing snps in both at and dt.
.
characterization and functional significance of snps and indels
using the annotation of the cotton reference genome , the distributions of snps and indels within various genomic features were revealed.
in general, a similar distribution pattern between snps and indels was observed in the five species .
the percentage of snps in genic regions ranged from  <dig>  to  <dig> %, while the percentage of indels ranged from  <dig>  to  <dig> %, among the five cotton species.
however, a substantial portion of snps  and indels  were identified in the intergenic regions.
a significant percentage of snps  and indels  were also detected in the 5 kb upstream  and 5 kb downstream regulatory regions among the five cotton species .fig.
<dig> annotation of snps and indels.
a, b distribution of snps  and indels  in different genomic regions for the five cotton species.
c, d distribution of snps  and indels  in different intergenic and genic regions among the five cotton species.
in cds region, the number of synonymous and non-synonymous snps detected has also been shown among the five cotton species
to further investigate the effect of snps in the coding sequence  region , amino acid substitution was analysed and a large proportion was detected to be non-synonymous among the five cotton species.
the ratio of non-synonymous to synonymous snps was approximately  <dig>  for gh_e <dig>   <dig>  for gb_3- <dig>   <dig>  for gt,  <dig>  for gm, and  <dig>  for gd .
further, these non-synonymous snps were present in  <dig>   <dig>   <dig>   <dig>  and  <dig> genes for gh_e <dig>  gb_3- <dig>  gt, gm, and gd, respectively .
the non-synonymous snps resided in the genes containing leucine-rich repeats , pentatricopeptide repeats , protein kinases, protein tyrosine kinases, wd domain repeats and nb-arc domains .
among the five cotton species,  <dig> large-effect snps and  <dig> large-effect indels were found in a total of  <dig> genes , which also contained many important functional domains, such as zinc-finger and p <dig> domains.fig.
<dig> the ratio of non-synonymous to synonymous snps, and frequency of top  <dig> pfam domains represented in the genes containing non-synonymous snp among the five cotton species.
a the ratio of non-synonymous to synonymous snps among five cotton species.
b-f frequency of top  <dig> pfam domains represents in the genes containing non-synonymous snp in gh_e <dig>  gb_3- <dig>  gt, gm and gd, respectively
gene ontology  enrichment analysis revealed putative functions of the genes containing non-synonymous snp variations among the five cotton species.
based on the reference genome annotation, a number of significantly enriched go terms in each of the three main categories  were identified.
the significantly enriched go terms, which were genes involved in biological processes, molecular function and cellular component of the go classification, revealed striking differences between the five cotton species .
for example, in gb_3- <dig>  the significant terms were related to reproduction, including multi-organism reproductive processes , recognition of pollen , pollen-pistil interaction , pollination , and reproductive processes .
in other cotton species, they also had their own unique characteristics; for instance, oxidoreductase activity  and oxidation-reduction process  were significantly enriched in gt, which was associated with stress resistance .
in gm, many significantly enriched categories were mainly associated with trna, which was involved in the synthesis of protein .
significant go terms of gd were highly concentrated in cellular components, including the organelle outer membrane , ribonucleoprotein complex  and intracellular parts .fig.
<dig> significantly enriched go terms and phylogenetic tree among five cotton species.
a the significantly enriched go terms are identified among five cotton species.
b phylogenetic tree was constructed among five cotton species.
c phylogenetic tree of a-homoeologs was constructed among five cotton species.
d phylogenetic tree of d-homoeologs was constructed among five cotton species
phylogenetic analysis of the cotton diploids and polyploids
understanding the phylogenetic relationships among the various genomes is not only significant for genetics research and breeding but also essential for conducting comprehensive studies of the plant genomes, especially in genome organization, function and evolution.
to further study the phylogenetic relationships in different cotton genomes, phylogenetic trees of the cotton diploids and polyploids were constructed.
specifically, and unsurprisingly, the five tetraploid cotton species were classified into three branches, one consisting of gh_e <dig>  and gt , one consisting of gm , and the third containing gb_3- <dig>  and gd .
additionally, two sets of snps were generated, and the phylogenetic trees for the at and dt homoeologs were constructed, independently .
for the a-subgenome, d <dig> clade as one outgroup species, it showed a high concordance in the phylogeny of the cotton diploids and polyploids.
the concordance analysis identified three main branches, which were the a-subgenomes of five tetraploid cotton clades,  clade and d <dig> clade .
compared with a-homoeologs, the  clade was one outgroup species, and analysis results of d-homoeologs showed that the clades were more consistent with the sister-species relationship.
the d <dig> clade had closer relationships with the d-subgenomes of the five tetraploid cotton clades than  .
interspecific variations for introgression genetics and breeding
to broaden the genetic basis of upland cotton, introgression line populations were developed in our laboratory, which used gh_e <dig> as the receptor parent and gb_3- <dig>  gm, gt, and gd as donor parents, respectively.
the snps/indels identified in this study facilitated the introgression genetics and breeding as molecular markers.
combined with our laboratory resources to investigate the extensive genetic diversity and the potential practicability, gh_e <dig> vs. gb_3- <dig>  gh_e <dig> vs. gt, gh_e <dig> vs. gm and gh_e <dig> vs. gd were analysed.
to gain insight into the differences in various combinations, the total number of homozygous snps/indels was identified between gh_e <dig> vs. gb_3- <dig>  gh_e <dig> vs. gt, gh_e <dig> vs. gm and gh_e <dig> vs. gd, which were  <dig> / <dig> ,  <dig> / <dig> ,  <dig> / <dig>  and  <dig> / <dig> , respectively , observing that gh_e <dig> vs. gb_3– <dig> had the higher dna polymorphisms than the others.
further, overlaps of the snps/indels among the four combinations were shown with a venn diagram .
the overlapped snps/indels between gh_e <dig> vs. gb_3– <dig> and gh_e <dig> vs. gd were larger than other two combinations.
the unique dna polymorphisms were maximum in gh_e <dig> vs. gm , and minimum in gh_e <dig> vs. gt .
intriguingly, the snps/indels were not evenly distributed cross the chromosomes.
the snps/indels in the at were more frequent than those in the dt .
the gh_e <dig> vs. gb_3– <dig> had the highest dna polymorphisms on each chromosome, followed in descending order by gh_e <dig> vs. gd, gh_e <dig> vs. gm and gh_e <dig> vs. gt, except for the chromosome a <dig>  where gh_e <dig> vs. gd had the maximum snps  and gh_e <dig> vs. gt had the maximum indels .
the snps/indels on a <dig> were the maximum in the four combinations, which were  <dig> / <dig>   <dig> / <dig>   <dig> / <dig>   <dig> / <dig> in gh_e <dig> vs. gb_3– <dig>  gh_e <dig> vs. gt, gh_e <dig> vs. gm and gh_e <dig> vs. gd, respectively.
in contrast, the minimum number of snps were on d <dig>  in gh_e <dig> vs. gb_3– <dig>  d <dig>  in gh_e <dig> vs. gt, d <dig>  in gh_e <dig> vs. gm and d <dig>  in gh_e <dig> vs. gd.
however, the minimum number of indels, were on d <dig>  in gh_e <dig> vs. gb_3– <dig>  d <dig>  in gh_e <dig> vs. gt, d <dig>  in gh_e <dig> vs. gm and d <dig>  in gh_e <dig> vs. gd, respectively .fig.
<dig> distribution of snps and indels between four gh_e <dig> vs. gb_3- <dig>  gh_e <dig> vs. gt, gh_e <dig> vs. gm and gh_e <dig> vs. gd.
a, b the total number of snps and indels between different combinations are illustrated in the bar graphs.
c, d the overlap of snps and indels.
e, f distribution of snps and indels on each chromosome are illustrated in the bar graphs in different combinations
the randomly selected  <dig> snps on chromosome a <dig> were used to design primers to check the gel-based polymorphism.
based on the sscp analysis, there were  <dig> polymorphic markers  between gh_e <dig> and gb_3– <dig>   <dig>  between gh_e <dig> and gt,  <dig>  between gh_e <dig> and gm, and  <dig>  between gh_e <dig> and gd.
discussion
in the present study, a high quality of  <dig> , <dig> 80-bp long paired-end reads and  <dig>  slafs were generated using high-throughput slaf-seq.
the average sequencing depth of each cotton species was  <dig> ×.
additionally, a total of  <dig>  homozygous snps, and  <dig>  homozygous indels were identified.
interspecific dna polymorphisms for gb_3– <dig>  gt, gm and gd occurred at a much higher rate compared to intraspecific dna polymorphisms for gh_e <dig> , especially in gb_3– <dig>  which may be that they belong to different cotton species and the distant phylogenetic relationships.
gb_3– <dig> and gd were found to have more overlap than the other three species, and a similar phenomenon also occurred in the following two species, gh_e <dig> and gt, which may be related to their genetic relationship .
as is known, the conventional five polyploid cotton species have been divided into three branches by previous cotton researchers; one branch is g. hirsutum and g. tomentosum, the other is g. mustelinum, and the third is g. barbadense and g. darwinii  <cit> .
our analysis, as a new resource, provides further evidence on the accuracy of the previous studies in their phylogenetic relationships.
for the snps, the non-synonymous/synonymous snp ratios in different cotton species were  <dig> – <dig>  , which was slight higher than ir64/pokkali  and ir64/n <dig>  in rice  <cit> , and js- <dig>  and upsm- <dig>  in soybean  <cit> .
in addition, the frequency of transitions was found to be significantly higher than transversions.
transition bias has been reported in rice, maize, chickpeas, and loquat .
the proportions of c/t were slightly lower than a/g transitions in gh_e <dig>  gt and gm, while the proportions of c/t were slightly higher than a/g transitions in gb_3– <dig> and gd, which may be associated with divergence times.
similar observations were previously reported in g. hirsutum and g. barbadense  <cit> .
however, the percentage of t/a transversions was relatively higher than others, namely, a/c, g/t and c/g.
similar observations were detected in rice and citrus  <cit> , which is remained to be explained as an unexpected observation.
the ratio of transitions and transversions  for gh_e <dig>  gb_3– <dig>  gt, gm and gd were  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig> , respectively, which were slightly higher than g. hirsutum  and g. barbadense  of a previous report in cotton  <cit> , ir64/pokkali  and ir64/n <dig>  in rice  <cit> , and js- <dig>  and upsm- <dig>  in soybean  <cit> .
further, at-to-gc conversion in gb_3– <dig> significantly outnumbered gh_e <dig>  followed in descending order by gd, gm and gt , which could provide the potential dna methylation sites, if occurring in tes  <cit> .
however, as for indels, the length distribution of indels was found to be larger than that in rice  <cit>  and less than that in soybean  <cit> .
intriguingly, indels in each of the individual five cottons also showed an upward bias towards single nucleotide insertions and deletions, which was similar to rice and soybean  <cit> .
the abundance of snps and indels were detected in 5 kb upstream and 5 kb downstream regulatory regions  <cit> , which may be due to the non-coding regulatory regions having lower sequence conservation and levels of purifying selection pressure than coding regions  <cit> .
in animals, the variations in the cis-regulatory regions have demonstrated the importance in the regulation of gene expression .
in plants, a few studies also have shown that such genetic variations play a major role in the regulation of gene expression and agronomic traits .
functional annotation of these genes containing variations revealed putative characteristics in different cotton species .
in gb_3– <dig>  a number of the significant terms were related to reproduction , which suggested that the variations may result in the disruption of regulatory genes determining reproductive fertility  <cit>  and may help to explain the hybrid breakdown and low yield in offspring between g. barbadense and g. hirsutum.
in previous studies, the genus gossypium has been researched broadly in phylogeny and a consensus phylogenetic tree was generated  <cit> .
the phylogenetic analysis has been reported in a previous study  <cit> , and it was based on the targeted genes sequence capture method.
here, we used high-throughput sequencing to reconstruct the phylogenetic trees that showed a high concordance with the historically hypothesized phylogeny, providing strong support for the phylogenetic relationships in diploid and allotetraploid cottons .
interspecific introgression breeding is a traditional method and requires labour, time and funds; however, it can be expedited by molecular marker-assisted selection.
in a recent study, the cotton breeding potential of introgression lines between g. hirsutum and g. barbadense was investigated for its yield and fibre quality improvement, including heterosis, combining ability and genetic effects  <cit> .
similarly, the introgression line populations were constructed in our laboratory with gh_e <dig> as the receptor and the other tetraploid species as donors.
currently, only  <dig>  gel-based loci were identified between gh_e <dig> and gb_3– <dig>  <cit> .
as is known, the gel-based markers are time-consuming and have low efficiency.
in this study, a huge number of snps and indels identified, plus the high throughput platform, will facilitate and speed up the genetics and breeding of the introgression lines .
although few polymorphisms were detected from the randomly selected  <dig> snp markers on a <dig>  because the sscp is still a gel-based genotyping method and has low resolution, we do believe that more polymorphisms will be detected with a high-resolution detecting platform, such as the snp chip or genotyping by resequencing.
CONCLUSIONS
a large number of high-quality dna polymorphisms through slaf-seq were identified in five tetraploid cotton species.
further comprehensive analysis of dna polymorphisms provided valuable insights into the different tetraploid cotton species, especially between cultivated and wild cottons.
to better utilize the novel characteristics of drought tolerance, defence responses and salt stress tolerance, etc.
in wild cottons, numerous interspecific variations were identified here, which will be very practical in future introgression breeding.
overall, the comprehensive data generated in this study provided insights into cotton evolution and a resource for future researches in high-density interspecific mapping, introgression breeding improvement in cotton.

 cotton is the most important natural textile fibre and is a significant oilseed crop, providing approximately 35% of the total fibre used in the world, profoundly affecting the world economy and changing human daily life. china is the world’s largest raw cotton producer and consumer, providing two-thirds of the world’s cotton together with india and the united states. the cotton genus consists of over species distributed in arid to sub-arid regions of the tropics and sub-tropics. among the cotton family, there are traditional ad tetraploid species that are inter-crossable to various degrees, and diploid species that are divided into diploid genome groups. there are four species have been domesticated for spinnable-fibre and collectively provide the world’s cotton production, including two diploids of g. herbaceum and g. arboreum, which are old world cottons from asia-africa, and two new world tetraploid cottons from the americas, which are g. hirsutum or upland cotton and g. barbadense or extra-long staple fibre cotton. tetraploid cottons are presumably derived from a single polyploidization event that occurred 1– million years ago between ancestors of g. arboreum and g. raimondii, accounting for 98% of the world’s cotton production. upland cotton is the primary source for cotton production, representing over 95% of the cotton fibre produced in the world; however, its genetic diversity is narrow. g. barbadense has been used to improve fibre quality in upland cotton, as the second domesticated tetraploid cotton. however, it has some limitations, such as low yield and limited adaption. on the contrary, wild cotton germplasms harbour extensive genetic diversity and potential practicability, and they have rich sources of novel traits and are currently being mined to improve many beneficial agronomic traits. for example, g. tomentosum contains many unique agronomic traits, including insect-pest resistance, salt tolerance, heat tolerance, drought tolerance, nectarilessness and lint colour. g. mustelinum has fine fibre quality and verticillium wilt resistantce. g. darwinii has many excellent traits, such as drought tolerance, finer fibre fineness, fusarium wilt and verticillium wilt resistance. in addition, with the rapid development of next-generation sequencing technologies, the recent publications on the genome sequence of the tetraploid ad1-genome and ad2-genome and the diploid a-genome and d-genome have improved the development of new analyses and comparative approaches for the genomics of both diploid and polyploid cottons. single nucleotide polymorphisms are known to be of considerable importance because they have a much higher abundance in the genome and are used to determine the population structure and for linkage disequilibrium analysis. similarly, insertions/deletions have been used in rice, cotton and chickpeas for fine mapping and marker-assisted selection. furthermore, gene expression and function could be affected by the position of snps and indels within a genome. for example, variations present in coding regions and regulatory sequences may change protein functions and induce/repress gene expression. therefore, the discovery of polymorphisms is very important in the study of genomic variation in crop species. to date, snp discoveries have been applied in other crops beyond cotton. snp discoveries in cotton have progressed using different methods, such as bac-end sequences, single copy sequences, transcriptome sequencing, reduced representation libraries techniques. nevertheless, based on these techniques, access to information on snps is still very limited. in addition, the resources of wild cotton have not yet been excavated due to the large and complex genome, which greatly hinders the genetic research and cotton breeding. fortunately, with the rapid developments and applications of ngs technologies, many complexity reduction approaches have been developed based on ngs platforms, such as genotype-by-sequencing, 2b-rad and rrls. therefore, with the successful application of these methods, a large number of sequence polymorphisms have been detected, including snps and indels. recently, based on double barcode genotyping systems and deep sequencing, slaf-seq was developed, which is an accurate and cost-effective high-throughput sequence-based technology. more importantly, the slaf-seq does not depend on the sequence of the reference genome, and it reduces the complexity of the reference genome, which is particularly important for species with an unknown genome. here, the slaf-seq method was used to identify dna polymorphisms between the reference and tetraploid cotton species, i.e., g. hirsutum cv. emian, g. barbadense acc. 3–, g. tomentosum, g. mustelinum and g. darwinii genomes. using the slaf-seq data, an integrated analysis was carried out. therefore, there were three main objectives in this study. first, interspecific variations were explored and characterized. second, the phylogenetic trees of tetraploid cotton species were reconstructed. third, the interspecific snp and indel markers were developed for future introgression breeding. in summary, the new dna polymorphisms presented in this study will dramatically increase the efficiency for future research into high-density interspecific mapping, introgression breeding, genetic dissection and gene utilization in cotton.  build variation resources of tetraploid cottons in the present study, a new variation resource was generated that had direct relevance to the current scenario of cotton breeding strategies. using an illumina high-throughput sequencing platform, a total of, 80-bp long paired-end reads were generated from the five tetraploid cotton species, including two cultivars and three wild species. after quality filtering,, high-quality reads were obtained, which varied from to million for the different materials. nearly 80% of the high-quality reads were mapped to the cotton reference genome, which covered approximately % of the total genome for each material. we also developed slaf tags among the five cotton species. slaf tags were compared, polymorphic slaf tags were identified, and it was observed that the number of polymorphic slaf tags were unevenly distributed between the a subgenome and d subgenome. the average depth for g. hirsutum cv. emian, g. barbadense acc. 3–, g. tomentosum, g. mustelinum and g. darwinii were,,,, and, respectively.table summary of sequence data and mapping statistics on the tm- genome mapq mapping quality of 30 identification and characteristics of snps and indels totally,, snps and indels were identified. additionally, after filtering out the low confidence, a total of polymorphisms, including snps and indels, were detected between the reference tm- genome and the five tetraploid cotton species. eventually, the total number of homozygous dna polymorphisms was for subsequent analysis, which included snps and indels. further, the snps and indels between each of the five tetraploid cotton species, and the reference tm- genome were also detected and showed that the gb_3- had the highest dna polymorphisms than other species, followed by gd; the dna polymorphisms of gh_e were the least. notably, the overlap of snps and indels revealed a moderate degree of overlap among the five tetraploid cotton species.fig. distribution and number of snps and indels detected in each of cotton species. a, b the number of snps and indels detected in each of cotton species are illustrated in the bar graphs by different colours as indicated. total number in each of cotton species is indicated by red bar graphs. the homozygous number in each cotton species is indicated by green bar graphs. c, d overlap of snps and indels across five species. e, f the distribution of snps and indels detected in five cotton species. i: the chromosomes; ii-vi: the total number of homozygous snps and indels of gd, gm, gt, gb_3– and gh_e respectively the snps and indels that were detected in each of the five cotton species were further analysed. the snps were classified as transitions and transversions based on nucleotide substitutions. among the transitions, the proportions of a/g were slightly higher than c/t transitions in gh_e gt, and gm, while the proportions of a/g was slightly lower than c/t transitions in gb_3- and gd. for transversions, the percentage of t/a transversions was relatively higher than others, namely, a/c, g/t and c/g. the ratio between transitions and transversions for gh_e gb_3- gt, gm and gd were,,, and, respectively. for indels, the length of the insertions ranged from 1 bp to 20 bp; however, the length distribution of deletions was 1 bp to 35 bp. the majority of indels were single nucleotides, accounting for %, di- to tetra-nucleotides were %, and the remaining % were ≥ 5 bp.fig. percentage of snps substitution types, the number of at-to-gc conversion in transposons, and length distribution of indels among five cotton species. a percentage of different snps substitution types in five cotton species. b the number of at-to-gc conversion in transposons in five cotton species. tes are indicated by blue bar graphs; non-tes are indicated by red bar graphs. c length distribution of indels among five cotton species. the number of insertions and deletions were shown on the y-axis in the bar graph. the numbers of various lengths are shown on the x-axis in the bar graph. gh_e22/ref_tm- gb/ref_tm- gt/ref_tm- gm/ref_tm- and gd/ref_tm- respectively indicates dna polymorphisms identified between gh_e gb_3- gt, gm, gd and the reference genome g. hirsutum acc. tm-1 the dna polymorphisms in transposable elements were also explored in the five cotton species. the percentages of snps/indels in tes for gh_e gb_3- gt, gm and gd were %/ %, %/ %, %/ %, %/ %, and %/ %, respectively. the number of snps and indels was calculated in te and non-te regions. further, through the analysis of snps and indels in tes, it was found that the dna polymorphisms mainly existed in long terminal repeat retrotransposons, particularly gypsy-type. additionally, at-to-gc conversion was larger in the te than non-te regions, which was significantly outnumbered in gb_3- than in other cotton species. genomic distribution of snps and indels across the reference tm- genome to survey the landscape of dna polymorphisms for each of the five cotton species, the distributions of dna polymorphisms were obtained. the number of homozygous snps and indels was detected for subsequent analysis, and their frequency varied across all chromosomes. herein, snps were observed to be unevenly distributed among the chromosomes and between the at and dt in each of cotton species, with more snps in the at than the dt. likewise, the distribution of indels among the chromosomes and between the at and dt was also not uniform. in addition, the total number of dna polymorphisms and their density varied across the cotton chromosomes in the five cotton species, and the total number of dna polymorphisms was found to be positively correlated with the chromosome length. among them, the maximum number of homozygous snps was found on chromosome a, and the minimum number was found on chromosome d. however, for indels, the maximum number was on chromosome a and the minimum number was on chromosome d. the distribution of gene densities across chromosomes from the tm- reference genome was also calculated, which showed the lowest gene density was in the middle of the chromosome. in addition, the genomic distribution of dna polymorphisms was determined by calculating the frequency of occurrence. the highest frequency of homozygous snps was on chromosome a, and the lowest frequency of snps was on chromosome d. similarly, the homozygous indel frequency was the highest on chromosome a, and the lowest on chromosome d. similarly, the total numbers of snps and indels were unevenly distributed amomg the chromosomes and between at and dt, with fewer dna polymorphisms in the d-subgenome than the a-subgenome. furthermore, a total of pairs of homeologous genes containing snps were identified between at and dt in the reference tm- genome, including pairs of homeologous genes with one containing snps in at, and the other gene did not contain snps in dt; pairs of homeologous genes with one containing snps in dt, and the other gene did not contain snps in at; and pairs of homeologous genes with genes containing snps in both at and dt.. characterization and functional significance of snps and indels using the annotation of the cotton reference genome, the distributions of snps and indels within various genomic features were revealed. in general, a similar distribution pattern between snps and indels was observed in the five species. the percentage of snps in genic regions ranged from to %, while the percentage of indels ranged from to %, among the five cotton species. however, a substantial portion of snps and indels were identified in the intergenic regions. a significant percentage of snps and indels were also detected in the 5 kb upstream and 5 kb downstream regulatory regions among the five cotton species.fig. annotation of snps and indels. a, b distribution of snps and indels in different genomic regions for the five cotton species. c, d distribution of snps and indels in different intergenic and genic regions among the five cotton species. in cds region, the number of synonymous and non-synonymous snps detected has also been shown among the five cotton species to further investigate the effect of snps in the coding sequence region, amino acid substitution was analysed and a large proportion was detected to be non-synonymous among the five cotton species. the ratio of non-synonymous to synonymous snps was approximately for gh_e for gb_3- for gt, for gm, and for gd. further, these non-synonymous snps were present in and genes for gh_e gb_3- gt, gm, and gd, respectively. the non-synonymous snps resided in the genes containing leucine-rich repeats, pentatricopeptide repeats, protein kinases, protein tyrosine kinases, wd domain repeats and nb-arc domains. among the five cotton species, large-effect snps and large-effect indels were found in a total of genes, which also contained many important functional domains, such as zinc-finger and p domains.fig. the ratio of non-synonymous to synonymous snps, and frequency of top pfam domains represented in the genes containing non-synonymous snp among the five cotton species. a the ratio of non-synonymous to synonymous snps among five cotton species. b-f frequency of top pfam domains represents in the genes containing non-synonymous snp in gh_e gb_3- gt, gm and gd, respectively gene ontology enrichment analysis revealed putative functions of the genes containing non-synonymous snp variations among the five cotton species. based on the reference genome annotation, a number of significantly enriched go terms in each of the three main categories were identified. the significantly enriched go terms, which were genes involved in biological processes, molecular function and cellular component of the go classification, revealed striking differences between the five cotton species. for example, in gb_3- the significant terms were related to reproduction, including multi-organism reproductive processes, recognition of pollen, pollen-pistil interaction, pollination, and reproductive processes. in other cotton species, they also had their own unique characteristics; for instance, oxidoreductase activity and oxidation-reduction process were significantly enriched in gt, which was associated with stress resistance. in gm, many significantly enriched categories were mainly associated with trna, which was involved in the synthesis of protein. significant go terms of gd were highly concentrated in cellular components, including the organelle outer membrane, ribonucleoprotein complex and intracellular parts.fig. significantly enriched go terms and phylogenetic tree among five cotton species. a the significantly enriched go terms are identified among five cotton species. b phylogenetic tree was constructed among five cotton species. c phylogenetic tree of a-homoeologs was constructed among five cotton species. d phylogenetic tree of d-homoeologs was constructed among five cotton species phylogenetic analysis of the cotton diploids and polyploids understanding the phylogenetic relationships among the various genomes is not only significant for genetics research and breeding but also essential for conducting comprehensive studies of the plant genomes, especially in genome organization, function and evolution. to further study the phylogenetic relationships in different cotton genomes, phylogenetic trees of the cotton diploids and polyploids were constructed. specifically, and unsurprisingly, the five tetraploid cotton species were classified into three branches, one consisting of gh_e and gt, one consisting of gm, and the third containing gb_3- and gd. additionally, two sets of snps were generated, and the phylogenetic trees for the at and dt homoeologs were constructed, independently. for the a-subgenome, d clade as one outgroup species, it showed a high concordance in the phylogeny of the cotton diploids and polyploids. the concordance analysis identified three main branches, which were the a-subgenomes of five tetraploid cotton clades, clade and d clade. compared with a-homoeologs, the clade was one outgroup species, and analysis results of d-homoeologs showed that the clades were more consistent with the sister-species relationship. the d clade had closer relationships with the d-subgenomes of the five tetraploid cotton clades than. interspecific variations for introgression genetics and breeding to broaden the genetic basis of upland cotton, introgression line populations were developed in our laboratory, which used gh_e as the receptor parent and gb_3- gm, gt, and gd as donor parents, respectively. the snps/indels identified in this study facilitated the introgression genetics and breeding as molecular markers. combined with our laboratory resources to investigate the extensive genetic diversity and the potential practicability, gh_e vs. gb_3- gh_e vs. gt, gh_e vs. gm and gh_e vs. gd were analysed. to gain insight into the differences in various combinations, the total number of homozygous snps/indels was identified between gh_e vs. gb_3- gh_e vs. gt, gh_e vs. gm and gh_e vs. gd, which were /, /, / and /, respectively, observing that gh_e vs. gb_3– had the higher dna polymorphisms than the others. further, overlaps of the snps/indels among the four combinations were shown with a venn diagram. the overlapped snps/indels between gh_e vs. gb_3– and gh_e vs. gd were larger than other two combinations. the unique dna polymorphisms were maximum in gh_e vs. gm, and minimum in gh_e vs. gt. intriguingly, the snps/indels were not evenly distributed cross the chromosomes. the snps/indels in the at were more frequent than those in the dt. the gh_e vs. gb_3– had the highest dna polymorphisms on each chromosome, followed in descending order by gh_e vs. gd, gh_e vs. gm and gh_e vs. gt, except for the chromosome a where gh_e vs. gd had the maximum snps and gh_e vs. gt had the maximum indels. the snps/indels on a were the maximum in the four combinations, which were / / / / in gh_e vs. gb_3– gh_e vs. gt, gh_e vs. gm and gh_e vs. gd, respectively. in contrast, the minimum number of snps were on d in gh_e vs. gb_3– d in gh_e vs. gt, d in gh_e vs. gm and d in gh_e vs. gd. however, the minimum number of indels, were on d in gh_e vs. gb_3– d in gh_e vs. gt, d in gh_e vs. gm and d in gh_e vs. gd, respectively.fig. distribution of snps and indels between four gh_e vs. gb_3- gh_e vs. gt, gh_e vs. gm and gh_e vs. gd. a, b the total number of snps and indels between different combinations are illustrated in the bar graphs. c, d the overlap of snps and indels. e, f distribution of snps and indels on each chromosome are illustrated in the bar graphs in different combinations the randomly selected snps on chromosome a were used to design primers to check the gel-based polymorphism. based on the sscp analysis, there were polymorphic markers between gh_e and gb_3– between gh_e and gt, between gh_e and gm, and between gh_e and gd. discussion in the present study, a high quality of, 80-bp long paired-end reads and slafs were generated using high-throughput slaf-seq. the average sequencing depth of each cotton species was ×. additionally, a total of homozygous snps, and homozygous indels were identified. interspecific dna polymorphisms for gb_3– gt, gm and gd occurred at a much higher rate compared to intraspecific dna polymorphisms for gh_e, especially in gb_3– which may be that they belong to different cotton species and the distant phylogenetic relationships. gb_3– and gd were found to have more overlap than the other three species, and a similar phenomenon also occurred in the following two species, gh_e and gt, which may be related to their genetic relationship. as is known, the conventional five polyploid cotton species have been divided into three branches by previous cotton researchers; one branch is g. hirsutum and g. tomentosum, the other is g. mustelinum, and the third is g. barbadense and g. darwinii. our analysis, as a new resource, provides further evidence on the accuracy of the previous studies in their phylogenetic relationships. for the snps, the non-synonymous/synonymous snp ratios in different cotton species were –, which was slight higher than ir64/pokkali and ir64/n in rice, and js- and upsm- in soybean. in addition, the frequency of transitions was found to be significantly higher than transversions. transition bias has been reported in rice, maize, chickpeas, and loquat. the proportions of c/t were slightly lower than a/g transitions in gh_e gt and gm, while the proportions of c/t were slightly higher than a/g transitions in gb_3– and gd, which may be associated with divergence times. similar observations were previously reported in g. hirsutum and g. barbadense. however, the percentage of t/a transversions was relatively higher than others, namely, a/c, g/t and c/g. similar observations were detected in rice and citrus, which is remained to be explained as an unexpected observation. the ratio of transitions and transversions for gh_e gb_3– gt, gm and gd were,,, and, respectively, which were slightly higher than g. hirsutum and g. barbadense of a previous report in cotton, ir64/pokkali and ir64/n in rice, and js- and upsm- in soybean. further, at-to-gc conversion in gb_3– significantly outnumbered gh_e followed in descending order by gd, gm and gt, which could provide the potential dna methylation sites, if occurring in tes. however, as for indels, the length distribution of indels was found to be larger than that in rice and less than that in soybean. intriguingly, indels in each of the individual five cottons also showed an upward bias towards single nucleotide insertions and deletions, which was similar to rice and soybean. the abundance of snps and indels were detected in 5 kb upstream and 5 kb downstream regulatory regions, which may be due to the non-coding regulatory regions having lower sequence conservation and levels of purifying selection pressure than coding regions. in animals, the variations in the cis-regulatory regions have demonstrated the importance in the regulation of gene expression. in plants, a few studies also have shown that such genetic variations play a major role in the regulation of gene expression and agronomic traits. functional annotation of these genes containing variations revealed putative characteristics in different cotton species. in gb_3– a number of the significant terms were related to reproduction, which suggested that the variations may result in the disruption of regulatory genes determining reproductive fertility and may help to explain the hybrid breakdown and low yield in offspring between g. barbadense and g. hirsutum. in previous studies, the genus gossypium has been researched broadly in phylogeny and a consensus phylogenetic tree was generated. the phylogenetic analysis has been reported in a previous study, and it was based on the targeted genes sequence capture method. here, we used high-throughput sequencing to reconstruct the phylogenetic trees that showed a high concordance with the historically hypothesized phylogeny, providing strong support for the phylogenetic relationships in diploid and allotetraploid cottons. interspecific introgression breeding is a traditional method and requires labour, time and funds; however, it can be expedited by molecular marker-assisted selection. in a recent study, the cotton breeding potential of introgression lines between g. hirsutum and g. barbadense was investigated for its yield and fibre quality improvement, including heterosis, combining ability and genetic effects. similarly, the introgression line populations were constructed in our laboratory with gh_e as the receptor and the other tetraploid species as donors. currently, only gel-based loci were identified between gh_e and gb_3–. as is known, the gel-based markers are time-consuming and have low efficiency. in this study, a huge number of snps and indels identified, plus the high throughput platform, will facilitate and speed up the genetics and breeding of the introgression lines. although few polymorphisms were detected from the randomly selected snp markers on a because the sscp is still a gel-based genotyping method and has low resolution, we do believe that more polymorphisms will be detected with a high-resolution detecting platform, such as the snp chip or genotyping by resequencing.  a large number of high-quality dna polymorphisms through slaf-seq were identified in five tetraploid cotton species. further comprehensive analysis of dna polymorphisms provided valuable insights into the different tetraploid cotton species, especially between cultivated and wild cottons. to better utilize the novel characteristics of drought tolerance, defence responses and salt stress tolerance, etc. in wild cottons, numerous interspecific variations were identified here, which will be very practical in future introgression breeding. overall, the comprehensive data generated in this study provided insights into cotton evolution and a resource for future researches in high-density interspecific mapping, introgression breeding improvement in cotton. 
more importantly, the identification of numerous interspecific genetic variations provides the basis and is very practical for future introgression breeding. the results presented herein provide a valuable genomic resource for new insights into the genetics and breeding of cotton. 3– based on the new data sets, we reconstructed phylogenetic trees that showed a high concordance to the phylogeny of diploid and polyploid cottons. the comprehensive analysis of dna polymorphisms provided invaluable insights into the different tetraploid cotton species. using slaf-seq, a large number of dna polymorphisms were identified. intriguingly, gene ontology enrichment analysis revealed that a number of significant terms were related to reproduction in g. barbadense acc. a large amount of interspecific genetic variations were identified, and some of them were validated by the single-strand conformation polymorphism method, which will be applied in introgression genetics and breeding with g. hirsutum cv. understanding the genomic variation of tetraploid cotton species is important for exploitation of the excellent characteristics of wild cotton and for improving the diversity of cotton in breeding. a total of, reads, specific length amplified fragment tags and high-quality dna polymorphisms were identified using specific length amplified fragment sequencing, including snps and indels between the reference genome and the five tetraploid cotton species. however, the discovery of dna polymorphisms in tetraploid cotton genomes has lagged behind other important crops.
 cotton, as the world’s most utilized textile fibre source, is an important, economically valuable crop worldwide. understanding the genomic variation of tetraploid cotton species is important for exploitation of the excellent characteristics of wild cotton and for improving the diversity of cotton in breeding. however, the discovery of dna polymorphisms in tetraploid cotton genomes has lagged behind other important crops.  a total of, reads, specific length amplified fragment tags and high-quality dna polymorphisms were identified using specific length amplified fragment sequencing, including snps and indels between the reference genome and the five tetraploid cotton species. intriguingly, gene ontology enrichment analysis revealed that a number of significant terms were related to reproduction in g. barbadense acc. 3– based on the new data sets, we reconstructed phylogenetic trees that showed a high concordance to the phylogeny of diploid and polyploid cottons. a large amount of interspecific genetic variations were identified, and some of them were validated by the single-strand conformation polymorphism method, which will be applied in introgression genetics and breeding with g. hirsutum cv. emian as the receptor and the other species as donors.  using slaf-seq, a large number of dna polymorphisms were identified. the comprehensive analysis of dna polymorphisms provided invaluable insights into the different tetraploid cotton species. more importantly, the identification of numerous interspecific genetic variations provides the basis and is very practical for future introgression breeding. the results presented herein provide a valuable genomic resource for new insights into the genetics and breeding of cotton. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords tetraploid cottonslaf-seqsingle-nucleotide polymorphismsinsertions/deletionsphylogenythe genetically modified organisms breeding major project of china2016zx08009lin zhongxu issue-copyright-statement© the author 2017 