{'fraction_non_alphanumeric': 0.025872093023255813, 'fraction_numerical': 0.05489341085271318, 'mean_word_length': 5.0601879036993545, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 36, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integrative genomics reveals paths to sex dimorphism in Salix purpurea L.

Sex dimorphism and gene expression were studied in developing catkins in 159 F2 individuals from the bioenergy crop Salix purpurea, and potential mechanisms and pathways for regulating sex development were explored. Differential expression, eQTL, bisulfite sequencing, and network analysis were used to characterize sex dimorphism, detect candidate master regulator genes, and identify pathways through which the sex determination region (SDR) may mediate sex dimorphism. Eleven genes are presented as candidates for master regulators of sex, supported by gene expression and network analyses. These include genes putatively involved in hormone signaling, epigenetic modification, and regulation of transcription. eQTL analysis revealed a suite of transcription factors and genes involved in secondary metabolism and floral development that were predicted to be under direct control of the sex determination region. Furthermore, data from bisulfite sequencing and small RNA sequencing revealed strong differences in expression between males and females that would implicate both of these processes in sex dimorphism pathways. These data indicate that the mechanism of sex determination in Salix purpurea is likely different from that observed in the related genus Populus. This further demonstrates the dynamic nature of SDRs in plants, which involves a multitude of mechanisms of sex determination and a high rate of turnover.


Introduction 51
Major progress has been made in recent years in identifying the master regulators of sex 52 determination in plants, but less is known about the transcriptional networks that account for sex 53 dimorphism. Transcription factors 1 (Asparagus), small RNAs 2 (Diospyros) and cytokinin 54 response regulators (Actinida and Populus) 3,4 have all been explicated as mechanisms of sex 55 determination in angiosperms. Development of separate sexes, however, requires that the genes 56 controlling sex determination regulate the transcription of many genes and metabolic pathways 57 in order to coordinate development of sex-specific characteristics, such as gametes and floral 58 morphology. With the notable exception of Diospyros 5 , little is known about the metabolic regulator genes in Salix purpurea is complicated by an SDR that comprises nearly 40% of 120 chr15W and contains 488 linked genes, including repetitive regions, and tandem duplications 14  genes showing at least four-fold expression differences between sexes (Tables S1-S4). Gene 154 Ontology enrichment for the male upregulated genes showed a significant overrepresentation of 155 76 terms, including pollen and anther development, male gamete development, the terpenoid 156 biosynthesis pathway, and cytokinin metabolism (Table S1). The male upregulated gene set 157 showed an underrepresentation for 71 terms, including terms relating to transcription and RNA 158 regulation, splicing, and modification (Table S2). Among the female upregulated genes, there 159 was an overrepresentation of RNA transcription and metabolism and shoot and organ 160 development (Table S3), and an underrepresentation of 47 terms including cell metabolism and 161 biosynthesis (Table S4)  A total of 604,150 methylated regions were identified using DMRfinder after combining nearby 191 methylated sites. Of these sites, 2,018 show differential methylation, including 1,465 with 192 greater methylation in males and 553 with greater methylation in females (Table S6). 170 genes 193 were identified with differentially methylated sites in their putative promoter regions (between 1 194 and 1000 bp ahead of transcription start site) (Table S7). 195

WGCNA Network Analysis and GO Enrichment 196
WGCNA network analysis was performed to explore pathways that may be involved in floral sex 197 dimorphism and development. Twenty-four modules were generated based on gene expression 198 similarities (Fig. 3, Table 2). 17,953 genes were assigned to modules, accounting for 49.8% of 199 all expressed genes and alternative transcripts. The remaining unassigned genes were placed in 200 the "grey" module. Three modules accounted for the majority of assigned genes: the "purple", 201 "cyan" and "brown" modules. The "purple" module was strongly correlated to the female sex 202 (r 2 =0.97) and captured 30.2% of all female expressed genes (Table 2), this module likely 203 represents the genes involved in primary and secondary female sex dimorphism. In the "purple" 204 module 17 GO terms were enriched, including RNA and nucleic acid metabolism and regulation, 205 photosynthesis, and phenylpropanoid metabolism (Table S8). This is consistent with the female 206 differential expression of the majority of genes in this module (Table S3) and supports the 207 contention that this module may be responsible for female sex-dimorphic traits. The cyan module 208 was the most male correlated (r 2 =0.92) and included 45.5% of all male upregulated genes likely 209 representing the genes in primary and secondary male sex dimorphism pathways. The "cyan" 210 module contained 230 significantly enriched GO terms, notably multiple terms related to pollen 211 development (Table S9). The "brown" module consisted mostly of genes not showing 212 differential expression, and probably represents gene pathways involved in basic cell and 213 biological processes. 214 215 216  Notably, there appears to be an eQTL "hotspot" on chr07, to which the 229 expression levels of 550 genes are associated (Fig. 2) The exceptional number of genomewide 230 eQTL at this locus suggests that it may have a major role in regulating sex dimorphic gene 231 expression in catkin tissue. Expression levels of 2,127 genes were found to be associated with 232 polymorphisms in the SDR, of which 1,686 were trans (Fig. 2). To identify top level regulatory 233 and intermediate pathway genes, 2,127 genes with eQTL in the SDR were subset for genes with 234 involvement in either secondary metabolism (22), hormone signaling (16), RNA splicing and 235 regulation (4) or transcription (96) (Table S10, Fig. S3). Of the 96 genome-wide transcription 236 factors found to have SDR eQTL, 70 were differentially expressed in either males or females. 237 Because these transcription factors include genes relating to floral development, phenylpropanoid 238 production, and cytokinin signaling, they are candidates for top-level regulatory genes that may 239 regulate further downstream expression. Fourteen MADS-box and floral development genes, 15 240 phenylpropanoid pathway genes, and five terpenoid pathway genes were also found to have 241 eQTL in the SDR (Table S10) coding hypothetical proteins met these criteria and were identified as candidate master regulator 253 genes (Table 3). 254 We have used a combination of differential expression, co-expression, and eQTL analyses to in tissue culture 4 . We found no evidence of an ARR17 RNA interference mechanism in S. 301 purpurea catkins. Salix purpurea has a similar truncated inverted repeat of ARR17 on chr15Z, 302 but we did not observe small RNAs mapping to the ARR17 genes and their proximal regions, nor 303 to the ARR17 homologs on Salix chr19. There were also no differential methylated regions in the 304 putative promoter regions of any of the ARR17 genes, and S. purpurea males show expression of 305 the ARR17 copies on chr19, whereas in P. trichocarpa males, there is no ARR17 expression. 306 Furthermore, Carlson et al. (2017) did not find that ARR17 was differentially expressed in shoot 307 tips containing floral primordia, indicating that this mechanism is not present at an earlier floral 308 development stage either 27 . Taken together, these results suggest that the RNA-interference 309 mechanism of ARR17 may be absent in S. purpurea. 310 The observation of ARR17 expression in both male (chr19) and female (chr19 and chr15W) S. 311 purpurea, combined with a lack of small RNA loci in these same regions, demonstrates that the 312 Salix sex determination mechanism is likely different from the model proposed by Müller et al. 313 (2020) 4 . Instead, our data suggest that if ARR17 is a master regulator in S. purpurea, it likely 314 involves a unique mechanism, possibly through gene dosage, such that a threshold of ARR17 315 expression must be reached to activate a switch from male to female development. Alternatively, 316 there may be another feature in the S. purpurea SDR that is suppressing this silencing 317 mechanism, one such candidate is the adjacent AGO4 homolog described below. 318

ARGONAUTE 4 (AGO4) 319
A single copy of an Arabidopsis AGO4 homolog, Sapur.15WG074400, is present within the 320 and silences mRNA 30 . In the bisulfite sequencing data, nearly three times as many regions 324 showed increased methylation in males (1465) compared to females (553), supporting that 325 methylation activity is downregulated in females and may have a role in mediating sex 326 dimorphisms (Fig. 2). The SDR AGO4 gene appears to be truncated to only 79 amino acids in 327 length compared to five other catkin-expressed AGO4 homologs in S. purpurea, which are 893 to 328 922 amino acids, and has multiple indels and substitutions when aligned (Fig. S4A). The most 329 similar AGO4 paralog to Sapur.15WG074400 by MUSCLE multiple sequence alignment is 330 Sapur.008G00580 which has a nearly seven-fold greater expression in males (Fig. S4B, Table  331 S11). We speculate it is possible that the truncated version of AGO4 is interfering with 332 expression of the full-length Sapur.008G00580 in males by a long non-coding RNA. This could 333 have wide-ranging effects on sexually dimorphic gene expression and could explain the 334 decreased genome wide methylation observed in females. The findings from the bisulfite-335 sequencing data indicate that methylation is globally reduced in females, suggesting the 336 possibility that the RDDM pathway may be compromised in females. We hypothesize that the 337 Sapur.15WG074400 could be competing for binding of siRNAs with a full-length AGO4 and 338 sequestering male-specific RDDM in females. Such a mechanism could also explain ARR17 339 expression levels, and why no small RNAs were observed mapping to ARR17 in Salix, despite 340 evidence for this mechanism in Populus. 341

Double Stranded RNA binding protein 1 (DRB1) 342
A copy of a DRB1 homolog, Sapur.15WG074300, is also located in the ARR17 inverted repeat 343 region of the chr15W SDR and is adjacent to AGO4. In Arabidopsis, DRB1 is involved in RNA-344 mediated post-transcriptional silencing, working directly with DCL1 in miRNA processing 31 . 345 RNA modification, RNA metabolic process, and regulation of transcription were all among the 346 enriched GO terms in the female differentially expressed genes and could be regulated directly or 347 indirectly by DRB1. 348

Transcription Factor GATA15 349
A female expressed homolog of GATA15, Sapur.15WG062800, is located in the W-specific 350 region of the SDR and shows a cis-eQTL association with polymorphisms on Chr15. GATA15 is 351 a transcriptional regulator that binds GAT or GATA motifs in gene promoters and is involved in 352 cell differentiation, morphogenesis, and development 32 . This is consistent with the GO 353 enrichment analysis of female expressed genes, which contains many significant terms related to 354 morphogenesis and development. Furthermore, chr15 GATA15 was found by Carlson et al. 355 (2017) to be differentially expressed in F1 S. purpurea shoot tips containing floral primordia 27 , 356 suggesting that this may be the earliest cue for floral sex differentiation, which would implicate it 357 as a master regulator gene. While functional genomics data is required to elucidate its precise 358 function, its expression in females both during floral differentiation and catkin emergence 359 suggests that it could be directly involved in gynoecium development. 360

Genes of Unknown Function 361
Four genes were identified that fit the criteria for candidate master regulator genes, but whose 362 functions are not known or whose annotations are insufficient for further analysis. These 363 included Sapur.15WG068800, a CCHC-type zinc finger, and three hypothetical proteins: 364 Sapur.15WG075300, Sapur.15WG074900, and Sapur.15WG075700. 365 While there is mounting evidence pointing towards ARR17 as the master regulator in Populus 366 spp. 4,24,25 the evidence for different expression profiles of ARR17 in Salix, as well as the presence 367 of additional candidate genes, suggests that the mechanism may be more complicated or 368 altogether different in Salix. Nevertheless, expression data from the ARR17 homologs in S. 369 purpurea does support a possible role in sex determination, either as a single gene master 370 regulator or part of a two-gene system in conjunction with another master regulator, and would 371 provide further evidence to support cytokinin response as a common mechanism for dioecy in 372  (Table S10). While differential expression of 380 MADS-box genes is expected in floral tissues, their association with the SDR through eQTL, 381 even after accounting for sex as a covariate, suggests that the SDR may have a direct role in 382 controlling expression of these genes and the resulting primary sex dimorphisms. These floral 383 development genes may have specific roles as "entry points" for the SDR into the floral 384 development pathway to regulate the development of a particular sex (Fig. 4). Cronk and Müller 385 (2020) proposed that ARR17 may act as a feminizing master regulator in Populus through the 386 suppression of PISTILLATA (PI) or APETALLA3 (AP3) MADS-box genes 25 . Importantly, there 387 was no association of either PI or AP3 expression with genomic variation in the SDR, despite the 388 fact that PI shows very high levels of expression in males (log2 M:F 13.21) (Table S11), further 389 suggesting that the mechanism of sex determination in Salix may be different from Populus. 390 eQTL analysis revealed several loci with an exceptionally high number of trans-eQTL (Fig. 2). 391 Intriguingly, The hotspot with the greatest number of trans-eQTLs is located in the region on 392 chr07 homologous to the S. nigra SDR 12 , which could either point towards its role as an ancestral 393 dimorphism genes in this region. Approximately 250 kb from this locus on chr07 is 395 Sapur.007G068100, a homolog of AGL32, a MADS-box gene involved in ovule development 33 . 396 The expression of this gene is associated in trans with the chromosome 15 SDR region, further 397 supporting its role as a top-level regulatory gene under direct regulation by the SDR. 398 Compounds produced from the phenylpropanoid and terpenoid pathways are well characterized 399 in Salicaceae, and there is evidence that floral volatile, terpenoid, and phenolic glycoside profiles 400 differ substantially between males and females, which affects both pollinator attraction and 401 herbivory, traits that are likely to be evolutionary drivers of dioecy and affect cultivar yield 22,23 . 402 In support of this, we identified five terpenoid pathway genes and 15 phenylpropanoid pathway 403 genes with eQTL in the SDR. These include both genes involved in the core phenylpropanoid 404 pathway, and biosynthesis of specific compounds, including naringenin, flavenol glucosides, and 405 sesquiterpenoids (Table S10). This provides evidence supporting a direct link between the SDR 406 and synthesis of these compounds, by an as-yet unknown mechanism. 407

Small RNA Regulation of Sex Dimorphism 408
Of the 45 identified miRNA loci, 18 were male differentially expressed and 7 were female 409 differentially expressed. Among the putative targets of these miRNAs were dicer-like genes, 410 squamosa promoter-like genes, and auxin response factors, and transcription factors (Table 1)  While all copies of both miR172 and miR156 show male biased expression in catkins, the overall 418 expression of miR156 is greater than miR172 in males. The majority of SPL genes that are 419 targeted are female upregulated in S. purpurea, including an SPL4 homolog on chr07 420 (Sapur.007G123900) that also shows increased methylation in males in its promoter region. 421 These data support that the miR156/172 pathway is upregulated in male catkins and may be 422 responsible for sex dimorphisms (Table S11) Genes with eQTL associations to the sex determination region are outlined in red boxes and 431 represent possible entry points for SDR regulation of the flowering pathway. 432 Intriguingly, one copy of miR156 is located in the SDR region unique to chr15Z. Alignment of 434 the chr15Z miR156 precursor sequence to chr15W reveals a single indel is responsible for this 435 chr15Z-specific mapping ( Figure S5), which may prevent transcription or processing of this 436 small RNA from chr15W. This could indicate a dosage dependent response in males, which have 437 four copies of this mature miR156 homolog, compared to only three in females. Female 438 upregulated miRNAs included miR403, which targets AGO2, and miR162, which targets a dicer-439 like gene. All of these are involved in small RNA signaling and DNA methylation, which is 440 consistent with the enrichment of transcription regulation terms in the female upregulated genes. 441 This may point towards a role of genome-wide DNA methylation or RNA silencing in regulating 442 sex dimorphism, possibly mediated by AGO4 or DRB1. 443

Conclusion 444
Taken together, RNA-Seq, small RNA-Seq, bisulfite sequencing, and eQTL mapping results 445 suggest that the sex determination mechanism in willow is different from the small RNA-446 mediated, single-gene sex determination pathway proposed in Populus, and reveals that many 447 miRNAs, transcription factors, floral development genes, and secondary metabolism genes are 448 involved in primary and secondary sex dimorphism pathways in Salix purpurea catkins. A major 449 factor inhibiting functional genomics research in Salix is the lack of a functional transformation 450 system, which precludes direct assessment of gene function in transgenic plants. We show here 451 that using a combination of genome-wide differential expression, co-expression, eQTL, small 452 RNA, and genomic variation can provide strong evidence supporting the involvement of these 453 proposed candidate master regulator genes in sex determination. 454

Plant material and growing conditions 456
The F1 family (Family 82) was generated from a cross between female S. purpurea clone 94006 457 and male clone 94001, both collected from different naturalized S. purpurea populations in 458 Upstate NY (Fig. S6)  values was generated using hclust, which showed that two male samples (10X-317-161 and 10X-495 317-020) displayed extreme differences in overall expression patterns from all of the other 496 samples (Fig. S8), and so were removed from downstream analyses. A topological overview 497 matrix was created to perform module identification. In order to correlate sex with module 498 expression, sex was coded numerically with values of -1 for females and 1 for males, thus a 499 module r 2 near one indicates strong correlation with males, and a module r 2 near negative one 500 indicates strong correlation with females. The minimum module size threshold was set at 30 501 genes per module. Identified modules were merged based on similar expression levels. A total of 502 32 unmerged modules were merged into 24 modules based on similar expression patterns. 503

Gene Ontology Enrichment and Pathway Analysis 504
Gene ontology (GO) enrichment at the lowest levels of the GO hierarchy was calculated in the 505 BiNGO app 50 in Cytoscape using the Arabidopsis homolog gene models. Arabidopsis gene 506 models were used as there is limited data available on gene function in Salix, and gene functional 507 annotations of non-model species are often ultimately based on experimental evidence from 508 Arabidopsis. Furthermore, using Arabidopsis gene models provided compatibility with existing 509 databases, such as TAIR. The background reference was created from the list of Arabidopsis 510 homologs for every transcript model in the S. purpurea 94006 v5.1 genome found in the gene 511 annotation information file on Phytozome (https://phytozome-512 next.jgi.doe.gov/info/Spurpurea_v5_1). The Arabidopsis homologs of the genes present in each 513 module or male:female differential expression cutoff (inclusive of all genes above or below 514 cutoff) were used as the query to determine overrepresented and underrepresented parent and