Seedlings of the resistant cultivar G. barbadense cv. Hai7124 and six cotton species (Gossypium anomalum, Gossypium trilobum, Gossypium aridum, Gossypium davidsonii, Gossypium thurberi, and Gossypium hirsutum) were planted in a 12-cm pot with sterilized soil at 28°C with a 14/10 light-dark photoperiod for 2 weeks, as previously described (Chen and Dai, 2010). The highly virulent V. dahliae strain Vd991 was cultured on potato dextrose agar (PDA) medium at 25°C for 10 days, the conidia were then washed with sterile water and the concentration was adjusted to 2 × 10⁷ conidia/mL. Six 2-week-old seedlings were each inoculated with 5 mL of conidial suspension by a soil drenching method. To prepare the samples for detection of GbaVd1 and GbaVd2 expression, the root tissues were collected 2, 6, 12, 24, and 48 h after inoculation.

RNA samples from inoculated plants were used for gene cloning. The RNA was extracted by using the cetyltrimethylammonium bromide (CTAB) method (Bekesiova et al., 1999). After tracing DNA in the RNA samples was removed by treatment with DNase I, cDNA was synthesized using a RevertAid™ First Strand cDNA Synthesis Kit from MBI (Fermentas, Glen Burnie, Maryland, USA). The primers specific to the cotton RLP gene homologs, GbaVd1 and GbaVd2 (cloned in our lab by homology-based cloning, the accession numbers is GU299533 and GU299534, respectively), were designed and used to amplify the cDNA samples. Primers used for gene cloning are listed in Table S1.

The open reading frames (ORFs) of GbaVd1 and GbaVd2 were determined by using ORF Finder (NCBI), and the protein sequences were deduced on the basis of gene sequences. The identities between GbaVd1 and GbaVd2 were determined by using the Blast program (http://www.ncbi.nlm.nih.gov/Blast.cgi) and Vector NTI software (Lu and Moriyama, 2004). Signal peptides and membrane-spanning structures were predicted by using SignalP 4.0 (Petersen et al., 2011) and TMHMM 2.0 (Sonnhammer et al., 1998), respectively. The primary protein structure was determined by protein sequence alignment with known RLPs, including Cf-2.2 (U42445), Cf-5 (AF053993), Cf-4 (AJ002235), Cf-9 (AJ002236), HcrVf1 (AJ297739), HcrVf2 (AJ297740), HcrVf3 (AJ297741), LeEix1 (AY359965), LeEix2 (AY359966), Ve1 (AF365929), and Ve2 (AF272366). ClustalX 1.83 software was used for multiple sequence alignment (Thompson et al., 1997), and data were exported with Boxshade 3.21 (http://www.ch.embnet.org/software/BOX_form.html, written by K. Hofmann and M. Baron). Phylogenetic trees were constructed in Mega 6.0 with the Maximum Parsimony method, using the Jones-Taylor-Thornton (JTT) model and performing 1000 bootstrap replicates (Tamura et al., 2013).

For the VIGS assays, GbaVd1 and GbaVd2 were integrated into a vector and introduced into A. tumefaciens GV3101. Agrobacterium the strains harboring the pTRV2-GbaVd1/GbaVd2 plasmid combined with strains harboring the pTRV1 vector were mixed in a 1:1 ratio and co-infiltrated into the cotyledon leaves of 2-week-old cotton plants. The effectiveness of the VIGS assay was evaluated by using the control cotton CLA1 gene as previously described (Gao and Shan, 2013). Approximately 14 days after the VIGS procedure, a visible phenotype of white-colored leaves were observed in plants in which the CLA1 gene had been targeted by VIGS, and all of the plantlets were subjected to V. dahliae inoculation with 5 mL of conidial suspension (2 × 10⁷ conidia/mL). The Verticillium wilt symptoms were investigated 3 weeks after inoculation. Fungal biomass in cotton were determined using a method described previously (Santhanam et al., 2013). qPCR was performed using a qPCR SYBR premix Ex Taq II kit (TaKaRa, Japan) with primers specific to the cotton 18S gene and V. dahliae elongation factor 1-α (EF-1α) (Table S1).

According to the sequences of GbaVd1 and GbaVd2, the ORF fragments were amplified with primers containing BamH I and BstE II enzyme sites. The fragments of GbaVd1 and GbaVd2 were integrated into the binary vector pCAMBIA1303 under control of the cauliflower mosaic virus (CaMV) 35S promoter. The recombinant vectors were transferred into A. tumefaciens strain LBA4404, and this was followed by genetic transformation of A. thaliana ecotype Col-0 via the floral dip transformation method (Clough and Bent, 1998). Homozygous transgenic Arabidopsis (T₃) plants were screened and used for this study. Transgenic plants were selected on MS medium containing 50 mg/L hygromycin and were identified by RT-PCR. The RT-PCR conditions consisted of an initial 94°C denaturation step for 10 min, which was followed by 35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s, and ubiquitin 4 (Ubi4, NM_122069.3) was used as a control. Two random transgenic plant lines carrying GbaVd1 and GbaVd2 were selected for the Verticillium wilt resistance assay. The inoculation method was performed as in cotton. In brief, 4-week-old seedlings were inoculated with 5 mL of conidial suspension (5 × 10⁶ conidia/ml). The inflorescence height was detected 4 weeks after inoculation. Likewise, the development of fungal biomass in transgenic plants was detected by qPCR, with primers to RNA-binding family proteins in Arabidopsis (Table S1).

Microarray analysis was performed using Arabidopsis Gene Expression Microarray Version 4.0 (43,803 probes) (Agilent Technologies, USA). Three independent transgenic lines of GbaVd1 and GbaVd2 overexpression were used for microarray analysis. Two-week-old seedling transgenic plants without inoculation were collected for RNA extraction, and the wild type (Col-0) was used as a control. Total RNA was extracted from the root and labeled with Cy-3 using a Low RNA Input Linear Amplification/Labeling Kit, One-color (Agilent Technologies), hybridized, and washed according to the manufacturer's instructions. Hybridized microarrays were then scanned with an Agilent Microarray Scanner (Agilent Technologies). Statistical data extraction processes were performed on three biological replicates for each treatment according to the instructions (Agilent, GeneSpring GX9 software package). Fold changes in expression levels in response to each transgenic plant were compared with those of the respective wild-type, P-values were adjusted for multiple comparisons using the Bonferroni correction. Sequence homologies of the deduced amino acids were determined using the Arabidopsis Information Resource (TAIR) database (http://www.Arabidopsis.org/). Functional annotation of GO and metabolic pathways were performed using the GeneSpring GX9 package (Agilent Technologies, USA), and statistically for overrepresented GO terms. The significance of GO catalog for differentially expressed genes were identified using a Fisher's Exact Test (filtered with FDR ≤ 0.05) as previously described (Chen et al., 2016).

The cDNA samples were used for qRT-PCR analysis with SYBR Green (Invitrogen, Carlsbad, California, USA), and the PCR was performed using a thermocycler (ABI PRISM 7500, Applied Biosystems) with the following program: 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s for 40 cycles. The cotton 18S gene was used as an endogenous control. Relative transcript levels of GbaVd1 and GbaVd2 were determined using the 2^−ΔΔCT method (Livak and Schmittgen, 2001), with three independent determinations.

Four-week-old seedlings of Arabidopsis transgenic lines were inoculated with 5 mL of conidial suspension (5 × 10⁶ conidia/ml). The internode of shoot tissue between the second and third nodes was used for cell wall reinforcement detection at 4 weeks after inoculation. The samples were stained with phloroglucinol solution (3% w/v phloroglucinol in 95% ethanol) for 2 min after slice treatment and were then immersed in concentrated HCl for a few seconds and covered with glycerol (Nakano and Meshitsuka, 1992). The tissue slices were observed with a Leitz Diaplan light microscope and photographed immediately.

The full-length of GbaVd1 and GbaVd2 were cloned from the cDNA sample of G. barbadense (cv. Hai7124) inoculated with V. dahliae by PCR (Figure S1A). Sequences analysis showed that GbaVd1 and GbaVd2 are 3,262 and 3,466 bp in length, respectively, without any introns, and each has a single open reading frame (ORF) of 1,020 or 1,082 amino acids (Figure S1B), respectively. BlastP analysis revealed that both contain an LRR ribonuclease inhibitor domain, a potential signal peptide and a transmembrane region (Figure S2), and both were predicted to be located in the plasma membrane (the testk used for kNN was 14, and the Nearest Neighbors numbers of GbaVd1 and GbaVd2 related to the plasma membrane were 9 and 8, respectively). Sequence alignment indicated that the protein-coding sequences of GbaVd1 and GbaVd2 are similar to those of known plant RLP genes, with GbaVd2 showing the highest identities to Ve genes, but with GbaVd1 displaying the highest identities to HcrVf genes (Figure 1A). Likewise, the alignment of encoded protein sequences also confirmed that GbaVd1 and GbaVd2 are most homologous to Ve and HcrVf genes (Figure 1A), respectively. Sequence alignment showed that GbaVd2 is a homolog to Ve-like proteins from G. barbadense, and is nearly identical to the previously cloned and published Gbvdr5 (Yang et al., 2015) with only one conservative amino acid change (Figure S3). GbaVd1 protein sequence showed high identity to a RLP from cotton (Figure S4). In addition, six allelic variants of the GbaVd1 and GbaVd2 were isolated from six different cotton species, and two of them are susceptible to V. dahliae. A DNA sequence alignment showed a total of 125 and 264 single nucleotide polymorphisms (SNPs) among the allelic variants compared with the reference sequence of GbaVd1 and GbaVd2 (Figures S5A,C), resulting in 74 and 177 amino acid changes (Figures S5B,D), respectively. However, the association of the amino acid substitutions with the Verticillium wilt resistance of cotton species was not found among these allelic variants (Figures S5B,D). Phylogenetic analysis showed a close relationship between GbaVd1 and HcrVfs, which clustered into the same branch, and between GbaVd2 and Ve genes (Figure 1B). These results suggested that both GbaVd1 and GbaVd2 display homology to the known RLPs but also show sequence divergence when compared with each other.

Although the sequence identities between GbaVd1 or GbaVd2 and known RLPs were relatively low (Figure 1A), the alignment results showed that the amino acid residues involved in the LRR structure are conserved (Figure S6). Sequence alignments between GbaVd1 or GbaVd2 and known RLPs indicated that both contain six conserved domains, Domains A through F (Figure 2). Domain A consists of the signal peptide and its cleavage site at the N-terminus of the protein. Domain B belongs to the N-terminus of the mature protein. Domain C contains LRR structures, with 32 and 33 LRRs in GbaVd1 and GbaVd2, respectively. Domain C is usually divided into two sections, domains C1 and C2, which are linked by a small peptide, and the number of LRRs (four) in domain C2 is identical between GbaVd1 and GbaVd2. Domains B and C also contain 22 and 26 putative N-glycosylation sites Nx(S/T) in GbaVd1 and GbaVd2, respectively. Domain D is a negatively charged extracytoplasmic domain predicted to play a role in orienting and anchoring the protein to the cell membrane and links domain C with the transmembrane structure of domain E. Domain F is a small cytoplasmic peptide with positively charged residues and is linked to domain E. In addition, the cytoplasmic domain of GbaVd1 possesses a YXXø signal sequence (YFRI), and GbaVd2 contains a KKX motif (KKH), both of which probably participate in the endocytosis process (Figure 2). A sequence similarity search using GbaVd1 showed high identity to a protein annotated as a LRR receptor-like serine/threonine-protein kinase FLS2 from cotton (Figure S4). However, there is no kinase domain in the small cytoplasmic region, suggesting improper annotation of the FLS2 homolog and that GbaVd1 encodes a RLP (not an RLK). Therefore, similar to the known RLPs, GbaVd1 and GbaVd2 encode extracytoplasmic glycoproteins in which the majority of the extracytoplasmic domain consists of multiple LRRs.

To identify the roles of GbaVd1 and GbaVd2 in Verticillium wilt resistance, tobacco rattle virus (TRV)-based virus-induced gene silencing was used in the island cotton cv. Hai7124. The cotton gene cloroplastos alterados 1 (CLA1), which is essential for chloroplast development, was used as a positive control to monitor VIGS efficiency. Plants in which the GbaVd1 or GbaVd2 expression was silenced by TRV-based VIGS were challenged with V. dahliae after all of the CLA1-silenced cotton plants presented an albino phenotype on their newly emerged leaves. At 21 days post-inoculation with the V. dahliae strain V991, plants treated with TRV-based VIGS of GbaVd1 or GbaVd2 but not the wild-type plants presented clear Verticillium wilt symptoms of leaf chlorosis, including withering and dwarfing (Figure 3A), thus indicating that the resistant cotton cv. Hai7124 displayed compromised resistance to V. dahliae after silencing of GbaVd1 or GbaVd2 by VIGS. Real-time quantitative polymerase chain reaction (qPCR) quantification of gene expression related to fungal biomass development demonstrated significantly increased V. dahliae strain Vd991 propagation in GbaVd1- or GbaVd2-silenced plants compared with the wild-type and empty vector-treated plants (Figure 3B). Moreover, GbaVd1 and GbaVd2 transcript accumulation was quickly and strongly induced 2–48 h after inoculation with V. dahliae and reached the highest expression levels (>10-fold change) 12 h after inoculation (Figure 3C). Thus, it can be concluded that the functions of GbaVd1 and GbaVd2 are closely associated with Verticillium wilt resistance in cotton.

To identify the Verticillium wilt resistance functions of GbaVd1 and GbaVd2, both were interfamily-transferred into the Arabidopsis genome by using the Agrobacterium tumefaciens-mediated transformation method. Ten GbaVd1 and nine GbaVd2 independent Arabidopsis transgenic lines were obtained by hygromycin-resistance selection and verified by reverse transcription PCR (RT-PCR), and the integrated genes were successfully expressed, except in one GbaVd2 Arabidopsis transgenic line (Figure S7). Two random transgenic plant lines of GbaVd1 and GbaVd2 were selected for Verticillium wilt resistance identification with the root dip method and were assessed for Verticillium wilt resistance by evaluation of the extent of leaf chlorosis and inflorescence heights. Compared with the wild-type Col-0 ecotype, overexpressed GbaVd1 or GbaVd2 in Arabidopsis improved the resistance of the transgenic plants to V. dahliae strain Vd991, as indicated by significant decreases in leaf chlorosis, withering and dwarfing (Figure 4A). More importantly, unlike the wild-type Arabidopsis, which lost the capacity to set seed, the transgenic plants presented with a normal seed set 4 weeks after inoculation with V. dahliae (Figure 4A), a result strongly suggesting that overexpression of GbaVd1 or GbaVd2 conferred Verticillium wilt resistance to Arabidopsis. Further examination by qPCR confirmed that the V. dahliae significantly developed less fungal biomass in GbaVd1 or GbaVd2 transgenic plants than in the wild-type plants (Figure 4B). In addition, plant dwarfing, a typical Verticillium wilt symptom, was also significantly alleviated in transgenic lines after inoculation with V. dahliae (Figure 4C). Together, these results strongly suggested that GbaVd1 and GbaVd2 function as resistance genes and that interfamily transfer of these genes confers Verticillium wilt resistance in Arabidopsis.

To identify the defense responses mediated by GbaVd1 and GbaVd2, the transcriptome of transgenic lines were assayed by using an Agilent 44 K custom oligo microarray system with three independent biological replicates. In total, 4,153 and 4,484 probes were detected in GbaVd1 and GbaVd2 transgenic lines (Supplementary Data S1), respectively. Of these detected genes, 257 (93 genes up-regulated and 164 genes down-regulated) showed differential expression in the GbaVd1 transgenic line (fold-change ≥1.5, p < 0.01), and 362 (226 genes up-regulated and 136 genes down-regulated) were differentially expressed in the GbaVd2 transgenic line (fold-change ≥1.5, p < 0.01) (Figure 5A; Table S2). A further comparison of the two differentially expressed gene sets revealed 45 common genes between the GbaVd1 and GbaVd2 transgenic lines, of which 35 displayed the same regulation patterns (nine genes up-regulated and 26 genes down-regulated), and the other 10 of which showed down-regulation in the GbaVd1 transgenic line but up-regulation in the GbaVd2 transgenic line (Figure 5A; Table S2). These results indicated that the regulation patterns mediated by GbaVd1 and GbaVd2 display some of the same aspects, but that each has unique characteristics. Gene ontology (GO) analysis indicated that subsets of the differentially expressed genes activated by GbaVd1 and GbaVd2 were involved in the categories cellular metabolic process (statistically in overrepresented, filtered with FDR ≤ 0.05), regulation of cellular process, response to stress, and transcription regulation (Figure 5B). These results suggested that GbaVd1 and GbaVd2 have similar functions in mediating the defense response after overexpression in Arabidopsis.

Generally, the endocytosis process and signal transduction are necessary for the disease resistance mediated by RLPs. GbaVd1 and GbaVd2 contain a conserved motif that is probably involved in the endocytosis process (Figure 2). As expected, several genes associated with endocytosis were differentially expressed after GbaVd1 and GbaVd2 overexpression in Arabidopsis, including genes encoding ubiquitin ligase family proteins and heat shock protein 70 (HSP70) (Figure 6A; Table S3). Specifically, there are three genes encoding ubiquitin ligase family proteins, two of which were commonly activated in the GbaVd1 and GbaVd2 transgenic lines. These proteins facilitate formation of clathrin-coated vesicles during endocytosis (Figure 6A; Table S3). In addition, two HSP70 genes (AT3G12580 and AT5G02490) involved in the formation of early endosomes were also differentially expressed in the GbaVd1 transgenic line (Figure 6A; Table S3). These results suggested that the endocytosis process associated with plant-pathogen interactions was affected after GbaVd1 and GbaVd2 overexpression in Arabidopsis. In addition, several genes encoding signal factors involved in plant-pathogen interactions were also differentially expressed in the GbaVd1 and GbaVd2 transgenic lines (Figure 6B; Table S3). In the GbaVd1 transgenic line, the genes encoding two heat shock protein 90 (HSP90) isoforms and the transcription factor WRKY33 were up-regulated, but two other genes encoding calmodulin-like (CaM/CML) protein 12 and the transcription factor WRKY22 were activated in the GbaVd2 transgenic line (Figure 6B). Moreover, 27 and 39 genes encoding transcription factors were significantly differentially expressed after overexpression of GbaVd1 and GbaVd2, respectively; these included members of the AP2 domain-containing protein family and the myeloblastosis (Myb) family (Figures 6C,D; Table S3). However, most transcription factor genes were specifically differentially expressed in either the GbaVd1 or GbaVd2 transgenic lines, except for five common genes (Figures 6C,D; Table S3). Together, these results indicated that the interfamily transfer of cotton GbaVd1 and GbaVd2 provided the capacity to differentially expressed defense response-related genes and mediate Verticillium wilt resistance in Arabidopsis.

Plants usually enhance disease resistance by activation of defense responses that are mediated by the plant-pathogen interaction pathway (PATHWAY: map04626, KEGG database). In this study, the plant-pathogen interaction pathway was significantly activated after GbaVd1 or GbaVd2 overexpression in Arabidopsis, thus resulting in the induction of defense responses, including the hypersensitive response, phytoalexin accumulation, and cell wall reinforcement (Figures 6A,B). Microarray analysis confirmed that several differentially expressed genes that were significantly enriched were involved in the phytoalexin (phenylpropanoids, alkaloids, terpenoids, etc.) biosynthesis pathway or encode cell wall-related proteins (Figure 7A; Table S4).

Functional analysis of differentially expressed genes revealed that 9 and 13 genes encoding cell wall-related proteins were differentially expressed by overexpression of GbaVd1 or GbaVd2, respectively, in Arabidopsis (Figure 7A). Of these genes, the most encoded were nodulin family proteins, glycine/proline-rich proteins, and lectin family proteins, and were mainly up-regulated in transgenic lines (Figure 7B). Specifically, for the nodulin family, four and six genes were significantly activated in the GbaVd1 and GbaVd2 transgenic lines, respectively (Figure 7B). Moreover, the identification of cell wall structures confirmed that the degree of lignification was significantly enhanced by the transcription of genes encoding cell wall-related proteins (Figure 7C), thus indicating that the cell wall is reinforced after GbaVd1 or GbaVd2 overexpression in Arabidopsis. In addition, of these cell wall-related genes, only two genes were commonly regulated between the GbaVd1 and GbaVd2 transgenic lines (Figure 7B), thus suggesting that the cell wall-related proteins were differentially activated in the two transgenic lines.

Although both GbaVd1 and GbaVd2 encode RLPs, the sequences and primary structures display divergence (Figures 1, 2), thus suggesting that the functions of GbaVd1 and GbaVd2 might diverge during plant-pathogen interactions. Microarray analysis revealed that at least four layers of defense responses were differentially activated after GbaVd1 or GbaVd2 overexpression in Arabidopsis (Figures 6, 7). Further analysis of the differentially expressed genes revealed that the defense responses mediated by GbaVd1 and GbaVd2 displayed significant differences that only five genes encoding transcription factors and two genes encoding cell wall-related proteins were common between the two transgenic lines (Figure 8). The activation of genes involved in endocytosis appeared to be similar between the GbaVd1 and GbaVd2 transgenic lines, with two common genes out of six differentially expressed genes, but more genes were activated in the GbaVd1 transgenic plants (Figure 8). Together, these results suggested that GbaVd1 and GbaVd2 mediate different resistance signaling and activation of defense responses to confer Verticillium wilt resistance after interfamily overexpression in Arabidopsis.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.