*2.3. Expression Analysis of GhTCP Genes in Different Tissues and under Various Stress Conditions*

To provide reliable information on the growth and developmental functions of TCP genes in *G. hirsutum*, their transcript accumulation patterns in mature leaves, stem, root, torus, petal, stamen, pistil, cylycle, ovules, and fibers of *G. hirsutum* was investigated (Figure 3). We obtained transcriptome

data from the NCBI Sequence Read Archive (accession number RJNA248163). Some TCP genes with close phylogenetic relationships showed similar or divergent expression patterns. For instance, the paralogous pair, *GhTCP2* and *GhTCP5*, was expressed highly in both the torus and petal, at moderate levels in the ovules at 20 days post anthesis (DPA), and at low levels in mature leaves. Most CYC/TB1-type genes were only weakly or not expressed in all tissues, suggesting that they were primarily expressed in other organs not tested or under special conditions. In contrast, *GhTCP31*, *GhTCP40*, and *GhTCP47* were constitutively expressed at very high levels in all tissues tested, indicating that these genes played regulatory roles during multiple development stages. Some TCP genes exhibited tissue-specific expression. *GhTCP31* and *GhTCP40* were highly expressed only in *Int. J. Mol. Sci.*  reproductive organs: Torus, petal, stamen, pistil, and calycles. **2018**, *19*, x FOR PEER REVIEW 6 of 14

**Figure 3.** Heat map representation of GhTCP gene expression in different tissues. The tissues used for expression profiling are indicated at the top. The genes are shown on the left of the expression bars and the phylogenetic relationship is shown. The −3 to 20 days post anthesis (DPA) indicate −3, −1, 0, 1, 3, 5, 10, and 20 days after pollination. **Figure 3.** Heat map representation of GhTCP gene expression in different tissues. The tissues used for expression profiling are indicated at the top. The genes are shown on the left of the expression bars and the phylogenetic relationship is shown. The −3 to 20 days post anthesis (DPA) indicate −3, −1, 0, 1, 3, 5, 10, and 20 days after pollination.

To predict possible functions of *TCP* genes in environmental adaptation, we investigated the transcriptional profile of *TCP* genes under various stress conditions, including heat, salt, and drought stresses. In total, 41 genes exhibited variations in expression (Figure 4). Of the three treatments, heat stress caused relatively more fluctuations in the transcript abundance of *TCPs* than did salt or drought stress. Under heat stress conditions, 18 *TCP* genes were downregulated and eight were upregulated. In response to salt treatment, the expression of five GhTCPs (*TCP7*, *14*, *25*, *33*, and *35*) increased instantly, and then decreased slowly during continued salt stress. Six GhTCP genes were selected at random for quantitative RT-PCR (qRT-PCR) analysis to determine the relative expression under salt and drought stresses (Figure 5). The qRT-PCR results indicated that these *GhTCP* genes showed similar expression patterns to the transcriptome sequencing results. To predict possible functions of *TCP* genes in environmental adaptation, we investigated the transcriptional profile of *TCP* genes under various stress conditions, including heat, salt, and drought stresses. In total, 41 genes exhibited variations in expression (Figure 4). Of the three treatments, heat stress caused relatively more fluctuations in the transcript abundance of *TCPs* than did salt or drought stress. Under heat stress conditions, 18 *TCP* genes were downregulated and eight were upregulated. In response to salt treatment, the expression of five GhTCPs (*TCP7*, *14*, *25*, *33*, and *35*) increased instantly, and then decreased slowly during continued salt stress. Six GhTCP genes were selected at random for quantitative RT-PCR (qRT-PCR) analysis to determine the relative expression under salt and drought stresses (Figure 5). The qRT-PCR results indicated that these *GhTCP* genes showed similar expression patterns to the transcriptome sequencing results.

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 7 of 14

**Figure 4.** Expression of GhTCP genes under heat, salt, and drought stresses. The genes are shown on the left of the expression bars and the phylogenetic relationship is shown. The abiotic stresses used for expression profiling are indicated at the top. The 1, 3, 6, and 12 h indicate hours after treatment. **Figure 4.** Expression of GhTCP genes under heat, salt, and drought stresses. The genes are shown on the left of the expression bars and the phylogenetic relationship is shown. The abiotic stresses used for expression profiling are indicated at the top. The 1, 3, 6, and 12 h indicate hours after treatment. the left of the expression bars and the phylogenetic relationship is shown. The abiotic stresses used for expression profiling are indicated at the top. The 1, 3, 6, and 12 h indicate hours after treatment.

**Figure 5.** Relative expression levels of six *GhTCP* genes under salt and drought stress. QRT-PCR **Figure 5.** Relative expression levels of six *GhTCP* genes under salt and drought stress. QRT-PCR analyses were performed using RNA generated from cotton leaves after NaCl and Polyethylene Glycol (PEG) treatment. Error bars represent standard error of the mean. **Figure 5.** Relative expression levels of six *GhTCP* genes under salt and drought stress. QRT-PCR analyses were performed using RNA generated from cotton leaves after NaCl and Polyethylene Glycol (PEG) treatment. Error bars represent standard error of the mean.

#### analyses were performed using RNA generated from cotton leaves after NaCl and Polyethylene Glycol (PEG) treatment. Error bars represent standard error of the mean. *2.4. Target Sites of miR319 in GhTCP Genes*

*2.4. Target Sites of miR319 in GhTCP Genes*  The miRNAs can cause endonucleolytic cleavage of mRNA by extension, which often perfect complementarity to mRNAs. In plants, miR319 was one of the first characterized and conserved miRNA families, which has been demonstrated to target *TCP* genes. In *G. hirsutum*, miR319 had only 1-nt mismatch compared with sequences in Arabidopsis (Figure 6). The predicted hairpin structures of the miR319 precursor had 191 nt. The miR319 sequence was located at the 3′-end of the pre-miRNAs and began with a 5′-uridine. Using a set of strict standards, *GhTCP21*, *GhTCP31*, and *GhTCP54* were predicted as targets of miR319. These three miR319 target sites were all located in the *2.4. Target Sites of miR319 in GhTCP Genes*  The miRNAs can cause endonucleolytic cleavage of mRNA by extension, which often perfect complementarity to mRNAs. In plants, miR319 was one of the first characterized and conserved miRNA families, which has been demonstrated to target *TCP* genes. In *G. hirsutum*, miR319 had only 1-nt mismatch compared with sequences in Arabidopsis (Figure 6). The predicted hairpin structures of the miR319 precursor had 191 nt. The miR319 sequence was located at the 3′-end of the pre-miRNAs and began with a 5′-uridine. Using a set of strict standards, *GhTCP21*, *GhTCP31*, and *GhTCP54* were predicted as targets of miR319. These three miR319 target sites were all located in the coding regions, and all miR319-targeted genes belonged to the CIN clade. Similarly, there were five The miRNAs can cause endonucleolytic cleavage of mRNA by extension, which often perfect complementarity to mRNAs. In plants, miR319 was one of the first characterized and conserved miRNA families, which has been demonstrated to target *TCP* genes. In *G. hirsutum*, miR319 had only 1-nt mismatch compared with sequences in Arabidopsis (Figure 6). The predicted hairpin structures of the miR319 precursor had 191 nt. The miR319 sequence was located at the 30 -end of the pre-miRNAs and began with a 50 -uridine. Using a set of strict standards, *GhTCP21*, *GhTCP31*, and *GhTCP54* were predicted as targets of miR319. These three miR319 target sites were all located in the coding regions, and all miR319-targeted genes belonged to the CIN clade. Similarly, there were five and three *TCP* genes containing miR319-binding sites in Arabidopsis and *P. mume*, respectively, and they also

coding regions, and all miR319-targeted genes belonged to the CIN clade. Similarly, there were five

belonged to the CIN group [26]. This suggests that miR319 held homologous target interactions during the evolution and diversification of plants. they also belonged to the CIN group [26]. This suggests that miR319 held homologous target interactions during the evolution and diversification of plants. and three *TCP* genes containing miR319-binding sites in Arabidopsis and *P. mume*, respectively, and they also belonged to the CIN group [26]. This suggests that miR319 held homologous target interactions during the evolution and diversification of plants.

and three *TCP* genes containing miR319-binding sites in Arabidopsis and *P. mume*, respectively, and

**Figure 6.** Mature and predicted fold-back structures of miR319 precursors in *G. hirsutum*. Sequences of mature miR319 are underlined. **Figure 6.** Mature and predicted fold-back structures of miR319 precursors in *G. hirsutum*. Sequences of mature miR319 are underlined. **Figure 6.** Mature and predicted fold-back structures of miR319 precursors in *G. hirsutum*. Sequences of mature miR319 are underlined.

Degradome sequencing had been widely used to identify plant miRNA cleavage sites. In this study, the *GhTCP* mRNA degradation sites were determined by BLASTing the sequenced degraded fragments against the *G. hirsutum* TCP genes. The degradome sequencing data are available at Gene Expression Omnibus (GEO accession number GSE69820). Using PairFinder software, miR319 was found to cleave *GhTCP21* and *GhTCP54* mRNA transcripts (Figure 7). The miR319–mRNA pair was at the cleavage site of the two TCP genes. There were both 67 raw reads at the position, with abundance at the position equal to the maximum on the transcript, and with only one maximum on the transcript. The 5′-ends of the mRNA fragments mapped to the nucleotide that paired to the tenth nucleotide of the miR319 sequence. To research the biological function of miR319, a negative-correlation expression test was undertaken for miR319 and its target *GhTCP21* and *GhTCP54* mRNAs using qRT-PCR. In response to salt and drought treatments, miR319 showed different degrees of upregulation. The transcriptome sequencing analysis showed that expression of Degradome sequencing had been widely used to identify plant miRNA cleavage sites. In this study, the *GhTCP* mRNA degradation sites were determined by BLASTing the sequenced degraded fragments against the *G. hirsutum* TCP genes. The degradome sequencing data are available at Gene Expression Omnibus (GEO accession number GSE69820). Using PairFinder software, miR319 was found to cleave *GhTCP21* and *GhTCP54* mRNA transcripts (Figure 7). The miR319–mRNA pair was at the cleavage site of the two TCP genes. There were both 67 raw reads at the position, with abundance at the position equal to the maximum on the transcript, and with only one maximum on the transcript. The 50 -ends of the mRNA fragments mapped to the nucleotide that paired to the tenth nucleotide of the miR319 sequence. To research the biological function of miR319, a negative-correlation expression test was undertaken for miR319 and its target *GhTCP21* and *GhTCP54* mRNAs using qRT-PCR. In response to salt and drought treatments, miR319 showed different degrees of upregulation. The transcriptome sequencing analysis showed that expression of *GhTCP21* and *GhTCP54* was downregulated. Degradome sequencing had been widely used to identify plant miRNA cleavage sites. In this study, the *GhTCP* mRNA degradation sites were determined by BLASTing the sequenced degraded fragments against the *G. hirsutum* TCP genes. The degradome sequencing data are available at Gene Expression Omnibus (GEO accession number GSE69820). Using PairFinder software, miR319 was found to cleave *GhTCP21* and *GhTCP54* mRNA transcripts (Figure 7). The miR319–mRNA pair was at the cleavage site of the two TCP genes. There were both 67 raw reads at the position, with abundance at the position equal to the maximum on the transcript, and with only one maximum on the transcript. The 5′-ends of the mRNA fragments mapped to the nucleotide that paired to the tenth nucleotide of the miR319 sequence. To research the biological function of miR319, a negative-correlation expression test was undertaken for miR319 and its target *GhTCP21* and *GhTCP54* mRNAs using qRT-PCR. In response to salt and drought treatments, miR319 showed different degrees of upregulation. The transcriptome sequencing analysis showed that expression of *GhTCP21* and *GhTCP54* was downregulated.

**Figure 7.** miR319 and its target genes, *GhTCP21* and *GhTCP54,* in *G. hirsutum*. (**A**) Target plot (t-plot) for *GhTCP21* and *GhTCP54*, which were targeted by miR319. Arrows indicate the signatures **Figure 7.** miR319 and its target genes, *GhTCP21* and *GhTCP54,* in *G. hirsutum*. (**A**) Target plot (t-plot) for *GhTCP21* and *GhTCP54*, which were targeted by miR319. Arrows indicate the signatures **Figure 7.** miR319 and its target genes, *GhTCP21* and *GhTCP54,* in *G. hirsutum*. (**A**) Target plot (t-plot) for *GhTCP21* and *GhTCP54*, which were targeted by miR319. Arrows indicate the signatures corresponding to the miRNA cleavage site. Partial mRNA sequences of target genes aligned with the miRNAs show perfect matches (straight lines) and G-U wobbles (circles). (**B**) Relative expression levels of miR319, *GhTCP21*, and *GhTCP54* under salt and drought stresses. QRT-PCR analyses were performed using RNA generated from cotton leaves after NaCl and Polyethylene Glycol (PEG) treatment. Error bars represent standard error of the mean.
