Structure, expression differentiation and evolution of duplicated fiber developmental genes in Gossypium barbadense and G. hirsutum

Background Both Gossypium hirsutum and G. barbadense probably originated from a common ancestor, but they have very different agronomic and fiber quality characters. Here we selected 17 fiber development-related genes to study their structures, tree topologies, chromosomal location and expression patterns to better understand the interspecific divergence of fiber development genes in the two cultivated tetraploid species. Results The sequence and structure of 70.59% genes were conserved with the same exon length and numbers in different species, while 29.41% genes showed diversity. There were 15 genes showing independent evolution between the A- and D-subgenomes after polyploid formation, while two evolved via different degrees of colonization. Chromosomal location showed that 22 duplicate genes were located in which at least one fiber quality QTL was detected. The molecular evolutionary rates suggested that the D-subgenome of the allotetraploid underwent rapid evolutionary differentiation, and selection had acted at the tetraploid level. Expression profiles at fiber initiation and early elongation showed that the transcripts levels of most genes were higher in Hai7124 than in TM-1. During the primary-secondary transition period, expression of most genes peaked earlier in TM-1 than in Hai7124. Homeolog expression profile showed that A-subgenome, or the combination of A- and D-subgenomes, played critical roles in fiber quality divergence of G. hirsutum and G. barbadense. However, the expression of D-subgenome alone also played an important role. Conclusion Integrating analysis of the structure and expression to fiber development genes, suggests selective breeding for certain desirable fiber qualities played an important role in divergence of G. hirsutum and G. barbadense.


Background
Cotton (Gossypium spp.) is the world's most important fiber crop plant. While most of the > 50 Gossypium species are diploid (n = 13), five are allopolyploids (n = 26), originating from an interspecific hybridization event between A-and D-genome diploid species. Humans have independently domesticated four different species for their fiber, two of which are diploids, Gossypium herbaceum and G. arboreum, and two are allopolyploids, G. hirsutum and G. barbadense [1].
Alhough G. hirsutum and G. barbadense probably originated from a single hybridization event between A-and D-diploid species, the two have very different agronomic and fiber quality characteristics. The high yield potential and diverse environmental and production system adaptability of G. hirsutum make it the most widely cultivated species, accounting for about 97% of the world's cotton fiber [2]. G. barbadense is a more modern species possessing superior fiber quality. Novel alleles are responsible for the improved fiber quality in G. barbadense. Despite its higher fiber quality, however, the narrow adaptation range and low yield of G. barbadense limit its cultivation. The two Gossypium species are sexually compatible, although partial sterility, longer maturity, and hybrid breakdown are often observed in later generations [3]. Nonetheless, the introgression of favorable alleles from G. barbadense to G. hirsutum would likely improve the fiber quality of G. hirsutum while simultaneously maintaining its high fiber yield [4].
The cotton fiber is a single cell without the complex cell division and multicellular development that develops from ovule's epidermal cells. Fiber development occurs in four distinct, but overlapping stages: initiation, elongation, secondary wall synthesis, and maturation [5]. To date, many of the genes predominantly expressed in cotton fiber development have been isolated and characterized. Gh14-3-3L was found to be predominantly expressed during early fiber development, and may be involved in regulating fiber elongation [6]. Yoder et al. [7] defined pectate lyase (PEL) as a cell wall modification enzyme. GhPel was found to play an essential role in fiber cell elongation by degradation of the deesterified pectin for cell wall loosening [8]. Ruan et al. [9] suggested the sucrose synthase gene (Sus) played an important role in the initiation and elongation of singlecelled fibers by influencing carbon partitioning to cellulose synthesis. GhBG (β-1,4-glucosidase), one of three cellulases, was specifically expressed in fiber cells and plays an important role in degradation of the primary cell wall and promotion of secondary cell wall synthesis [10]. Cotton CelA1 and CelA2 genes, encoding the catalytic subunit of cellulose synthase, are expressed at high levels during active secondary wall cellulose synthesis in developing cotton fibers [11]. Two cotton Rac genes, GhRacA and GhRacB, expressed in the fibers at the initiation and elongation stages, might play an important role in early fiber development [12]. In addition, several genes are expressed specifically or preferentially in fibers [13][14][15][16][17][18], although their exact functional roles remain unclear.
In theory, there are two homologs in tetraploid cotton species, representing descendants from the A-genome and D-genome donors at the time of polyploidy formation. The goals of this study were to: 1) better understand the genetic basis of cotton fiber development, 2) identify the structural difference of duplicated genes, and 3) reveal the expression and evolution of fiber quality differences between upland and sea-island cotton. To complete this study, we selected 17 fiber development genes accessioned in National Center for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov) to study structure and expression differences of the two cultivated tetraploid species. To investigate their frame and sequence divergence, we initially cloned these genes in the genome DNA of the G. hirsutum accession TM-1, the G. barbadense cultivar Hai7124, and their two putative diploid progenitors. The chromosomal locations of each homeolog of several studied genes, having effective single nucleotide polymorphism (SNP) or amplification polymorphism loci between TM-1 and Hai7124 were determined by linkage analysis in allotetraploid cotton using an interspecific BC 1 mapping population (TM-1×Hai7124)×TM-1] [19][20][21][22]. Finally, expression patterns of each gene and each homeolog were explored. A more thorough understanding of interspecific divergence of cotton will provide a solid foundation from which key fiber quality genes may be exploited in cotton molecular breeding.

Results
Sequence and structure analysis of fiber development genes The orthologs of each of the 17 genes were cloned and sequenced (GenBank accession numbers (GQ340731-GQ340736, HQ142989-HQ143048, and HQ143055-HQ143090; Table 1). Phylogenetic groupings and sequence comparisons allowed the copy number for all genes, except Exp1 in the Hai7124 cultivar, to be independently isolated with a single copy from the diploid species and two homeologs from each tetraploid species for each gene. There was a single copy of Exp1 in the diploid species and two distinct copies in TM-1, however, the Exp1 sequence from Hai7124 was of only one type, though more than 10 clones were selected randomly to sequence. This result was further validated by a different primer pair of this gene (see Additional file 1: Supplemental Table S1 for list of primer pairs). The sequence from Hai7124 has a closer relationship with G. raimondii than with G. herbaceum. Further, southern blotting of Exp1 performed on the four species, showed two distinct hybridizing bands after digestion with EcoRI and HindI in Hai7124 and TM-1, and one hybridizing band in G. herbaceum and G. raimondii (Figure not shown), which indicated that Exp1 had two copies in both TM-1 and Hai7124. Combining sequence and southern blot analysis, the homeolog of Exp1 in the A-subgenome of Hai7124 was colonized to a type resembling that of the D-subgenome via nonreciprocal homoeologous exchange [23].
The lengths of the genomic DNA sequences isolated from the four species varied from 1020 bp (Exp1) to 6126 bp (CelA3) ( Table 2). Based on the alignments between the genomic DNA and the cDNA sequences in orthologs, twelve genes (70.59%) shared the same intron/exon structures in different genomes, and variations in length were mainly caused by insertion/deletion events within introns (Additional file 2: Supplemental Figure S1A). The remaining five genes (29.41%), CIPK1, CAP, BG, ManA2 and CelA3, produced some structure differences caused by different exon length or numbers (Additional file 2: Supplemental Figure S1B).
Seventeen gene trees were constructed using the NJ method to distinguish the duplicated genes independent of evolution or local interlocus recombination after tetraploid formation. Two major clades, one including G. herbaceum and the A-subgenomes of TM-1 and Hai7124, the other including G. raimondii and the D-subgenomes of TM-1 and Hai7124, were formed for 15 genes (Additional file 3: Supplemental Figure S2A). High bootstrap values supported duplicated genes independent of evolution after tetraploid formation. Two genes were determined to have local interlocus recombination or colonization after tetraploid formation (Additional file 3: Supplemental Figure S2B). ACT1 from G. raimondii was more closely related to ACT1s from G. herbaceum and the A-subgenomes than with ACT1s from the D-subgenomes. This relationship suggests that  ACT1s from the D-subgenomes evolved at an accelerated rate, relative to ACT1s from the A-subgenomes. The Exp1 sequence in Hai7124 was closer to that found in G. raimondii.
To locate all 17 homeolog gene pairs on our backbone genetic map [23], the subgenome-specific PCR primers (Additional file 4: Supplemental Table S2) and singlenucleotide amplified polymorphisms (SNAP) primers (Additional file 5: Supplemental Table S3) were used to detect polymorphisms between TM-1 and Hai7124. Polymorphic primer pairs were also used to survey 138 individuals of the ([(TM-1×Hai7124)×TM-1]) BC 1 mapping population (Table 3). Eight gene pairs were located on their corresponding homeologous chromosomes, and each of six pairs were located on one of their homeologous chromosomes, while three pair of genes, Exp, Exp1 and CelA1 could not be mapped because no available polymorphic loci were found between TM-1 and Hai7124. A large body of data was compiled by integrating previously reported cotton fiber quality quantitative trait locus (QTL) [24][25][26][27][28][29][30][31][32][33][34] with the 22 identified fiber quality-related genes within 20 cM. Most genes had at least one fiber quality QTL; some had several (Table 3), indicating important fiber quality roles.

Rates of sequence evolution
With only a few exceptions, purifying selection, as indicated by Ka/Ks < 1, appears to be in place for most of pairwise comparisons ( Table 4). The exceptions include G. herbaceum (A 1 genome, A 1 ) vs G. hirsutum A-subgenome (Ath) of POD2 and D 5 vs Dth of RacA, which demonstrated positive selection (Ka/Ks > 1), as well as D 5 vs Dtb of RacA, A 1 vs Atb of 14-3-3L, A 1 vs Ath and A 1 vs Atb of CAP, which had exceptionally strong positive selection (Ka/Ks >> 1).
Because there were no nucleotide substitutions in A 1 vs Atb of RacA, A 1 vs Atb of EXP, D 5 vs Dth and D 5 vs Dtb of LTP3, and all RacBs pairs, it was not possible to compare their evolutionary rates. Although those particular "not applicable" comparisons (NA vs NA) were Table 3 Integration analysis of chromosomal locations of genes and fiber quality QTL reported in other studies

Gene
Subgenome Chromosomal location QTL associated with specific chromosome a  excluded from our analysis, the "NAs" were considered zero when they were compared with others whose Ks were not NA. In all pairwise comparisons of nucleotide diversity for each gene between subgenomes within a species [32 pairs, 16 in TM-1 and 16 in Hai7124 (RacB was excluded)], 62.5% (10 in TM-1 and 10 in Hai7124) had a higher evolutionary rate in the D-subgenome than in the A-subgenome. Furthermore, in the 16 gene pairs (excluding RacB) from the A-subgenomes of TM-1 and Hai7124, 62.5% (10 out of 16) had a higher evolutionary rate in TM-1 than in Hai7124, 31.25% (5 of 15) and were reversed and 6.25% (1 of 16, ACT1) showed an equivalent evolutionary rate between TM-1 and Hai7124. Similarly, in the 15 gene pairs (RacB and LTP3 were excluded) from the D-subgenomes of TM-1 and Hai7124, 60% (9 of 15) had a higher evolutionary rate in TM-1 than in Hai7124, 26.67% (4 of 15) and were reversed and 13.33% (2 of 15, Pel and Exp) showed an equivalent evolutionary rate between TM-1 and Hai7124. Phylogenetic relationships are reflected in the nucleotide substitution results (Additional file 3: Supplemental Figure S2). Based on branch length, all of homeologs from the two tetraploid species had unequal rates of sequence evolution following allopolyploid formation. The rates at which the deviations occurred in allopolyploids are sufficient to generate branch length inequality between the A-and D-subgenomes [35].
Ka/Ks ratio comparisons showed that selection had altered the molecular evolutionary rate of some genes due to allopolyploid formation. Four genes, Pel, RacA, Exp and Sus1, in TM-1, and five genes, Pel, CIPK1, 14-3-3L, CAP and CelA3, in Hai7124, yielded higher Ka/Ks ratios in A-At, D-Dt and At-Dt comparisons than in the A-D comparison, indicating that selection for some genes related with fiber development had acted at the tetraploid level.

Differential expression fiber development genes
After the specificity of homeolog-specific primer pairs were confirmed by PCR amplification of genomic DNA from G. herbaceum (A-genome), G. raimondii (D-genome), TM-1 and Hai7124 (Figure 1), their homeolog transcripts in young tetraploid cotton fiber were further detected by qPCR analysis. The relative expression values at 10 different fiber development stages were obtained by combining the homeolog transcripts of each gene at the same stage. Expression for the 17 genes could be broken down into five categories (Additional file 6: Supplemental Figure S3 Figure S3D); 5) the whole fiber developmental period, such as CelA3, CAP, ACT1, RacA and CIPK1 (Additional file 6: Supplemental Figure S3E). In the last category, however, transcript preference was shown at some stages. For example, CelA3 and CAP were expressed preferentially at the fiber elongation and secondary cell wall thickening stages (8-23 DPA), but had moderate expression at 0-5 DPA. Gene expression differences in TM-1 and Hai7124 were further clarified by statistical analysis of least signification difference (LSD). Greater expression in Hai7124 than in TM-1 was observed for 14-3-3L except at 20 DPA, and for CelA3 except at 5, 17 and 23 DPA. Other gene transcripts showed different expression advantages in the two cotton species at various fiber developmental stages.
At fiber initiation and early elongation (0-8 DPA), most genes, including Exp, ManA2, Sus1, RacB, CelA3, CAP and RacA, had significantly higher expression levels Figure 1 Amplification products in four cotton species using subgenome-specific qPCR primer pairs. First line includes amplified results from A-genome specific primers; second line includes amplified results from D-genome specific primers. "M" represents marker, "A" represents G. herbaceum var. africanum, "D" represents G. raimondii, "T" represents G. hirsutum acc. TM-1, "H" represents G. barbadense cv. Hai7124. Numbers represent the sizes of the makers (bp). in Hai7124 than in TM-1. During fiber elongation (5-14 DPA), the expression profiles of genes preferentially expressed during that period were either biased to TM-1 or Hai7124 or were equally expressed between the two. Five genes, Exp1, Pel, CAP, CIPK1 and RacA, were expressed preferentially in TM-1 or equally between TM-1 and Hai7124, except at 8 DPA, that these same genes showed significantly greater expression levels in Hai7124; LTP3 and ACT1 showed significantly higher expression levels in TM-1 than in Hai7124; expression of CelA3 was higher in Hai7124.
During primary-secondary cell wall transition (17-23 DPA), peak expression occurred earlier in TM-1 than in Hai7124 for most genes (CAP, CelA3, CIPK1, CEL, BG, RacB, Sus1 and 14-3-3L). ACT1 and RacA expressed equally in TM-1 and Hai7124 at 17 and 20 DPA, but significantly greater in Hai7124 at 23 DPA. The extended fiber development period, as indicated by higher expression at a later DPA, may help explain why G. barbadense has an extra long staple cotton. One gene, CelA1, showed no significant expression difference between TM-1 and Hai7124.
The 17 fiber development genes were clustered into three comparison patterns between TM-1 and Hai7124. First, homeologs for CelA3, Exp, Exp1 and CIPK1 in both TM-1 and Hai7124 were equally expressed between the A-and D-subgenomes in the preferentially-expressed stages (Figure 2). Of these, Exp1 had equal transcript levels from the two homeologs in TM-1 and Hai7124, with two distinguishable copies in TM-1 and two undistinguishable copies in Hai7124. These data were consistent with the fact that the duplicated loci for Exp1 in Hai7124 had the same sequence as the D-subgenome ( Figure 2).
Second, the transcripts of 11 genes, CEL, Pel, Sus1, 14-3-3L, RacA, CelA1, ManA2, RacB, CAP, LTP3 and POD2, were A-or D-subgenome biased (Figure 3). Among these, CEL, Pel, Sus1, 14-3-3L and RacA were A-subgenome biased and CelA1, ManA2 and RacB were D-subgenome biased in both TM-1 and Hai7124 at all stages. The transcripts of the homeologs of CAP, LTP3 and POD2 were significantly altered in the preferentially expressed stages in TM-1 and Hai7124. In TM-1, the transcripts of CAP and LTP3 were significantly A-subgenome biased. However, the transcripts of the two genes in Hai7124 were equivalently expressed at most stages, only D-subgenome bias in LTP3 in the primary-secondary cell wall transition period detected. Expression of POD2 was A-subgenome biased at 0, 3 and 10 DPA and D-subgenome biased at 1 DPA in TM-1. In Hai7124, POD2 expression was A-subgenome biased at 0, 3 and 10 DPA and D-subgenome biased at 1, 5 and 8 DPA.
Third, BG was significantly (P < 0.001) affected only from the A-subgenome, and ACT1 was significantly (P < Figure 2 Q-PCR analysis for homeologous expression of genes expressed equally between A-and D-subgenomes. Significant values were obtained by comparison between the two subgenomes. * P < 0.05, ** P < 0.01. See Table 2 for abbreviation designations. Vertical bars represented standard deviation (STD).  (Figure 4). Based on the comparison patterns and the structural analysis of the two genes, we proposed that the homeolog of BG from the D-subgenome might be silenced and that of ACT1 from the A-subgenome may have novel roles in other species (neofunctionalization).
Differences between TM-1 and Hai7124 in transcriptome contributions of the subgenome at key fiber developmental stages were detected. During initiation and early elongation of the fiber, 10 gene transcriptions showed greater expression levels in Hai7124 than in TM-1. Of those, the D-subgenome contributed higher amounts of ACT1, RacB and Man2, while the A- subgenome contributed higher amounts of Sus1, CIPK1 and 14-3-3L. CelA3, Exp, CAP and RacA were equally supplied by both subgenomes.
At the primary-secondary transition period, the expression of 10 genes, CelA3, CAP, ACT1, RacA, CIPK1, 14-3-3L, Sus1, RacB, CEL and BG, occurred earlier in TM-1 than in Hai7124. The transcripts of ACT1 and RacB were mainly from D-subgenomes; those of BG, CEL and Sus1 were mainly from A-subgenomes, and other genes were from both A-and D-subgenomes. Based on these data we inferred that the expression accumulation of the A-subgenome, or the combination of A-and D-subgenomes, played critical roles in fiber quality divergence of G. hirsutum and G. barbadense. However, the expression of D-subgenome alone also played an important role.

Evolutionary fate of duplicated genes
For each gene that was studied, allopolyploid species should have two homelogs, representing descendants from the A-genome and D-genome donors at the time of polyploidy formation. Cronn et al. [35] indicated that most duplicated genes in allopolyploid cotton evolved independently of each other. Our phylogenetic analyses support this hypothesis, and the independent evolution of several genes was distinctively evident in their structure, in our study. For example, CAPs and RacBs had the same structure between each diploid and its counterpart in allopolyploid cotton (A-and At-subgenome, D-and Dt-subgenome), but the different structures were apparent in the A-D comparison (Additional file 2: Supplemental Figure S1B). Though expression of ManA2 from the At-subgenome of Hai7124 ceased rather early in the growth process, the structure difference between the A-, At-subgenomes and D-, Dt-subgenomes was also distinct. The fact that the structure of the At-and Dt-subgenomes mirrored their putative ancestral diploid species suggested the difference may have occurred before allopolyploid formation and evolved independently in allopolyploid cotton. CelA3s from At-subgenome of TM-1 and Hai7124 displayed the same mutation, which altered their coding regions, indicating not only independent evolution, but also parallel evolution between TM-1 and Hai7124. This change, however, was not detected in their putative ancestral diploid species, suggesting accelerated evolution of CelA3 in the At-subgenome after allopolyploid formation. Though most genes independently evolved in allopolyploid cotton, there were some exceptions. For example, Exp1 from At-subgenome were colonized in Hai7124 by Dt-subgenomes.
Relative to expression, duplicate genes can follow one of three evolutionary paths. First, one copy may evolve into a nonfunctional pseudogene [37][38][39][40][41]. Second, the multiple copies can contribute to an increase in the gene expression level [42,43] or both copies can suffer mutations but the combined action of both gene copies is necessary to maintain original function and expression levels (subfunctionalization) [40,44,45]. Third, one copy may gain a novel beneficial function (neofunctionalization) that is selectively maintained within the genome [40,[46][47][48]. We measured homeolog-specific contributions to the transcriptome in allopolyploid cotton fiber by Q-PCR analysis. Because the majority (64.70%) of diagnostic genes exhibited subgenome-specific bias to the A or D-subgenome, subgenomebiased expression in cotton fiber developmental stages was considered commonplace. This result was consistent with previous studies [49][50][51][52][53][54][55]. Most of genes in our study exhibited the same expression bias in the two cultivated cotton species, TM-1 and Hai7124. However, some inconsistencies were detected in three genes (CAP, LTP3 and POD2), suggesting that these genes may have had different roles in the interspecific divergence between G. hirsutum and G. barbadense. Artificial selection by humans of certain desirable fiber traits may have also influenced G. hirsutum and G. barbadense genetic structure [55].
Synthesizing structures and expression profiles of the duplicates, their possible fates are inferred. BG accumulated solely in A-subgenome transcripts (D-subgenome silenced), in both TM-1 and Hai7124 (Figure 4). BGs obtained from the D-subgenomes of TM-1 and Hai7124 had a nucleotide deletion and a nonsense mutation, respectively, which altered the ORFs ( Table 2). The structure difference suggests that BG in the D-subgenomes of TM-1 and Hai7124 may be pseudogenes. On the other hand, while CAP and CelA3 had different A-and D-subgenome structures in both TM-1 and Hai7124 (Additional file 2: Supplemental Figure S1B), their A-and D-subgenome expression profiles were active (Figure 2, 3). Therefore, the duplicated genes of CAP and CelA3 may be subfunctionalized. Similarly, the functions of duplicated genes from CEL, Sus1, 14-3-3L, RacA, RacB, Exp, Exp1, CIPK1, CelA1, Pel, ManA2, LTP3 and POD2 were also subfunctionalized ( Figure 3). Because ACT1 transcripts of A subgenomes could not be detected at all stages of fiber development (Figure 4), they may have evolved new functions.

Domestication of allopolyploid cotton
Numerous plant species have been selectively bred over the course of human social evolution [56]. Allopolyploid cotton species are believed to have formed about 1-2 million years ago, by hybridization between a maternal Old World diploid A-genome G. herbaceum [57] and paternal New World diploid D-genome G. raimondii [57][58][59]. The allotetraploid lineage gave rise to five extant tetraploid species, including G. barbadense and G. hirsutum, known for their superior fiber quality and high yield, respectively. In the present study, the Ka/Ks ratios among four cotton species indicated that selection of fiber development genes occurred at the tetraploid level. By comparing the nucleotide diversity between TM-1 and Hai7124 within the same subgenome, most genes (62.5% in A-subgenome and 60% in D-subgenome) had a higher evolutionary rate in TM-1 than in Hai7124, which may be associated with longer and more frequent cultivation of TM-1. Given these data, we propose that diversity evolution between A-and D-subgenomes within a species or between TM-1 and Hai7124 within the same subgenome was due to both natural and artificial selection pressure [55].

Gene expression differences between TM-1 and Hai7124
G. hirsutum and G. barbadense are two domesticated cotton species possessing very different agronomic and fiber quality characteristics with G. barbadense having superior fiber quality. Rapp et al. (2010) studied the transcriptomes of cotton fibers from wild and domesticated accessions (G. hirsutum) and found that human selection during the initial domestication and subsequent crop improvement had resulted in a biased upregulation of components of the transcriptional network during fiber development [60]. In this study, of the 17 fiber development-related genes, 14 had the similar expression pattern and three that did not, in TM-1 and Hai7124 (Figure 2,  3, 4). Of three genes, the transcripts of homeologs were significantly A-or D-subgenome biased in TM-1. However, in Hai7124, homeolog transcripts were equally expressed between the two subgenomes or D-subgenome biased. Though 14 genes had the same expression patterns between TM-1 and Hai7124, the relative expression levels were different at most stages. While the same Aor D-biased or equal expression profile in the two cultivated cotton species might be related to functional partitioning of genomic contributions during cellular development after allopolyploid formation, significant alternation of homoelog A/D ratio and expression difference at the same fiber developmental time points between G. hirsutum and G. barbadense indicated that domestication for different fiber qualities may play an important role in fiber quality divergence of G. hirsutum and G. barbadense.
In previous study, fiber growth curves have shown longer fiber elongation phases in domesticated G. hirsutum than that in wild G. hirsutum, and further comparative gene expression profiling of isolated cotton fibers over a developmental time course of fiber differentiation indicated that domesticated TM-1 displayed a much higher level of transcriptional variation between the sampled time points than the wild accession did [60]. In the study, the expression peak of transcripts in most genes was earlier in TM-1 than in Hai7124, especially at the primary-secondary transition period, which indicated that most genes related to fiber development expressed longer and more intensely in Hai7124. Resulting differences in mRNA levels may lead to changes in enzyme activity, further contributing to phenotypic differences between the two cotton species. Several genes that are differentially expressed in TM-1 and Hai7124 should be further mined. The 14-3-3 protein is an important regulatory protein. Shi et al. [6] proposed that the Gh14-3-3L transcripts are highly accumulated during early cotton fiber development, suggesting that Gh14-3-3L may be involved in regulating fiber elongation. Our data showed that, although 14-3-3L is expressed preferentially in the early development stages of cotton fibers in both TM-1 and Hai7124, the relative expression values were significantly different. The expression of 14-3-3L was significantly higher in most stages in Hai7124 (Additional file 6: Supplemental Figure S3E), than in TM-1. Furthermore in the primary-secondary transition period of fiber development, a secondary expression peak of 14-3-3L was also detected, indicating 14-3-3L may also be involved in thickening of the secondary cell wall, a function that has been reported in other species [61]. At 23 DPA, the expression of 14-3-3L decreased in TM-1, but was still increasing in Hai7124, suggesting a longer transcriptional period.
CelA1 was solely expressed during secondary wall cellulose synthesis, with no significant expression difference between TM-1 and Hai7124 (Additional file 6: Supplemental Figure S3C). The importance of CelA genes in secondary wall cellulose synthesis in developing cotton fibers has been reported by other researchers [16]. Expression of CelA3 at 17 DPA was significantly higher in TM-1 than in Hai7124. However, strength dropped much faster in TM-1 (e.g., at 20 DPA, Additional file 6: Supplemental Figure S3E) than in Hai7124 after that, which more stable expression was detected over a longer period of time in Hai7124.
Rac genes (coding for GTP-binding Rac proteins) are involved in various physiological processes, including cell polar growth [62][63][64][65], synthesis of the secondary wall [66][67][68], resistance response [69,70] and signal transduction [71][72][73]. In the present study, RacAs were predominantly expressed at the fiber elongation and secondary wall synthesis stages (Additional file 6: Supplemental Figure S3E), indicating that the role of RacA may be similar to that of Rac13, which is involved in the signal transduction pathway for cytoskeleton organization [66]. Similarly, RacB was predominantly expressed at the secondary wall synthesis stage in TM-1, and preferentially expressed at the initiation of fiber development and secondary wall synthesis in Hai7124 (Additional file 6: Supplemental Figure S3D), suggesting multiple roles for RacB. Because RacA and RacB expressions were consistently higher at most stages in Hai7124 relative to TM-1, the two genes may be contributing to fiber quality diversity.
Sucrose synthase (Sus) affects initiation and elongation of the single-celled fibers, and also has a major role in cell wall cellulose synthesis [9,[74][75][76][77]. Cotton fiber cell initiation and elongation are very sensitive to changes in Sus activity [9]. The expression peak of cotton sucrose synthase genes transcripts was earlier in wild cotton than in TM-1 [60]. Our study showed that greater Sus1 expression in Hai7124 than in TM-1 at initiation (0 and 1 DPA) and 8 DPA and at 3 and 5 DPA was reversed. At 20 DPA, Sus1 expression was significantly higher in TM-1, possibly as a result of earlier termination of fiber elongation and earlier initiation of cell wall synthesis in TM-1.
Recent studies of comparative gene expression profiling of isolated cotton fibers [60,78] have identified genes controlling differences in fiber growth between wild and domesticated cotton [e.g., genes encoding tubulin isoforms, endotransglycosylase/hydrolases, cytochrome P(450) monooxygenase, and antioxidant enzymes]. These studies will provide new clues in the future studies to better understand the interspecific divergence of fiber development in the two cultivated tetraploid cotton species, G. barbadense and G. hirsutum.

Putative role of D-subgenome in interspecific divergence
Most genes (62.5% in TM-1 and 62.5% in Hai7124) exhibited a higher evolutionary rate in the D-subgenome than in the A-subgenome, indicating that the D-subgenome of the allotetraploid evolved faster than the A-subgenome [59,[79][80][81][82]. Considering the weight of evidence provided by genome-specific expression of homeologs, chromosome location of genes and QTL distribution for fiber qualities (Table 3), we also determined the subgenome transcriptome contribution in interspecific divergence between TM-1 and Hai7124. In the integration interval of genes and QTLs, the gene transcripts were significantly A-subgenome biased, and the QTLs associated with fiber qualities were detected in corresponding A-subgenome chromosome bins. Similarly, the gene transcripts were significantly D-subgenome biased and relative fiber QTLs were detected in corresponding D-subgenome chromosome bins. While the transcripts of CAP was significantly A-subgenome biased in TM-1, but equally expressed in Hai7124 in the fiber developmental period, implying an important D-subgenome contribution for elite fiber quality traits in Hai7124. Integration of CAP duplicated loci with fiber length and strength QTLs in D13, indicates domestication and artificial selection of tetraploid plants for superior fiber quality resulted in more significant evolutionary dynamics in the D-subgenome than in the A-subgenome.

Conclusion
The study provided us the systematic report on the interspecific divergence for fiber development between Gossypium barbadense and G. hirsutum by analyzing structures, molecular evolution and transcripts levels of fiber development-related genes. The results indicated that selective breeding for certain desirable fiber qualities may have played an important role in divergence of the two cultivated tetraploid cotton species.

Plant materials
Of the 17 genes selected for this study (Table 1) were described in previous reports as being expressed in developing cotton fibers; their sequence information was deposited in GenBank. The orthologous loci of each of the 17 genes were isolated from two allotetraploid cottons, G. hirsutum acc. TM-1, G. barbadense cv. Hai7124, and living models of their two ancestral genome donors, G. herbaceum var. africanum (A-genome) and G. raimondii Ulbrich (D-genome). For expression analysis, developing cotton ovules and attached fibers were harvested from three replicates of TM-1 and Hai7124 on 0, 1, 3, 5, 8, 10, 14, 17, 20 and 23 days postanthesis (DPA). The mapping population was comprised of 138 BC 1 individuals generated from the cross [(TM-1×Hai7124) ×TM-1] [19].

PCR amplification, cloning and sequencing
Based on GenBank deposited sequences, gene-specific PCR primer pairs were individually designed for PCRamplification of full-length sequences of the 17 genes. For the smaller genes, we designed one primer pair; at least two nested primer pairs were designed for the larger genes (Additional file 1: Supplemental Table S1).
Standard polymerase chain reaction (PCR) was completed using High-fidelity ExTaq DNA polymerase (TaKaRa Biotechnology (Dalian) Co., Ltd., China). The PCR products were cloned into pMD18-T Vector (TaKaRa) according to the manufacturer's instructions, and sequenced from plasmid DNA templates. In order to obtain the sequence from both the A-subgenome and D-subgenome, at least 10 clones for each gene from each of the tetraploid species, TM-1 and Hai7124, were picked randomly and sequenced. Homeologs were identified by comparison to sequences from their diploid progenitors. To avoid possible complications originating from PCR recombination in allopolyploid cotton [83], a minimum of three clones was used to determine the gene sequence in each duplicated copy. In cases where PCR recombination occurred, 10 additional clones were sequenced to verify the corresponding correct sequence [83]. When the copy number of genes was difficult to identify by sequencing, it was confirmed by southern blot analysis.
Lengths were determined using alignment analysis, which was applied to each relevant sequence of applicable gene. For genes with one primer pair, the A-and D-subgenome sequences were clustered using the Neighbor-Joining method in MEGA3.1 http://www. megasoftware.net/. For those genes having more than one primer pair, sequence segments from the same genomes and subgenomes were contiged by CAP3 software.

Gene structure and phylogenetic analyses
Six DNA sequences for each gene, including homeologs and orthologs, isolated from the four species (A-, D-genome and the A and D subgenome from the allotetraploid), were aligned by ClustalX http://www.ebi.ac.uk/ Tools/msa/clustalw2/. Based on the sequences of genomic DNA and the cDNA, the structure of each gene was illustrated by the Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn/chinese.php). Phylogenetic analyses of the orthologs for each gene were performed using the Neighbor-Joining method in MEGA3.1 http://www.megasoftware.net/.

Chromosomal location
In order to locate all 17 pairs of homeologs on our genetic map, PCR primers were designed according to the sequence size difference between the two allotetraploid cottons (Additional file 4: Supplemental Table S2). For homeologs with no size polymorphism, singlenucleotide amplified polymorphism (SNAP) primers were designed based on a SNP between the TM-1 and Hai7124 sequences at a putative locus (Additional file 5: Supplemental Table S3). SNAP is a modified allele-specific amplified method, it introduce an additional base pair change in the primer to increase the specificity of the primer [84]. All SNAP primers was designed following Drenkard et al. [84] and using SNAPER http://ausubellab.mgh.harvard.edu.
All the primer pairs exhibiting polymorphisms between TM-1 and Hai7124 were used to survey 138 individuals of the BC 1 mapping population. The polymorphic loci were integrated in our backbone map [22] using Joinmap 3.0 software [85].

Evolution and ratios of sequence forms
For each gene, synonymous substitution rates (Ks), nonsynonymous substitution rates (Ka) and their ratios (Ka/Ks) in various comparisons among the four species (A 1 VS D 5 , A 1 VS At, D 5 VS Dt, At VS Dt) were calculated for the coding region using KaKs_Calculator http://sourceforge.net/projects/kakscalculator2/. Ks is the rate of evolution and the Ka/Ks ratio indicates the selective force acting on the protein [86]. A Ka/Ks ratio of 1 represents neutral evolution, that is, the number of nonsynonymous changes at each possible nonsynonymous site is the same as the number of synonymous changes per synonymous site. Ka/Ks < 1 indicates a purifying selection (selection generally eliminates deleterious mutations and maintains protein status quo); Ka/Ks > 1 indicates positive selection, where selection may be modifying the protein. The higher the ratio is, the stronger the evidence that selection is occurring.

Development of homeolog-specific PCR primer pairs
Homeolog-specific primers were designed and applied based on sequence differences between duplicated loci (Additional file 7: Supplemental Table S4). Specificity of those primers was detected by PCR amplification of genomic DNA from G. herbaceum (A-genome), G. raimondii (D-genome), TM-1 and Hai7124. Primers specific for the A-genome resulted in amplifications from G. herbaceum, TM-1 and Hai7124, but not from G. raimondii, which indicated that the amplification products of TM-1 and Hai7124 were from their A-subgenomes. Conversely, the D-genome specific primers only amplified from G. raimondii and D-subgenome of the allopolyploid.

RNA isolation and Q-PCR analysis
Total RNA was extracted from young fibers of G. hirsutum acc. TM-1 and G. barbadense cv. Hai7124 at 0, 1, 3, 5, 8, 10, 14, 17, 20 and 23 DPA using the CTAB-acid phenol extraction method [87]. RNA samples were treated with DNase I (Ambion, Austin, TX, USA) according to manufacturer's instructions to remove trace contaminants of genomic DNA. Total RNA samples (1 μg per reaction) were reversely transcribed into cDNAs by avian myeloblastosis virus (AMV) reverse transcriptase, and the cDNAs were used as templates in Q-PCR reactions.
For each gene, the homeolog-specific primers were designed based on the SNP between A-and D-subgenome homeologs (Additional file 7: Supplemental Table  S4). A cotton elongation factor (EF1a) gene was used as a standard control in the Q-PCR reactions (Additional file 7: Supplemental Table S4). The Q-PCR reaction and calculated relative value for expression level of each homeolog were determined as described by Wang et al. [8]. The information of the Q-PCR analysis based on the MIQE checklist was submitted in Additional file 8.