Co-expression network analysis of the genes associated with pistillody-stamen development in wheat CURRENT

Background: Crop male sterility has great values in both theoretical research and breeding application. Wheat pistillody-stamen is an important male sterility phenomenon, and HTS-1 is an important pistillody-stamen material. However the molecular mechanism of HTS-1 stamens transformed into pistils or pistil-like structures remains a mystery. Weighted gene co-expression network analysis (WGCNA) are widely used to explore hub genes and gene interaction networks from high throughput data in various plants. Results: In the present study, for exploring gene networks associated with wheat pistillody-stamen development, WGCNA was employed to analyze 11 RNA-sequencing (RNA-seq) data of wheat tissues, including stamens of CSTP, pistils and pistillody-stamen of HTS-1. 19 out of 25 merged modules were highly associated with specific wheat tissues, and the MEdarkseagreen1 module was highly related to wheat pistillody-stamen (correlation with weight r =0.7, correlation p-value p =0.02). Then 180 genes about wheat flower development were identified from the MEdarkseagreen1 module by GO term analysis. Among 180 genes, the hub gene number associated with anther, filament, style, and ovary development were 12, 3, 3, and 10, respectively. We compared the published pistillody related proteins with proteins of HTS-1 by BLAST. A total of 58 pistillody-stamen development associated proteins were validated by BLAST. MADS-box and YABBY transcription factor about pistillody-stamen development were also analyzed in wheat flower. There were 47 of MADS-box and 17 of YABBY transcription factors were identified. BLAST program was used to align the published pistillody associated MADS-box and YABBY transcription factors with transcription factors identified in wheat flower. Totally, 36 of 47 MADS-box and 14 of 17 YABBY transcription factors were considered to regulate the development of pistillody-stamen, which had never been reported yet. Conclusion: These results have systematically identified the key candidate genes about the development of HTS-1 substructures flower. The tissue-specific correlation network analyses provide important insights into the molecular interactions underlying psitillody-stamen development. RNA Miniprep Kit (Biomiga, China). DNase treatment was performed before proceeding with complementary DNA (cDNA) synthesis. The quality of the cDNA was determined by gel electrophoresis. cDNA concentration was calculated and adjusted to 50 ng/µl by spectrophotometer (Thermo, Nanodrop 2000).

[1]. According to the report of United States department of agriculture, the global output of wheat in 2017, 2018, and July 2019 were accounted for 761.88, 730.90, and 771.46 million tons, respectively [2]. Although the global output of wheat are steadily increased in recent years, some studies show that wheat yields are hard to meet demand of increasing human population [3]. Thus many researchers strive to improve grain productivity and qualities. The development of wheat flowers affect grain yields and qualities. There are various wheat flower mutant materials have been reported in the past decade. Mutant (cr)-N26 and (cr)-CSdt7BSare male sterility material, whose stamen of floret will transform into pistil or pistil-like structure (called pistillody) [4][5][6]. Both (cr)-N26 and (cr)-CSdt7BS mutants contain exogenous genes of goat weed, so the pistillody-stamen phenomena of (cr)-N26 and (cr)-CSdt7BS are caused by heterologous cytoplasmic male sterility [5]. In 2013, Peng et al.
Since HTS-1 mutants do not contain exogenous cytoplasm, the pistillody-stamen phenomena of HTS-1 mutant is controlled by nuclear genes [8]. Previous reports show that the molecular mechanisms of wheat stamens transformed into pistils of floret are very complex. In 2002, Murai et al. found that reduced WAP3 expression could be associated with the induction of pistillody in (cr)-CSdt7BS [5]. The mitochondrial gene orf25 of Ae. crassa cytoplasm had been considered to be a candidate as the cytoplasmic factor related to pistillody in the alloplasmic wheat [9]. Zhu et al. reported a mitochondrial gene called orf260 cra , which causing the stamen of (cr)-CSdt7BS transformed into pistillody [10]. Situ hybridization analysis showed that a calmodulin-binding protein gene called WCBP1 was highly expressed in the pistil-like stamens at early to late developmental stages. This suggested that WCBP1 gene played an important role in pistil-like stamen development [11].The homologous transformation of stamens into pistils of HTS − 1 floret was determined by the interaction of two recessive karyogenes hts1 and hts2 [8]. In 2015, Yang et al. identified 206 differentially expressed genes (DEG) correlated to HTS-1 stamen and pistil development using RNA-seq analysis.
123 of 206 genes were down-regulated genes and the remaining of 83 genes were up-regulated genes [12]. In 2019, sun et al. reported a TaEPFL1 gene, which was highly expressed in pistillody-4 stamens than in pistils and stamens. Heterologous expression of the TaEPFL1 gene in Arabidopsis caused shortened filaments and pedicels, which might induce stamens transformation into pistils or pistil-like structures in wheat [13]. Furthermore, some transcription factors of MADS-box and DROOPING LEAF gene of YABBY family have been considered to regulate in the development of wheat pistillody-stamen [6,9,[14][15][16][17][18]. The expression pattern analysis indicate that TaDL1, TaDL2, and TaDL3 genes might cause the stamen of HTS-1 complete or partial transforming into pistils [16].
These researches let us to know the development of wheat pistillody stamen. Whereas the genes mechanisms of stamen transformation into pistil are still unclear.
With rapid improvement of high throughput sequencing and computer analysis technology, we can systematically study complicated biological problems using variety of high throughput sequencing data. Various network analysis methods have been constructed and widely used to reveal the interaction of molecular in high throughput sequencing data analysis, such as neural network analysis [19], protein-protein interaction network analysis [20] and WGCNA [21]. WGCNA is an important method for studying gene regulation network. According to gene expression levels in different tissues, genes with similar expression patterns can be clustered into the same module by WGCNA [21]. Therefore, WGCNA is widely used to explore hub genes, relationship between modules and traits, and gene interaction networks of flower development in various plants [22][23][24]. Particularly, Ramírez-González et al. using 209 RNA-seq samples from 22 tissue types analyzed gene expression pattern of wheat Azhurnaya cultivar, and established tissue-and stress-specific co-expression networks that revealed extensive coordination of homoeolog expression throughout wheat development [1]. These networks, including spikelet development networks, and detailed gene expression atlases were deposited in wheat eFP browser (http://bar.utoronto.ca) [1]. However networks in wheat eFP browser don't reveal and provide the information about mechanism of genes for pistillody-stamen development. Because Azhurnaya cultivar is not a pistillody-stamen wheat. Furthermore, systematic identifications of genes associated pistillody-stamen development by WGCNA have never been reported yet.
In the present study, to explore genes related to wheat pistillody-stamen development and their 5 interaction mechanisms, we employed WGCNA to analyze 11 RNA-sequencing (RNA-seq) data of wheat tissues, including stamens of CSTP, pistils and pistillody-stamen of HTS-1 [12]. We surveyed the relationship between modules and specific flower tissue from wheat RNA-seq data. Then we identified genes associated with pistillody-stamen development from the modules highly by KOBAS. In order to validate the results explored by WGCNA, expression pattern of 20 genes related to the development of pistillody-stamen in HTS-1 were performed by Quantitative real-time PCR (qRT-PCR).
This work will provide further insights into the molecular mechanisms underlying wheat pistillodystamen development.

Results
A total of 58,576 genes obtained by Salmon were used to construct WGCNA networks. There were 47,492 genes remained after filtering out the genes with low expression level (TPM < 1). Once the coexpression network was constructed, the next step was module identification. Modules in coexpression network were clusters of highly interconnected genes. Figure 1B showed the gene dendrogram clustering on Tom-based dissimilarity and modules obtained by the dynamic tree cut.
There were 224 modules were clustered by the similarity of their eigengenes, and each module described by a branch of different color. If a module was highly correlated to other modules, they should be merged into one module. Finally, Fig. 1B depicted the major tree branches constitute 25 merged modules labeled by different colors (cutHeight = 0.25).
WGCNA enabled us to reveal the relationships between the modules and tissue types (Fig. 2). When the p-value was less than 0.05, it indicated that this module was significant correlation with the tissue type. Figure 2 depicted that 19 out of 25 modules were highly associated with specific wheat tissue types, excluding MEburlywood, MEdarksalmon, MEcoral2, MEcoral4, MEgrey, and MElavenderblush.
Each tissues had at least associated with a specific modules, excluding pistil and late ovary. Figure 2 displayed that the MEdarkseagreen1 module was highly related to wheat pistillody-stamen (correlation with weight r = 0.7, correlation p-value p = 0.02). The MEmediumpurple3 (r = 0.67, p = 0.02) and MEturquoise (r = 0.99, p = 3e-09) modules were significantly correlated with wheat stamen.
To investigate genes related to HTS-1 flower development in MEdarkseagreen1 module, we performed a Gene Ontology (GO) analysis by Kobas 3.0 (Additional file 1: Table S1). In Additional file 1: Table S1, there were 691 GO term records related to wheat flower substructure development. Further analysis revealed that the 691 GO terms were identified from 180 genes (Additional file 1: Table S1). Table 1 depicted genes number corresponding to different GO terms of wheat flower development, which were accounting from Additional file 1: Table S1. Table 1 showed that there were 51 genes about anther development, which were assigned to 8 different GO terms. There were 5 genes related to wheat style development in GO term of GO:0048479. There were 3 wheat filament development associated genes in GO term of GO:0080086. A total of 40 genes related to ovary development were identified in GO term of GO:0035670, and there were 20 genes about ovule development in GO term of GO:0048481. Wheat stamen and pistil had 3 and 2 types of GO terms, respectively. Moreover, there were 423 genes involved in wheat flower development, which were belonging to 22 types of GO terms (Table 1). Table 1 The genes number of different GO terms related to wheat flower development in pistillody stamen.    (Table 2). Table 2 Transcription factor identified from different wheat flower sub-structure tissue associated modules.  Figure 3 depicted the genes networks related to wheat flower development in MEdarkseagreen1 module. In Fig. 3, each node represents a gene.
The lines between nodes represent the correlations. The networks showed the genes networks of wheat anther (Fig. 3A), style (Fig. 3B), filament (Fig. 3C), and ovary (Fig. 3D), respectively. In Fig. 3, the larger the node, the more connectivity its input gene had in the networks. This indicated that they had more important roles involved in wheat flower development. There were twelve hub genes in blue color associated anther development, such as TraesCS4A02G059000.1, TraesCS3D02G425800.1, TraesCS4D02G235800.1 (Fig. 3A). Genes of TraesCS4A02G058800.3, TraesCS4A02G058800.1, and TraesCS4A02G058800.5 in pink color had greater connectivity than other genes,demonstrating that they had more regulating role with style development (Fig. 3B). Figure 3C displayed that genes of TraesCS5B02G202000.1, TraesCS5A02G203300.1, and TraesCS5D02G209700.1 played a significant role in filament development. There were 10 hub genes affected ovary development showing in Fig. 3D, such as TraesCS2A02G514200.1, TraesCS2B02G542400.2, TraesCS4A02G422700.1.

Discussion
WGCNA can identify hub genes from tissue-specific modules, explore gene interaction networks, and carry out functional annotation about unknown genes, according to genes with similar functions always have similar expression levels [21]. Therefore it has been widely used to reveal genes interaction networks of flower development in different plants [22,23,25].  Table S1). We have queried for the information of 180 wheat flower development associated genes in wheat eFP browser. The query results in wheat eFP browser has proved that all of the 180 genes are expressed in wheat spikelet, which suggest that they are affected with the development of wheat flower. Hence the query results prove that WGCNA is an efficient and accurate analysis method to cluster genes from high throughput data. However networks in wheat eFP browser don't reveal and provide information about mechanism of genes for pistillodystamen development. Because Azhurnaya cultivar is not a pistillody-stamen wheat. WGCNA results of this research can not only reveal the gene modules associated pistillody-stamen development, but also identify specific gene modules involving in other wheat tissues, except pistil and late ovary (Fig. 2). The MEdarkseagreen1 module is highly related to wheat pistillody-stamen (Fig. 2). Moreover, gene number related HTS-1 style and ovary development is 5 and 40, respectively (Additional file 1: Table S1). In order to verify the accurate of candidate genes related to pistillody-stamen development, we selected 20 candidate genes to verify their expression pattern by qRT-PCR (Fig. 4).
In Fig. 4, 14 proteins of 20 sequences are related to ovary development of pistillody-stamen, and 5 of 20 proteins are associated with the development of style in pistillody-stamen (Additional file 1: Table   S1). Among 14 ovary genes mentioned above, 7 gene are hub genes showing in Fig. 3D, which IDs are as follows: TraesCS2A02G514200.1, TraesCS6B02G012800.1, TraesCS4D02G101700.1, TraesCS7D02G060400.2, TraesCS3D02G442000.5, TraesCS2B02G542400.2, TraesCS5A02G155100.1. The pistillody and pistil expression level of the 7 hub genes of ovary development are apparently higher than that of in stamen (Fig. 4). This result infer that the 7 hub genes must participate to regulate the development of ovary. Among 5 style associated genes, 3 genes are hub gene related the development of style; Their IDs are as follows: TraesCS4A02G058800.3, TraesCS4A02G058800.1, TraesCS4A02G058800.5 (Fig. 3B). The expression level of the 3 style related hub genes are barely expressed in stamen, and their pistillody-stamen expression are greater than in stamen (Fig. 4). This indicate that TraesCS4A02G058800.3, TraesCS4A02G058800.1, TraesCS4A02G058800.5 must regulate the development of pistillody-stamen. qRT-PCR results show that the prediction genes associated with pistillody-stamen development by WGCNA are very reliable (Fig. 4).
The development of flower features directly affect crop yield and other agronomic traits. Male sterility has been widely studied in flower development, and has very important application value in breeding [8]. Stamen transformation into pistil or pistil-like structure (called pistillody-stamen) is rare phenomenon of male sterility. Previous studies have reported some pistillody-stamen material in wheat, such as (cr)-n26, (cr)-csdt7bs, and HTS-1 [4,5,8], which let us have some insight into understanding the molecular interaction mechanism of wheat pistillody-stamen genes [5,[10][11][12]. We compare the published pistillody related proteins with proteins of HTS-1 by BLAST. The alignment results are deposited in Additional file 2: Table S2. The BLAST results display that protein of WAP3 gene (NCBI Accession: BAA33459.1) [5] have two highly similar proteins in HTS-1, which IDs are TraesCS7B02G286600.2 and TraesCS7A02G383800.1 (Additional file 2: Table S2). Both of the identity value and the ratio of alignment sequence cover the total length between WAP3 gene and HTS-1 genes are greater than 95% (Additional file 2: Table S2). Murai et al. reported that the expression level of WAP3 decrease could induce the pistillody of (cr)-csdt7bs [5]. TraesCS7B02G286600.2 and TraesCS7A02G383800.1 of HTS-1 are confirmed to be expressed in stamen of HTS-1 (Additional file 2: will maintain a certain expression level in stamen. When its amount of expression decrease below a certain level, the stamen will transform into pistil. The mitochondrial gene orf25 protein (BAA82046.1) [9] are highly similar with TraesCS4D02G236200.1 and TraesCS7A02G099300.1 of HTS-1. This result suggest that TraesCS4D02G236200.1 and TraesCS7A02G099300.1 can induce pistillody of HTS-1.
TraesCS4D02G140600.1 and TraesCS7A02G237800.1 will be used as candidate gene for further study. Protein of WCNP1 gene (BAM65843.1) was significantly up-regulated in young spikes of the pistillody-stamens in (cr)-csdt7bs induced by mitochondrial retrograde signaling [11]. BAM65843.1 shows high similarity with HTS-1 proteins of TraesCS5D02G102100.1 and TraesCS5B02G095800.1, whose identity and ratio are greater than 95% (Additional file 2: Table S2). Whereas the function of  Table S2). The identities and ratios between 123 DEG and 214 records are above 95% and 90%, respectively (Additional file 2: Table S2). stamens [13]. The alignment identity of TaEPFL1 gene against TraesCS6D02G296500.1 and TraesCS6B02G347500.1 were 100%. This result is in accord with the result of Sun et al. [13]. Finally, a total of 58 pistillody-stamen development associated proteins are validated by BLAST in HTS-1. These 58 sequences will be considered as pistillody-stamen development candidate genes for future studies.
Transcription factors play an important role in regulating the development of plants. Some previous studies have reported that MADS-box [9,14,16,17,26,27,28] and DROOPING LEAF gene of YABBY family [6,15,16] associated with the development of wheat pistillody-stamen. We identify MADS-box and YABBY transcription factors from wheat flower substructure tissues (Table 2). Totally, there are 47 MADS-box identified in early ovar, pistillody-stamen, and stamen (Table 2) Table S1 by BLAST. The alignment results are shown in Additional file 3: Table S3. In Additional file 3: Table S3, DROOPING LEAF gene of TaDL3 sequence is complete the same as proteins of TraesCS4B02G245900.4 and TraesCS4B02G245900.5 of HTS-1.
Furthermore, TraesCS4A02G058800.3, TraesCS4B02G245900.4 and TraesCS4B02G245900.5 have been verified highly expressed in pistil or pistillosy-stamen of HTS-1 by qRT-PCR (Fig. 4)  Analysis of gene quantification Salmon [29] was utilized to quantify gene expression of 11 wheat tissues mentioned above, and 133,346 CDS were used as reference sequences. We filtered out the gene quantification of transcripts per million (TPM) values less than 0.5 in pistillody stamen, and there were 58,576 genes remained. A total of 58,576 genes with TPM values in 11 wheat tissues were obtained.

WGCNA analysis
The WGCNA R package [30] was applied to construct the co-expression network of 11 wheat tissues (58,576 genes). Firstly, genes with mean TPM value less than 1 were filter out, and there were 47,492 genes remained. In order to make the genes in the network conform to the scale-free topology distribution, a soft threshold (power) must be set. When scale-free topology fitted index of 0.9, the first soft threshold (power) of 18 was used (Fig. 1A). Next the adjacency matrix was converted into topological overlap matrix (TOM), and the function dissTom = 1-TOM was used to invert the TOM into their invert matrix. Then Dynamic Tree Cut was adopted for module identification, and the minimum module size was 30. The function hclust was utilized to merge the similar modules (mergeCutHeight = 0.25) (Fig. 1B). To identify the tissue specific modules, we analyzed the correlation between the tissue and module eigengene (ME) (Fig. 2). The function corPvalueStudent() was employed to do the significance test between the trait and module eigengene.

GO analysis and visualization
The result of WGCNA identified that the MEdarkseagreen1 module was significantly related to the development of pistillody-stamen of HTS-1. Gene Ontology (GO) analysis was planned to find out genes related to pistillody-stamen development in MEdarkseagreen1 module. A total of 4,702protein sequences were obtained from the Ensembl Plants database (version release-41), which edge values were more than 0.3 in MEdarkseagreen1 module. Kobas 3.0 was applied to do GO annotation of the 4,702 proteins [31]. Additional file 1: Table S1exhibited genes related wheat pistillody-stamen development by Kobas 3.0. According to GO analysis results showed in Additional file 1: Table S1, genes involved in wheat anther, stigmas, style, filament, and ovary development were collected from MEdarkseagreen1 module. Table 1 displayed the genes number of different GO terms associated with wheat flower development in pistillody-stamen, which were accounting from Additional file 1: Table   S1. Then, for identifying the unpublished proteins about pistillody-stamen development of HTS-1, we aligned the published pistillody related proteins against HTS-1 proteins by BLAST. The alignment results were deposited in Additional file 2: Table S2.
Cytoscape V3.3 [32] was used to generate genes network visualization of different flower tissues (Fig. 3). Genes with higher connectivity were considered as hub genes [1].

Transcription factors analysis
To seek for the transcription factors about pistillody-stamen development, we identified MADS-box and YABBY transcription factors from specific modules highly correlated to different flower tissues using the method of Chen et al. [33]. The MADS-box and YABBY transcription factors for different flower tissues were shown in Table 2. In order to find out unpublished transcription factors related to pistillody-stamen development, we utilized BLAST program to compare the published wheat MADSbox and YABBY transcription factors related to pistillody with HTS-1 transcription factors in Table 2.
The BLAST results were displayed in Additional file 3: Table S3.

Plant materials and RNA extraction
Both CSTP and HTS-1 are common wheat maintained in our laboratory. CSTP is a near-isogenic line of Chinese Spring, which has three pistils in each floret [8]. HTS-1 is a pistillody wheat mutant derived from CSTP [8]. CSTP and its pistillody mutant of HTS-1 were planted in a field in the China West  Table S4. Each reaction contained 5 µl of SSofast (Bio-Rad), 1 µl of gene-specific primers, 0.5 µl of cDNA, 3.5 µl of ddH 2 O, and in a final volume of 10 µl.
The wheat genes (GenBank No. AB181911) and wheat housekeeping genes Ubiq (DQ086482) were used as a reference gene [34]. Fold-changes of RNA transcripts were calculated via the 2-ΔΔCt method [35]. Expression pattern of 20 genes associated pistillody-stamen development were involved in study design, data collection and analysis, or preparation of the manuscript.

Availability of data and materials
The RNA-seq datasets can be obtained from the NCBI SRA database, which have been mentioned in materials and methods of this study. All datasets generated during this study are available from the corresponding author on reasonable request.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interest.  Expression patterns of twenty pistillody-stamen development genes verified by real-time PCR. P represents fully developed wheat pistill of HTS-1; PS represents completely developed wheat pistiilodystamen of HTS-1; S represents fully developed wheat stamen of CSTP.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. Table S3.xlsx   Table S4.xlsx   Table S1.xlsx Table S2.xlsx