Gene Ontology enrichment was performed on the male and female differentially expressed genes showing at least twofold expression differences between sexes. Gene Ontology enrichment for the male-upregulated genes showed a significant over representation of 76 terms, including pollen and anther development, male gamete development, the terpenoid biosynthesis pathway, and cytokinin metabolism . The male-upregulated gene set showed an under representation for 71 terms, including terms relating to transcription and RNA regulation, splicing, and modification . Among the female-upregulated genes, there was an over representation of RNA transcription and metabolism and shoot and organ development , and an under representation of 47 terms, including cell metabolism and biosynthesis .WGCNA network analysis was performed to explore pathways that may be involved in floral sex dimorphism and development.In total, 17,953 genes were assigned to modules, accounting for 49.8% of all expressed genes and alternative transcripts. The remaining unassigned genes were placed in the “gray” module. Three modules accounted for the majority of assigned genes: the “purple”, “cyan”, and “brown” modules. The “purple” module was strongly correlated to the female sex and captured 30.2% of all female-expressed genes , this module likely represents the genes involved in primary and secondary female sex dimorphism. In the “purple” module, 17 GO terms were enriched,fodder system including RNA and nucleic acid metabolism and regulation, photosynthesis, and phenylpropanoid metabolism . This is consistent with the female differential expression of the majority of genes in this module and supports the contention that this module may be responsible for female sex-dimorphic traits.

The cyan module was the most male-correlated and included 45.5% of all male-upregulated genes likely representing the genes in primary and secondary male sex dimorphism pathways. The “cyan” module contained 230 significantly enriched GO terms, notably multiple terms related to pollen development . The “brown” module consisted mostly of genes not showing differential expression, and probably represents gene pathways involved in basic cell and biological processes.A total of 1,381,813 cis-eQTL and 811,499 trans-eQTL were identified after accounting forcovariates. Notably, there appears to be an eQTL “hotspot” on chr07, to which the expression levels of 550 genes are associated The exceptional number of genome-wide eQTL at this locus suggests that it may have a major role in regulating sex dimorphic gene expression in catkin tissue. Expression levels of 2127 genes were found to be associated with polymorphisms in the SDR, of which 1686 were trans . To identify top-level regulatory and intermediate pathway genes, 2127 genes with eQTL in the SDR were a subset for genes with involvement in secondary metabolism , hormone signaling , RNA splicing and regulation , or transcription . Of the 96 genome-wide transcription factors found to have SDR eQTL, 70 were differentially expressed in either males or females. Because these transcription factors include genes relating to floral development, phenylpropanoid production, and cytokinin signaling, they are candidates for top-level regulatory genes that may regulate further downstream expression. However, confirmation of such roles will require further investigation using methods such as ChIPSeq and DAP-Seq. Fourteen MADS-box and floral development genes, 15 phenylpropanoid pathway genes, and five terpenoid pathway genes were also found to have eQTL in the SDR , representing candidates for intermediate pathway genes directly responsible for dimorphisms in floral morphology, pigmentation, and volatile and secondary metabolite profiles.We identified eleven genes that are strong candidates as master regulators of sex using the following criteria: presence on chr15W and absence from chr15Z, a significant log2 M:F < − 1, presence in the femalecorrelated “purple” WGCNA module, and gene annotation either consistent with a possible floral sex dimorphism pathway or of unknown function.

Genes meeting all four of these criteria are expected to be present only in females and have expression levels and module membership that would implicate them in sex dimorphism. Four copies of ARR17, a truncated AGO4 gene, DRB1, GATA15, a CCHC zinc finger nuclease, and three genes coding hypothetical proteins met these criteria and were identified as candidate master regulator genes .Total RNA-Seq and small RNA-Seq captured the unique sex-specific transcriptomic profiles during catkin development, after floral meristem differentiation and prior to maturation of any stamens or pistils. Within a single maturing catkin, there are hundreds of individual flowers across a range of developmental stages, resulting in pooled expression data from across floral development time points as well as tissue types . In addition to the primary sex dimorphism genes responsible for anther and carpel development, this enables the identification of secondary sex dimorphisms, such as genes involved in pigmentation, volatile emission, and differences in catkin phenology, which can also inform differences in vegetative emergence and secondary metabolites. By using network analysis and incorporating genomic data through eQTL, we can hypothesize how the SDR may regulate differential gene expression in catkins. Nearly two-thirds of all expressed genes in the floral tissue exhibited differential expression between males and females. This number is due in part to the large sample size of 159 individuals, whereby there is enough statistical power to detect even slight differences in expression. Nevertheless, over 21% of the expressed genes showed at least twofold expression differences between sexes, providing evidence of global expression differences, which would require robust transcriptional regulation, ultimately leading back to the sex-determinant genes in the SDR. These genes provide important clues about the regulation of sex determination in this species and the molecular mechanism responsible for diecy and floral sex dimorphism, as described in more detail below.Four copies of ARR17, a type A cytokinin-response regulator, in the SDR, show high levels of expression in female S. purpurea: Sapur.15WG073500, Sapur.15WG073900, Sapur.15WG074000, and Sapur.15WG075200.

Two additional copies of ARR17 are present on chr19 but are not differentially expressed. The cytokinin signaling pathway has been proposed as a common pathway for sex determination in angiosperms. Cytokinin-response regulators serve as feminizing factors in Actinidia, where they are master regulators, and Diospyros, where they act as top-level regulators downstream of the SDR. There is recent evidence implicating ARR17 as the master regulator of sex in the closely related genus Populus, where it may function as a feminizing factor whose expression is suppressed in males by small RNAs. The presence of two complete copies of ARR17 on S. purpurea chr19, expressed in both males and females, suggests that the dosage of ARR17 may play a role in sex determination in willow. Interestingly, these findings suggest a different mechanism for ARR17 than the leading model for sex determination in Populus proposed by Müller et al.. They proposed that functional ARR17 in P. alba is a feminizing factor, and in XY species, ARR17 is silenced by inverted repeats on the Y chromosome through the RNA-directed DNA methylation pathway, leading to a male phenotype. To confirm this, they silenced the ARR17 gene in an early-flowering female line and observed male flowers in tissue culture. We found no evidence of an ARR17 RNA interference mechanism in S. purpurea catkins. Salix purpurea has a similar truncated inverted repeat of ARR17 on chr15Z, but we did not observe small RNAs mapping to the ARR17 genes and their proximal regions, nor to the ARR17 homologs on Salix chr19. There were also no differential methylated regions in the putative promoter regions of any of the ARR17 genes, and S. purpurea males show expression of the ARR17 copies on chr19, whereas in P. trichocarpa males, there is no ARR17 expression. Furthermore, Carlson et al. did not find that ARR17 was differentially expressed in shoot tips containing floral primordia, indicating that this mechanism is not present at an earlier floral development stage either. Taken together, these results suggest that the RNAinterference mechanism of ARR17 may be absent in S. purpurea. The observation of ARR17 expression in both male and female S. purpurea, combined with a lack of small RNA loci in these same regions,fodder system for sale demonstrates that the Salix sex-determination mechanism is likely different from the model proposed by Müller et al.. Instead, our data suggest that if ARR17 is a master regulator in S. purpurea, it likely involves a unique mechanism, possibly through gene dosage, such that a threshold of ARR17 expression must be reached to activate a switch from male-to-female development. Alternatively, there may be another feature in the S. purpurea SDR that is suppressing this silencing mechanism, one such candidate is the adjacent AGO4 homolog described below.A single copy of an Arabidopsis AGO4 homolog, Sapur.15WG074400, is present within the ARR17-inverted repeat region of the chr15W SDR that exhibits a log2 M:F expression of −7.94, and has a cis-eQTL in the SDR. AGO4 is a component of the RNA-induced silencing complex in the RNA-dependent DNA methylation pathway, where it binds small RNAs and silences mRNA. In the bisulfite sequencing data, nearly three times as many regions showed increased methylation in males compared with females , supporting that methylation activity is down regulated in females and may have a role in mediating sex dimorphisms . The SDR AGO4 gene appears to be truncated to only 79 amino acids in length compared with five other catkin expressed AGO4 homologs in S. purpurea, which are 893–922 amino acids, and has multiple indels and substitutions when aligned .

The most similar AGO4 paralog to Sapur.15WG074400 by MUSCLE multiple-sequence alignment is Sapur.008G00580 which has a nearly seven fold greater expression in males . We speculate it is possible that the truncated version of AGO4 is interfering with expression of the full-length Sapur.008G00580 in males by a long noncoding RNA. This could have wide-ranging effects on sexually dimorphic gene expression and could explain the decreased genome-wide methylation observed in females. The findings from the bisulfite-sequencing data indicate that methylation is globally reduced in females. We hypothesize that the Sapur.15WG074400 could be competing for binding of siRNAs with a full-length AGO4 and sequestering male-specific RDDM in females. These global methylation differences could be responsible for sex determination, such as in Melandrium album where demethylation of male plants results in monoecy, with no effect on female plants. Such a mechanism could also explain ARR17 expression levels, and why no small RNAs were observed mapping to ARR17 in Salix, despite evidence for this mechanism in Populus.A female-expressed homolog of GATA15, Sapur.15WG 062800, is located in the W-specific region of the SDR and shows a cis-eQTL association with polymorphisms on Chr15. GATA15 is a transcriptional regulator that binds GAT or GATA motifs in gene promoters and is involved in cell differentiation, morphogenesis, and development. This is consistent with the GO enrichment analysis of female-expressed genes, which contains many significant terms related to morphogenesis and development. Furthermore, chr15 GATA15 was found by Carlson et al. to be differentially expressed in F1 S. purpurea shoot tips containing floral primordia, suggesting that this may be the earliest cue for floral sex differentiation, which would implicate it as a master regulator gene. While functional genomics data are required to elucidate its precise function, its expression in females both during floral differentiation and catkin emergence suggests that it could be directly involved in gynecium development.Four genes were identified that fit the criteria for candidate master regulator genes, but whose functions are not known or whose annotations are insuffificient for further analysis. These included Sapur.15WG068800, a CCHC-type zinc finger, and three hypothetical proteins: Sapur.15WG075300, Sapur.15WG074900, and Sapur.15 WG075700. While there is mounting evidence pointing toward ARR17 as the master regulator in Populus spp, the evidence for different expression profiles of ARR17 in Salix, as well as the presence of additional candidate genes, suggests that the mechanism may be more complicated or altogether different in Salix. Nevertheless, expression data from the ARR17 homologs in S. purpurea do support a possible role in sex determination, either as a single gene master regulator or part of a two-gene system in conjunction with another master regulator, and would provide further evidence to support cytokinin response as a common mechanism for diecy in angiosperms, assuggested by Montalvão et al. Further functional genomics studies will be necessary to elucidate the precise functions of candidate master regulators and their role in sex determination.Among the floral development genes with eQTL in the SDR were homologs of AGL11 and AGL32, AGL29 and AGL30, and AGL6, as well as TOC1, WOX1, RGA, and CONSTANS. While differential expression of MADS-box genes is expected in floral tissues, their association with the SDR through eQTL, even after accounting for sex as a covariate, suggests that the SDR may have a direct role in controlling expression of these genes and subsequent primary sex dimorphisms.