Genome-Wide Characterization and Expression of Two-Component System Genes in Cytokinin-Regulated Gall Formation in Zizania latifolia

The thickening of Zizania latifolia shoots, referred to as gall formation, depends on infection with the fungal endophyte Ustilago esculenta. The swollen and juicy shoots are a popular vegetable in Asia. A key role for cytokinin action in this process was postulated. Here, trans-zeatin stimulated swelling in fungi-infected Z. latifolia. A two-component system (TCS) linked cytokinin binding to receptors with transcriptional regulation in the nucleus and played important roles in diverse biological processes. We characterized 69 TCS genes encoding for 25 histidine kinase/histidine-kinase-like (HK(L)) (21 HKs and 4 HKLs), 8 histidine phosphotransfer proteins (HP) (5 authentic and 3 pseudo), and 36 response regulators (RR; 14 type A, 14 type B, 2 type C, and 6 pseudo) in the genome of Z. latifolia. These TCS genes have a close phylogenetic relationship with their rice counterparts. Nineteen duplicated TCS gene pairs were found and the ratio of nonsynonymous to synonymous mutations indicated that a strong purifying selection acted on these duplicated genes, leading to few mutations during evolution. Finally, ZlCHK1, ZlRRA5, ZIRRA9, ZlRRA10, ZlPRR1, and ZlPHYA expression was associated with gall formation. Among them, ARR5, ARR9, and ZlPHYA are quickly induced by trans-zeatin, suggesting a role for cytokinin signaling in shoot swelling of Z. latifolia.


Introduction
Zizania latifolia belongs to the wild rice genus Zizania and was an important grain crop in ancient China. Around 2000 years ago, Z. latifolia was infected by the fungal endophyte Ustilago esculenta which causes swelling in shoot apical tissues, resulting in juicy galls. Nowadays, the swollen shoots of Z. latifolia are a popular vegetable in China (called "Jiaobai"), East Asian countries, and Japan [1,2].
In nature, plant tumors or "galls" can be induced by different pathogens such as bacteria, viruses, fungi, protists, and insects. Despite different pathogenic mechanism to induce galls, most of these plant pathogens demonstrate the ability to change the phytohormone level in the host plants, especially cytokinins. For example, the well-studied gall-forming Agrobacterium tumefaciens carries the Ti plasmid, which can integrate a DNA segment called T-DNA into the host's genome [3]. The T-DNA carries several genes including isopentenyl transferase (IPT), a key enzyme of cytokinin biosynthesis, which instigates increased cytokinin production in the transformed cells of the host [4].

Identification of TCS Genes in Z. latifolia
To identify members of the TCS family in Z. latifolia, HMM searches and BLASTP analysis were performed by employing 415 TCS protein sequences from Arabidopsis [14], rice [21], maize [25], soybean [24], Chinese cabbage [26], and tomato [27] as queries. After repetitive sequences were removed, all identified sequences were reserved and submitted to NCBI CDD, Pfam, and SMART to confirm their typical domains. A total of 69 non-redundant sequences, including 25 HK(L)s, 8 HPs, and 36 RRs, were identified in the Z. latifolia genome. All TCS members in Z. latifolia were named according to the homologous genes in Arabidopsis [27,28,32]. This nomenclature has been used in tomato [27], cucumber, and watermelon [28]. The TCS genes have been intensively studied in model plants and some important crops. The number of TCS genes in Z. latifolia (n = 69) is similar to that in tomato (n = 65) but larger than that of all reported species except soybean (n = 98) and Chinese cabbage (n = 85; Table S1).

HK Proteins in Z. latifolia
Sequence analysis of the entire Z. latifolia genome revealed 25 distinct putative histidine protein kinases genes, including 21 histidine protein kinases (ZlHKs) and 4 related proteins (ZlHKLs), based on whether their HK domains were conserved or divergent (Table 1). These histidine protein kinase homologs can be further divided into distinct subfamilies: Eight cytokinin receptors (CHKs), one CKI-like protein, seven ethylene receptors (ERSs and ETRs), five phytochrome photoreceptors (PHYs), and four pyruvate dehydrogenase kinases (PDKs; Table 1). A typical HK domain (corresponding to motifs 1 and 7 in Figure S2) has five conserved signature motifs-H, N, G1, F, and G2-with the conserved His as the key feature in the H motif. The other four motifs were defined as the nucleotide-binding cleft [33]. Through multiple-sequence alignment, all HKs and HKLs proteins were found to be conserved in the H motif with His residue, except ZlCHK2 and ZlCHK3 ( Figure S1A). All proteins had complete N motifs. However, the G1, F, and G2 motifs occurred only in ZlCHK3 and ZlCHK4, suggesting that other HK/HKL proteins may not have full HK function or that these motifs may not be necessary for HK function ( Figure S1A). Notably, only ZlCHK4 had a full HK domain, suggesting that it has a critical function in Z. latifolia. Table 1 and Figure S1B also show that both CHASE domains (motif 11 in Figure  S2) and 1-4 transmembrane sequences (TM) occurred in these cytokinin receptors, except for ZlCHK2 and ZlCHK3. Both CHASE and TM were demonstrated to be crucial for membrane-associated cytokinin recognizing and binding [17]. In addition, the receiver domain (motifs 9 and 13 in Figure S2) with Asp residue was conserved in all ZlCHKs except for ZlCHK7 ( Figure S1B). Two AHK2-like genes (ZlCHK1 and ZlCHK2) were 65% and 46% similar, respectively, to their Arabidopsis counterparts of AHK2, whereas two AHK3-like members (ZlCHK3 and ZlCHK4) were 61% and 63% similar, respectively, to Arabidopsis AHK3. Another four AHK4-like cytokinin receptor genes (ZlCHK5, ZlCHK6, ZlCHK7, and ZlCHK8) were highly similar (55%-72%) to Arabidopsis CRE1/AHK4/WOL (Table 1). AHK2, AHK3, and CRE1/AHK4/WOL are cytokinin receptors and can regulate cell division in Arabidopsis [15,[34][35][36][37][38][39][40][41]. Some of the cell cycle and meristem control master genes, such as CYCD3s and STM1, respond to cytokinins and function downstream of CRE1/AHK4/WOL to regulate shoot and root development [42].
The ethylene receptors, including ZlETR1-4 and ZlERS1-3, contain HK, HATPase, and cyclic GMP adenylyl cyclase FhlA (GAF) domains (Table 1). Although histidine-kinase activity can subtly modulate the ethylene response, no major role has yet been identified in ethylene signal transduction [43][44][45]. However, histidine kinase activity could allow cross-talk between ethylene perception and other TCS pathways such as cytokinin signal transduction. All five PHY subfamily members (ZlPHYA, B, C, D, E) possess PHY (chromophore-binding), GAF, and PAS (signal sensor) domains (Table 1). These domains were found to be essential for responding to red and far-red light signals during plant development in Arabidopsis [46]. The HATPase domain was identified in all ZlPDKs proteins (ZlPDK1-4). A newly published paper has revealed that PDK1 could regulate auxin transport and vascular development through phosphorylation of AGC1 kinase PAX in Arabidopsis [47]. To further understand the structural diversity of Z. latifolia HK genes, we analyzed the exon-intron organization and conserved protein motifs ( Figure S2). These analyses revealed that HK genes in the same group usually had a similar gene structure and motif composition, which strongly support the reliability of the phylogenetic classification.

HP Proteins in Z. latifolia
Eight ZlHP genes were identified: Five authentic HPs and three pseudo-HPs (Table 2). All ZlHPs were small proteins with 149-276 amino acids and had two conserved motifs (motifs 1 and 2 in Figure  S3) that encode the Hpt domain. Five authentic HPs (ZlHP1, ZlHP2, ZlHP4, ZlHP7, ZlHP8) contained conserved His (H) residues in the HP signature sequence of XHQXKGSSXS, whereas in the other three pseudo-HPs (ZlHP3, ZlHP5, and ZlHP6), the histidine phosphorylation site was replaced by Gln (Q) residue ( Figure S4). Notably, all HPs except ZlHP1 and ZlHP3 were localized in the nucleus (Table 2), which might be essential for their phosphorelay, during which HPs translocate from the cytoplasm to the nucleus [48,49]. In addition, all ZlHPs had a very similar gene structure containing 5-7 introns and motif compositions ( Figure S3).

RR Proteins in Z. latifolia
Thirty-six genes were identified as ZlRRs: 14 type A RRs, 14 type B RRs, 2 type C RRs, and 6 pseudo-RRs (Table 3). The type A RRs (ZlRRA1-14) were a group of small proteins with 76-269 amino acids. All type A RRs contained a receiver domain corresponding to three motifs (motifs 1, 3, and 4 in Figure S5) along with short N-and C-terminal extensions (Table 3 and Figure S5). Most type A RRs contained four introns. The 14 ZlRRAs were highly similar (52%-81%) to their Arabidopsis homologs (Table 3).
In addition to having a conserved receiver domain, all type B RRs (ZlRRB1-14) contained long C-terminal extensions with a Myb-like DNA binding domain (motif 2 in Figure S5), suggesting their function as transcription factors. Except for ZlRRB6, all type B RRs were predicted to be localized in the nucleus (Table 3). Most type B RRs had four or five introns ( Figure S5). Type B RRs were 30%-66% similar to their counterparts in Arabidopsis.
Different from type A RRs, two identified type C RRs (ZlRRC1, 2) contained a receiver domain (motifs 1 and 4 in Figure S5) but lacked a long C-terminal extension. Both proteins had three introns and were 40% and 41% similar to their Arabidopsis counterpart, ARR24, respectively. In addition, six proteins identified as pseudo-RRs (ZlPRR1-6) all had a conserved receiver domain (motifs 1 and 4 in Figure S5) and a CCT domain (motif 5), which was found to play an important role in regulating circadian rhythms [54]. These pseudo-RRs were 39%-54% similar to their counterparts in Arabidopsis. Except ZlPRR2, all these type C and pseudo-RRs were localized in the nucleus (Table 3, Figure S5). Figure S5 also shows that the ZlRRs genes in the same category usually had very similar gene structure and motif composition, thereby supporting the evolutionary conservation and reliability of the phylogenetic classification.

Phylogenetic Relationship of TCS Members in Z. latifolia
To investigate the evolutionary relationships of TCS genes, we used all amino acid sequences of 124 HK(L)s, 53 HPs, and 275 RRs from Arabidopsis [14], rice [21], maize [25], soybean [24], Chinese cabbage [26], tomato [27], and Z. latifolia to construct maximum likelihood (ML) trees [55]. On the basis of the bootstrap values and the topology of the tree, all HK(L)s proteins in the seven species were divided into seven distinct subfamilies, including cytokinin receptor, ethylene receptor, PHY-like, CKI1-like, CKI2/AHK5-like, AHK1-like, and PDK-like subfamilies (Figure 2). Similar phylogenetic structures were viewed in previous studies [14,17,21]. The HK(L)s in Z. latifolia usually have much closer relationships to those of rice and maize than other species, as they are all members of the monocotyledon. This indicates that the HK(L)s gene expansion occurred after the divergence of monocot and dicot plants. In addition, HK(L)s genes of monocots and eudicots showed an alternative distribution pattern in each subfamily, but no monocots occurred in the AHK1-like subgroup.  All HP proteins from the seven species were divided into three clades-clades I, II, and III ( Figure 3). Clade II contained HPs only from monocots, whereas all HPs in clade III were from dicots. These results suggest that the gene expansion events for HPs occurred after the divergence of monocots and dicots. Four authentic ZlHPs except ZlHP4 were grouped into clade II, whereas clade I mainly contained pseudo-HPs from both monocots and dicots. All ZlHPs showed the closest phylogenetic relationship to rice OsHpts (Figure 3).  The 275 RR proteins from the seven plant species were classified into type A, type B, type C, and pseudo-RR groups ( Figure 4). Generally, the ZlRRs in Z. latifolia were phylogenetically closer to rice OsRRs in all these groups. Phylogenetic analyses indicated that all type A RRs shared a close evolutionary relationship and an alternating distribution between monocots and dicots plants, suggesting that type A RR genes might already have expanded before the monocot-dicot divergence. Type B RRs from these species could be divided into three subgroups (I, II, and III), consistent with previous studies [17,26]. In detail, type B II subgroups contain RRs only from dicots of Arabidopsis, Chinese cabbage, and tomato; the RRs from monocots might have been lost during evolution. All type B ZlRRs from Z. latifolia except for ZlRRB3 and ZlRRB14 were classified into the type B I subfamily. A similar phenomenon was also viewed in soybean and cucumber [24,28]. All pseudo-RRs from seven plant species were subclassified into type B PRRs or clock PRRs (Figure 4).

Strong Purifying Selection for TCS Genes in Z. latifolia
To further understand how the TCS gene family in Z. latifolia expanded during evolution, the gene duplication events were investigated. Nineteen duplicated gene pairs in Z. latifolia TCS genes were identified including 10 ZlHK pairs, 3 ZlHP pairs, and 6 ZlRR pairs (Table S2). To study the selection pressures on the TCS gene family, nonsynonymous (Ka), synonymous (Ks), and Ka/Ks ratios were calculated. In this study, the Ka/Ks values from the 19 pairs of Z. latifolia TCS genes were much less than 1, except ZlPDK1/ZlPDK4, ZlPDK2/ZlPDK4, and ZlPDK3/ZlPDK4 (Table S2, Figure 5). This suggests that the TCS gene family in Z. latifolia underwent strong purifying selection, with a slow rate of accumulation of missense mutations during evolution.

Strong Purifying Selection for TCS Genes in Z. latifolia
To further understand how the TCS gene family in Z. latifolia expanded during evolution, the gene duplication events were investigated. Nineteen duplicated gene pairs in Z. latifolia TCS genes were identified including 10 ZlHK pairs, 3 ZlHP pairs, and 6 ZlRR pairs (Table S2). To study the selection pressures on the TCS gene family, nonsynonymous (Ka), synonymous (Ks), and Ka/Ks ratios were calculated. In this study, the Ka/Ks values from the 19 pairs of Z. latifolia TCS genes were much less than 1, except ZlPDK1/ZlPDK4, ZlPDK2/ZlPDK4, and ZlPDK3/ZlPDK4 (Table S2, Figure 5). This suggests that the TCS gene family in Z. latifolia underwent strong purifying selection, with a slow rate of accumulation of missense mutations during evolution.

TCS Gene Expression During Z. latifolia Shoots Swelling
To gain more insight into the TCS gene function in cytokinin-regulated gall formation, the transcriptome measured by RNA sequencing was mined to investigate the expression of TCS genes in shoot apices during shoot swelling in Z. latifolia. Levels of transcripts of TCS genes in wild-type (collected from uninfected shoots) were compared to samples from young but infected unswollen shoots (collected from fungi-infected plant with eight leaves) and samples from developing galls ( Figure 1A).
We noted the transient kinetics of the ZICHK1 and ZIETR2 receptor transcripts, the transcripts of the gene encoding for the red-light receptor phytochrome A (ZIPHA) transiently accumulated during gall formation and the activation of ZIPDK2 in young infected shoots ( Figure 6A). Elevated expression of the phosphotransfer protein ZIAHP1, ZIAHP2, and ZIAHP6 was associated with infected shoots ( Figure 6B). We noted that the homologous gene pairs of ZlHP2/ZlHP7 showed different expression patterns, suggesting their functional divergence after the duplication event ( Figure 6B). The response regulator encoding ZlRRA5 and ZIARR9 was activated early, whereas ZIAAR7 and ZlRRA10 type A and ZlPRR1 response regulator genes were activated later in more developed galls. Overall, levels of ARRB type encoding transcripts were lower, with ZIRRB12 displaying a higher level of transient expression kinetics ( Figure 6C). To validate the transcriptome data, the transcript levels of five genes with distinct kinetics including an HK, HPs, and an RR and PHYA were monitored by qPCR. The expression trends agreed with the RNA seq data, suggesting the RNA seq results were reliable ( Figure S6).

TCS Gene Expression During Z. latifolia Shoots Swelling
To gain more insight into the TCS gene function in cytokinin-regulated gall formation, the transcriptome measured by RNA sequencing was mined to investigate the expression of TCS genes in shoot apices during shoot swelling in Z. latifolia. Levels of transcripts of TCS genes in wild-type (collected from uninfected shoots) were compared to samples from young but infected unswollen shoots (collected from fungi-infected plant with eight leaves) and samples from developing galls ( Figure 1A).
We noted the transient kinetics of the ZICHK1 and ZIETR2 receptor transcripts, the transcripts of the gene encoding for the red-light receptor phytochrome A (ZIPHA) transiently accumulated during gall formation and the activation of ZIPDK2 in young infected shoots ( Figure 6A). Elevated expression of the phosphotransfer protein ZIAHP1, ZIAHP2, and ZIAHP6 was associated with infected shoots ( Figure 6B). We noted that the homologous gene pairs of ZlHP2/ZlHP7 showed different expression patterns, suggesting their functional divergence after the duplication event ( Figure 6B). The response regulator encoding ZlRRA5 and ZIARR9 was activated early, whereas ZIAAR7 and ZlRRA10 type A and ZlPRR1 response regulator genes were activated later in more developed galls. Overall, levels of ARRB type encoding transcripts were lower, with ZIRRB12 displaying a higher level of transient expression kinetics ( Figure 6C). To validate the transcriptome data, the transcript levels of five genes with distinct kinetics including an HK, HPs, and an RR and PHYA were monitored by qPCR. The expression trends agreed with the RNA seq data, suggesting the RNA seq results were reliable ( Figure S6).

TCS Gene Response to Exogenous Cytokinin
Our observations indicated that exogenous cytokinin can increase the frequency of gall formation and promote the growth of galls in Z. latifolia shoots. To verify whether the gall-associated TCS genes were responsive to exogenous cytokinin, the expression of 9 selected genes was measured through real-time PCR upon trans-zeatin treatment (Figure 7). Transcript levels of ZlCHK1, ZIRRA5, ZlRRA7, ZIRRA8, ZIRRA9, ZlRRA10, ZIPHYA, ZIRRB7, and pseudo-RR gene of ZlPRR1 in samples of 150 mg/L trans-zeatin sprayed shoots were compared with mock-sprayed samples after 0 (9:00 am in the morning), 2, 4, 8, and 24 h. The mock-sprayed samples displayed a similar trend for all genes measured. In these mock samples, we observed a transient decline in transcript levels followed by an increase after 8 h. In contrast, transcript levels of ZIRRA5, ZIRRA9, and ZIPHYA were elevated by zeatin treatment after two to four hours (Figure 7), and this response was transient. These results indicate that the expression of only a subset of A-type RR genes is triggered quite rapidly upon cytokinin treatment, indicating the complexity of cytokinin signaling-regulated Z. latifolia gall formation. Interestingly, the expression of red and far-red light sensor genes of ZlPHYA was quickly induced after 4 h trans-zeatin treatment (Figure 7). In Arabidopsis, PHYA, the counterpart gene of ZlPHYA played various roles in light-regulated plant development such as seed germination [56,57], internode elongation [58], and root hair growth [59]. The function of ZlPHYA in cytokinin-regulated gall swelling in Z. latifolia needs further elucidation.

TCS Gene Response to Exogenous Cytokinin
Our observations indicated that exogenous cytokinin can increase the frequency of gall formation and promote the growth of galls in Z. latifolia shoots. To verify whether the gall-associated TCS genes were responsive to exogenous cytokinin, the expression of 9 selected genes was measured through real-time PCR upon trans-zeatin treatment (Figure 7). Transcript levels of ZlCHK1, ZIRRA5, ZlRRA7, ZIRRA8, ZIRRA9, ZlRRA10, ZIPHYA, ZIRRB7, and pseudo-RR gene of ZlPRR1 in samples of 150 mg/L trans-zeatin sprayed shoots were compared with mock-sprayed samples after 0 (9:00 am in the morning), 2, 4, 8, and 24 h. The mock-sprayed samples displayed a similar trend for all genes measured. In these mock samples, we observed a transient decline in transcript levels followed by an increase after 8 h. In contrast, transcript levels of ZIRRA5, ZIRRA9, and ZIPHYA were elevated by zeatin treatment after two to four hours (Figure 7), and this response was transient. These results indicate that the expression of only a subset of A-type RR genes is triggered quite rapidly upon cytokinin treatment, indicating the complexity of cytokinin signaling-regulated Z. latifolia gall formation. Interestingly, the expression of red and far-red light sensor genes of ZlPHYA was quickly induced after 4 h trans-zeatin treatment (Figure 7). In Arabidopsis, PHYA, the counterpart gene of ZlPHYA played various roles in light-regulated plant development such as seed germination [56,57], internode elongation [58], and root hair growth [59]. The function of ZlPHYA in cytokinin-regulated gall swelling in Z. latifolia needs further elucidation. latifolia. Data are the means ± standard error (n = 3). Two percent EtOH was used as mock treatment. * Indicates the significant differences (p < 0.05) and ** indicates significant differences (p < 0.01).

Plant Growth and Treatments
Z. latifolia cv. Ivory, was provided by Zibing Ge, Shucheng Agricultural Science Research Institute, Demonstration Park of Agricultural Science Institute of Taoxi Town, Shucheng, China. The plants were grown in the field of the Anhui Agricultural University Farm (Hefei, China). To investigate the morphological changes of Z. latifolia in response to cytokinin, a 150 mg/L trans-zeatin (HENUO WEIYE Co., China) solution in 2% EtOH was sprayed onto the leaves of fungi-infected plants two weeks before their shoot swelling began. Two percent EtOH was used as mock treatment. The frequency of swelling was recorded daily over a period of 18 days. As the swelling ceased, shoots were harvested, and weight and size were recorded.
To examine the expression profiles of TCS genes responding to cytokinin, 150 mg/L trans-zeatin solution was sprayed until the fully expanded leaves were all covered with the solution. The shoots' apical tissues were then collected at 0, 1, 2, 4, 8, and 24 h after treatments. All samples were collected in three biological replicates with 3 shoots each and then stored at −80 • C until RNA extraction.

Identification of TCS Genes in Z. latifolia
The genome sequences of Z. latifolia were downloaded from the Ricerelativesb GD (http://ibi.zju. edu.cn/ricerelativesgd/download.php). The protein sequences of other TCS genes from Arabidopsis, rice, maize, tomato, soybean, and Chinese cabbage were downloaded from Phytozome (http:// phytozome.jgi.doe.gov/pz/portal.html) or NCBI (www.ncbi.nlm.nih.gov). All TCS sequences were used to build a local blast database and then used as queries to perform BLASTP (Basic Local Alignment Search Tool: Protein BLAST) searches with an E-value of 1e −5 as the threshold [24,60]. However, the putative TCS proteins of Z. latifolia were searched with the hidden Markov model (HMM) of TCS characteristic domains from Pfam (http://pfam.janelia.org/). On the basis of the above two independent methods, the redundancy removal sequences were obtained. To confirm the reliability of the searched results, all candidate protein sequences were then confirmed with SMART (http: //smart.embl-heidelberg.de/) [61], Pfam (http://pfam.xfam.org/) [62,63], and NCBI Conserved Domain Database (CDD) (https://www.ncbi.nlm.nih.gov/cdd/) [64] according to whether these sequences possess the specific structural characteristics and conserved domains of TCS elements, including HisKA, HATPase, CHASE, His-containing phosphotransfer (HPt), and receiver (Rec) domains.

Genes Structure Prediction, Protein Sequence Analysis, and Phylogenetic Analysis
Analysis of exons and introns were conducted using Gene Structure Display Server 2.04 (http: //gsds.cbi.pku.edu.cn/) [65] by comparing the coding sequences with their corresponding gene sequences. The isoelectric point (PI) and protein molecular weight in kDa were obtained using ExPASy proteomics server (http://www.expasy.org/tools/) [66]. Protein subcellular localizations were predicted using Target P online software (http://www.cbs.dtu.dk/services/TargetP/). Multiple Expectation Maximization for Motif Elicitation (MEME; http://meme.sdsc.edu/meme4_3_0/intro.html) [67] was used for motif analysis to confirm the presence of the conserved motifs in these TCS proteins. Multiple-sequence alignment for the predicted peptide sequences of conserved domains was generated using default parameters [68]. Next, phylogenetic analysis was performed with the MEGA 7 program by using the maximum likelihood (ML) method with 1000 replicates of the bootstrap based on the full-length protein sequences [69].

Calculation of Nonsynonymous (Ka) to Synonymous (Ks) Substitutions
The TCS paralogous genes in Z. latifolia were subjected to a manual BLAST search, and two genes with a similarity of >75% were considered paralogous gene pairs. However, the tandem or segmental duplication genes could not be determined due to the lack of chromosome information. The selection pressure of paralogs was then determined by analyzing the synonymous and nonsynonymous nucleotide substitutions rates (Ka/Ks ratio) of each paralogous gene pair. In theory, Ka/Ks < 1 indicates purifying or negative selection, Ka/Ks = 1 indicates neutral selection, and Ka/Ks > 1 indicates positive selection [70]. The Ka/Ks ratio was calculated using DnaSP v 5.0 software [71].

Transcriptome Analysis and Expression of TCS Genes
To obtain a better understanding of the functions of the TCS genes, we investigated their transcriptional levels in the different gall formation stages in Z. latifolia. Samples from shoots apical tissue (0.5 cm in length) were collected from wild-type Z. latifolia (no fungi infected and never swollen), and galls in different size including non-swollen, 5, 10, and 20 cm, respectively. All samples were collected in three biological replicates with 3 shoots each, frozen immediately in liquid nitrogen, and then stored at −80 • C until used. RNA was extracted and the cDNA library was constructed using Illumina TruseqTM RNA sample prep Kit. The Illumina Novaseq 6000 platform was used for transcriptome sequencing, which generated a total of 254.5 Gb of sequence data. Upon filtering, the sequencing data to eliminate poor quality sequences, the sequences were mapped to the Z. latifolia genome using TopHat2 (http;//tophat.cbcb.umd.edu/) [72]. The Transcripts Per Million reads (TPM) algorithm was used for quantification of gene expression. The heatmap of the expression patterns of TCS genes was generated by Heml 1.0 (Heatmap Illustrator) software (CUCKOO, Wuhan, China). All the transcriptomic data have been deposited at GenBank under accession PRJNA664353.

RNA Isolation and Real-Time PCR
Total RNA was extracted using an RNAprep Pure Plant Plus Kit (polysaccharides and polyphenolics rich; TIANGEN, Beijing, China). The first cDNA strand was synthesized with 2 µg of total RNA by using a FastKing RT Kit (with gDNase) (TIANGEN, Beijing, China), as per the manufacturer's instructions. qRT-PCR was performed with the CFX96 Touch C1000 Real-Time PCR System (Bio-Rad, Hercules, CA, USA) by using specific primers designed by Primer 5 Software (Table S3). Three biological and three technical replicates for each sample were performed with 20 µL of reaction volume using the SYBR ® Green Realtime PCR Master Mix (Toyobo, Osaka, Japan). The ACTIN gene was used as the housekeeping loading reference. The relative gene expression was calculated using the 2 −∆∆Ct method [73].

Conclusions
A role for cytokinin in gall formation was postulated. Here we observed a positive effect of cytokinin treatment on the gall formation in Z. latifolia in line with earlier reports [7,8], which raised the question of how cytokinin signaling regulates the gall formation. To address this problem, 69 two-component system genes, including 25 HK(L)s (21 HKs and 4 HKLs), 8 HPs (5 authentic and 3 pseudo-HPs), and 36 RRs (14 type A, 14 type B, 2 type C, and 6 pseudo-RRs), were identified in the Z. latifolia genome; the classification was supported by conserved motifs and domains, exon and intron structure, and phylogeny. Our phylogenetic analysis revealed the closest relationships between rice and Z. latifolia in seven species comparisons. Gene duplication events contributed to the TCS gene expansion in the Z. latifolia genome, and Ka/Ks analysis suggested that these duplicated gene pairs experienced strong positive selection during evolution.
We observed the activation of the expression of several of these TCS genes during gall formation, either in the very early stages or during the development of the gall. Among these, two genes encoding for RRA-type response regulators were also rapidly induced by cytokinin treatment. This transient upregulation of these RRA type was in line with RRA genes being a direct target of cytokinin signaling in other plant systems. Furthermore, the far-red light receptor phytochrome A was revealed to be a putative target.
In conclusion, these observations suggest an important role for cytokinin signaling in the process of gall formation and they highlight an important role for cytokinins in the ancient relationship between the endophyte Ustilago esculenta and its host Z. latifolia.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2223-7747/9/11/1409/s1, Figure S1: Amino acid sequence alignment of ZlHK(L)s in Z. Latifolia; Figure S2: Phylogenetic relationship, gene structures, and conserved motif of all HK(L) genes in Z. latifolia; Figure S3: Phylogenetic relationship, gene structures, and conserved motif of the HP family members in Z. Latifolia; Figure S4: Amino acid sequence alignment of ZlHP proteins in Z. latifolia; Figure S5: Phylogenetic relationship, gene structures, and conserved motif of RR genes in Z. latifolia; Figure S6: The validation of TCS gene expression by qPCR; Table S1: Summary of the two-component system (TCS) gene numbers identified in plants; Table S2: The Ka/Ks ratios of duplicated TCS genes in Z. latifolia; Table S3: The primer sequences.