BACKGROUND
cotton  is the leading fiber crop worldwide.
there are four cultivated cotton species, two diploids from africa-asia, g. herbaceum l.  and gossypium arboreum l. , and two tetraploids from americas, g. hirsutum l.  and g. barbadense l. .
at present, g. hirsutum is the most widely cultivated cotton species, accounting for more than 95% of the world cotton production , followed by g. barbadense  and g. arboreum , while g. herbaceum is rarely cultivated.
g. hirsutum has the characteristics of high yield, broad adaptation, and medium fiber quality.
g. barbadense typically has a longer growing period, matures later, and produces smaller bolls that give a yield significantly lower than that of g. hirsutum.
in spite of these drawbacks, however, g. barbadense possesses superior fiber properties, which makes it an important raw material for high-grade or special cotton fiber textiles  <cit> .
recently, a number of genome resources have been developed from the genus gossypium including the construction of high-density tetraploid cotton genetic linkage maps  <cit> , construction of large-insert bac libraries  <cit> , and analyses of expressed sequence tags  related to fiber development  <cit> .
of these, est analysis is not only the most efficient approach for gene discovery, but also an effective approach for the development of polymorphic dna markers.
as of jan.  <dig>   <dig>  approximately  <dig>  cotton est sequences are available in genbank ests database .
among them,  <dig>  ests were from g. hirsutum,  <dig>  from g. raimondii,  <dig>  from g. arboreum, while only  <dig>  and  <dig> were from g. barbadense and g. herbaceum, respectively.
compared with great amount of est resources from g. hirsutum, ests for g. barbadense are relatively scarce and have hindered the exploration of its economic importance.
large scale transcript sequences offer an efficient resource for targeted marker development.
in cotton, ssr markers have been mined widely based on existing sequences data from different cotton species  <cit> , and applied widely in characterizing variations of genes, genome-wide mapping  <cit> , and as a tool for marker-assisted selection  <cit> .
in addition to ssr loci, the distribution of insertion-deletion  and single nucleotide polymorphisms  variations are more widespread in the whole genome.
recently, indel and snps are increasing in their application in studies of cotton genomics  <cit> .
as an application,  <dig> snp markers were mapped in  <dig> recombinant inbred lines from a cross between tm- <dig> and 3– <dig>  <cit> .
information from these large-scale markers will serve as a foundation for constructing high-density genetic maps, cloning and mapping of important genes, improving gene prediction and annotation, and elucidating the interspecific divergence in different cotton species.
in this study, we constructed two non-normalized cdna libraries from fibers and ovules mixtures at − <dig> to  <dig> days post-anthesis  and fibers at  <dig> to  <dig> dpa, for efficient generation of unique ests from g. barbadense cv.
hai <dig>  in total,  <dig>  ests were generated by a large-scale 5′ end single-pass sequencing of randomly picked cdna clones from the two libraries.
upon adding the released  <dig>  g. barbadense est sequences in the national center for biotechnology information , est assembly for g. barbadense was performed and functional categories of unigenes were assigned.
further, by meta-analysis of fiber ests between g. barbadense and g. hirsutum, some putative g. barbadense-specific expressed fiber genes and indel loci of orthologous or homoeologous variation between/within g. barbadense and g. hirsutum were developed, which put the foundation for promoting the identification of genes related to superior agronomic traits, genetic differentiation and comparative genomic studies between g. hirsutum and g. barbadense.
RESULTS
generation of g. barbadense fiber ests
two high-quality g. barbadense cotton fiber cdna libraries were constructed using hai <dig> developing fiber tissues, one for − <dig> to  <dig> dpa  and the other for  <dig> to  <dig> dpa  respectively.
once each library was constructed, the lengths of inserts were identified by amplifying randomly  <dig> clones with universal t <dig> primers.
the results showed that most of the inserts were among  <dig> - <dig>  kb, with an average length of  <dig>  kb for the gb <dig> library and  <dig>  kb for the gb <dig> library.
based on this,  <dig>  cdna clones  were randomly isolated and sequenced to generate ests from the 5′ end.
after the removal of low quality sequences below the q <dig> threshold, short sequences less than  <dig> bp in length, and those with vector sequences and poly-a tails,  <dig>  high-quality sequences with scores ranging from  <dig> to  <dig> were obtained for further analysis.
most of the sequences ranged from  <dig> to  <dig> bp, with an average sequence length of  <dig> bp with the longest est at  <dig> bp.
all high quality sequences have been deposited in the dbest division of genbank under the accession numbers jk <dig> to jk <dig>
assembly and function annotation of g. barbadense fiber ests
to reduce the redundancy and facilitate the process of gene annotation and mapping, all g. barbadense est sequences available in genbank and herein were combined for est processing and assembly analysis.
of  <dig>  ests,  <dig>  were singletons and the remaining were assembled into  <dig>  contigs, yielding  <dig>  unigenes for further analysis.
the average sequence length of unigenes was  <dig> bp, and the longest unigene was  <dig>  bp.
we subjected  <dig>  unigenes to blastx against the non-redundant  database to identify their similarities with known proteins .
in total,  <dig>   sequences shared the high similarity with public protein sequences, of these,  <dig>   matched known proteins and  <dig>  were unknown or hypothetical proteins.
the remaining  <dig>  sequences had no similarities in the tested databases and might be either 3′ or 5′ untranslated regions  of genes with too short of a coding sequence, or novel genes  <cit> .
function classification and metabolic pathway analysis of unigenes
of the  <dig>  unigenes,  <dig>  were mapped to the go hierarchy with characterized biochemical and physiological functions using gene ontology annotation .
of these,  <dig>  were mapped to go categories, and were involved with biological processes, molecular functions, and cellular components.
at a second level, the majority of the go terms were grouped into cellular process  and metabolic process  categories within biological processes, binding  and catalytic activity  categories within molecular functions, and cell part , cell , and organelle  categories within cellular components.
at a third level, the different categories were further mined.
in biological processes , the metabolic process was the most enriched as “cell metabolic process” ranked first, accounting for  <dig> % of the unigenes involved with biological processes.
“primary metabolic processes” was second with  <dig> %, and third was the “macromolecular metabolic processes” with  <dig> % of the unigenes.
within the molecular function category , the highly enriched go terms in “binding” were distributed in the following way:  <dig> % represented “protein binding”,  <dig> % were “nucleic acid binding”,  <dig> % were “ion binding”, and  <dig> % were involved with “nucleotide binding”.
the top three highly enriched go terms in “catalytic activity” were  <dig> % belonging to “hydrolase activity”,  <dig> % on transferase activity, and  <dig> % on oxidoreductase activity.
furthermore, the cellular component mostly took place in “intracellular” with  <dig> %, “intracellular part” with  <dig> % and “intracellular organelle” with  <dig> % .
in addition, by a comparison with the kyoto encyclopedia of genes and genomes database , the metabolic-related enzymes encoded by  <dig>  unigenes were also located in metabolic maps based on the kegg pathway classification .
the metabolism category was the richest , followed by the categories of environmental information processing , cellular processes , genetic information processing , organismal systems  and metabolic pathways related to human diseases .
in the category of metabolism, the mapped enzymes were mostly involved in carbohydrate metabolism , amino acid metabolism , and energy metabolism .
in the category of cellular process,  <dig> % were related to cell growth and death .
in the category of organismal systems, the endocrine system , and immune system  were major contributors.
in the category of environmental information processing,  <dig> % of unigenes were involved in signal transduction.
in the category of genetic information processing, folding, sorting and degradation was the majority , followed by transcription  and replication and repair .
at a third level, the mappings of metabolism associated with oxidative phosphorylation , glycolysis/gluconeogenesis , and starch and sucrose metabolism  were further tagged .
transcript abundance and highly enriched genes during fiber development
the level of representation in a cdna library generally correlates with transcript abundance in non-normalized conditions.
we assessed the distribution of contigs based on the number of assembled ests.
compared with singletons, most of the contigs consisted of  <dig> to  <dig> ests and a few were comprised of more than  <dig> ests .
the top  <dig> most highly expressed unigenes and their annotations can be found in the table  <dig>  including non-specific lipid-transfer protein family, copper-zinc superoxide dismutase, e <dig>  fiber protein , which was responsible for fiber development and growth in previous studies  <cit> .
comparisons with the g. hirsutum ests
ests data from g. hirsutum were downloaded from genbank.
poly a/t sequences occurring at the end of ests and vectors were removed from original sequences.
as a result, a set of cleaned ests generated from g. hirsutum was clustered and assembled using the cap <dig> assembly program.
the clustering of ests from g. hirsutum yielded  <dig>  unigenes, with  <dig>  contigs and  <dig>  singletons, respectively.
to find putative specific expression genes in g. barbadense,  <dig>  unigenes from g. barbadense were used to detect similarity with g. hirsutum unigenes using the tblastx program with a significant similarity  threshold of 10- <dig>  as a result,  <dig>  unigenes from g. barbadense were detected the similarity with g. hirsutum.
the remaining  <dig> sequences had no hits against the g. hirsutum unigenes database.
the putative  <dig> specific expression genes in g. barbadense were assigned functions by blastx against the nr protein database.
of these,  <dig> sequences had one or more similarities with proteins in the non-redundant protein database, however,  <dig> sequences had no significant similarity.
according to ests abundance, we also confirmed the expression specificity/predominance of these genes by selecting seven putative genes at random to analyze their expression patterns at  <dig> and  <dig> dpa in the fiber developmental process between g. hirsutum acc.
tm- <dig> and g. barbadense cv.
hai <dig> using rt-pcr method .
these putative g. barbadense specific expressed genes encoded mainly for the vacuolar h+-atpase c subunit, profiling, anthocyanidin synthase, proline-rich protein, kinetochore protein, etc.
the detailed annotation information and abundance of specific expressed genes in g. barbadense were listed in supplementary additional file 4: table s <dig>
large-scale discovery and confirmation of indels between/within g. barbadense and g. hirsutum ests
indels can lead to differentiation of species, and small indels were widely distributed in the different species.
in the study,  <dig>  g. hirsutum and  <dig>  g. barbadense est sequences were assigned as two distinct data sets for all possible indels loci based on in silico pcr strategy.
a set of primer pairs from g. barbadense est sequences were run against the g. hirsutum est sequences dataset based on in silico pcr analysis with a threshold less than  <dig> mismatch bases.
in total,  <dig>  indel loci derived from  <dig>  g. barbadense ests were identified using a custom indel_pipeline.pl script.
after redundancy analysis,  <dig>  unique est indel loci were detected for further study.
of these,  <dig>  contained only  <dig> bp indel size,  <dig>  with  <dig> to  <dig> bp,  <dig>  with  <dig> to  <dig> bp, and  <dig>  had more than  <dig> bp.
the average length of the indels was  <dig> bp, with the longest indels up to  <dig> bp .
based on in silico pcr analysis,  <dig>  indel markers were developed with indel lengths > = <dig> bp .
to verify the accuracy and efficiency of indel markers,  <dig> randomly selected indel primer pairs were synthesized for the polymorphisms detection between/within g. barbadense and g. hirsutum using two tetraploids, g. hirsutum acc.
tm- <dig> and g. barbadense cv.
hai <dig>  and two diploids, g. herbaceum and g. raimondii.
of  <dig> randomly selected indel loci,  <dig> primer pairs showed distinguishable orthologous and/or homoeologous polymorphisms , with  <dig> %  orthologous polymorphisms  and  <dig> %  homoeologous polymorphisms .
in addition,  <dig> indel primer pairs were unable to produce detectable polymorphic electrophoresis bands between g. hirsutum and g. barbadense , and  <dig> indel primer pairs failed to amplify pcr products in tested cotton accessions.
amplification information for  <dig> indels markers is included in the electronic supplementary additional file 7: table s <dig>
for confirming the accuracy of indel loci in the transcriptional regions, the amplification products from  <dig> indel primer pairs with different types  were individually recovered from polyacrylamide gels and sequenced.
the length of all amplicons was equal to or greater than that expected from g. barbadense.
sequences alignment analysis from orthologous and/or homoeologous loci amplified by  <dig> indel primer pairs , showed that, except for size difference in intron region, indels from the corresponding exon regions were all confirmed the existence, just like the result of in silico analysis.
to undetectable polymorphic electrophoresis bands from indel <dig>  sequencing analysis also showed the complete consistency of orthologous loci from four cotton species, indicating that the indels from the transcriptional level might be due to alternative splicing events with no difference at the genome level.
as for failure pcr amplification between g. hirsutum and g. barbadense, this might be related with the reason that the primer sequence from est spanned the intron–exon boundary of the genome and need to be confirmed in the future.
insertion-deletion est variation in tetraploid cultivated cotton species might be responsible for fiber quality
during fiber developmental stage, cellulose synthesis and secondary cell wall  thickening is one of the most important events for fiber quality.
to further confirm whether insertion-deletion est variation in tetraploid cultivated cotton species might be related to fiber quality difference between g. hirsutum and g. barbadense,  <dig>  candidate indel loci with indel lengths > = <dig> bp were further assessed by blastx against tair  database  to identify genes involved in secondary wall synthesis.
as a result,  <dig> indel loci with significant function similarity with known genes associated with secondary cell wall synthesis process of arabidopsis thaliana trichomes .
of them,  <dig> were targeted to encode the cellulose synthase related genes, such as cellulose synthase, ras-related gtp-binding family, germin-like protein, beta- <dig> -glucanase, glycosyl transferase , fasciclin-like arabinogalactan , chitinase-like protein , arabinogalactan , and sucrose synthase .
other nine were targeted to encode the transcriptional factor related genes, such as myb-related protein, homeodomain protein, and zinc finger protein.
in these predicted genes, some family genes have been identified as a major role in cotton fiber secondary wall synthesis, such as ghcesa   <cit> , ghfla  <cit> , ghagp   <cit> , rac <dig> <cit> , ghctl   <cit> , and sus   <cit> .
the insertion-deletion est variation of these key genes related to fiber secondary cell wall synthesis might play the important roles in fiber quality in tetraploid cultivated cotton species.
discussion
g. barbadense ests - indispensable cotton breeding sequences resource
g. barbadense has been used more commonly as a gene donor for the high-quality cotton fiber trait.
the hybrids produced with g. hirsutum are expected to have the desirable characteristics including a high yield, exceptional fiber length, fineness and strength.
however, the lack of g. barbadense genomic resources has seriously hampered the exploration of modern functional genomic approaches for selective breeding purposes.
expressed sequence tag  sequencing strategies are efficient in identifying a large number of genes expressed in a given tissue, which are particularly relevant when no genomic data are available  <cit> .
prior to this publication, only one g. barbadense dataset was publicly released, which included  <dig>  ests from g. barbadense cv.
3– <dig> fiber tissues  <cit> .
in our study,  <dig>  high quality sequences were generated from two non-normalized cdna libraries prepared by using a mixture of hai <dig> fibers and ovules at − <dig> to  <dig> days post-anthesis  and fibers at  <dig> to  <dig> dpa.
the assembly resulted in approximately  <dig>  unigenes, involving  <dig>  contigs and  <dig>  singletons.
about 83% of contigs and singletons reported here were matched with high confidence to the nr protein database.
these sets of ests are a large contribution to enhance the number of g. barbadense ests available in public databases.
in addition, a high level of representation in a cdna library generally correlates with transcript abundance in non-normalized conditions.
abundance analysis revealed that large clusters existed , such as lipid-transfer protein precursors, copper-zinc superoxide dismutase, e <dig>  and arabinogalactan protein, which was in agreement with many previously reported cotton fiber key genes.
orford and timmis  suggested lipid transfer protein gene family might play an important role in fiber development  <cit> .
john et al.
described the isolation of the e <dig> gene that was expressed preferentially in cotton fiber tissues on the 15th and 24th days after flowering  <cit> .
superoxide dismutases, commonly found in all oxygen-consuming organisms, are a class of metal enzymes, which can remove superoxide anion radicals.
hu et al.
reported that the copper/zinc superoxide dismutase mrna expressing level showed regular changing in the fiber developmental process  <cit> .
other highly expressed genes such as anthocyanidin reductase  <cit> , pectate lyase  <cit> , arabinogalactan protein  <cit>  and fiber protein fb <dig> also were found to play an important role in the cotton fiber development.
comparative analysis between g. hirsutum and g. barbadense
comparative analysis of ests from fiber development between g. hirsutum and g. barbadense will facilitate the identification of genes related to the important agronomic traits such as fiber growth and development.
tblastx comparisons between g. barbadense fiber sequences and those from g. hirsutum released publicly in genbank revealed that  <dig>  unigenes  from g. barbadense were detected as similarities with those from g. hirsutum.
the remaining  <dig> sequences had no hits against the g. hirsutum unigenes database.
though the e-value setting, utrs regions' interference, and insufficient coverage of g. hirsutum ests database might inevitably cause certain false positives in the process of similarity analyses, the expression predominance of selected seven unigenes in g. barbadense was confirmed experimentally by comparing  <dig> and  <dig> dpa fiber tissues between g. hirsutum acc.
tm- <dig> and g. barbadense cv.
hai <dig>  based on the fact that g. barbadense possesses superior fiber quality, key genes related to fiber quality in g. barbadense may be exploited in cotton molecular breeding programs for improving fiber quality.
the predominantly expressed sequences found here will provide an insight into exploring interspecific fiber quality divergence between g. barbadense and g. hirsutum.
further function confirmation of these sequences will be verified by experimental analyses.
elite genes can be mined effectively by large-scale discovery of indels between g. barbadense and g. hirsutum
currently, there is an increasing focus on polymorphisms of short insertions and deletions  types in genomic research of model species such as humans  <cit> , mice  <cit>  and rice  <cit> .
indels have been recognized as an abundant source of genetic markers that are widely spread across the genome, and can be genotyped with simple size separation based on pcr amplification and gel electrophoresis analyses.
thus, the indel molecular markers were highly informative sequence-based markers suitable for high-density map construction, genome-wide association studies, genomic selection, and alignment of the whole genome sequence information.
in allopolyploids, mutational variation may arise between homologous sequences within individual subgenomes and between homoeologous sequences among subgenomes, in addition to paralogous variation between duplicated gene copies  <cit> .
when information regarding genomic sequences is not rich, indels polymorphism of tetraploid cotton species obtained by computational analysis might include both homologous and homoeologous loci variation.
nevertheless, mining these differences will be useful for verifying the function of the target genes and elucidating evolution in polyploid.
in the study, we first reported the feasibility of using large-scale ests sequence data from different cotton species to extract putative insertion and deletion polymorphisms.
<dig>  indel loci were identified between g. barbadense and g. hirsutum, and  <dig>  potential gh-gb indel markers were further developed for genotyping and evolutionary research.
further, the indel-derived ests/genes related to fiber growth and development was detected.
in the paper, we found some genes involved in secondary cell wall synthesis process with potential est indel loci, such as those encoding cellulose synthase, sucrose synthase, beta- <dig> -glucanase, glycosyl transferase, fasciclin-like arabinogalactan protein, arabinogalactan protein, chitinase-like protein, and other transcriptional factor genes involved in secondary cell wall synthesis.
in arabidopsis, indel mutations of cellulose synthase genes in secondary cell wall synthesis had been known to cause structural weaknesses in vascular bundles presumably due to cellulose deficiency  <cit> .
so, the deep research of these genes will accelerate our understanding to fiber growth and development process in cotton.
due to the fact that all indel loci were discovered based on transcriptional level and all g. barbadense ests were collected from different stages of fiber development, gh-gb indels enriched many genes related to cotton fiber development, which could greatly promote their deep research for structural and expressional differences related to the fiber development processes in g. barbadense and g. hirsutum.
by randomly selecting  <dig> indels primer pairs to conduct pcr analysis of genomic dna, with the sequencing confirmation from  <dig> indel loci in four cotton species, more than 70% indel polymorphic loci  were detected effectively.
from the result, not only homologous and homoeologous loci difference were found, but also some evolutionary events could be inferred.
the diversity among the allotetraploid cottons such as orthologous and/or homoeologous polymorphisms could be traced back to ancient diploid ancestors by banding patterns analysis, which indicated the independent evolution or different degrees of colonization by comparing g. barbadense and g. hirsutum, with two diploids, g. herbaceum and g. raimondii, as controls  <cit> .
taken together, indel molecular detection provides a new tool for effectively mining genes related to superior agronomic traits, and selecting appropriate g. barbadense or g. hirsutum germplasm in cotton breeding.
it also offers a novel model for the study of the origin, evolution, and genetic differentiation of g. barbadense and g. hirsutum and their adaptation to various environmental changes.
CONCLUSIONS
transcriptome analysis from different tissues and organs provides the basis for functional genomics research.
in the present work, we constructed two g. barbadense fiber cdna libraries and obtained  <dig>  high-quality sequences.
a systematic analysis and utilization of ests was further performed including assembling, annotation, go classification and comparative analysis with g. hirsutum.
the resulting dataset yielded nearly  <dig>  putative unigenes, from which over 80% had similarities with publicly available proteins.
furthermore, putative ests indels loci involved in the orthologous and/or homoeologous difference between/within g. barbadense and g. hirsutum were discovered by in silico analysis and confirmed by experimental analysis.
the large-scale g. barbadense ests in the study were a significant contribution for public g. barbadense ests databases, and either expression or candidate est indel difference will provide a new tool for effectively mining genes related to superior agronomic traits.
these data will provide a solid foundation for molecular breeding, functional genomics studies, and comparative genomics analysis in gossypium.
with the continuous growth of sequence information from non-model organisms such as cotton, we suggest that indels will be a crucial source for next-generation mapping of key genes in these accessions.

 cotton is the leading fiber crop worldwide. there are four cultivated cotton species, two diploids from africa-asia, g. herbaceum l. and gossypium arboreum l., and two tetraploids from americas, g. hirsutum l. and g. barbadense l.. at present, g. hirsutum is the most widely cultivated cotton species, accounting for more than 95% of the world cotton production, followed by g. barbadense and g. arboreum, while g. herbaceum is rarely cultivated. g. hirsutum has the characteristics of high yield, broad adaptation, and medium fiber quality. g. barbadense typically has a longer growing period, matures later, and produces smaller bolls that give a yield significantly lower than that of g. hirsutum. in spite of these drawbacks, however, g. barbadense possesses superior fiber properties, which makes it an important raw material for high-grade or special cotton fiber textiles. recently, a number of genome resources have been developed from the genus gossypium including the construction of high-density tetraploid cotton genetic linkage maps, construction of large-insert bac libraries, and analyses of expressed sequence tags related to fiber development. of these, est analysis is not only the most efficient approach for gene discovery, but also an effective approach for the development of polymorphic dna markers. as of jan. approximately cotton est sequences are available in genbank ests database. among them, ests were from g. hirsutum, from g. raimondii, from g. arboreum, while only and were from g. barbadense and g. herbaceum, respectively. compared with great amount of est resources from g. hirsutum, ests for g. barbadense are relatively scarce and have hindered the exploration of its economic importance. large scale transcript sequences offer an efficient resource for targeted marker development. in cotton, ssr markers have been mined widely based on existing sequences data from different cotton species, and applied widely in characterizing variations of genes, genome-wide mapping, and as a tool for marker-assisted selection. in addition to ssr loci, the distribution of insertion-deletion and single nucleotide polymorphisms variations are more widespread in the whole genome. recently, indel and snps are increasing in their application in studies of cotton genomics. as an application, snp markers were mapped in recombinant inbred lines from a cross between tm- and 3–. information from these large-scale markers will serve as a foundation for constructing high-density genetic maps, cloning and mapping of important genes, improving gene prediction and annotation, and elucidating the interspecific divergence in different cotton species. in this study, we constructed two non-normalized cdna libraries from fibers and ovules mixtures at − to days post-anthesis and fibers at to dpa, for efficient generation of unique ests from g. barbadense cv. hai in total, ests were generated by a large-scale 5′ end single-pass sequencing of randomly picked cdna clones from the two libraries. upon adding the released g. barbadense est sequences in the national center for biotechnology information, est assembly for g. barbadense was performed and functional categories of unigenes were assigned. further, by meta-analysis of fiber ests between g. barbadense and g. hirsutum, some putative g. barbadense-specific expressed fiber genes and indel loci of orthologous or homoeologous variation between/within g. barbadense and g. hirsutum were developed, which put the foundation for promoting the identification of genes related to superior agronomic traits, genetic differentiation and comparative genomic studies between g. hirsutum and g. barbadense.  generation of g. barbadense fiber ests two high-quality g. barbadense cotton fiber cdna libraries were constructed using hai developing fiber tissues, one for − to dpa and the other for to dpa respectively. once each library was constructed, the lengths of inserts were identified by amplifying randomly clones with universal t primers. the results showed that most of the inserts were among - kb, with an average length of kb for the gb library and kb for the gb library. based on this, cdna clones were randomly isolated and sequenced to generate ests from the 5′ end. after the removal of low quality sequences below the q threshold, short sequences less than bp in length, and those with vector sequences and poly-a tails, high-quality sequences with scores ranging from to were obtained for further analysis. most of the sequences ranged from to bp, with an average sequence length of bp with the longest est at bp. all high quality sequences have been deposited in the dbest division of genbank under the accession numbers jk to jk assembly and function annotation of g. barbadense fiber ests to reduce the redundancy and facilitate the process of gene annotation and mapping, all g. barbadense est sequences available in genbank and herein were combined for est processing and assembly analysis. of ests, were singletons and the remaining were assembled into contigs, yielding unigenes for further analysis. the average sequence length of unigenes was bp, and the longest unigene was bp. we subjected unigenes to blastx against the non-redundant database to identify their similarities with known proteins. in total, sequences shared the high similarity with public protein sequences, of these, matched known proteins and were unknown or hypothetical proteins. the remaining sequences had no similarities in the tested databases and might be either 3′ or 5′ untranslated regions of genes with too short of a coding sequence, or novel genes. function classification and metabolic pathway analysis of unigenes of the unigenes, were mapped to the go hierarchy with characterized biochemical and physiological functions using gene ontology annotation. of these, were mapped to go categories, and were involved with biological processes, molecular functions, and cellular components. at a second level, the majority of the go terms were grouped into cellular process and metabolic process categories within biological processes, binding and catalytic activity categories within molecular functions, and cell part, cell, and organelle categories within cellular components. at a third level, the different categories were further mined. in biological processes, the metabolic process was the most enriched as “cell metabolic process” ranked first, accounting for % of the unigenes involved with biological processes. “primary metabolic processes” was second with %, and third was the “macromolecular metabolic processes” with % of the unigenes. within the molecular function category, the highly enriched go terms in “binding” were distributed in the following way: % represented “protein binding”, % were “nucleic acid binding”, % were “ion binding”, and % were involved with “nucleotide binding”. the top three highly enriched go terms in “catalytic activity” were % belonging to “hydrolase activity”, % on transferase activity, and % on oxidoreductase activity. furthermore, the cellular component mostly took place in “intracellular” with %, “intracellular part” with % and “intracellular organelle” with %. in addition, by a comparison with the kyoto encyclopedia of genes and genomes database, the metabolic-related enzymes encoded by unigenes were also located in metabolic maps based on the kegg pathway classification. the metabolism category was the richest, followed by the categories of environmental information processing, cellular processes, genetic information processing, organismal systems and metabolic pathways related to human diseases. in the category of metabolism, the mapped enzymes were mostly involved in carbohydrate metabolism, amino acid metabolism, and energy metabolism. in the category of cellular process, % were related to cell growth and death. in the category of organismal systems, the endocrine system, and immune system were major contributors. in the category of environmental information processing, % of unigenes were involved in signal transduction. in the category of genetic information processing, folding, sorting and degradation was the majority, followed by transcription and replication and repair. at a third level, the mappings of metabolism associated with oxidative phosphorylation, glycolysis/gluconeogenesis, and starch and sucrose metabolism were further tagged. transcript abundance and highly enriched genes during fiber development the level of representation in a cdna library generally correlates with transcript abundance in non-normalized conditions. we assessed the distribution of contigs based on the number of assembled ests. compared with singletons, most of the contigs consisted of to ests and a few were comprised of more than ests. the top most highly expressed unigenes and their annotations can be found in the table  including non-specific lipid-transfer protein family, copper-zinc superoxide dismutase, e fiber protein, which was responsible for fiber development and growth in previous studies. comparisons with the g. hirsutum ests ests data from g. hirsutum were downloaded from genbank. poly a/t sequences occurring at the end of ests and vectors were removed from original sequences. as a result, a set of cleaned ests generated from g. hirsutum was clustered and assembled using the cap assembly program. the clustering of ests from g. hirsutum yielded unigenes, with contigs and singletons, respectively. to find putative specific expression genes in g. barbadense, unigenes from g. barbadense were used to detect similarity with g. hirsutum unigenes using the tblastx program with a significant similarity threshold of 10- as a result, unigenes from g. barbadense were detected the similarity with g. hirsutum. the remaining sequences had no hits against the g. hirsutum unigenes database. the putative specific expression genes in g. barbadense were assigned functions by blastx against the nr protein database. of these, sequences had one or more similarities with proteins in the non-redundant protein database, however, sequences had no significant similarity. according to ests abundance, we also confirmed the expression specificity/predominance of these genes by selecting seven putative genes at random to analyze their expression patterns at and dpa in the fiber developmental process between g. hirsutum acc. tm- and g. barbadense cv. hai using rt-pcr method. these putative g. barbadense specific expressed genes encoded mainly for the vacuolar h+-atpase c subunit, profiling, anthocyanidin synthase, proline-rich protein, kinetochore protein, etc. the detailed annotation information and abundance of specific expressed genes in g. barbadense were listed in supplementary additional file 4: table s large-scale discovery and confirmation of indels between/within g. barbadense and g. hirsutum ests indels can lead to differentiation of species, and small indels were widely distributed in the different species. in the study, g. hirsutum and g. barbadense est sequences were assigned as two distinct data sets for all possible indels loci based on in silico pcr strategy. a set of primer pairs from g. barbadense est sequences were run against the g. hirsutum est sequences dataset based on in silico pcr analysis with a threshold less than mismatch bases. in total, indel loci derived from g. barbadense ests were identified using a custom indel_pipeline.pl script. after redundancy analysis, unique est indel loci were detected for further study. of these, contained only bp indel size, with to bp, with to bp, and had more than bp. the average length of the indels was bp, with the longest indels up to bp. based on in silico pcr analysis, indel markers were developed with indel lengths > = bp. to verify the accuracy and efficiency of indel markers, randomly selected indel primer pairs were synthesized for the polymorphisms detection between/within g. barbadense and g. hirsutum using two tetraploids, g. hirsutum acc. tm- and g. barbadense cv. hai and two diploids, g. herbaceum and g. raimondii. of randomly selected indel loci, primer pairs showed distinguishable orthologous and/or homoeologous polymorphisms, with % orthologous polymorphisms and % homoeologous polymorphisms. in addition, indel primer pairs were unable to produce detectable polymorphic electrophoresis bands between g. hirsutum and g. barbadense, and indel primer pairs failed to amplify pcr products in tested cotton accessions. amplification information for indels markers is included in the electronic supplementary additional file 7: table s for confirming the accuracy of indel loci in the transcriptional regions, the amplification products from indel primer pairs with different types were individually recovered from polyacrylamide gels and sequenced. the length of all amplicons was equal to or greater than that expected from g. barbadense. sequences alignment analysis from orthologous and/or homoeologous loci amplified by indel primer pairs, showed that, except for size difference in intron region, indels from the corresponding exon regions were all confirmed the existence, just like the result of in silico analysis. to undetectable polymorphic electrophoresis bands from indel sequencing analysis also showed the complete consistency of orthologous loci from four cotton species, indicating that the indels from the transcriptional level might be due to alternative splicing events with no difference at the genome level. as for failure pcr amplification between g. hirsutum and g. barbadense, this might be related with the reason that the primer sequence from est spanned the intron–exon boundary of the genome and need to be confirmed in the future. insertion-deletion est variation in tetraploid cultivated cotton species might be responsible for fiber quality during fiber developmental stage, cellulose synthesis and secondary cell wall thickening is one of the most important events for fiber quality. to further confirm whether insertion-deletion est variation in tetraploid cultivated cotton species might be related to fiber quality difference between g. hirsutum and g. barbadense, candidate indel loci with indel lengths > = bp were further assessed by blastx against tair database to identify genes involved in secondary wall synthesis. as a result, indel loci with significant function similarity with known genes associated with secondary cell wall synthesis process of arabidopsis thaliana trichomes. of them, were targeted to encode the cellulose synthase related genes, such as cellulose synthase, ras-related gtp-binding family, germin-like protein, beta- -glucanase, glycosyl transferase, fasciclin-like arabinogalactan, chitinase-like protein, arabinogalactan, and sucrose synthase. other nine were targeted to encode the transcriptional factor related genes, such as myb-related protein, homeodomain protein, and zinc finger protein. in these predicted genes, some family genes have been identified as a major role in cotton fiber secondary wall synthesis, such as ghcesa, ghfla, ghagp, rac, ghctl, and sus. the insertion-deletion est variation of these key genes related to fiber secondary cell wall synthesis might play the important roles in fiber quality in tetraploid cultivated cotton species. discussion g. barbadense ests - indispensable cotton breeding sequences resource g. barbadense has been used more commonly as a gene donor for the high-quality cotton fiber trait. the hybrids produced with g. hirsutum are expected to have the desirable characteristics including a high yield, exceptional fiber length, fineness and strength. however, the lack of g. barbadense genomic resources has seriously hampered the exploration of modern functional genomic approaches for selective breeding purposes. expressed sequence tag sequencing strategies are efficient in identifying a large number of genes expressed in a given tissue, which are particularly relevant when no genomic data are available. prior to this publication, only one g. barbadense dataset was publicly released, which included ests from g. barbadense cv. 3– fiber tissues. in our study, high quality sequences were generated from two non-normalized cdna libraries prepared by using a mixture of hai fibers and ovules at − to days post-anthesis and fibers at to dpa. the assembly resulted in approximately unigenes, involving contigs and singletons. about 83% of contigs and singletons reported here were matched with high confidence to the nr protein database. these sets of ests are a large contribution to enhance the number of g. barbadense ests available in public databases. in addition, a high level of representation in a cdna library generally correlates with transcript abundance in non-normalized conditions. abundance analysis revealed that large clusters existed, such as lipid-transfer protein precursors, copper-zinc superoxide dismutase, e and arabinogalactan protein, which was in agreement with many previously reported cotton fiber key genes. orford and timmis suggested lipid transfer protein gene family might play an important role in fiber development. john et al. described the isolation of the e gene that was expressed preferentially in cotton fiber tissues on the 15th and 24th days after flowering. superoxide dismutases, commonly found in all oxygen-consuming organisms, are a class of metal enzymes, which can remove superoxide anion radicals. hu et al. reported that the copper/zinc superoxide dismutase mrna expressing level showed regular changing in the fiber developmental process. other highly expressed genes such as anthocyanidin reductase, pectate lyase, arabinogalactan protein and fiber protein fb also were found to play an important role in the cotton fiber development. comparative analysis between g. hirsutum and g. barbadense comparative analysis of ests from fiber development between g. hirsutum and g. barbadense will facilitate the identification of genes related to the important agronomic traits such as fiber growth and development. tblastx comparisons between g. barbadense fiber sequences and those from g. hirsutum released publicly in genbank revealed that unigenes from g. barbadense were detected as similarities with those from g. hirsutum. the remaining sequences had no hits against the g. hirsutum unigenes database. though the e-value setting, utrs regions' interference, and insufficient coverage of g. hirsutum ests database might inevitably cause certain false positives in the process of similarity analyses, the expression predominance of selected seven unigenes in g. barbadense was confirmed experimentally by comparing and dpa fiber tissues between g. hirsutum acc. tm- and g. barbadense cv. hai based on the fact that g. barbadense possesses superior fiber quality, key genes related to fiber quality in g. barbadense may be exploited in cotton molecular breeding programs for improving fiber quality. the predominantly expressed sequences found here will provide an insight into exploring interspecific fiber quality divergence between g. barbadense and g. hirsutum. further function confirmation of these sequences will be verified by experimental analyses. elite genes can be mined effectively by large-scale discovery of indels between g. barbadense and g. hirsutum currently, there is an increasing focus on polymorphisms of short insertions and deletions types in genomic research of model species such as humans, mice and rice. indels have been recognized as an abundant source of genetic markers that are widely spread across the genome, and can be genotyped with simple size separation based on pcr amplification and gel electrophoresis analyses. thus, the indel molecular markers were highly informative sequence-based markers suitable for high-density map construction, genome-wide association studies, genomic selection, and alignment of the whole genome sequence information. in allopolyploids, mutational variation may arise between homologous sequences within individual subgenomes and between homoeologous sequences among subgenomes, in addition to paralogous variation between duplicated gene copies. when information regarding genomic sequences is not rich, indels polymorphism of tetraploid cotton species obtained by computational analysis might include both homologous and homoeologous loci variation. nevertheless, mining these differences will be useful for verifying the function of the target genes and elucidating evolution in polyploid. in the study, we first reported the feasibility of using large-scale ests sequence data from different cotton species to extract putative insertion and deletion polymorphisms. indel loci were identified between g. barbadense and g. hirsutum, and potential gh-gb indel markers were further developed for genotyping and evolutionary research. further, the indel-derived ests/genes related to fiber growth and development was detected. in the paper, we found some genes involved in secondary cell wall synthesis process with potential est indel loci, such as those encoding cellulose synthase, sucrose synthase, beta- -glucanase, glycosyl transferase, fasciclin-like arabinogalactan protein, arabinogalactan protein, chitinase-like protein, and other transcriptional factor genes involved in secondary cell wall synthesis. in arabidopsis, indel mutations of cellulose synthase genes in secondary cell wall synthesis had been known to cause structural weaknesses in vascular bundles presumably due to cellulose deficiency. so, the deep research of these genes will accelerate our understanding to fiber growth and development process in cotton. due to the fact that all indel loci were discovered based on transcriptional level and all g. barbadense ests were collected from different stages of fiber development, gh-gb indels enriched many genes related to cotton fiber development, which could greatly promote their deep research for structural and expressional differences related to the fiber development processes in g. barbadense and g. hirsutum. by randomly selecting indels primer pairs to conduct pcr analysis of genomic dna, with the sequencing confirmation from indel loci in four cotton species, more than 70% indel polymorphic loci were detected effectively. from the result, not only homologous and homoeologous loci difference were found, but also some evolutionary events could be inferred. the diversity among the allotetraploid cottons such as orthologous and/or homoeologous polymorphisms could be traced back to ancient diploid ancestors by banding patterns analysis, which indicated the independent evolution or different degrees of colonization by comparing g. barbadense and g. hirsutum, with two diploids, g. herbaceum and g. raimondii, as controls. taken together, indel molecular detection provides a new tool for effectively mining genes related to superior agronomic traits, and selecting appropriate g. barbadense or g. hirsutum germplasm in cotton breeding. it also offers a novel model for the study of the origin, evolution, and genetic differentiation of g. barbadense and g. hirsutum and their adaptation to various environmental changes.  transcriptome analysis from different tissues and organs provides the basis for functional genomics research. in the present work, we constructed two g. barbadense fiber cdna libraries and obtained high-quality sequences. a systematic analysis and utilization of ests was further performed including assembling, annotation, go classification and comparative analysis with g. hirsutum. the resulting dataset yielded nearly putative unigenes, from which over 80% had similarities with publicly available proteins. furthermore, putative ests indels loci involved in the orthologous and/or homoeologous difference between/within g. barbadense and g. hirsutum were discovered by in silico analysis and confirmed by experimental analysis. the large-scale g. barbadense ests in the study were a significant contribution for public g. barbadense ests databases, and either expression or candidate est indel difference will provide a new tool for effectively mining genes related to superior agronomic traits. these data will provide a solid foundation for molecular breeding, functional genomics studies, and comparative genomics analysis in gossypium. with the continuous growth of sequence information from non-model organisms such as cotton, we suggest that indels will be a crucial source for next-generation mapping of key genes in these accessions. 
by gel electrophoresis combined with sequencing verification, % candidate indel loci were reconfirmed orthologous and/or homoeologous loci polymorphisms using g. hirsutum acc tm- and g. barbadense cv hai blastx result showed among indel loci, with significant function similarity with known genes associated with secondary wall synthesis process, indicating the important roles in fiber quality in tetraploid cultivated cotton species. sequence comparisons and indel markers development will lay the groundwork for promoting the identification of genes related to superior agronomic traits, genetic differentiation and comparative genomic studies between g. hirsutum and g. barbadense. however, a systematic analysis and utilization of cdna sequences from g. barbadense fiber development remains understudied. a total of high quality sequences were generated from two non-normalized cdna libraries prepared by using a mixture of g. barbadense hai fibers and ovules. sequences comparison between g. barbadense and g. hirsutum revealed that unigenes from g. barbadense were detected the similarity with those released publicly in g. hirsutum, however, the remaining sequences had no hits against g. hirsutum unigenes database. furthermore, putative ests indels loci involved in the orthologous and/or homoeologous differences between/within g. barbadense and g. hirsutum were discovered by in silico analyses, and indel markers were developed by ests with more than five insertions or deletions. the molecular functions of these unigenes were mostly related to binding and catalytic activity, and carbohydrate, amino acid, and energy metabolisms were major contributors among the subsets of metabolism. of those, were matched to known proteins and were assigned to functional categories. cotton is the leading fiber crop worldwide.
 cotton is the leading fiber crop worldwide. gossypium barbadense is an important species of cotton because of its extra-long staple fibers with superior luster and silkiness. however, a systematic analysis and utilization of cdna sequences from g. barbadense fiber development remains understudied.  a total of high quality sequences were generated from two non-normalized cdna libraries prepared by using a mixture of g. barbadense hai fibers and ovules. after assembly processing, a set of unigenes were obtained. of those, were matched to known proteins and were assigned to functional categories. the molecular functions of these unigenes were mostly related to binding and catalytic activity, and carbohydrate, amino acid, and energy metabolisms were major contributors among the subsets of metabolism. sequences comparison between g. barbadense and g. hirsutum revealed that unigenes from g. barbadense were detected the similarity with those released publicly in g. hirsutum, however, the remaining sequences had no hits against g. hirsutum unigenes database. furthermore, putative ests indels loci involved in the orthologous and/or homoeologous differences between/within g. barbadense and g. hirsutum were discovered by in silico analyses, and indel markers were developed by ests with more than five insertions or deletions. by gel electrophoresis combined with sequencing verification, % candidate indel loci were reconfirmed orthologous and/or homoeologous loci polymorphisms using g. hirsutum acc tm- and g. barbadense cv hai blastx result showed among indel loci, with significant function similarity with known genes associated with secondary wall synthesis process, indicating the important roles in fiber quality in tetraploid cultivated cotton species.  sequence comparisons and indel markers development will lay the groundwork for promoting the identification of genes related to superior agronomic traits, genetic differentiation and comparative genomic studies between g. hirsutum and g. barbadense. 