Genome-Wide Identification of Cotton MicroRNAs Predicted for Targeting Cotton Leaf Curl Kokhran Virus-Lucknow

1. Introduction

The allotetraploid upland cotton (Gossypium hirsutum L.) is an important essential fiber-producing industrial cash crop which is grown on several continents of the world attributed 40% world’s agriculture production. Cotton is an important natural textile fiber made up of 90% cellulose and is a renewable natural source for global textiles industry. The quality parameters of cotton fibers have a direct influence on agro-based textile industry [1,2,3]. Being Allotetraploid nature, upland cotton genome contains 52 chromosomes (2n = 4x = 52)[4,5]. The first draft genome sequences of model polyploidy upland allotetraploid cotton species was accomplished and physical map was released in 2015 [6]. Cotton leaf curl disease (CLCuD) is a major biotic and economic constraint of cotton and is caused by at least four species in the genus Begomovirus (Geminiviridae) [7,8]. CLCuD complex is caused by monopartite begomoviruses (BGVs) [9,10]. CLCuD outbreaks caused by CLCuKoV significantly constrain cotton production in Pakistan [11,12,13]. BGVs are known to be spread efficiently by polyphagous whitelyfly vector, i.e., Bemisia tabaci [14,15,16,17].

Cotton leaf curl Kokhran virus-Lucknow (CLCuKoV-Lu) has recently emerged, and has been reported as new ‘strain’ of CLCuKoV in the begomovirus genus of the Geminiviridae family [18,19]. The architecture of the CLCuKoV genome inside virion is made up of a single copy of circular, single-stranded (ss) DNA molecule. The monopartite genome of a CLCuKoV encodes six proteins. The transcription and regulation of these proteins are governed from a large intergenic region (LIR) using bidirectional mode of transcription [12,20,21]. The plus (+) virion-sense (VS) and negative (-) complementary-sense (CS) strands encode ORFs V1 and V2, and ORFs C1, C2, C3 and C4, respectively [15,18,22]. CLCuKoV-Lu was first noticed in Lucknow in 2010 on guar (Cyamopsis tetragonoloba) plant, exhibiting leaf curl symptoms [19]. Rolling circle amplification (RCA) is widely used molecular technique to identify and isolate BGVs. Loop-mediated isothermal amplification (LAMP) assay, quantitative real-time polymerase chain reaction (qPCR) assay, multiplex PCR and immunofluorescence assay are also optimized and standardized as molecular diagnostics methods for CLCuD [23,24,25,26,27,28,29].

Although efforts have been redirected to control CLCuD, there are no resistant cultivars available yet in Pakistan. Considering the DNA nature of CLCuKoV, RNA interference (RNAi) has emerged as a robust and efficient tool for targeting microRNA-induced silencing complex (miRISC)-mediated gene silencing in eukaryotes [30,31,32,33]. RNAi-based silencing is a sequence-specific double-stranded (ds) RNA-mediated gene regulatory innate antiviral, natural mechanism to inhibit viral gene replication and transcription. The RNAi machinery is composed of two core central components, Dicer and Argonaute. These are involved for loading of 20-30 nucleotides long RNA molecule for processing of miRNA/miRNA* duplex and further was incorporated into RISCs. The dsRNA was processed and are broken down to short 21–24 nucleotides siRNA that was degraded. [34,35,36,37,38]. The plant microRNAs (miRNAs) are the smallest, widespread, highly conserved, non-coding, single stranded (ss) RNA molecules ranging from 18-24 nt encoded by MIR genes. The endogenous miRNAs play important role in regulating plant gene expression and key biological process by recruiting RISC complex [39,40,41].

Upland Cotton plant has been investigated for potential diverse molecular mechanism and was explored mature miRNAs that are critical for normal growth [42,43,44]. Cotton miRNAs are good source of immunity to control biotic and abiotic stress response networks [45,46,47,48,49]. The artificial microRNA (amiRNA)-based gene silencing impart resistance against invading viruses [50]. The amiRNA construct was transformed to induce gene silencing of target plant virus in 2006 [51]. Experimentally verified locus-derived mature miRNAs in the cotton genome was predicted to regulate gene expression. A subset of high-confidence mature miRNAs in cotton is designed to have predicted target sites in the CLCuKoV-Lu genome.

Our present an algorithmic in silico approach relied on prediction of upland allotetraploid cotton genome-encoded miRNAs targeting CLCuKoV-Lu. The research work summarized key strategies for analyzing most effective target sites of ghr-miRNAs in the CLCuKoV-Lu genome. To explore host-virus complex mechanism, miRNA-mRNA target site interactions were also explored. The study aims to elucidate the predicted ghr-miRNAs for the generation of CLCuKoV-Lu-resistant cotton plants in future.

2. Materials and Methods

2.1. Biological Data Retrieval

There are currently 80 mature cotton locus-derived ghr-microRNAs (commonly called Gossypium hirsutum -microRNAs (ghr-miR156-ghr-miR7514) (Accession IDs: MIMAT0005806-MIMAT0029164) (Table S1 (supplementary Materials)) identified. In a previous study, 78 stem-loop (precursor) cotton locus-derived ghr-miRNAs (ghr-MIR156-ghr-MIR7514) (Accession IDs: MI0005638-MI0024206) have been identified (Table S2). The cotton locus-derived ghr-miRNAs were retrieved from the miRBase latest version (v22.1) database (http://mirbase.org/) (accessed on 26 December 2018). miRBase offers biological information and annotation on miRNAs [52]. The full-length genome (2761 nucleotides) of CLCuKoV-Lu (accession number GU385879) was downloaded from the NCBI GenBank database (http://ncbi.nlm.nih.gov) (accessed on 19 December 2018) [53].

2.2. Target Prediction

The current in silico approach was based on the most widely used, publicly available miRNA prediction algorithms such as miRanda, RNA22, psRNATarget, and RNAhybrid. These algorithms were used to identify the CLCuKoV-Lu genome predicting the ‘most effective’ miRNA binding sites of the cotton miRNAs (Table 1). The G.hirsutum locus-derived ghr-miRNA sequences and the CLCuKoV-Lu genome-predicted transcripts (in FASTA format) were analyzed.

2.3. miRanda

miRanda is considered seed-matching, scoring scheme, dynamic programming and conservation consisting of three major steps to identify target sites since its release in 2003 [54]. miRanda has been used to predict the candidate target sites (CTSs) based on minimum free energy hypothesis. It works on RNA duplex dimerization and target sequence complementarity score [55]. The miRanda algorithm (written in C programming language) is a user-friendly web tool. The miRanda software was downloaded using the online source website (http://www.microrna.org/) (accessed on 23 March 2019). Prediction analysis was performed at default parameters. The miRanda program was set at MFE threshold: −20 Kcal/mol, score threshold: 140.00, gap open penalty: −9.000, gap extend penalty: −4.000 and scaling parameter: 4.00.

2.4. RNA22

The RNA22 algorithm is considered patter-recognition approach based on non-seed interaction for the prediction of target binding sites.[56]. The RNA22 is non-seed based algorithm and can be accessed web server (http://cm.jefferson.edu/rna22v1.0/) (accessed on January 26, 2019). Non-seed-based biological significant miRNA-mRNA interaction was calculated using target patterns and maximum folding energy (MFE)[57]. Site complementarity is also unique feature of RNA22 algorithms. Default parameters were set after loading cotton miRNA sequences and CLCuKoV-Lu genomic template: (63%), specificity (61%) and MFE (−15.00 Kcal/mol).

2.5. psRNATarget

The psRNATarget algorithm is considered small RNA target prediction using plant specific features coupled with seed region matching and scoring scheme. The web server is used to predict binding sites of the plant miRNAs listed in the web server on the basis of complementary scoring schema [58]. The psRNATarget web server has been updated and developed to report the inhibition pattern [59]. The fasta sequence of CLCuKoV-Lu genome was set in the server after selecting G. hirsutum miRNAs in the psRNATarget web server (http://plantgrn.noble.org/psRNATarget) (accessed on 26 October 2020). Target sites of cotton miRNAs were identified using default criteria of prediction: expectation cut-off value 6.5 and mode of inhibition ‘Cleavage’ was selected.

2.6. RNAhybrid

RNAhybrid is a seed-based flexible online web-based computational algorithm that is based on intermolecular hybridization predicting accurate miRNA binding sites in the target sequence. The RNAhybrid algorithm is MFE-based hybridization model. Site complementarity, MFE and seed match key features [60]. The fasta sequence of CLCuKoV-Lu genome and miRNAs were loaded after setting defaults parameter in the RNAhybrid web server http://bibiserv.techfak.uni-bielefeld.de/rnahybrid (accessed on 2 19 December 2020). MFE is a model feature and set at (−20.00 Kcal/mol).

2.7. RNAfold and RNAcofold

RNAfold is a web-based algorithm to identify accurate secondary structures of the target single-stranded (ss) miRNA precursor [61]. Precursor sequences was uploaded in the web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) (accessed on September 6, 2022) under user-defined default settings. RNAcofold is a web-based algorithm and specifically intended to estimate the cofolding free energy (ΔG) of RNA duplex sequences. MFE and base-pairing pattern of miRNA-mRNA target duplex were estimated to evaluate miRNA-mRNA interaction [62]. FASTA sequences of unique duplex pair were set in the RNAcofold web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAcofold.cgi) (accessed May 6, 2022) under user-defined default settings.

2.8. Mapping of Network-based miRNA-Target Interactions

A Circos plot was generated using CIRCOS algorithm. We used the circos package v0.69-9R-Language [63].

2.9. Identification of miRNA Binding Site Distribution

The genomic transcript of monopartite CLCuKoV strains were downloaded from NCBI Genbank nucleotide database (accessed on 26 March 2018). The accession numbers of CLCuKoV strains under study are classified: GU385879.2, AM421522.1, AJ496286.1, HF549182.1, and FN5520001.1[18]. MEGA X (version 10.0.5) was used to align conservation of binding affinity of the predicted consensus miRNAs [64]. MUSCLE, a multiple sequence alignment algorithm was used to align a set of CLCuKoV genomic sequences [65]. The miRNA binding site sequences were mapped CLUSTALW algorithm[66].

2.9. Statistical Analysis

miRNA target prediction data were further processed for graphical representations using R-language (version 3.1.1, software version 3.5.1 [67].

2.10. Genome Annotation

pDRAW32 DNA analysis (AcaClone 1.1.147) is DNA annotation tool and was downloaded from the source web. It is used for ssDNA genomic annotation of CLCuKoV-Lu and further editing. The genome analysis and production of graphical output of CLCuKoV-Lu ORFs were generated.

3. Results

3.1. Cotton miRNAs-mRNA Interaction Pairs in CLCuKoV-Lu Genome

The genetic map of CLCuKoV-Lu has a single copy of circular ssDNA biomolecule that is composed of 2750 nucleotides (Figure 1).

Figure 1. The genetic makeup of CLCuKoV-Lu is composed of six ORFs. The colored arrows are showing ORFs. The plus (+)-strand comprises virion-sense (encapsidated strand) ORFs (V1 and V2). The negative (−)-strand is composed of complementary-sense strand ORFs (C1, C2, C3 and C4). Replication of viral genome and transcription of the viral coding strands are governed by a non-coding large intergenic region (LIR).

Figure 1. The genetic makeup of CLCuKoV-Lu is composed of six ORFs. The colored arrows are showing ORFs. The plus (+)-strand comprises virion-sense (encapsidated strand) ORFs (V1 and V2). The negative (−)-strand is composed of complementary-sense strand ORFs (C1, C2, C3 and C4). Replication of viral genome and transcription of the viral coding strands are governed by a non-coding large intergenic region (LIR).

The CLCuKoV-Lu genome possesses six overlapping open reading frames (ORFs) that are shown to encode four major proteins. Using the miRBase web-based tool for miRNA genomics, and 3 distinct prediction algorithms (miRanda, RNA22 and psRNATarget), we search out cotton genome encoded miRNAs possess the capacity to target the CLCuKoV-Lu genome. For this we first accessed the newly emerged novel CLCuKoV-Lu genome from GenBank and the experimentally verified mature cotton locus-derived ghr-miRNAs were retrieved from the miRBase (version 22) to evaluate the binding strength of miRNA-mRNA target interactions. CLCuKoV-Lu genomic sequence potentially targeted by cotton miRNAs was searched by the miRanda algorithm which predicted 11 miRNA-mRNA target pairings. RNA22 predicted effective target binding sites of 11 cotton miRNAs at 11 genomic loci in the CLCuKoV-Lu genome. psRNATarget predicted cleavable target candidates: 20 cotton miRNAs and 26 loci. RNAhybrid 76 miRNA-mRNA target pairs (Figure 2) (Table S3-S4 and File S1).

Figure 2. The Venn diagram showing miRNA-mRNA target pairs in the CLCuKoV-Lu genome. Four widely used computational databases (miRanda, RNA22, psRNATarget and RNAhybrid) were mined to estimate host-virus interaction pairs. Binding sites are key feature to determine degree of overlap exhibited by in silico tools used in this study. The intersection of four computational tools’ graph concludes four common ghr-miRNAs.

Figure 2. The Venn diagram showing miRNA-mRNA target pairs in the CLCuKoV-Lu genome. Four widely used computational databases (miRanda, RNA22, psRNATarget and RNAhybrid) were mined to estimate host-virus interaction pairs. Binding sites are key feature to determine degree of overlap exhibited by in silico tools used in this study. The intersection of four computational tools’ graph concludes four common ghr-miRNAs.

3.2. V1 encoding coat protein (CP)

The begomoviral V1 ORF (291-1061 (770 nt) encodes the coat protein (CP), which is required for encapsidation of the begomoviral ssDNA genome into virions, vector-mediated transmission, and virus movement [68,69,70]. V1 was targeted by four cotton locus-derived ghr-miRNAs: ghr-miR7486 (a, b) (locus 846), ghr-miR7497 (locus 349) and ghr-miR7506 (locus 509), as indicated by miRanda (Figure 3A). RNA22 identified two miRNAs: ghr-miR169a (locus 691) and ghr-miR7512 (locus 917) (Figure 3B). The psRNATarget algorithm predicted ten miRNAs: ghr-miR827 (a, b and c) (locus 740), ghr-miR3476-5p (locus 765), ghr-miR7492 (a, b and c) (locus 901), ghr-miR7497 (locus 459 ghr-miR7500 (locus 674) and ghr-miR7510a (locus 805) (Figure 3C). RNAhybrid identified nine ghr-miRNAs: ghr-miR393, ghr-miR482 (a, b), ghr-miR7486 (a, b), ghr-miR7490, ghr-miR7504a, ghr-miR7510a and ghr-miR7512 at nucleotide positions 611, 581, 849, 670, 694 and 917, respectively (Figure 3D and Table 2) and (Table S3-S4 and File S1).

3.3. V2 encoding pre-coat protein

V2 ORF (131-487 nt) is composed of 356 nucleotides encoding a pre-coat protein. The V2 protein of the viral genome is involved in symptom development, movement, transmission and regulation [70,71,72,73]. The ghr-miR7497 targeted the V2 gene at (locus 349) and (loci 155, 459), as identified by the miRanda and psRNATarget algorithms, respectively (Figure 3A and 3C). RNAhybrid identified six ghr-miRNAs to target overlapping region of V1 and V2 ORF: ghr-miR164, ghr-miR479, ghr-miR3476-5p, ghr-miR7497, ghr-miR7498 and ghr-miR7507 at nucleotide positions 375, 371, 421, 362, 370 and 364, respectively (Figure 3D and Table 2) and (Table S3-S4 and File S1).

3.4. C1 encoding replication-associated protein (Rep)

The C1 ORF (1505-2581 (1076 bases) encodes a replication-associated protein (Rep), that is essential for ssDNA replication, transcription, regulation [18,74,75,76]. miRanda predicted six miRNAs: ghr-miR390 (a, b and c) (locus 2278), ghr-miR7486 (a, b) (locus 2488) and ghr-miR7503 (locus 2214) (Figure 3A). C1 was also targeted by six miRNAs: ghr-miR390 (a, b and c) (locus 2281), ghr-miR393 (locus 1735) and ghr-miR7486 (a, b) (locus 2488) by RNA22 (Figure 3B). Five potential miRNA candidates were detected for silencing the C1 gene by psRNATarget: ghr-miR7486 (a, b) (locus 2488), ghr-miR7505 (locus 2049), ghr-miR7510b (locus 1581) and ghr-miR7513 (locus 2540) (Figure 3C). The RNAhybrid algorithm predicted twenty one ghr-miRNAs: ghr-miR156 (a, b, c, d), ghr-miR162a, ghr-miR166b, ghr-miR169a, ghr-miR398, ghr-miR827 (a, b, c), ghr-miR2949-3p, ghr-miR3476-3p, ghr-miR7491, ghr-miR7492 (a, b, c), ghr-miR7500, ghr-miR7501, ghr-miR7505, and ghr-miR7506 at nucleotide positions 2500, 1665, 2038, 1820, 2469, 2473, 1655, 2470, 1743, 2566, 2467, 2428, 2036 and 1795, respectively (Figure 3D and Table 2) and (Table S3-S4 and File S1).

3.5. C2 encoding transcription activator protein (TrAP)

The C2 ORF of begomoviruses (1153-1599) (446 nt) encoded transcriptional activator protein (TrAP) that is essential for symptom development in infected plants [15,77,78,79]. Among the targeted genes of CLCuKoV-Lu, C2 has few binding sites of cotton miRNAs. miRanda predicted hybridization of ghr-miR7513 at locus 1350 in the overlapping region of C2 and C3 ORFs (Figure 3A). RNA22 predicted two miRNAs: ghr-miR7489 (locus 1408) and ghr-miR7513 (2540) in the C2 and C3 ORFs (Figure 3B). Three miRNA were predicted to have binding affinity with C2 by psRNATarget: ghr-miR396 (a, b) at locus 1250 in the C2 and C3 overlapping ORFs and ghr-miR7510b at locus 1581 in the C2 and C1 overlapping ORFs (Figure 3C). RNAhybrid predicted cotton ghr-miRNAs: ghr-miR394 (a, b), ghr-miR7504b in the C2 ORF (Figure 3D) and (Table 2, Table S3-S4 and File S1).

3.6. C3 encoding replication enhancer protein (REn)

The C3 ORF of begomoviruses (1058-1459) (401nt) encode a replication enhancer protein (REn) that is involved in replication of the viral ssDNA genome[80,81]. The psRNATarget algorithm identified five miRNAs in the C3 ORF: (locus 1250), ghr-miR7484 (a, b) (1081) and ghr-miR7492 (a, b and c) (1094). In addition psRNATarget predicted in the overlapping region of C2 and C3: ghr-miR396 (a, b) (Figure 3C). RNAhybrid predicted two cotton miRNAs targeting C3 ORF: ghr-miR7484 (a, b) at nucleotide position 1077. Further, RNAhybrid identified cotton miRNAs in the C2 and C3 overlapping region: ghr-miR167 (a, b), ghr-miR396 (a, b), ghr-miR2949 (a-5p, b, c), ghr-miR7489, ghr-miR7493, ghr-miR7494, ghr-miR7511 and ghr-miR7513 at genomic positions 1400, 1226, 1242, 1447, 1352, 1399 and 1350 respectively (Figure 4D) and (Table 2, Table S3-S4 and File S1).

3.7. C4 encoding transcription regulator protein

The C4 ORF of begomoviruses (2091-2429) (338 bases) encodes a transcription regulator protein that work as viral effector [76,82,83]. miRanda and RNA22 predicted binding of ghr-miR390 (a, b and c) at consensus genomic locus 2281. The ghr-miR7503 targeted the C1 gene to have binding affinity at genomic locus 2214, as indicated by miRanda (Figure 3A-3B). RNAhybrid predicted cotton miRNAs in the C4 and C1 overlapping region: ghr-miR160, ghr-miR172, ghr-miR390 (a, b and c), ghr-miR399d, ghr-miR7488, ghr-miR7495 (a, b), ghr-miR7503, ghr-miR7508, ghr-miR7509 and ghr-miR7510b at nucleotide positions 2044, 2177, 2196, 2047, 2258, 2276, 2214, 2282, 2162 and 2175, respectively (Figure 3D) and (Table 2, Table S3-S4 and File S1).

3.8. Large Intergenic Region

Large intergenic region (LIR) drives the transcriptional regulation of the V1 and C1 ORFs of CLCuKoV [20,21,84]. Four distinct prediction algorithms (miRanda, RNA22, psRNATarget and RNAhybrid) predicted hybridization binding site of ghr-miR2950 at consensus genomic locus 82 of CLCuKoV targeting LIR (Figure 3 A-D). In addition, LIR was targeted by three miRNAs: ghr-miR7484 (a, b) and ghr-miR7497, as indicated by psRNATarget (Figure 3C). RNAhybrid predicted cotton ghr-miRNAs in the LIR: ghr-miR399 (a, b, c), ghr-miR2948-5p, ghr-miR2950, at nucleotide positions 2633, 2616, 2611, and 82, respectively (Figure 3D) and (Table 2, Table S3-S4, File S1).

3.9. Predicting Common Cotton miRNAs

Among all the targeting miRNAs of cotton to silence the CLCuKoV-Lu genome, four miRNAs (ghr-miR2950, ghr-miR7486 (a, b) and ghr-miR7513) were predicted by all the prediction tools employed in this study (Figure 2-4) and (Table S3 and File S1).

3.10. Predicting Consensual Cotton miRNAs

Of the 80 targeting mature G.hirsutum locus-derived ghr-miRNAs investigated, 7 G.hirsutum ghr-miRNAs: ghr-miR390 (a, b, c), ghr-miR7484 (a, b) and ghr-miR7503, ghr-miR7512 at nucleotide positions 2281, 1081, 2214 and 917, respectively were predicted to have potential binding sites within the CLCuKoV-Lu genome based on consensus genomic loci by at least two online miRNA prediction tools ( Table 2- 4). Out of 80 cotton miRNAs, three consensual G. hirsutum ghr-miRNAs: ghr-miR7486 (a, b) and ghr-miR7513 were predicted to have consensus genomic binding sites at nucleotide positions 2488 and 1350, respectively, was identified by union of consensus between the multiple algorithms used in this study. In this current study, only one G.hirsutum ghr-miRNA (ghr-miR2950) was predicted to have target binding site at common genomic position 82 by all the algorithms used (Table 2- 4, Figure 5). In terms of the CLCuKoV-Lu genome, ghr-miR2950 was predicted to target the non-coding LIR, ghr-miR2488 (a, b) were predicted to target the coding region C1 gene while ghr-miR7513 targets the overlapping region of ORFs C2/C3 (Figure 5 and Table 3).

Out of eleven consensual G. hirsutum ghr-miRNAs investigated, only one ghr-miRNA of G.hirsutum (ghr-miR2950 at nucleotide position 78-97), with a MFE of −27.38 Kcal/mol, was detected as top effective candidate (Table 4-5). The ‘cleavage’ efficacy of the ghr-miR2950 was verified against CLCuKoV-Lu by the RNAi-mediated suppression as concluded by Brodersen [85].

Figure 5. Intersection plot of predicted consensus binding sites of cotton ghr-miRNAs by multiple algorithms. Binding sites were predicted at least three algorithms at consensus genomic positions targeting specific genes of CLCuKoV-Lu.

Figure 5. Intersection plot of predicted consensus binding sites of cotton ghr-miRNAs by multiple algorithms. Binding sites were predicted at least three algorithms at consensus genomic positions targeting specific genes of CLCuKoV-Lu.

Table 3. Target binding sites of predicted consensus cotton ghr-miRNAs in the CLCuKoV-Lu genome were detected by different algorithms. .

Cotton miRNA	Target Site miRanda	Target Site RNA22	Target Site psRNATarget	Target Site RNAhybrid	MFE * miRanda	MFE ** RNA22	Expectation psRNATarget	MFE * RNAhybrid
ghr-miR390 (a, b, c)	2278	2281			−21.48	−18.00
ghr–miR2950	78	78	78	82	−27.38	−23.70	6.5	−30.20
ghr–miR7484 (a, b)			1081	1077			6.5	−20.90
ghr-miR7486 (a, b)	2488/846	2488	2488	849	−23.15/−29.28	−21.48	5.0	−30.70
ghr-miR7503	2214			2214	−23.35			−27.00
ghr-miR7512		917		917		−16.70		−23.50
ghr-miR7513	1350	1350		1351	−21.45	−17.50		−26.80

Table 3. Target binding sites of predicted consensus cotton ghr-miRNAs in the CLCuKoV-Lu genome were detected by different algorithms. .

Cotton miRNA	Target Site miRanda	Target Site RNA22	Target Site psRNATarget	Target Site RNAhybrid	MFE * miRanda	MFE ** RNA22	Expectation psRNATarget	MFE * RNAhybrid
ghr-miR390 (a, b, c)	2278	2281			−21.48	−18.00
ghr–miR2950	78	78	78	82	−27.38	−23.70	6.5	−30.20
ghr–miR7484 (a, b)			1081	1077			6.5	−20.90
ghr-miR7486 (a, b)	2488/846	2488	2488	849	−23.15/−29.28	−21.48	5.0	−30.70
ghr-miR7503	2214			2214	−23.35			−27.00
ghr-miR7512		917		917		−16.70		−23.50
ghr-miR7513	1350	1350		1351	−21.45	−17.50		−26.80

Table 4. Binding sites of predicted consensus cotton ghr-miRNAs targeting different genes.

miRNA ID	Accession ID	Mature Sequence (5′–3′)	Target Genes ORF(s)	Target Binding Locus Position
ghr-miR390a	MIMAT0005815	AAGCUCAGGAGGGAUAGCGCC	C1/C4	2278–2298
ghr-miR390b	MIMAT0005816	AAGCUCAGGAGGGAUAGCGCC	C1/C4	2278–2298
ghr-miR390c	MIMAT0005817	AAGCUCAGGAGGGAUAGCGCC	C1/C4	2278–2298
ghr-miR2950	MIMAT0014348	UGGUGUGCAGGGGGUGGAAUA	LIR	78–97
ghr-miR7484a	MIMAT0029124	UUUGUAUAUUAGAUCAAAGAGCAA	C3	1081–1105
ghr-miR7484b	MIMAT0029125	UUUGUAUAUUAGAUCAAAGAGCAA	C3	1081–1105
ghr-miR7486a	MIMAT0029127	AAGGAAGCGCUUUGUCCACGUGGA	C1/V1	2488–2510/846-871
ghr-miR7486b	MIMAT0029128	AAGGAAGCGCUUUGUCCACGUGGA	C1/V1	2488–2510/871
ghr-miR7503	MIMAT0029150	AGAUCGAUGGCUGAACAAGUUAGA	C4/C1	2214–2237
ghr-miR7512	MIMAT0029161	UGCUACUUGUAGUUAUGCAUG	V1	917–938
ghr-miR7513	MIMAT0029162	AAUCAGCCAGGAAUCGUUUGA	C2/C3	1350–1372

Table 4. Binding sites of predicted consensus cotton ghr-miRNAs targeting different genes.

miRNA ID	Accession ID	Mature Sequence (5′–3′)	Target Genes ORF(s)	Target Binding Locus Position
ghr-miR390a	MIMAT0005815	AAGCUCAGGAGGGAUAGCGCC	C1/C4	2278–2298
ghr-miR390b	MIMAT0005816	AAGCUCAGGAGGGAUAGCGCC	C1/C4	2278–2298
ghr-miR390c	MIMAT0005817	AAGCUCAGGAGGGAUAGCGCC	C1/C4	2278–2298
ghr-miR2950	MIMAT0014348	UGGUGUGCAGGGGGUGGAAUA	LIR	78–97
ghr-miR7484a	MIMAT0029124	UUUGUAUAUUAGAUCAAAGAGCAA	C3	1081–1105
ghr-miR7484b	MIMAT0029125	UUUGUAUAUUAGAUCAAAGAGCAA	C3	1081–1105
ghr-miR7486a	MIMAT0029127	AAGGAAGCGCUUUGUCCACGUGGA	C1/V1	2488–2510/846-871
ghr-miR7486b	MIMAT0029128	AAGGAAGCGCUUUGUCCACGUGGA	C1/V1	2488–2510/871
ghr-miR7503	MIMAT0029150	AGAUCGAUGGCUGAACAAGUUAGA	C4/C1	2214–2237
ghr-miR7512	MIMAT0029161	UGCUACUUGUAGUUAUGCAUG	V1	917–938
ghr-miR7513	MIMAT0029162	AAUCAGCCAGGAAUCGUUUGA	C2/C3	1350–1372

Table 5. Features of predicted consensus G.hirsutum ghr-miRNA-target pairs.

Cassava miRNA	miRNA-Target Pair	Locus Position		MFE (Kcal/mol)	Score	Complementarity (%)	Mode of Inhibition
ghr-miR2950	Query: 3' auaagGUGGGGGACGUGUGGu 5' | : | : | | | | | | | | | : Ref: 5' aataaCGCTCCC-GCACACTa 3'		78-97	−27.38	142	93.33	Cleavage
ghr-miR7486 (a, b)	Query: 3' aggUGCACCUGUUUCGCGAAGGAa 5' | : | | | : | | | | : | | | | | | | Ref: 5' tgaATTTGGG-AAAGTGCTTCCTc 3'		2488-2510	−23.15	171	90.00	Cleavage
ghr-miR7513	Query: 3' agUUUGCUAA--GGACCGACUAa 5' : : | | | : | | | | | | | | | | | Ref: 5' atGGACGGTTGACGTGGCTGATg 3'		1350-1372	−21.45	162	85.00	Cleavage

Table 5. Features of predicted consensus G.hirsutum ghr-miRNA-target pairs.

Cassava miRNA	miRNA-Target Pair	Locus Position		MFE (Kcal/mol)	Score	Complementarity (%)	Mode of Inhibition
ghr-miR2950	Query: 3' auaagGUGGGGGACGUGUGGu 5' | : | : | | | | | | | | | : Ref: 5' aataaCGCTCCC-GCACACTa 3'		78-97	−27.38	142	93.33	Cleavage
ghr-miR7486 (a, b)	Query: 3' aggUGCACCUGUUUCGCGAAGGAa 5' | : | | | : | | | | : | | | | | | | Ref: 5' tgaATTTGGG-AAAGTGCTTCCTc 3'		2488-2510	−23.15	171	90.00	Cleavage
ghr-miR7513	Query: 3' agUUUGCUAA--GGACCGACUAa 5' : : | | | : | | | | | | | | | | | Ref: 5' atGGACGGTTGACGTGGCTGATg 3'		1350-1372	−21.45	162	85.00	Cleavage

3.10.1. Visualization of miRNA Target

In order to integrate miRNA-mRNA target gene interaction analysis, we constructed Circos plots to identify novel host miRNA targets. The mapped mature cotton locus-derived ghr-miRNAs are depicted in the CLCuKoV-Lu genome.

3.10.2. Secondary Structure Analysis

In this study, we validated the efficacy of predicted consensus miRNAs by producing their secondary structures. The secondary structures were developed using manually curated cotton precursors. The MFE is a standard feature for stability of the map of the CLCuKoV-Lu RNA structurome. The current study validated and determined characteristic features of seven consensus precursors (Table 4).

Table 6. Characteristics features of precursor of cotton locus-derived ghr-miRNAs were predicted.

miRNA ID	Accession ID	Length Precursor	MFE*/Kcal/mol	AMFE**	MFEI***	(G+C)%
ghr-MIR2950	MI0013555	108 nt	−48.10	−44.53	−1.002	44.44
ghr-MIR7486a	MI0024169	105 nt	−81.90	−78.00	−1.436	54.29
ghr-MIR7486b	MI0024170	101 nt	−69.50	−68.81	−1.336	51.49
ghr-MIR7513	MI0024204	103 nt	−36.70	−35.63	−0.965	36.89

Table 6. Characteristics features of precursor of cotton locus-derived ghr-miRNAs were predicted.

miRNA ID	Accession ID	Length Precursor	MFE*/Kcal/mol	AMFE**	MFEI***	(G+C)%
ghr-MIR2950	MI0013555	108 nt	−48.10	−44.53	−1.002	44.44
ghr-MIR7486a	MI0024169	105 nt	−81.90	−78.00	−1.436	54.29
ghr-MIR7486b	MI0024170	101 nt	−69.50	−68.81	−1.336	51.49
ghr-MIR7513	MI0024204	103 nt	−36.70	−35.63	−0.965	36.89

3.10.3. Free Energy (ΔG) Computation

In order to validate the predicted miRNAs we evaluated the Free Energy (ΔG) of duplex. The free energies (ΔG) of four seven consensual cotton locus-derived ghr-miRNAs were computed (Table 7).

3.10.4. Conserved Genomic Binding Sites Analysis

Among all the predicted genomic binding sites, highest level of conservation was obserevd in ghr-miR2950 (78-97) in different strains of CLCuKoV-Lu (Figure 6).

Figure 6. A Circos plot showing cotton locus-derived ghr-miRNAs targeting CLCuKoV-Lu ORFs. CLCuKoV-Lu ORFs are represented with colored lines. (A) Showing interaction Circos map by seed-based algorithm miRanda. (B) Representing interaction Circos map by non-seed based algorithm psRNATarget.

Figure 6. A Circos plot showing cotton locus-derived ghr-miRNAs targeting CLCuKoV-Lu ORFs. CLCuKoV-Lu ORFs are represented with colored lines. (A) Showing interaction Circos map by seed-based algorithm miRanda. (B) Representing interaction Circos map by non-seed based algorithm psRNATarget.

4. Discussion

CLCuD is a multifaceted disease. It has very complex molecular mechanism of development and progression. The current study has identified, by computational approach, and validated, by different algorithms, a combination of 3 cotton miRNAs in 80 mature miRNAs with potential to target CLCuKoV-Lu genome. For the last three decades, as a primitive virus, CLCuKoV can severely affect cotton production in Pakistan [9,11,12,86,87]. CLCuKoV-Bu infection triggered widespread gene silencing as an adaptive defence using cotton miRNAs. We filtered the false positive prediction data and considered to validate using the computational algorithms. We presented an equitable and integrative approach for the evaluation of predicted miRNA biological data at three different levels. Performance of data-driven algorithms is validated at union of intersections levels for interpretable predicted biological data. Assessing the potential cotton miRNA–mRNA targets interactions could be consensually validated (Figure 5). The present study predicts mature cotton genome encoded miRNAs (ghr-miR2950, ghr-miR7486 (a, b) targeting CLCuKoV-Lu to develop resistance in cotton cultivar. The predicted miRNAs have their interactions with the C1 and LIR of CLCuKoV-Lu. Our data indicated that evolutionary conserved cotton miRNA (ghr-miR2950) was selectively employed by CLCuKoV-Lu. Previous studies has reported host-virus interaction using online computational tools to identify binding affinity of genome encoded miRNAs in RTV1 [88], SCBV [89], SCYLV [90], ZYMV [91], SCBGAV [92], ToBRFV and PhCMoV [93], RYMV [94], MCMV [95] and ICMV-Ker[96]. We have reported similar findings using online computational tools for best target prediction against sugarcane and rubber tree viruses [88,89,90,92]. In this study, we employed in-silico algorithms for computational prediction, i.e. consensual target binding sites of ghr-miR2950 at locus (locus 78), ghr-miR7484 (a, b) at locus (1081) and ghr-miR7486 (a, b) at locus (2488), while TAPIR was predicted no binding site. Host delivered plant miRNAs are primary source of inducing the degradation in the viral target and are responsible for typical pattern of base pairing. This study demonstrated that CLCuKoV-Lu genomic components (C1 and LIR) are highly susceptible to be targeted by consensual miRNAs. Among the 80 cotton miRNAs investigated, ghr-miR2950 was predicted to have consensus genomic binding site within the LIR of the CLCuKoV-Lu genome (Figure 7). LIR is the key component of CLCuKoV genome to govern the bidirectional mode of transcription of C1 and V1 genes and works as a bidirectional promoter [20,21].

Union and intersection approaches of prediction are important to control false-positive prediction. Union level of prediction approach relies on combination of target prediction algorithms to find true targets. In this case, sensitivity level of predicted data was increased at the cost of lower level of specificity. Whereas an intersectional approach was predicted a different mode and combines two or more algorithmic tools at specific level at the cost of lower sensitivity [97,98,99,100]. Here, we present overview computational methods for miRNA target interaction prediction achieved the best outcomes (Figure 2 and 5). Previous studies on plant host-delivered miRNAs had also predicted the silencing viral genome using list of in silico tools. Experimental evaluation of potential plant host genome encoded miRNAs targeting different plant viruses have been reported successfully [31,51,101].

This study was designed an equal computational approach to predict cotton miRNAs in the CLCuKoV-Lu genome to combat of Begomovirus in cotton cultivars. Application of RNAi discovery to varietal improvement in cotton to control CLCuD infection is the highly recommended strategy to decrease yield losses [102,103,104,105]. However, gene pyramiding for enhanced resistance to CLCuKoV-Lu for sustainable upland cotton is highly complicated due to the complex allotetraploid nature of the cotton genome. Lower regeneration efficiency of cotton callus is another constraint for the development of CLCuD-resistant allotetraploid upland cotton. The differential expression profile of ghr-miR2950 was reported against at early stages of V.dahlliae infection [48,106,107]. The cotton miRNA (ghr-miR2950) exhibited a differential expression in PHYA1 RNAi cotton [108]. The ghr-miR2950 is a key role to encoding a gibberellin 3 hydroxylase [109]. This enzyme was experimentally validated to accumulate in fibers at high level and is responsible for fiber cell elongation using GA signaling in the PHYA1 RNAi cotton [108,110]. The ghr-miR2950 was characterized and identified during the growth and development of ovule and fiber in cotton. The ghr-miR2950 is also involved in inducing Root-knot nematodes (RKN) infection [111,112]. RNAi is a technique used to screen host-delivered factors for identifying various cellular functions against viruses [113,114,115]. In this study, we employed 80 experimentally validated mature cotton miRNAs with annotated targets in the CLCuKoV-Lu genome. In the current study, we cover bioinformatics workflows for CLCuKoV-Lu genome silencing provides evidence for the generation of antiviral agents. The design, construction and validation of the amiRNA-based construct was presented to harbor a modified miRNA/miRNA* duplex of the precursor (ghr-MIR-2950) (Figure 8). Furthermore, we summarized current knowledge to demonstrate a novel methodology to minimize the antiviral effects of host genome encoded miRNAs against CLCuKoV-Lu.

5. Conclusions and Recommendations

CLCuKoV-Lu is an emerging pathogen associated with the CLCuD pandemic. CLCuD reduces yield and quality of all cotton cultivars in Pakistan. Overall this study provides an optimized prediction tools and parameters to identify the best-candidate miRNAs against CLCuKoV-Lu begomovirus. Prior to molecular cloning, we employed in silico tools and approaches for predicting effective binding affinity of mature candidate cotton miRNAs in the CLCuKoV-Lu genome. Among the 80 cotton miRNAs investigated, we predicted 3 consensual cotton locus-derived ghr-miRNAs that have miRNA-mRNA hybridization pairings of allotetraploid upland cotton miRNAs within the CLCuKoV-Lu. We identified ghr-miR2950 (out of 80 miRNAs) as the miRNA with highest affinity for the CLCuKoV-Lu genome. By way of comparison to other CLCuKoV strains, we have also established a reliable microRNA-binding region in the CLCuKoV genome that we predicted in course of an analysis of multiple sequence alignment. Hence, our efforts in mapping the cotton miRNA–mRNA target interaction may be aided in untangling molecular underpinnings of hereditary and CLCuD.