Translational Regulation of Duplicated Gene Expression Evolution in Allopolyploid Cotton

1. Introduction

Polyploidy, the whole genome duplication (WGD) event, has been a fundamental and recurrent force in the evolutionary history of angiosperms, profoundly influencing genome structure, function, and organismal complexity [1,2]. By doubling or multiplying genome content, polyploidy provides a reservoir for genetic innovation, leading to phenotypic diversification and adaptation [3,4]. One key realization in the past decade is the profound impact of polyploidy on transcriptomic landscapes [5,6,7]. Studies have reported substantial transcriptomic rewiring in polyploid species, encompassing biased homoeolog expression, condition-specific gene regulation, and transgressive expression levels. These transcriptional changes often underlie phenotypic and adaptive variations upon polyploidy formation. However, despite extensive research on transcriptome and proteome [8], the consequences of polyploidy at the translational level remain largely unexplored.

Ribosome profiling, known as Ribo-seq, is a high-throughput sequencing technique that provides a comprehensive snapshot of active translation within a cell by mapping the positions of ribosomes on mRNA molecules [9]. This approach involves the isolation and deep sequencing of mono-ribosome-protected mRNA fragments, often referred to as ribosome footprints. By analyzing the distribution and abundance of these footprints, researchers can quantitatively assess translation rates, identify translation start sites, and uncover the full repertoire of translated open reading frames (ORFs) [10]. Beyond its utility in identifying translated regions, Ribo-seq enables the study of translational regulation, including the impact of cis-regulatory elements, ribosome pausing, proteins, and environmental cues on protein synthesis [11,12,13]. Moreover, by comparing Ribo-seq data with transcriptomic data, translational efficiency (TE), a metric that reflects the coupling of transcription and translation could be measured [14,15].

The application of Ribo-seq has rapidly expanded across various plants systems. In Arabidopsis, for example, global translatome profiling was performed in conjunction with RNA sequencing, revealing tightly regulated translation with poor correlation with transcriptional changes during immune induction, leading to the discovery of novel regulators of the immune response [16]. Similarly, translational reprogramming was demonstrated a crucial layer of stress responses in crop species including rice (Oryza sativa) [17,18], potato (Solanum tuberosum) [14], citrus (Citrus reticulata) [19], and maize (Zea mays) [15]. Beyond its role in stress responses, Ribo-seq has been instrumental in investigating allele-specific translation in hybrids, revealing the intricate interplay between genetic variation and translational regulation [13,20]. These studies have highlighted the potential of translational regulation to contribute to heterosis, where hybrid offspring exhibit superior traits compared to their parents. While some progress has been made in understanding the impact of polyploidy on translational regulation in plants, as exemplified by studies in soybean (Glycine max) [21], the role of translational control on gene expression evolution accompanying allopolyploidization remains largely underexplored.

Cotton, a globally significant fiber crop, represents a model system for studying allopolyploidization [22]. The allopolyploid cotton genome was formed approximately 1-2 million years ago, through hybridization and subsequent genome doubling of A- and D-genome diploid progenitors. Previous research on cotton gene expression evolution accompanying allopolyploidization have primarily focused on transcriptomics [23,24] and proteomics [25,26], with limited insights from translational regulation. To address this gap and elucidate the impact of allopolyploidization on translational regulation, we conducted a comprehensive analysis of gene expression at both the transcriptional and translational levels in allopolyploid cotton G. hirsutum (AD₁), and the representatives of its progenitor diploid species G. arboreum (A₂) and G. raimondii (D₅). To expand our comparative analysis, we incorporated public RNA-seq and Ribo-seq data from another allopolyploid species G. barbadense (AD₂) [27] for comparison. Initially, we compared the correlation and expression variances between transcriptomic and translational datasets, revealing a dominant correlation and larger expression variances in transcriptomes, demonstrating greater differential ranges in transcriptomic data. Subsequently, we compared differential gene expression patterns between transcript and translation levels, finding nearly no genes with completely opposite regulatory trends between these two levels, with only 5.46% showing consistent regulatory trends, primarily enriched in terpene functional genes associated with plant stress resistance. We then compared TE and found that the TE of genes showed a significant negative correlation with transcript levels and that genes with high TE were enriched for specific functions, which was also affected by allopolyploidy. Finally, comparison of expression between homologous A and D subgenomes, allopolyploidy was found to reduce homoeologs expression bias genes at the translational level, and a significant D subgenomic expression dominance at the translational level was detected in AD₂. In conclusion, the results of this study provide insights into the effects of allopolyploidy on gene expression and translation in cotton.

2. Materials and Methods

2.1. Plant Materials

Three Gossypium species were used in the present study: one allopolyploid species Gossypium hirsutum cv. TM-1 (AD₁) and two diploid species representing the model diploid progenitors, Gossypium arboreum cv. Zhongya 1 (A₂) and Gossypium raimondii PI 530899 (D₅). Seeds were germinated in a growth chamber in darkness under controlled conditions (28±1°C and 60±5% humidity). After germination, seedlings were cultivated with a 16-hour light/8-hour dark photoperiod cycle under the same controlled condition for six weeks. Leaf tissue from the fourth to fifth true leaves were collected from three biological replicates (5 plants per replicate) of each genotype, immediately frozen in liquid nitrogen, and stored at -80°C.

2.2. Ribosome Profiling (Ribo-Seq)

Ribosome-protected fragments (RPFs) were isolated using a modified protocol [28]. Briefly, 1 gram of leaf tissue was ground in liquid nitrogen with 10% [w/w] polyvinylpolypyrrolidone poly(vinylpolypyrrolidone) in liquid nitrogen and resuspended in to 5ml of lysis buffer (150 mM Tris-HCl pH 8.0, 40 mM KCl, 20 mM MgCl₂, 2% [v/v] polyoxyethylene, 0.4% [w/v] sodium deoxycholate, 1.5mM dithiothreitol, 50 μg/ml chloramphenicol, 50 μg/ml cycloheximide, 2% [v/v] polyvinylpyrrolidone). Adding 25 μL DNase I to a final concentration of 10 units/mL, the mixture was incubated on ice for 30 minutes. After centrifugation at 4°C for 10 minutes at 8000×g, the supernatant was transferred to a new centrifuge tube and further centrifuged at 4°C for 10 minutes at 16000×g. The supernatant was layered onto an 8 ml cushion (150 mM Tris-HCl pH 8.0, 40 mM KCl, 20 mM MgCl₂, 5 mM EGTA, 1.5 mM dithiothreitol, 50 μg/ml chloramphenicol, 50 μg/ml cycloheximide and 1.75 M sucrose) and subjected to ultracentrifugation at 4°C, 257,000×g for 3 hours in a Beckman ultracentrifuge (Optima XE-100, United States) using a Type 70 Ti Rotor. The pellet was washed twice with pre-cooled ddH2O and resuspended in 300 μl RNase I digestion buffer (10 mM Tris-HCl pH 8.0, 100 mM NaCl, 10 mM EDTA, 5 mM EGTA,1% [w/v] SDS). The resuspended mixture was transferred to a new tube and centrifuge at 1,000×g at 4 °C for 30s. Six units of RNase I was added to the supernatant and incubated at room temperature for 1.5 hours. For rRNA depletion from RPFs, the Ribo-pool rRNA Depletion Kit (siTOOLs p-k024-000031) was used. The obtained RNA fragments were recovered using the RNA Clean & Concentrator-5 kit (Zymo R1014) and purified using the 15% denatured PAGE gel (7M Urea, 0.5×TBE, 15% [v/v] Acrylamide-Bisacrylamide, 0.08% [v/v] Ammonium persulfate substitute, and 0.05% [v/v] TEMED) targeting 27-34 nt in size, according to the protocol [29].

To construct Ribo-seq libraries from purified and size selected RPFs, the phosphate group of RPFs was treated with T4 PNK, followed by 3′ adaptor (5′-rAppGATCGGAAGAGCACACGTCT-NH2) ligation using truncated T4 RNA ligase 2 (NEB). The RNA strand with added adapters was denatured at 65°C along with the designed RNA primers. Subsequently, reverse transcription was carried out using SuperScript II (ThermoFisherScientific, Waltham, MA, USA) and a RT primer (5′-GATCGTCGGACTGTAGAACTCTGAACGTGTAGATCTCGGTGGTCGCCGTATCATT/iSp18/CACTCA/iSp18/CAGACGTGTGCTCTTCCGATCT) at 42°C for 1.5 hours. RPFs were circularized with Circligase ssDNA ligase (Epicentre CL4111K) and amplified by PCR. Finally, PCR amplification to generate a 150 bp cDNA library was performed using E×Taq (14 cycles, 60 °C annealing, and primer sequences are listed in Table S1), Which was sequenced using Illumina NovaSeq 6000 by Pasenno Biotechnology Co. (Shanghai, China).

2.3. RNA Sequencing (RNA-Seq)

Total RNA was extracted using the TIANGEN Total RNA extraction kit (TIANGEN X0817). RNA-seq libraries were prepared using the Illumina TruSeq Stranded mRNA-seq library preparation kit (refer to online protocol accessed on March 26, 2023) for paired-end sequencing (150 bp×2) on an Illumina HiSeq 2500 platform (Novogene Technology Co. Beijing, China).

2.4. Ribo-Seq and RNA-Seq Data Processing

In addition to data generated from A₂, D₅, and AD₁, public cotton leaf RNA-seq and Ribo-seq data for another allopolyploid species, G. barbadense (AD₂), were downloaded from the National Center for Biotechnology Information (NCBI) SRA depository (details in Supplemental Table S2) to enable comparisons across different polyploid genotypes. Public cotton ribosomal RNA (rRNA) sequences were retrieved from NCBI and subjected to de-redundancy using CD-HIT (v4.8.1) (https://github.com/weizhongli/cdhit) to generate a reference set for rRNA filtering.

Quality control, adapter trimming, and low-quality read filtering were performed using FASTP (0.23.4) [30]. Bowtie2 (v2.3.5.1) [31] was used to remove rRNA reads, utilizing the non-redundant cotton rRNA sequence set as reference. To facilitate direct comparison between allopolyploid and diploid datasets, the number of reads in the matched cleaned data from A₂ and D₅ diploids was adjusted to equal amounts and then combined to generate a “Combined Diploid” dataset (A₂D₅). The processed data from AD₁, AD₂, and A₂D₅ were then subjected to further analysis.

Processed reads from RNA-seq were aligned to the G. hirsutum cv. TM-1 UTX (v2.0) reference genome [32]. Similarly, using Hisat2 (v2.2.1) [33] with default parameters (details provided). Read counts and transcripts per million (TPM) values were quantified for each gene using the featureCounts function from the Rsubread package (v2.16.1) [34].

Processed reads from Ribo-seq were aligned to the same reference genome using STAR (v2.7.9a) [35] with parameters ‘--outFilterMismatchNmax 1 --outFilterMultimapNmax 1′ to retain only uniquely mapped reads with at most one mismatch. The R function featureCounts was used to obtain read counts and TPM values. BAM files from STAR alignment were subjected to determine read length distribution, and three-nucleotide (3-nt) periodicity using RiboCode (v1.2.15) [36].

2.5. Analysis of Gene Expression at Transcriptional and Translational Levels

Differential gene expression analysis was performed using the DESeq2 package in R (v1.44.0) [37] to identify genes with significant expression changes between allopolyploids (AD₁ and AD₂) and the combined diploid dataset (A₂D₅). Genes with adjusted p-values < 0.05 by Benjamini-Hochberg correction [38] were considered differentially expressed. Pairwise comparisons were conducted between AD₁ vs. A₂D₅ and AD₂ vs. A₂D₅ at both transcriptional and translational levels. To investigate the coordinated regulation of gene expression, genes were categorized into eight groups by cross tabulating their differential expression patterns (i.e., increased, decreased, and unchanged relative to A₂D₅) between two levels. Except for the category with no change in expression at both levels, functional enrichment analysis of Gene Ontology (GO) terms was performed for genes in each of the other eight groups using R-clusterProfiler (v4.8.3) [39], with p-values adjusted for multiple testing using the Benjamini-Hochberg method. A threshold of adjusted p-value < 0.05 was applied to identify significantly enriched GO terms. Functional protein annotation is performed by Revigo was employed to summarize and visualize significant GO terms (v 3.4.4) [40].

2.6. Calculation of Translational Efficiency

TE was calculated for each gene as the log2 ratio of RPF to mRNA expression levels: TE = log2(RPF / mRNA). To account for potential library size differences, both RPF and mRNA expression data were normalized using the median of ratios normalization method built-in DESeq2. Pairwise comparisons between matched Ribo-seq and RNA-seq samples for each genotype were conducted using DESeq2 to obtain log2FoldChange values as TE estimates. Genes with adjusted p-value < 0.05 were considered to have significantly different TE. A positive TE value indicates higher translation efficiency relative to transcription, while a negative value suggests lower translation efficiency. Genes exhibiting extreme TE changes (|TE| > 6) were subjected to Gene Ontology (GO) enrichment analysis using R-clusterProfiler (v4.8.3) [39]. The R package pheatmap (v1.0.12) was utilized for visualization [41].

2.7. Analysis of Homoeolog Expression Bias (HEB)

HEB was defined as the unequal expression of homoeologous gene copies in the allopolyploid [5]. In allopolyploid cottons (AD genome), the direction of bias was determined by the subgenome (At or Dt) exhibiting higher expression. The At and Dt homoeologous relationships within the allopolyploid G. hirsutum reference genome was previously inferred [42].

To assess HEB at the transcriptional and translational levels, expression levels (mRNA and RPF) were compared between At and Dt homoeologs using DESeq2. Genes with adjusted p-values < 0.05 (Benjamini-Hochberg correction) were considered to exhibit significant HEB. The magnitude of HEB was quantified as the log2fold-change in expression between homoeologs. At the genomic scale, the prevalence of bias towards a specific subgenome was evaluated to assess genome-wide HEB.

To assess HEB at the translational level, TE values were compared between At and Dt homoeologs using Student’s t-test. Significant differences in TE between homoeologs were determined using a Fisher’s combined p-value < 0.05.

3. Results

3.1. Transcriptional and Translational Profiling of Diploid and Allopolyploid Cottons

A total of 183 million, 130 million, and 226 million Ribo-seq reads were generated for G. arboreum (A₂), G. raimondii (D₅), and G. hirsutum (AD₁) seedling leaves, respectively (Table S2). After removal of low-quality reads and rRNA, an average of 18 million clean reads remained per sample, representing 30.8% of the total. The average unique mapping rate ranged from11.2% to 39.6%, tetraploid (AD₁) having the lowest rate most likely due to its higher level of sequence homology compared to diploids A₂ and D₅. To facilitate direct comparisons between parental diploids and allopolyploids, we combined the clean reads of A2 and D5 into a ‘Combined Diploid’ dataset (A₂D₅) and used it for subsequent analyses, which exhibited 14.5%-16.5% unique mapping rates, comparable to those of AD₁.

To validate Ribo-seq data quality, we examined RPF (Ribosome Protected Fragment) sizes, which ranged from 26 to 32 nucleotides (Figure S1), consistent with known cotton Ribo-seq data [27]. The RPFs exhibited a periodicity of 3 nucleotides (Figure 1A), a typical characteristics of ribosome profiling [43]. Pairwise Pearson correlation coefficients were calculated between all samples to examine data reproducibility between replicated experiments, including a previously published dataset from G. barbadense (AD₂) leaves [27]. One AD₁ sample (AD₁-2) with relatively low correlation coefficients (R=0.85-0.86) (Figure 1B) and poor clustering in principal component analysis (PCA) (Figure 1C) with other replicates, was excluded from subsequent analyses.

For RNA-seq, a total of 131 million, 143 million, and 155 million read pairs were obtained from A₂, D₅, and AD₁ cotton leaf samples, respectively (Table S3). The clean read pairs obtained per sample ranged 39 to 57 million, with an average unique mapping read rate of 87.4% (Table S3). Diploid data from A₂ and D₅ were combined into the A₂D₅ dataset as for the Ribo-seq data. Pearson correlation coefficients ranging from 0.79 to 0.99 were observed between transcriptome and translatome data (Figure 1B). Notably, expression variances were higher in transcriptomes than in translatomes (Figure 1D; variance of log10(TPM): 0.914 vs 0.686), potentially indicating a wider expression range in mRNA levels than translated proteins.

3.2. Concordantly Decreased Transcription and Translation Accompanying Allopolyploidization

To study gene expression changes during allopolyploidization, we conducted differential gene expression analyses between diploid and tetraploid cotton using RNA-seq and Ribo-seq data (criterion: |log2fold change| ≥ 1 and FDR < 0.05). We identified significant differences as follows (Figure 2A): comparing AD₁ to A₂D₅, 11,707 (15.6%, 5,711 up-regulated and 5,996 down-regulated) transcriptional changes and 1,622 (2.1%, 472 up-regulated and 1150 down-regulated) translational changes were identified. Comparing AD2 to A2D5, more DEGs were identified at both levels, with 21,506 (28.7%, 10,773 up-regulated and 10,733 down-regulated) transcriptional changes and 6,803 (9.1%, 2,293 genes up-regulated and 4510 down-regulated) translational changes. Between AD₁ and AD₂, 14,370 (19.2%, 6919 up-regulated and 7451 down-regulated) transcriptional changes and 2148 (2.9%, 1259 up-regulated and 889 down-regulated) translational changes were identified. Concordantly, fewer DEGs were identified at the translational level by Ribo-seq. Regarding the four inter-ploidy contrasts (i.e., AD₁ vs A₂D₅ and AD₂ vs A₂D₅ at either transcriptional or translational levels), we identified a total of 7,686 genes exhibiting allopolyploidy changes common to both tetraploids AD1 and AD2 (Figure 2B, numbers in red). Of these, 189 genes showed common changes only at the translational level, 7,078 genes showed changes only at the transcriptional level, and 419 genes showed changes at both levels. Considering only those exhibiting conserved up- or down-regulated changes relative to A₂D₅, allopolyploid-specific regulation was inferred for 7,393 genes, which were further categorized into eight groups (A-H) based on their expression changes in Ribo-seq and RNA-seq (each being increased, decreased, or no changes) (Figure 2C).

We specifically focus on the categories exhibiting allopolyploid-specific regulation at the translational level, including categories A, B, C, F, G, and H (582 genes, 7.9%). Among these, category C (53, 0.7%) and F (352, 4.6%) exhibited concordant regulatory directions between transcription and translation in allopolyploids relative diploids, with a prominent bias toward decreased expression (more F than C genes). Category F was enriched for GO terms associated with terpenoid synthesis and metabolism, indicating a negative impact of allopolyploidization on terpenoids at both transcriptional and translational levels. Categories C was not enriched for any GO terms.

Few genes were found exhibiting opposite regulatory directions between transcription and translation, as classified to category A (0 genes) or H (2 genes), nor enriched for GO terms.

Category B (21 genes, 0.3%) and G (155 genes, 2%) include genes with unchanged transcription levels but higher or lower translation in allopolyploids, respectively. Only category G was enriched for GO terms related to ribosome function and protein import into chloroplasts.

3.3. Impact of Allopolyploidy on Translational Efficiency

To isolate the effect of translational change independent of transcriptional changes, we calculated the translation efficiency (TE) values for genes, defined as the ratio of translational to transcriptional expression levels (i.e., log2(RPF/mRNA)). A gene with TE = 0 is considered to have consistent transcriptional and translational products; TE ≠ 0 indicates regulation at the translational level, with TE > 0 suggesting higher and TE < 0 suggesting lower translation relative to transcription.

We detected a total of 65,001 transcriptionally expressed genes (RNA-seq read counts above 1 in at least one biological replicate) in A₂D₅, with 9,810 (15.1%) having TE < 0 and 7,274 (11.2%) having TE > 0. In AD₁, there were 62,343 transcriptionally expressed genes, with 8,664 (13.9%) having TE < 0 and 5,845 (9.4%) having TE < 0 and. In AD₂, there were 63,871 transcriptionally expressed genes, with 9,231 (14.5%) having TE < 0 and 4,645 (7.3%) having TE > 0 (Figure 3A). These results indicate that at least one-fifth of genes exhibit translational level regulation, with a general trend of more down-regulation than up-regulation of TE. Interestingly, the magnitude of translational level regulation was slightly reduced in allopolyploids compared to diploids.

To explore the relationship between TE and mRNA levels, we analyzed their correlations in each genotype. TE was significantly negatively correlated with transcriptional expression in A₂D₅, AD₁, and AD₂, with Pearson’s R of -0.2, -0.2, and -0.4, respectively (Figure 3B). This suggested that, similar to observations in maize and potato [14,15], the mRNA translation efficiency and ultimately protein synthesis capacities are limited when mRNA levels are high.

We further analyzed the genes with the highest TE (> 6) and the lowest TE (< -6) in in A₂D₅ (1,062 and 801 genes), AD₁ (753 and 978) and AD₂ (808 and 1,936) (Figure 3C). GO enrichment analysis revealed that genes with TE > 6 in A₂D₅ were primarily enriched for energy metabolism-related terms, including ATPase activity, proton-transporting ATP synthase complex, ribonucleoside triphosphate metabolic process, and photosynthesis (Figure 3D). Similar GO terms were enriched in AD₁ and AD₂, with additional enrichment for RNA polymerase activity-related terms in AD₂. This suggests that energy metabolism-related genes tend to be translated with high efficiency to increase protein production, a trend consistent before and after allopolyploidy.

For genes with TE < -6 (Figure 3D), A₂D₅ showed enriched GO terms related to cellular component biogenesis and ciliary were predominantly. In contrast, both AD₁ and AD₂ were predominantly enriched for stomatal opening, cyclin-dependent protein, and GO terms related to serine/threonine protein kinase. This suggests that allopolyploidy has a profound impact on genes with low TE.

3.4. Reduced Amount of Homoeolog Expression Biases at the Translational Level

To investigate the subgenomic contribution to gene expression, 22,889 ortho-homoeolog groups (OGs) were previously inferred between the At and Dt homoeologs in the AD₁ reference genome. At the transcriptional level, we observed 23%-37% of OGs exhibiting homeolog expression bias (HEB): A₂D₅ had 5,285 (23%, 2,574 A-bias vs 2,711 D-bias), AD₁ had 5,237 (22.9%, 2,576 vs 2,661), and AD₂ had 8,407 (36.7%, 4,165 vs 4,242), with no significant imbalance towards either the A- or D-subgenome. At the translational level, the number of OGs exhibiting HEB decreased to 3%-9%: A₂D₅ had 2,224 (9.7%, 1,115 vs 1,109), AD₁ had 1,013 (4.4%, 484 vs 529), and AD₂ had 562 (2.5%, 244 vs 318), with only AD₂ showing significantly more D-biases of translational expression.

When comparing TE, there were 1,076 (4.7%, 391 vs 685), 625 (2.7%, 320 vs 305), and 547 (2.4%, 237 vs 310) significant differences between At and Dt in A₂D₅, AD₁, and AD₂ respectively. Both A₂D₅ and AD₂ showed significant D -biased translational level regulation.

These results indicated that allopolyploidization may have reduced the subgenomic differences at the translational level, leading to fewer biases at this level. Additionally, a D-genome dominant translational regulation was observed in AD₂, may reflecting species-specific translational regulation with respect to subgenomic contribution.

Table 1. Homoeolog expression bias (HEB) genes on transcriptional and translational levels.

		A2D5		AD1		AD2
		Number	Proportion	Number	Proportion	Number	Proportion
RNA	Total	5285	23.09%	5237	22.88%	8407	36.73%
	A>D	2574	11.25%	2576	11.25%	4165	18.20%
	A<D	2711	11.84%	2661	11.63%	4242	18.53%
	p-value*	0.0595		0.2402		0.401
Ribo	Total	2224	9.72%	1013	4.43%	562	2.46%
	A>D	1115	4.87%	484	2.11%	244	1.07%
	A<D	1109	4.85%	529	2.31%	318	1.39%
	p-value*	0.8988		0.1574		0.001799
TE	Total	1076	4.70%	625	2.73%	547	2.39%
	A>D	391	1.71%	320	1.40%	237	1.04%
	A<D	685	2.99%	305	1.33%	310	1.35%
	p-value*	2.20E-16		0.5485		0.001801

*Chi-square test was performed to compare the numbers of A>D and A<D.

Table 1. Homoeolog expression bias (HEB) genes on transcriptional and translational levels.

		A2D5		AD1		AD2
		Number	Proportion	Number	Proportion	Number	Proportion
RNA	Total	5285	23.09%	5237	22.88%	8407	36.73%
	A>D	2574	11.25%	2576	11.25%	4165	18.20%
	A<D	2711	11.84%	2661	11.63%	4242	18.53%
	p-value*	0.0595		0.2402		0.401
Ribo	Total	2224	9.72%	1013	4.43%	562	2.46%
	A>D	1115	4.87%	484	2.11%	244	1.07%
	A<D	1109	4.85%	529	2.31%	318	1.39%
	p-value*	0.8988		0.1574		0.001799
TE	Total	1076	4.70%	625	2.73%	547	2.39%
	A>D	391	1.71%	320	1.40%	237	1.04%
	A<D	685	2.99%	305	1.33%	310	1.35%
	p-value*	2.20E-16		0.5485		0.001801

*Chi-square test was performed to compare the numbers of A>D and A<D.

4. Discussion

This project primarily conducted transcriptome and translatome sequencing on an allotetraploid upland cotton AD₁ and its two diploid progenitors, A₂ and D₅, aiming to investigate translational changes associated with cotton polyploidization.

4.1. Ribo-Seq Technical Difficulties

Ribo-seq technology involves the degradation of non-ribosome-protected RNA fragments using RNase I, allowing the recovery of ribosome-protected mRNA fragments typically sized between 27 to 34 nucleotides [43]. These fragments, known as ribosome footprints (RFs), are used to infer various aspects of translation such as translation start site identification, codon usage bias, upstream open reading frames, and translation efficiency [44]. However, Ribo-seq has its limitations. The enzymatic cleavage conditions are challenging to control, and constructing libraries can be difficult due to the small size of the target fragments. Moreover, these small fragments may map to multiple locations, leading to potential false positives [45]. Additionally, the presence of polysaccharides and polyphenols in cotton can interfere with nucleic acid extraction. To address these issues, we first used PVPP and PVP-40 to adsorb excess polysaccharides and polyphenols from cotton materials, facilitating nucleic acid extraction. Subsequently, we employed polysome profiling to optimize enzyme digestion conditions for the recovery of single ribosomes, thereby increasing the yield of initial samples to ensure sufficient mRNA fragment libraries. To mitigate false positives in subsequent analyses, we retained only uniquely mapped reads during mapping. Finally, to obtain ample data for analysis, we increased sequencing depth, ensuring adequate coverage for comprehensive translation analysis. In 2023, Ghulam Qanmber similarly utilized ribo-seq technology to analyze the translatome data of the allopolyploid island cotton G. barbadense (AD₂) [27]. Given that the diploid ancestors of G. barbadense (AD₂) are also G. arboreum (A₂) and G. raimondii (D₅), we analyzed the transcriptome and translatome data of its leaf samples together (Figure 1B and C).

4.2. Cotton Leaf Gene Express Patterns in Transcriptional and Translational Levels

According to our results, the number of genes expressed at the transcriptional level is more than that at the translational level, which is consistent with the conclusions from other translational studies [43]. This discrepancy arises because RPFs captured by translationomics represent only those mRNA segments occupying ribosomes. The abundance of detected RPFs reflects the biological state at that moment [46]. More than one-fifth of transcriptionally expressed genes show some degree of translational regulation, similar to observations in allopolyploids soybean leaves [21]. Based on this information, it can be hypothesized that the discrepancy between transcriptional and translational gene expression levels highlights a dynamic regulatory process in cellular biology. This suggests that the presence of mRNA does not necessarily equate to active translation, as indicated by the selective detection of RPFs in translationomics. Understanding the extent and mechanisms of translational regulation is crucial for deciphering how cells allocate resources and respond to environmental cues, which warrants further investigation into the temporal and spatial dynamics of protein synthesis in biological systems.

4.3. Expression Bias of Homologous Genes in Cotton

At the transcriptional level, more than 20% of homologous genes exhibiting expression bias in allopolyploid cotton leaves, consistent withthe results in other studies [47,48,49]. While D-biased genes slightly outnumber A-biased genes across three samples, the global genomic expression is unbiased [49]. At the protein level, despite examining protein expression only in cotton seeds compared to root tips, D-biased gene numbers exceed those of A-biased genes among detected proteins [25,26]. Correspondingly, we observed more D-biased genes than A-biased genes at the translational level in AD₂. Although the homozygosis transcripts from different subgenome exhibit various TE, the regulation mechanism is elusive. Some well-known translational regulation manners may give a clue to understand translation bias in allopolyploid species, such as codon usage, uORF, and secondary structures in the downstream of translation sites. Polyploidy species could increase protein abundance, however, some crucial biology processes need preciously control in dosage effect. Translational regulation as one of the gene expression control manners shows translational bias in different subgenome transcripts.

4.4. Translational Regulation that May Be Caused by Polyploidy

In studies of translational regulation of homologous genes during polyploid evolution, greater differences were found between diploids than post-polyploidization, where differences in AD₁ did not show significant bias. This phenomenon, noted in the limited studies on polyploid translational regulation [21], suggests a tendency for translational regulation to reduce disparities in gene expression between combined allopolyploid transcriptomes and their diploid progenitors. This adaptive response of the genome likely ensures basic expression balance as it adjusts over time to new gene copy numbers. Such expression adjustments may involve various mechanisms like DNA methylation [50], histone modifications [51], and transcription factor activities [52,53]. Therefore, despite significant expression differences between parental diploids, these disparities may diminish post-polyploidization, indicating a re-establishment of expression equilibrium to maintain cellular functionality and homeostasis [54].

In summary, polyploidization influences gene expression far beyond mere changes in gene copy numbers, encompassing intricate regulatory adjustments. These adaptive mechanisms likely provide essential flexibility and stability for organisms to thrive in dynamically changing environments. Future research will continue to explore how these mechanisms interact and contribute to shaping and maintaining species adaptability and survival in diverse ecological niches.