Functional Characterization of Candidate Genes, Gohir.D05G103700 and Gohir.D12G153600, Identified Through the Expression QTL Analysis using Virus Induced Gene Silencing in Upland Cotton (Gossypium hirsutum L.)

1. Introduction

The genus Gossypium emerged almost 12.5 million years ago (mya), and about 5-10 mya diverged into diploid A and D sub-genome species [1]. Subsequently, around 1-2 mya, the two-diploid species underwent trans-oceanic hybridization and chromosome doubling, resulting in an allopolyploid AD clade that radiated from the Americas and diversified into multiple species [1]. The allopolyploid cotton, Gossypium hirsutum L. (Upland cotton), and Gossypium barbadense L. (Egyptian cotton) have originated from two independent domestication events [2]. Today, G. hirsutum dominates the world’s cotton fiber production (90%). Whereas G. barbadense, G. arboretum, and G. herbaceum together contribute only 10% to it [3]. Cotton production generates an annual economic impact of around 600 billion US dollars globally, making it a major contributor to the world’s economy [4].

The plant architecture plays a crucial role in determining crop production and its management. It is determined by the meristematic cells and could be indeterminate, leading to continuous vegetative growth, to determinate, where meristematic cells are consumed in producing a terminal structure. Primary meristems are dynamically added in cotton-producing monopodial or sympodial branches, inflorescences, and terminal shoots, i.e., flowers. Similarly, the secondary meristems (cambium) add continuously and contribute to wood production. The position of organ-specific meristems called primordia contributes to the final plant architecture [5]. Flowering in cotton depends on the differentiation of flower primordia that depends on axillary bud primordia, which in turn depends on the shoot apical meristem (SAM) differentiation. Axillary bud primordia initiate the transition from vegetative to the reproductive phase in cotton by differentiating from leaf primordia to flower primordia at a specific time point during plant development.

Transcription factors (TF) regulate tissue-specific gene expression by binding to the promoter region of the target genes. Subsequently, the DNA-TF complex promotes the recruitment of the RNA polymerase onto the target sequence and promotes gene expression. The GATA TFs (one of the important family of transcription factors) in plants express at different developmental stages, especially the light-induced processes, such as embryo development, flowering, petal differentiation and expansion, and maturation [6,7,8]. The NTL1 was the first GATA TF gene identified in Nicotiana tabacum [9]. So far, over 30 GATA TFs have been identified and characterized in Arabidopsis thaliana [10]. Additionally, 64 GATA TFs have been identified in Glycine max, 30 GATA TFs in Gossypium spp., and 28 GATA TFs in Oryza sativa [10,11,12]. These studies highlight the importance of GATA TFs in plant growth and development.

Another diverse family of transcription factors, SQUMOSA PROMOTER BINDING PROTEIN-LIKE (SPL), plays a crucial role in plant growth and development [13,14]. The SPL transcription factors play various roles in the growth and development of plants, such as controlling vegetative to reproductive phase change in A. thaliana [15,16], ovary and fruit development in tomato [17], yield-related traits in wheat [18], and floral organ size and ovule production in cotton [19].

Virus-induced gene silencing (VIGS) is an important reverse genetics tool for in vivo gene function study in plants [20]. VIGS depends on posttranscriptional gene silencing (PTGS) machinery in a sequence-specific manner. Briefly, in this machinery, double-stranded RNAs (dsRNAs), such as the replication intermediates of plant RNA viruses and highly structured genomic RNA, are processed by Dicer-like proteins (DCLs) to form small interfering RNAs (siRNAs). The guide strand of siRNA gets loaded onto the RNA-induced silencing complex (RISC) and cleaves complementary messenger RNA (mRNA) [21]. For VIGS, fragments of the gene of interest are cloned into a viral vector and applied to plants where the endogenous RNA silencing machinery in the host has degraded the viral RNA and, in this process, produced complementary siRNAs that reduce the target gene expression. VIGS has been demonstrated using numerous plant-virus combinations [22,23]. In the past decade, several viral genomes have been modified as a powerful reverse genetic tool for the functional characterization of plant genes. Some examples of VIGS constructs are the Tobacco rattle virus (TRV) [24], Apple latent spherical virus (ALSV) [25,26], African cassava mosaic virus (ASMV) [27], Cucumber mosaic virus (CMV) [28], and Barley streak mosaic virus (BSMV) [29].

TRV consists of a bipartite genome, TRV1, and TRV2 [30], where TRV1 determines viral movement [31] and TRV2 codes for the coat protein and other nonstructural proteins. Also, it shows a great deal of variability among different viral isolates [30]. These nonstructural proteins are involved in the nematode transmission of this virus [32], but these proteins are not essential for plant infection. Therefore, TRV2’s nonstructural protein-coding genes were replaced by multiple cloning site where fragments of the target gene could be introduced [33]. TRV1 and TRV2 were later cloned into a binary vector for Agrobacterium-mediated plant delivery [24,33].

TRV-VIGS has been effectively used to manipulate gene expression in different plant species such as A. thaliana, Solanum lycopersicum, Nicotiana spp., G. arboreum, and Petunia hybrida [34]. TRV-based VIGS is one of the most effective methods to study gene function in plants due to several advantages, such as mild symptoms, spread to large cell patches, ability to migrate to the growing meristems, and systemic spread to new tissues and all plant parts [35,36,37].

Given the desirable attributes of TRV and its successful implementation in VIGS in different plant species, in this study, we undertook the TRV-based VIGS of two candidate genes, Gohir.D05G103700 and Gohir.D12G153600. Using upland cotton mini-core collection, we identified these genes earlier by expression QTL (eQTL) analysis for five floral induction and meristem identity genes (FT, SOC1, AP1, FUL, and LFY) using the upland cotton mini-core collection (Naveed et al. unpublished data). In this analysis, Gohir.D05G103700 showed an association with the APETALA 1 (AP1) gene expression trait and annotated as GATA transcription factor 11-like, and Gohir.D12G153600 showed an association with FRUITFUL (FUL) and SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) gene expression traits and annotated as SQUAMOSA promoter binding-like transcription factor. This work is focused on characterizing the molecular function of these eQTLs and determining their role in cotton architecture development. Unveiling the molecular functions of these genes might have implications in breeding cotton for reduced regrowth after defoliation or annual growth habit with improved lint yield.

2. Materials and Methods

2.1. Plant husbandry

Cotton cultivar Coker 201 was selected for the VIGS experiment, as this cultivar was used in the past for genetic transformation experiments using Agrobacterium, which is also an intended method for delivering the TRV constructs here [38]. Seeds of cotton cultivar Coker 201 were germinated in 0.5-gallon plastic pots filled with an autoclaved 3:1 mixture of Fafard 3B and playground sand. The pots were placed in a walk-in growth chamber (GR96, Conviron, Pembina, ND) set at 30°C daytime and 20°C nighttime temperature with a 16h photoperiod. Three sets of six plants were placed in the growth chamber (one set each for two genes of interest and one as a mock control). All plants were treated alike regarding watering and fertilization for the entire length of the study, and the plant ruminants were discarded after autoclaving at the end of the study. The whole experiment was repeated once.

2.2. Construction of the gene silencing constructs

We selected two candidate genes, Gohir.D05G103700 and Gohir.D12G153600 (Table 1), which showed association with the cotton AP1 and FUL/SOC1 gene expression traits, respectively, for the functional characterization using the TRV-based VIGS system. Target sites for PCR primer designing from the cotton Gohir.D05G103700 and Gohir.D12G153600 genes were selected based on the conservation of sequences between the homoeologous genes and were determined using multiple sequence alignment (MSA) on the Clustal Omega software. Primers flanking the conserved regions in these genes were designed using the Primer 3 software (https://primer3.ut.ee/). The primers were used to amplify a 175 bp and a 166 bp fragment for Gohir.D05G103700 and Gohir.D12G153600 genes from cotton cultivar Coker 201 (for the list of primers, see Supplementary Table S1). For PCR cloning of these fragments in TRV2, primers were supplemented with the XbaI and KpnI restriction sites. After PCR, the product was digested with XbaI and KpnI and ligated with the TRV2 plasmid pYL156 linearized formerly using the same restriction enzymes. Subsequently, the construct was transformed into E. coli DH5α cells by heat shock method. The plasmid DNA was isolated and validated by PCR, restriction digestion, and DNA sequencing (McLAB, San Francisco, CA), and later whole plasmid sequence was obtained from Plasmidsaurus (Eugene, OR) (Supplemental Figure S1).

2.3. Agrobacterium tumefaciens transfection with TRV1 and TRV2 constructs and plant inoculation

The detailed structure of TRV1 and TRV2 constructs amenable for Agrobacterium-mediated gene silencing has been described previously [24,33,34]. The TRV2 plasmid containing the Gohir.D05G103700 gene sequence (pYL156-Gohir.D05G103700), hereafter referred to as TRV2-GATA-TF and pTRV1 were transfected into Agrobacterium tumefaciens strain GV3101 by the heat shock method and selected on plates supplemented with kanamycin (50μg/mL), gentamycin (50μg/mL), and rifampicin (10μg/mL). Likewise, a TRV2-SPL-TF and an empty viral control TRV2-MCS (mock control) were developed and transfected into A. tumefaciens strain GV3101. Subsequently, the leaves were agroinfiltrated with a 1:1 mixture of TRV1 and recombinant TRV2 constructs [34]. A schematic representation of these steps is presented in Figure 1. The leaf tissues were collected from the infiltrated plants 14 days after inoculation in liquid nitrogen and stored at -80°C until used.

2.4. RNA extraction and cDNA conversion

Total RNA was extracted from leaf tissue using the Spectrum™ Plant Total RNA Kit (Sigma-Aldrich, USA) following the manufacturer’s instructions with minor modifications. RNA concentrations were determined using NanoDrop™ (Thermo Scientific, USA). The total RNA (1μg) was converted to cDNA using the Revert Aid first strand cDNA synthesis kit (Thermo Scientific, USA) following the manufacturer’s instructions using the oligo d(T) primer.

2.5. Quantitative Real-Time PCR (qRT-PCR)

The quantitative real-time PCR (qRT-PCR) analysis was performed using the iTaq™ Universal SYBR® Green Supermix (BioRad, USA) following the manufacturer’s instructions on iCycler iQ™ (BioRad, USA). The following PCR profile was used for the gene expression analysis: initial denaturation at 95°C for 3 min followed by 50 cycles at 95°C for 10s and 56°C (Gohir.D05G103700 and ACTIN1) or 59.5°C (Gohir.D12G153600) for 30 s and 50 cycles of melt curve analysis at 56°C to 93.5°C for 10 s and final incubation at 4°C. The cotton ACTIN1 gene was used as an internal control for data normalization. The gene expression data were analyzed using the 2^–∆∆Ct method and presented as log 2 of the fold change relative to the reference gene [39].

2.6. Phenotypic data analysis

The data on plant height, number of leaves per plant, number of squares per plant, and number of bracts per square were collected from the inoculated and mock control plants in the walk-in growth chamber (Supplementary Table S2). The statistical analysis of the phenotypic data was performed using the SAS package.

3. Results

3.1. Sequence analysis for the antisense construct development

The DNA sequence of the cotton Gohir.D05G103700 and Gohir.D12G153600 genes identified as eQTLs on chromosomes D05 (SNP markers i09222Gh and i00443Gh at 8730.32 Mb) and D12 (SNP marker i08185Gh at 48443.48 Mb) for the cotton AP1 and FUL/SOC1 gene expression traits respectively were retrieved from NCBI (National Center for Biotechnology Information) website. Both genes existed as homoeologous gene pairs. The former gene mapped to upland cotton chromosomes A05 (LOC107906238; location 8729.803-8733.656 Mb) and D05 (LOC107959021; location 9724.645-9728.563 Mb) and coded for a GATA transcription factor and the latter gene mapped to chromosomes A12 (LOC107934496; location 91532.20-91536.64 Mb) and D12 (LOC107946275; location 48443.14-48445.05 Mb) and coded for a SQUAMOSA PROMOTER BINDING LIKE protein (SPL). When the homoeologous gene sequences were aligned, these showed 95-96% similarity. A 175 bp antisense construct for GATA-TF and a 166 bp antisense construct for the SPL-TF genes were designed from sites showing 99-100% sequence similarity or conservation (Figure 2). The primer pairs were designed to amplify and PCR-clone products of expected size from the upland cotton cultivar ‘Coker-201’ in the TRV2 vector. The resulting plasmids were confirmed via DNA sequencing, and the maps were constructed. After sequence confirmation, the construct was transformed to the Agrobacterium tumefaciens strain GV3101. Likewise, the TRV1 construct was also transformed to the Agrobacterium tumefaciens GV3101 and bacteria with TRV1 and TRV2 were mixed at 1:1 proportion, and the plants were inoculated with the medley [34].

3.2. Agrobacterium mediated-VIGS of Gohir.D05G103700 and Gohir.D12G153600 genes in the upland Cotton Cultivar ‘Coker-201’

Three-week-old plants were inoculated with TRV2-GATA-TF plus TRV1 and in a separate experiment with TRV2-SPL-TF plus TRV1. Leaf samples were collected from inoculated and control plants two weeks after the treatment for RNA extraction and qRT-PCR analysis. The qRT-PCR was performed for the cotton Gohir.D05G103700 and Gohir.D12G153600 genes, and the cotton housekeeping ACTIN1 gene was used for the normalization of the gene expression data. The TRV2-GATA-TF showed an about 1.4-fold reduced abundance of the cotton Gohir.D05G103700 transcript in treated plants, G3 and G4 relative to the uninoculated control, hence, suggested downregulation of the gene (Figure 3A). Other inoculated plants, G2, G5, and G6 except G1 showed a negligible reduction in Gohir.D05G103700 transcript abundance, whereas, no silencing was observed in G1.

Likewise, qRT-PCR analysis of TRV2-SPL-TF inoculated plants showed up to a 2.3-fold reduction in Gohir.D12G153600 transcript abundance in P05 relative to uninoculated control (Figure 3B). Whereas, about 0.8-fold reduction in transcript abundance was observed in rest of the treated plants, P01, P04, and P06, except P02 that exhibited wild type expression level.

3.3. Expression profiling of the cotton AP1 genes in TRV2-GATA-TF inoculated plants and the FUL and SOC1 genes in TRV2-SPL-TF inoculated plants

Since the cotton Gohir.D05G103700 and Gohir.D12G153600 genes were identified as eQTLs for the cotton AP1 and FUL/SOC1 genes, respectively, we studied the effect of silencing these genes on the transcription levels of the cotton AP1 gene in TRV2-GATA-TF inoculated plants and the cotton FUL and SOC1 genes in TRV2-SPL-TF inoculated plants. Among the TRV2-GATA-TF inoculated plants, G3 exhibited a ~1.4-fold reduction in the transcript level of the cotton AP1 gene. In contrast, no reduction was observed in the case of G1, G4, G5, and G6, and a slight over-accumulation was observed in the case of G2 (Figure 3C).

Similarly, among the TRV2-SPL-TF inoculated plants, P01 and P03 showed up to a 1-fold reduction in the cotton FUL transcript level, P04 and P05 exhibited negligible reductions, and P02 showed no reduction. On the other hand, P06 showed over-accumulation (Figure 3D). The transcription level of the cotton SOC1, TRV2-SPL-TF inoculated plants showed up to a 1.2-fold reduction in P05 and ~ 1-fold reduction in P03. A negligible reduction in transcript abundance of cotton SOC1 was observed in P04. On the other hand, P01, P02, and P06 did not show any reduction in gene expression (Figure 3E).

3.4. Phenotypic observations

Data on various traits such as plant height, leaf and flower number per plant, and number and arrangement of floral organs (including the number of floral bracts per square) were recorded from the inoculated and mock control plants grown in the walk-in growth chamber (Supplementary Table S2). Interestingly, we observed loss of the supernumerary (fourth) floral bract in the squares of the TRV2-GATA-TF inoculated plants, G3 and G4 with a relatively higher level of Gohir.D05G103700 suppression (Figure 4). None of the other TRV2-GATA-TF inoculated plants and mock-treated plants exhibited suppression of the fourth bract in the floral buds of cotton cv. Coker-201. On the other hand, the TRV2-SPL-TF inoculated plants showed no observable phenotypic change relative to the mock-treated plants.

4. Discussion

Upland cotton an allopolyploid has a large and complex genome. Genetic transformation in this crop was successful and was the source of several herbicide tolerant and insect pest resistant genotypes available commercially [41]. However, genetic transformation is still not commonplace in cotton and is highly genotype dependent and time consuming [42]; hence cannot be applied to test gene functions. Virus-induced gene silencing offers a solution to this problem, as the experiments are relatively less time consuming. However, the gene-silencing effects are broadly transient, which has both pros and cons [43]. There are two viral systems available for gene silencing and gene delivery in upland cotton, Cotton leaf crumple virus (CLCrV) and Tobacco rattle virus (TRV) [43]. Both systems were shown to work well; however, TRV-based system is more commonly used and hence was adopted in the present study for functional characterization of Gohir.D05G103700 and Gohir.D12G153600.

Despite the fact that we used an upland cotton cultivar ‘Coker 201’ that earlier showed susceptibility to Agroinfiltration [38], hence consistent delivery and systemic spread of the recombinant virus, we observed plant-to-plant variation in the silencing efficiency in our experiments. For instance, five of the six plants showed some level of suppression; however, the level of suppression ranged from about 0.4-fold to 1.4-fold for Gohir.D05G103700 and from about 0.8-fold to a little over 2.2-fold for Gohir.D12G153600 in two and a single plant, respectively. A review of the literature suggested that variation in the level of gene suppression is not uncommon in VIGS experiments, as was documented in several earlier studies performed in different crop plants [44]. Several reasons were proposed to contribute to the observed variation in the level of gene suppression, from the effect of abiotic factors to genetic factors [45,46]. Further, in this study, we observed incomplete gene silencing for either of the genes. This incomplete silencing could be attributed to multiple factors, including the plant genotype, ploidy level, the plant development stages, temperature, and inoculation method [47,48,49]. Further, it was demonstrated that TRV-based gene silencing was more efficient in the diploid (G. arboreum and G. herbaceum) than in tetraploid cotton (G. hirsutum) [49]. Furthermore, the inoculation method is a critical factor for the successful delivery of the VIGS construct and subsequent steps involving systemic spread and gene silencing [29].

In addition, the two genes Gohir.D05G103700 and Gohir.D12G153600, differed in the level of suppression. There could be several possible reasons, including VIGS is a random process and may not lead to a similar level of suppression in two different attempts, as the long antisense RNA after complementation with the target RNAs process is somewhat random fashion. Also, the two homoeologous copies of a gene do not carry perfect similarity; hence, they do not exhibit similar levels of suppression. Unfortunately, it is quite difficult to develop homoeologue-specific primers to study the level of suppression of individual homoeologues. Hence, a cumulative level of gene suppression was studied in most of the studies [50,51]. In some cases, the antisense was able to silence both copies equally well, and some only one of the two copies. Additionally, often times a compensatory increase in the expression level of one of the homoeologues was observed upon silencing the other homoeologue [45].

Phenotypic analysis revealed that Coker 201 tends to produce a supernumerary bract (fourth bract) in each flower (visible in the mock control), which was absent from GATA-VIGS inoculated plants (Figure 4A). The earlier research suggested that cotton tends to produce a fourth bract in certain environmental conditions when the meristem in the axil of the true leaves is not consumed fully in flower production see Figure 4B [40]. It was proposed that the three bracts in the cotton bud are modifications of one true leaf and two stipules, and the fourth bract represents the second leaf (indicative of indeterminate growth). These results suggested that the GATA transcription factor identified as an eQTL for the cotton AP1 gene has a target site in the AP1 promoter, where it binds and occupies a site not letting the FT/FD assembly bind and induce/promote its expression, which triggers flowering and hence consumption of the meristem in floral development [12,52,53,54]. We hypothesized that VIGS of the GATA transcription factor released the site from suppression, allowing FT/FD-induced activation of the AP1 gene. The qRT-PCR analysis, however, did not exhibit AP1 over accumulation. Furthermore, the fourth bract in cotton is an undesirable trait [40], as its enlargement creates a gap allowing thrips and other insects to enter and feed on the developing fruit. Functional characterization of the Gohir.D05G103700 gene identified as the cause of the undesirable fourth floral bract and made available the target to select for this trait. However, we did not see any observable phenotypic changes in the TRV2-SPL-TF plants, which suggested the functional redundancy of this transcription factor in cotton or a need to completely silence the gene to have a notable phenotypic effect. Quantitative RT-PCR analysis of the AP1 and the SOC1/FUL genes established a connection between these genes and eQTLs, Gohir.D05G103700 and Gohir.D12G153600. In sum, this study reconfirmed the utility of TRV-based VIGS in cotton by functional characterization of an eQTL target for the AP1 expression trait. Also, the study reflected on the transcriptional regulation of the cotton AP1 gene and provided a target to breed for thrip resistance in cotton cultivars with a supernumerary bract (fourth bract).