An F₁ individuals of Salix matsudana Koidz was produced by cross-breeding the male parent “9901” and the female parent “Yanjiang” in 2014 (Liu et al., 2021c). The branches of 30 female F₁ progenies and 30 male F₁ progenies were clipped at 10 cm lengths and 1 cm thicknesses and hydroponically cultured in water under two conditions, with water only (i.e., the normal condition, marked as “CK”) and with 0.5% NaCl solution (g/v) (i.e., the salinity stress condition, marked as “T”). The clipped branches were grown under each condition in three biological replications for RNA-seq at Nantong University in March 2022. After 30 days, the newly sprouted roots and shoots were collected from each replicate. The excised roots were immediately frozen in liquid nitrogen and stored at -80°C until use (Liu et al., 2021a).

The excised roots were ground in liquid nitrogen. The Plant RNA Reagent Kit (Tiangen, China) was used to extract total RNA from three replications. A Nanodrop ND 2000 spectrophotometer (NanoDrop, Thermo, Waltham, MA, USA) was then used to quantify the RNAs. The RNAs were then stored at -80°C before performing RNA sequencing. Finally, Illumina sequencing technology (Illumina, San Diego, CA, USA) was employed to perform RNA sequencing by Biomarker (Beijing, China) (Chen et al., 2020).

The transcriptome reads were processed into clean, full-length reads by removing the low-quality and adapter reads (Chen et al., 2022). The assembled Salix matsudana Koidz. (“Yanjiang”) genome sequence was selected as the reference for paired-end read mapping (Zhang et al., 2020). The clean reads were aligned to genes of the reference genome using HiSAT2 software with default parameters (Kim et al., 2019). Then, StringTie2 was used to detect new transcripts (Pertea et al., 2015). RSEM was chosen to calculate the fragments per kilobase transcriptome per million mapped reads (FPKM) by normalizing to the length of the gene and to the number of mapped reads. To identify the differentially expressed genes (DEGs), DESeq2 was selected. Two standards were used to detect the DEGs: (1) the fold change should be no less than 2 between different libraries, and (2) the adjusted false discovery rate (FDR) should be less than 0.05 (Chen et al., 2020).

The identified genes were annotated by using the BLASTx search in the NCBI nonredundant protein database. Then, Gene Ontology (GO) categorization, clusters of eukaryotic orthologous groups (KOG), Cluster of Orthologous Groups of proteins (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using BMKCloud (www.biocloud.net, version 2.0) (Liu et al., 2021a). The data that supported the findings of this study were deposited into the CNGB Sequence Archive (CNSA) of the China National GeneBank Database (CNGBdb) with accession number CNP0003818.

The root depth, root width, and root cap area of 30 male and 30 female progenies under normal and saline conditions in three biological replications were measured by using Win RHIZO TRON (Xu et al., 2020). Then, all seedlings were harvested and oven-dried to measure the biomass of the whole plant by using the student t test.

The phenotype of 30 male F₁ progenies and 30 female F₁ progenies under normal and salinity conditions were measured. Under normal conditions, the root depth and width of male and female plants did not differ significantly ( Figure 1A ). In contrast, under salinity stress, the roots of male and female plants subjected to intrasexual competition were significantly different. Under salinity stress, male plants showed deeper and wider roots than female plants ( Figure 1B ). The root cap area of males was also higher than that of females ( Figure 1C ). This result means that Salix matsudana showed sexual dimorphism under salinity stress. Unlike in Populus deltoides (Li et al., 2016), male Salix matsudana have much more developed roots than females. We then measured the whole-plant dry mass, and the average biomass accumulation of males was higher than that of females, significantly ( Figure 1D ).

To characterize the role of the response of active genes to salinity stress in males and females, deep sequencing libraries were generated using total RNA extracted from roots under normal and salinity stress conditions. After trimming off the adapter sequences and removing the low-quality reads, we obtained 19,468,575–27,771,058 clean reads for the 12 libraries; these libraries had a single read length of 90 bp and a Q30 percentage (percentage of the sequences with sequencing error rates lower than 0.1%) over 90% ( Supplementary Table S1 ). The clean reads were then mapped onto the reference genome of S. matsudana using HISAT2. In total, 44,906 (77.64% of the 57,841 gene models in the reference genome) genes were identified as being expressed in at least one library.

The DEGs under normal and salinity stress conditions in males and females were identified using a threshold FDR ≤ 0.05 and an absolute value of log2-fold change ≥ 1. A total of 4906 DEGs (2604 in males and 3227 in females) were identified, ( Figure 2A ). More genes were identified in the response to salinity stress in females than in males. However, only approximately one-fifth (925 in 4906) of the DEGs were identified as differentially expressed in both males and females. Compared to those under the normal condition, most genes under salinity stress conditions showed low expression ( Figure 2B ). Only 327 and 546 DEGs showed similar expression trends in males and females, respectively. Moreover, 52 DEGs were identified as showing opposite expression trends in males and females. Thirty-seven DEGs were up-regulated in females but down-regulated in males, and 15 DEGs were up-regulated in males but down-regulated in females. The results indicate that the mechanisms of the response to salinity stress in males and females may differ.

Weighted gene co-expression network analysis was performed to further detect common and unique salinity stress response genes in males and females. Among the 4906 DEGs, 969 could be classed into 3 modules ( Figure 3A ). In “MEbrown”, genes showed similar expression trends in males and females. These genes were up-regulated in both males and females under salinity stress, indicating that the regulation of these genes in males and females is similar ( Figure 3B ). Genes in “MEblue” showed a high correlation with female plants under normal conditions, indicating that these genes were only highly expressed in females under normal conditions and were reduced under salinity stress. In males, these genes showed similar expression levels between normal and salinity stress conditions. The genes in “MEturquoise” showed a high correlation with male plants responding to salinity stress. These genes were highly expressed in males under salinity stress and did not change much in females ( Figure 3C ).

A total of 48 genes were detected in “MEbrown”, and COG enrichment analysis was performed. These genes were found to be enriched in carbohydrate transport and metabolism and secondary metabolite biosynthesis, transport and catabolism ( Figure 4A ). KEGG enrichment analysis showed that selenocompound metabolism was enriched at high levels ( Figure 4B ). GO enrichment analysis was further performed. For biological processes, the DEGs were mostly enriched in systemic acquired resistance and carbohydrate metabolic processes ( Figure 4C ). In contrast, for molecular functions, the DEGs were mostly enriched in nutrient reservoir activity ( Figure 4E ). cell wall and apoplast were detected enriched in cellular components ( Figure 4D ). According to the common salinity stress-induced protein interaction network, several genes encoding alcohol dehydrogenase (ADH) proteins, oxygenase-related genes, glycosyl hydrolase and ethylene synthesis genes were detected as hub genes ( Figure 5A , Table S2 ). The up-regulation of these genes could help resist salinity stress (Zhao et al., 2019).

The genes responding to salinity stress that were unique to males or females were further detected. A total of 718 genes in “MEturquoise” were up-regulated in males under salinity stress. The COG enrichment analysis showed that these genes were enriched in carbohydrate transport and metabolism, secondary metabolite biosynthesis, transport and catabolism, signal transduction mechanisms and defense mechanisms ( Figure 6A ). The KEGG enrichment analysis showed that plant−pathogen interaction and starch and sucrose metabolism were significantly enriched ( Figure 6B ). According to the GO analysis, carbohydrate metabolic process and response to oxidative stress were enriched in biological processes ( Figure 6C ); membrane and extracellular region were detected in cellular components ( Figure 6D ); and many oxygen-related terms, i.e., peroxidase activity and oxidoreductase activity, were significantly enriched in molecular functions ( Figure 6E ). These results indicate that males could activate some pathways to defend against salinity stress, such as peroxidase activity, defense mechanisms, and sucrose metabolism (Jia et al., 2019). We then constructed the salinity stress-induced protein interaction network for males, and 16 hub genes were identified ( Figure 5B , Table S2 ). These hub genes were annotated as dihydrofolate reductase (DHFR), ubiquitin, glutathione transferase (GST), and DNAJ, among others ( Table S3 ). These genes were reported to reduce the damage caused by active oxygen.

In females, 203 DEGs in “MEblue” were down-regulated under salinity stress, which means that the expression of these genes may work to the disadvantage of females when resisting salinity stress. Secondary metabolite biosynthesis, transport and catabolism, and signal transduction mechanisms were detected in the COG enrichment analysis ( Figure 7A ). In the KEGG enrichment analysis, nitrogen metabolism and plant−pathogen interaction were significantly enriched ( Figure 7B ). The GO analysis also detected several enriched nitrogen-related terms ( Figures 7C–E ). These results indicate that female plants could reduce the metabolism of nitrogen to adjust to salinity stress. These genes were then used to construct the protein interaction network of females. As a result, genes encoding nitrate reductase (NR) and nitrate transporter (NRT) were detected as hub genes ( Figure 5C , Table S2 ). Nitrate is one of the most important sources of nitrogen in plants and can affect the growth and development of plants. Researchers have already found that NRT could regulate the distribution of lateral roots. The down-regulation of nitrate-related genes reduced the metabolism of nitrogen and plant growth. These results showed that female plants could reduce their metabolism to adjust to salinity stress.

Female and male plants experience different selection pressures and have different evolutionary directions. The difference in their evolution could induce a series of morphological, physiological and ecological differences between males and females, which is called sexual dimorphism. In this research, Salix matsudana of both sexes were selected to identify their sexual dimorphisms under salinity stress. To suppress the interference of different genotypes, an F₁ population was selected to perform the experiments. Males and females did not differ significantly in root depth and width under normal conditions. However, sexual dimorphism was evident under salinity stress. Male plants showed a higher tolerance to salinity stress than female plants, with deeper and wider roots and larger cap areas. Usually, males have a higher tolerance for abiotic stress because females may invest more energy into reproductive growth. Previous studies have found that under salinity stress, Populus cathayana and Populus yunnanensis males showed higher osmotic regulation ability, water use efficiency and antioxidant enzyme activity than females. This study also showed that Salix matsudana males gained an advantage under salinity stress. However, in Populus deltoides and Ginkgo, females are less sensitive to salinity stress than males. The results indicate that even within the same family or genus, sexual dimorphism may also differ. Interestingly, the sexual dimorphism under different abiotic stresses is not fixed. For example, male Populus cathayana showed an advantage under salinity (Chen et al., 2011), deficient nitrogen (Wu et al., 2021), phosphorus deficiency (Xia et al., 2020) and drought stress (Chen et al., 2014), while female Populus cathayana showed an advantage under high phosphorus supply (Xia et al., 2020). A probable reason for the changing sexual dimorphism could be that males and females have different selection pressures.

To further understand the molecular response to salinity stress in Salix matsudana males and females, we classified the DEGs into 3 module types: male-unique modules, female-unique modules, and common modules. In the common modules, most genes were annotated into the systemic acquired resistance term. ADH and oxygenase-related genes were identified as hub genes. These genes were reported to regulate plant resistance to abiotic stress. ADH has been reported to take part in abiotic stress responses, such as those to cold stress (Davik et al., 2013; Song et al., 2017), drought stress (Senthil-Kumar et al., 2010), salinity stress (Shi et al., 2017) and flooding stress (Liu and Adams, 2007; Komatsu and Ahsan, 2009; Komatsu et al., 2012). In Arabidopsis, overexpression of AtADH1 could increase the accumulation of soluble sugar and produce a stronger salt tolerance phenotype than the wild type (Shi et al., 2017).

The unique modules in males and females showed different expression trends. In males, unique response genes were up-regulated under salinity stress, while in females, most genes showed down-regulated expression trends. In males, many abiotic stress-related pathways were activated, such as starch and sucrose metabolism, peroxidase activity, and oxidoreductase activity. The hub genes, i.e., DHFR, GST, DNAJ, and ubiquitin, were annotated as being able to reduce the harm to the roots caused by oxidation. Previous studies found that DHFR could affect the content of chlorophyll (Van Wilder et al., 2009), and the suppression of DHFR could down-regulate cell proliferation and lead to cell death (Assaraf, 2007; Zheng, 2009). GST can expel oxygen free radicals from cells and reduce the damage caused by stress (Csiszar et al., 2014). In Arabidopsis, the expression of GST could increase salt tolerance and accelerate plant growth (Qi et al., 2010). In Solanum lycopersicum, DNAJ could reduce reactive oxygen species accumulation and enhance the tolerance to cold stress and heat stress (Kong et al., 2014). The overexpression of DNAJ could increase the tolerance to salt in Arabidopsis (Bekh-Ochir et al., 2013). In tobacco, the overexpression of ubiquitin could improve the resistance to cold, high salt and drought (Guo et al., 2008).

In contrast, unique genes in females were down-regulated under salinity stress. Most of these genes were enriched in nitrogen metabolism-related pathways. The hub genes were also annotated to encode nitrate reductase (NR) and nitrate transporter (NRT). The low expression of nitrogen metabolism-related genes could decrease the cell viability and growth of plants (Tabata et al., 2014). Previous studies also found that nitrate-related proteins could be involved in the plant response to salt stress. Under salt stress, the expression of NRT1.5 decreased to reduce the transport of NO₃ ^- to the bud and prevent harmful Na⁺ from entering the bud and causing injury to plants (Lin et al., 2008). By comparing the unique modules in males and females, we could develop a hypothesis: males preferred to activate abiotic stress response genes to adjust to salinity stress, while females preferred to reduce their basic nitrogen metabolism and regulate the transport of NO₃ ^- to decrease the harm caused by salinity stress. This hypothesis could explain the phenomenon of male Salix matsudana plants showing higher tolerance to salinity stress than female plants.