In silico analysis of key regulatory networks related to microfibril angle in Populus trichocarpa Hook.

Dissection of regulatory network that control wood structure is highly challenging in functional genomics. Nevertheless, due to the availability of genomic, transcriptomic and proteomic sequences, a large amount of information is available for use in achieving this goal. MicroRNAs, which compose a class of small non-coding RNA molecules that inhibit protein translation by targeting mRNA cleavage sites and thus regulate a wide variety of developmental and physiological processes in plants, are important parts of this regulatory network. These findings and the availability of sequence information have made it possible to carry out an in silico analysis to predict and annotate miRNAs and their target genes associated with an important factor affecting wood rigidity, microfibril angle (MFA), throughout the Populus trichocarpa Hook. genome. Our computational approach revealed miRNAs and their targets via ESTs, sequences putatively associated with microfibril angle. In total, 250 miRNAs were identified as RNA molecules with roles in the silencing and post-transcriptional regulation of the expression of nine genes. We found SHY2, IAA4 (ATAUX2–11), BZIP60, AP2, MYB15, ABI3, MYB17, LAF1 and MYB28 as important nodes in a network with possible role in MFA determination. Other co-expressed genes putatively involved in this regulatory system were also identified by construction of a co-expression network. The candidate genes from this study may help unravel the regulatory networks putatively linked to microfibril angle.


Introduction
Wood is a biological material made by heterogeneous cells, and its structure evolved to provide rigidity, to withstand pressure from the mass of the tree and to respond to different environmental factors (Plomion et al. 2001). These characteristics that are beneficial for woody plants also have important commercial consequences for the lumber industry. For this reason, during breeding of different tree species, efforts have been made to improve tree quality.
The mechanical properties of wood are influenced by a number of factors, including anatomical structure, the walls of constituent cells, the adhesion between them and their turgidity (Ryden et al. 2003). In turn, the mechanical characteristics of plant cell walls are largely controlled by the structure of the cytoskeleton. In terms of mechanical properties, while microtubules and actin filaments are of paramount importance (Fletcher and Mullins 2010), within the cell wall, the orientation of cellulose microfibrils is the most important factor (Bidhendi and Geitmann 2016). A number of studies have shown that cellulose microfibril orientation is the predominant mechanical property of both primary and secondary cell walls of xylem (Barnett and Bonham 2004;Courchene et al. 2006); the angle between cellulose microfibrils and the long axis of the cells is estimated to account for merely 86% of the variation in wood rigidity within the same species (Evans and Ilic 2001).
By the mid-1980s, a general agreement had emerged that the orientation of microfibrils and microtubules is mediated by auxin (Taiz 1984;Bergfeld et al. 1988). Investigations supporting the hypothesis that auxin controls the orientation of newly deposited microfibrils have been conducted on the basis of electron micrographs (Bergfeld et al. 1988) and via approaches based on the analysis of a series of auxinresistant mutations (Leyser et al. 1993;Del Pozo et al. 1998;Del Pozo and Estelle 1999). It is clear that auxin induces rapid changes in the orientation of microtubules in growing cells (Heisler and Byrne 2020;Winnicki 2020). Although it is already known, that the orientation of plant microtubules follows geometrical and mechanical cues (Durand-Smet et al. 2020), the following question remains unanswered: is the orientation of microtubules determined by the effector signal directly or by changes in cell growth? Despite the answer given by Chen et al. (2014) that the inhibition of cell expansion is mediated by the rapid ABP1-dependent action of auxin on microtubules, the molecular mechanism underlying the microfibril angle is still unclear. Since there is evidence for a causal relationship between wall-relaxing processes and the rapid regulation of cell growth by auxin (Cosgrove 2005;Perrot-Rechenmann 2010), data concerning its direct effector signal are still scarce.
MiRNAs regulate the expression of multiple genes underlying agriculturally important traits in plants. Therefore, miRNAs are potential targets for molecular plant breeding programmes (Mondal and Ganie 2014;Min et al. 2017). miRNAs are able to directly silence complexes because of their complementary target sites, regulating the expression of target genes via translational repression or endonucleolytic mRNA cleavage (Bartel 2004). In addition to their functions in organ development, plant miRNAs might be involved in various stress responses, such as dehydration, mineral nutrient deficiency and mechanical stress (Yang et al. 2007), and may function in defence systems associated with structural and mechanical fitness (Lu et al. 2005).
Although miRNA sequences are not conserved between animals and plants, they are highly conserved within different kingdoms, and this evolutionary conservation is one of their most important characteristics (Ambros et al. 2003). Mature plant miRNAs pair with their target sites via complementarity, recognizing the target site in the coding region and guiding the mRNA for cleavage (Bartel 2004). Based on this high degree of complementarity, to confirm the diverse roles of miRNAs in plants, many methods have been developed. By synthesizing bioinformatics, genetics and molecular biology applications, researchers have characterized some miRNAs and validated their targets, and their functions have been confirmed (Dugas and Bartel 2004;Yang et al. 2007).
The initial objective of this study was to identify by in silico analysis the mature miRNAs and their target genes that have possible roles in determining microfibril angle in P. trichocarpa using available transcriptomic (ESTs, miR-NAs), genomic (P. trichocarpa) and plant protein resources. To provide support for the use of de novo markers targeting conserved regions putatively linked to microfibril angle across a gymnosperm and an angiosperm species, as input data, expression profiles related to juvenile wood from a gymnosperm (Pinus radiata L.) were used. During the analyses, while the regulatory functions of miRNA target genes were predicted, our attention focused on the interaction enrichment of auxin-response target genes, indicating a biologically associated network. Therefore, our study also provides bioinformatic data concerning the role of auxin in the reorientation of microfibrils and microtubules and the regulation by miRNAs of genes mediated by auxin, facilitating an improved understanding of the regulatory networks inherent to the microfibril angle.

Data sources
A list of identified candidate genes for low and high microfibril angles (MFAs) was compiled based on the work of Li et al. (2012) and is summarized in Table S1. Based on this list, an EST library containing parts of putative candidate genes likely influencing juvenile wood microfibril angle was constructed with publicly available ESTs (103,038 sequences of eudicots), accessed and retrieved from the NCBI database (https:// www. ncbi. nlm. nih. gov/). Since Li et al. performed their experiment on Pinus radiata L., and because the taxonomic distance may be related to structural differences in these genes between the two species, to determine the largest number of those in the Populus trichocarpa genome, our EST database was compiled from all publicly available sequences identified from eudicots. The latest P. trichocarpa reference genome assembly (RefSeq assembly accession: GCF_000002775.4) was downloaded from the National Center for Biotechnology Information database (Tuskan et al. 2006). A mature miRNA library was created with all currently available mature miRNA sequences (10,414 sequences) from all plant species within from the miRBase data repository (Kozomara et al. 2019; http:// www. mirba se. org/). A protein BLASTX database comprising 760,913 sequences from all plant species was also generated with protein sequences accessed and retrieved from the NCBI database.

Bioinformatics data mining
The complete scheme followed is shown as a flowchart, presented in Fig. 1.

EST sequence processing and assembly with EGassembler
The EST library was subjected to sequence cleaning and assembly using EGassembler (Masoudi-Nejad et al. 2006; https:// www. genome. jp/ tools/ egass embler/). The sequence cleaning was accomplished enforcing the following thresholds: 96% minimum identity for an alignment with a contaminant and a minimum length of 11 of a terminal vector hit to be considered. The repeat masking process was carried out with a cut-off score for masking repeats of 225 using the Arabidopsis Repbase repeats library. By this method, all repetitive sequences were masked, trimmed and discarded. The vector sequences were masked by setting the minimum length of matching words to 10, the minimum alignment score to 20, the number of potential vector bases at the beginning of each read to 80, the maximum length of matching words to 20 and the level of masking for reporting a match to 80 using the NCBI vector library. Sequences having similarity to those of plastids and mitochondria were masked and discarded by setting the same parameters using the Arabidopsis thaliana (L.) Heynh. plastids and plant mitochondrial library. The Fig. 1 Diagrammatic representation of the most important bioinformatic steps of the identification of key regulatory networks in Populus trichocarpa Hook. As entry, the EST library was filled with eud-icot EST sequences based on the P. radiata genome information from Li et al. (2012) 1 3 remaining sequences were considered high-quality ESTs and were assembled into contigs and singletons with a base quality cut-off for clipping of N > 12, a max gap length in any overlap of N > 20, an overlap length cut-off of N > 40, an overlap percentage identity cut-off of N > 80, an overlap similarity score cut-off of N > 900 and a clipping range of N > 250.
Searching for sequence similarities in the P. trichocarpa genome A search of contigs and singletons generated from highquality sequences of eudicot ESTs against the P. trichocarpa genome was performed by the NCBI BLAST+ BLASTN package in the Linux environment; the contigs and singletons were searched as nucleotide queries and against the genome as a nucleotide database. For this purpose, with the P. trichocarpa contigs (FASTA sequences served as inputs) used as the first step, a BLAST database was generated using makeblastdb, applied with no additional masking data. The BLASTN application was performed by setting shorter word sizes and using the best-hit-filtering algorithm, searching for only the best matches for each query reporting a match. The BLASTN search results were displayed in tabular format.

Extraction of BLAST output sequences
Following the BLAST activity, sequences using the tabular BLAST files (as a result of the similarity search of the contigs and singletons) were extracted with a basic shell script executable in the Linux environment SeqExtractor version 1.1 (Pereira et al. 2017), producing FASTA files through the tab files containing the IDs of the BLAST results.

MiRNAs and their target predictions
The previously best-matched BLAST contig and singleton sequences were searched against all mature miRNAs from plant species available in the miRBase data repository. MiR-NAs and their potential target sites in P. trichocarpa were identified using a two-step strategy. First, miRanda (Enright et al. 2003) was used to provide an algorithm for identifying genomic targets for microRNAs. The second phase constituted the validation of the miRanda results, which were identified by a subsequent BLASTX search of coding sites for proteins involved in the molecular control of microfibril angle. The validation process was accomplished by using another plant small RNA target analysis server, psRNATarget (Dai et al. 2018), with all the parameters set to the defaults.

Prediction of proteins involved in cellulose microfibril angle determination
Before the aforementioned analysis performed by psRNATarget, the putative proteins involved in the molecular control of microfibril angle and whose biosynthesis is regulated by miRNAs were identified and selected by mining the reference protein library of all plant species using the basic local alignment search tool BLASTX algorithm, with the previous miRanda-detected miRNA target genes used as queries. For this purpose, with all plant protein libraries (FASTA sequences as inputs) compiled serving as the first step, a BLAST database was generated using makeblastdb, performed with no additional masking data. The BLASTX activity was performed by the Galaxy/Europe NCBI BLAST +s BLASTX tool (Camacho et al. 2009;Cock et al. 2015), with the following parameters: standard query genetic code, E-value of 1e-3, filtering of lowcomplexity regions, 0 maximum hits to consider, 1 maximum HSP to retain for any single query-subject pair and 99-100 minimum query coverage per HSP. The BLASTX search results were displayed in tabular format.

Gene ontology annotation
GenBank identifiers of the protein sequences were converted into UniProt identifiers using g:Profiler (Raudvere et al. 2019). To understand the biological meanings of the output of the BLASTX search result list, GO annotation was performed using PANTHER version 14 (Mi et al. 2019), the String web server (Snel et al. 2000) and the ShinyGO (Ge et al. 2020) web application. The enriched biological themes, particularly GO terms, the enriched functional-related gene groups, the redundant cluster annotation terms pointing out the interactive proteins and the highlighted protein functional domains and motifs were identified. The identified targeted genes were formed into a network using the GeneMANIA (Warde-Farley et al. 2010) prediction server, revealing functionally similar genes with predictive value from our dataset. Enriched regulatory interactions in the network analysis were determined and visualized with Cytoscape version 3.8.0 (Shannon et al. 2003). Identification and annotation of genetically mobile domains of proteins and analysis of domain architectures were accomplished via Simple Modular Architecture Research Tool (SMART) version 8.0 (Letunic and Bork 2018).

EST sequence processing and assembly, sequence similarities in the P. trichocarpa genome
Out of the 103,172 eudicot ESTs assembled in our library, cleaning resulted in 103,079 high-quality sequences, while the remaining 93 sequences were discarded. EGassembler generated 12,947 contigs and 16,459 singletons. Among all the ESTs, 83% were grouped into contigs. The average contig length detected was 907 bp -the shortest being 101 bp, and the longest being 3275 bp. The NCBI BLAST+ BLASTN search of all the contigs and singletons generated from high-quality sequences of the eudicot ESTs in the P. trichocarpa genome revealed 2870 sequences (from singletons) and 3612 sequences (from contigs).

MiRNAs and their target genes
The identified P. trichocarpa sequences that were compared with all mature miRNA sequences from all the plant species in miRanda resulted in 1,538,816 miRNAs vs. P. trichocarpa target site combinations. Although many sequences occurred multiple times in different combinations, sequences from all the combinations were selected as inputs for the next step of analysis, considering that, according to the literature (Moss and Tang 2003;Wu et al. 2010), one miRNA can exert a regulatory effect on several targets and that the expression of a specific gene can also be regulated by multiple miRNA sequences.

Proteins putatively involved in cellulose microfibril angle determination
The P. trichocarpa sequences from all miRNA-target site combinations were used to mine all the plant protein references in the library. The BLASTX search with the minimum query coverage per HSP set to 99-100 revealed nine proteins [auxin-responsive protein IAA3 (SHY2), BZIP transcription factor 60, an AP2-like ethylene-responsive transcription factor, transcription factor MYB15, transcription factor LAF1, transcription factor MYB17, B3 domain-containing transcription factor ABI3, auxin-responsive protein IAA4 (ATAUX2-11), transcription factor MYB28] putatively involved in the molecular control of microfibril angle and whose biosynthesis is regulated by miRNAs (Table 1). All these proteins were identified in A. thaliana.

MiRNA results validation
From all the miRanda results, by filtering and selecting the previous nine protein coding sequences and their regulatory miRNAs, we identified 7436 miRNAs as RNA molecules with roles in silencing and post-transcriptional regulation of the expression of these genes (Table S2). However, considering the questionable suitability of miRanda for predicting targets in species other than Arabidopsis, we checked all the data with psRNATarget. Stringent searching with psR-NATarget yielded congruent but fewer microRNAs, which were accepted as the final results for this analysis (Table S3).

Gene ontology annotation results
As a result of the ontology association performed by PAN-THER, we determined the molecular functions, biological processes and the cellular components in the case of each protein. According to these results, auxin responsive protein IAA3 (SHY2) participates in the auxin-related pathway, plays a role in the regulation of transcription in response to auxin, and is auxin templated. The transcription factor MYB17 (a MYB41 orthologue) has a role in cell differentiation and the regulation of flower development. BZIP transcription factor 60 is involved in immune system processes in response to different stresses. The transcription factor LAF1 is part of the far-red light signalling pathway, the transcription factor MYB15 has a role in cell differentiation and in stress responses, and the auxin-responsive protein IAA4 (ATAUX2-11) is part of the auxin signalling pathway. The B3 domain-containing transcription factor ABI3 also regulates transcription in response to auxin. The transcription factor MYB28 is also Table 1 Identified proteins putatively involved in the molecular control of microfibril angle and whose biosynthesis is regulated by miR-NAs. The nine contigs were selected from the assembled EST library, following the BLASTN search of all contigs and singletons against the P. trichocarpa genome. The proteins resulted from mining the reference protein library of all plant species by BLASTX algorithm, using these P. trichocarpa sequences as query. All proteins were identified in A. thaliana

Contig number
GenBank a transcriptional regulator that acts in response to different biotic and abiotic stresses, and AP2 is part of the ethyleneactivated signalling pathway (Table S4). By performing a functional classification of the resulting BLASTX protein list using PANTHER and by determining the list of supported organisms in Arabidopsis, we identified 13 GO classes (Table 2). Kappa statistics quantitatively measuring the degree of agreement regarding how genes share similar annotation terms resulted in 13 biological modules, of which 11 presented a kappa value of 1 (Table 2); their functions were related to the regulation of transcription, DNA-dependent regulation of transcription, the nucleus, transcriptional regulator activity, transcription factor activity, DNA binding, biological regulation, regulation of RNA metabolic processes, responses to stimuli and responses to organic substances. There were three biological modules with a kappa value of 0.4210, and they were related to responses to endogenous stimuli, cellular processes and DNA binding (Table 2).
To provide stronger support for the results of the previous Gene Ontology analysis, our protein list was annotated by STRING (Table S5). This annotation produced additional data, and according to the results, SHY2, in addition to its regulation of transcription, plays a central role in auxininduced regulation of root growth, gravitropism, and lateral root formation. BZIP60 transcription factor 60 can activate the promoters of chaperones in the endoplasmic reticulum, AP2 binds to the GCC-box pathogenesis-related promoter, and MYB15 may increase the expression levels of genes involved in abscisic acid (ABA) biosynthesis and signalling. ABI3 is a member of the AP2/B3-like transcription factor family, MYB17 contributes to meristem identity transition, and LAF1 activates the expression of light-induced genes. IAA4 (ATAUX2-11) functions by interacting with auxinresponse factors (ARFs) that bind to auxin-responsive promoter elements, while MYB28 acts in stress responses by promoting aliphatic glucosinolate biosynthesis and repressing indolic glucosinolate biosynthesis.

Regulatory networks putatively inherent to microfibril angle
STRING analysis previously predicted an association between (IAA4) ATAUX2-11 and SHY2 (auxin-responsive protein IAA3) (both of which are AUX/IAA transcriptional regulators) based on protein homology, co-expression, text mining and experimental evidence. Similarly, a connection between ABI3 and MYB15 (both of which are transcription factors) was predicted by text mining and co-expression evidence. As shown in Fig. 2a, among all nine proteins identified, only the aforementioned four were associated. For each of the proteins detected as fitting into a network (Fig. 2b), we increased the number of additional proteins by 51.
To summarize the correlations among significant pathways, our proteins were hierarchically clustered by ShinyGO v0.61 into two major clusters (Fig. S1), with those associated with the response to endogenous stimulus and to hormonerelated terms being grouped into the branch at the bottom. The lower branch was characterized by a greater number of terms, the most significant number of which were related to regulation.
Given the wealth of genomic and proteomic data provided by GeneMANIA, the generation of hypotheses about our genes resulted in a functional genomic dataset, which is summarized in Table 3. Notably, among the FDR values of all the functions, the highest FDR value was related to the response to auxin. On the basis of their FDR values, the other functions were as follows (in descending order): negative regulation of nitrogen compound metabolic processes; DNA-templated negative regulation of transcription; negative regulation of nucleobase-containing compound metabolic processes; responses to gravity, tropism, and gravitropism; root system function; and root development. Considering these results and assuming that our selected genes are part of a protein complex, a Gene-MANIA network of interactions was constructed (Fig. 3). Figure 3 deserves special attention, since 21 of the 29 proteins from this protein complex are related to the response to auxin, six are related to root system development, four are related to the response to tropism, and three are related to the negative regulation of transcription and biosynthesis and metabolic processes.

Domain architectures of the proteins putatively involved in cellulose microfibril angle determination
Using SMART, we identified nine select low-complexity regions of the proteins. BZIP60, AP2 and MYB17  have low-complexity regions positioned towards protein sequence extremities. In the case of LAF1, MYB15, and MYB28, these regions were centrally located, while SHY2 and IAA4 (ATAUX2-11) presented no regions with low complexity. Interestingly, ABI3 presented regions with low complexity along the entire length of the sequence. Our results concerning LCR positions in each protein sequence are presented in Fig. 4.

Discussion
Our assembly process resulted in different numbers of contigs and singletons, all of which were then searched within the P. trichocarpa genome. The sequences remaining as singletons could be considered genes expressed at low levels, while contigs with different lengths (between 101 and 3275 bp) presumably indicate expressed genes with different sizes. Taking into consideration that all our contigs and singletons are derived from EST sequences in conjunction with the microfibril angle and that the reorganization of the microfibrils is tissue specific and is affected by different environmental stimuli (Kutschera 1989), these differences in length can be deduced as being associated with different expression levels.
The use of miRanda resulted in P. trichocarpa nucleotide sequences revealing more than 1.5 million miRNA-P. trichocarpa target site combinations. This is not unexpected, since many P. trichocarpa targets and many miRNA sequences occurred multiple times in these combinations in each contig. We also found that not only miRNAs linked to multiple target sites but also target sequences were subjected to several miRNAs, all under the minimum free energy threshold. Of these, BLASTX searches revealed nine proteins involved exclusively in the molecular control of the microfibril angle. The filtering and selection of the regulatory miRNAs of these nine proteins revealed 7436 miRNAs (identified by miRanda) and 395 (identified by psRNATarget) as RNA molecules with roles in silencing and post-transcriptional regulation of the expression of these genes. These findings lead us to draw two main conclusions: first, one miRNA can exert its regulatory effect on several targets, and second, the expression of a specific gene can be regulated by multiple miRNA sequences. Cooperative translational control, namely, the phenomenon by which multiple miRNAs can target the same gene, was first predicted by Krek et al. (2005) and Lewis et al. (2005) and then confirmed experimentally by Wu et al. (2010), suggesting that neither a specific miRNA nor the combination of miRNAs Fig. 4 Results of sequence and structure analysis of the identified nine proteins by SMART (Letunic and Bork 2018). The low-complexity regions for each protein are marked in grey. The structural homology models are presented on the right 1 3 is a determinant of target gene expression. Target multiplicity, read as the ability of one miRNA to have several target genes, has also been described (Enright et al. 2003) and relies on the possession of conserved miRNA binding sites (Moss and Tang 2003).
Functional classification of the aforementioned proteins revealed 13 biological modules, of which 10 presented a kappa value of 1, while three modules presented a kappa value of 0.42 (Table 2). Given that Cohen (1960) suggested the kappa results from 0.81-1.00 could be interpreted as meaning near-perfect agreement and that those from 0.41-0.60 could be interpreted as meaning moderate but acceptable agreement, we consider the inter-rater reliability of all 13 detected modules as meaning adequate agreement. As shown in Table 2, the functions of the whole set of modules are related to cellular processes involving transcriptional regulation, depending on the DNA, in response to different stimuli, in the nucleus. Consequently, the detected proteins (and their coding sequences) are part of a process that results in a change in the state or activity of the cell in terms of gene expression as a result of an organic substance stimuli. These results were also confirmed by the hierarchical tree summarizing the significant GO terms (Fig. S1). Tables S3 and S4 present the collection of data that provide more detailed information about the molecular function, biological process and cellular response of the detected proteins. By a thorough examination of these tables, we can conclude that the biological functions of all these proteins depict a frame of a network in which the nodes have an important role in transcriptional regulation during cell development and differentiation in relation to auxin-activated, mitochondrial nucleus and abscisic acid-activated signalling pathways in response to abiotic (cold, water, red or far-red light, sulfur, salt) and biotic (insect) factors. Since, according to the literature (McFarlane et al. 2014;Baskin 2015), during primary growth, cell expansion largely depends on microfibril reorientation in the primary cell wall and before secondary cell wall formation, the pattern of microfibril orientation is already fixed, and the transcriptional regulation of the development and differentiation of the cell may also have a role in microfibril angle determination. As we mentioned before, the nodes of the resulting network are related to auxin, which not only plays a central role in plant growth and development but also is an important stress signalling mediator. IAA3 regulates cell division, cell elongation, embryo polarity, vascular differentiation, apical dominance and tropic responses to light and gravity (Woodward and Bartel 2005) and has been identified as an important mediator of auxin-dependent regulation of cambial proliferation activity (Nilsson et al. 2008). However, data from the literature also confirm the direct role of auxin in determining microfibril angle. Members of the plasma membrane protein family PIN-FORMED (PIN) act in auxin efflux and correlate with cortical microtubule array orientation in Arabidopsis shoot apical meristems (Heisler et al. 2010). The fact that SHY2/IAA3 acts as a negative regulator of auxin signalling was suggested by Tian et al. (2002), and the gravitational stimulus interaction with the IAA signal transduction pathway and its effect on MFA was suggested by Hellgren et al. (2004) and Tobias et al. (2020). During these interactions, due to their position and structure, arabinogalactans act as mechanical sensors, activating important signal mediators in the cascade of events that leads to arabinogalactan gene transcription. A link between arabinogalactans and microtubules was suggested by Andeme-Onzighi et al. (2002) and Sardar et al. (2006). Iwata et al. (2008) suggested that BZIP60 transcription factor 60 activates the promoters of chaperones in the endoplasmic reticulum, and findings indicating that cortical microtubules play an important role in abiotic stressinduced signalling pathways can be found in the work of Wang et al. (2011). AP2 has conserved functions that are required for the regulation of disease resistance pathways (Gutterson andReuber 2004), andSakamoto et al. (2018) demonstrated that the expression of AP2 modifies the characteristics of primary cell walls. MYB15 may enhance the expression levels of genes involved in abscisic acid (ABA) biosynthesis and signalling, and ABI3 is part of the AP2/ B3-like transcription factor family, whose members play a role in ABA biosynthesis. A link between plant phytohormones such as ABA and turgor-mediated cell expansion, synthesis, deposition and remodelling of plant cell walls was demonstrated by Bai et al. (2012), and the fact that ethylenedependent genes may be related to cell wall composition and microtubule orientation was suggested by Seyfferth et al. (2019). Several members of the MYB family control the secondary wall synthesis regulatory network, so they may directly regulate the expression of cell wall synthesis genes . MYB17 contributes to meristem identity transition. The molecular regulatory network underlying this process interferes with the structural elements of the cell, in particular the cell wall, to induce specific morphogenic events (Traas 2019). LAF1 activates the expression of light-induced genes, and as we mentioned before, auxin regulates responses to light. According to the STRING annotation, IAA4 (ATAUX2-11) exerts its role by interacting with auxin-response factors, and MYB28 is considered a common target of PtrWNDs in the transcriptional networks of secondary cell wall formation (Zhong et al. 2010). All of the above-mentioned data from the literature suggest that the following genes identified in the present study may represent important nodes in a network, with a role in MFA determination: SHY2 and IAA4 (ATAUX2-11) (AUX/IAA transcriptional regulation); BZIP60, AP2, and MYB15 (stress signal transduction); ABI3 (responses to auxin and abscisic acid); MYB17 (activator protein 1 regulation); LAF1 (photomorphogenesis); and MYB28 (cell differentiation in response to biotic and abiotic stresses). In the network of predicted associations performed by the STRING database (Fig. 2), between ABI3 and MYB15, we found co-occurrence evidence. These results indicate that plant development can adapt in response to environmental stresses and that environmental stress-related signals are integrated into auxin regulatory networks (Tognetti et al. 2012).
However, the network of predicted associations created with the detected proteins only (Fig. 2) did not indicate clear relationships among all nine proteins. For integration into a single connected network, we have to add nodes that are connected to our proteins that have scores better than or equal to the confidence cut-off. This raises the question of why, starting from genes that can be linked to MFA according to the literature, our entire study shed light on only nine proteins partially connected to the network (shown in Fig. 2b). In connection with this, we would like to draw attention to two aspects. First, our EST sequences were from eudicots, and many of these sequences are presumably not found in the P. trichocarpa genome. Second, our search was performed only on genes (and their coded proteins), which are subject to negative regulation of miRNAs. From this point of view, our results may indicate the most important nodes at which combinations of miRNAs intervene to regulate the network determining microfibril angle.
A number of studies (Tompa 2002;DePristo et al. 2006;Coletta et al. 2010) have revealed that low-complexity regions within protein sequences have position-dependent roles. Those preferentially positioned towards the protein sequence extremities compared with those centrally located show correlations between their lengths and degrees of connectivity. Centrally located LCRs were enriched with transcription-related GO terms, and terminal LCRs were enriched with translation-and stress response-related terms. By undertaking an investigation of low-complexity regions to explore the possible functional significance of the detected nodes (proteins), we showed that LAF1, MYB28 and IAA4 (ATAUX2-11) may have transcription-related roles, since the LCRs of these proteins were centrally located. BZIP60, AP2 and MYB17, with LCRs in the sequence extremities, may be translation and stress-response related, while both ABI3 with LCRs, with LCRs across their entire sequence length, might participate in transcription, translation and stress responses. Taken together, all these findings perfectly correlate with our previous GO annotations obtained by PANTHER (Table S4). According to Coletta et al. (2010), terminal LCRs play major roles in low-specificity biological events that involve large protein complexes. Protein chaperones play a major role in the stress response and have low-specificity binding properties due to the large variety of partners they bind (Tompa and Csermely 2004), and as we found, BZIP60 can activate the promoters of the ER chaperones BiP1 and BiP. AP2 may also be involved in the regulation of gene expression by stress factors and by components of stress signal transduction pathways (Table S5). As shown in Table S5, LAF1, MYB15, MYB28 and IAA4 (ATAUX2-11) are all transcription factors with centrally located LCRs (Fig. 4.). As they have multiple binding partners, proteins with central LCRs need to be highly specific (Coletta et al. 2010), and transcription regulatory events, as in the case of transcription factors, limit the number of proteins present simultaneously (Reményi et al. 2004). Altogether, our results shed light on networks with a large spectrum of components with possible role in regulation of an important wood determinant, microfibril angle. As potential targets for molecular plant breeding programmes, the identified miRNAs, proteins and genes, as well as the interpretation and analysis of networks generated by these components, warrant further investigation.

Conclusions
Our study demonstrates that gene expression data can be applied to define a corresponding regulatory network, to identify known and new regulatory elements and to predict co-expression networks. We revealed a set of miRNA sequences regulating genes with a putative role in MFA determination in poplar and a possible role in the response to abiotic and biotic stresses. Our data are consistent with the notion of cooperative translational control, namely, that multiple miRNAs can target the same gene. Based on the specified miRNA target genes, we identified the proteins they encode and their biological functions, and we outlined a network in which the nodes may have an important role in transcriptional MFA regulation during cell differentiation in relation to auxin-, mitochondrial nucleus and abscisic acid-activated signalling pathways in response to abiotic and biotic factors. Thus, our results indicated that plant development is adaptively modulated under conditions of environmental stress and that environmental stress-related signals are integrated into auxin regulatory networks. Moreover, the results of our investigation concerning low-complexity regions to explore the possible functional significance of the detected nodes also demonstrate that the regulation of microfibril angle may involve components of stress signal transduction pathways. We established a correlative relationship between expression profiles and MFA, overall, our findings can lead to more complex insights into MFA regulation, with data for molecular plant breeding programmes in poplar.
Authors' contribution ZAK designed the research, collected and analysed the data and wrote the article. EGT made significant contributions to the data analysis, the interpretation of the results, the figures and the manuscript editing. ABE, KC and ABO contributed by providing feedback on the manuscript drafts. All authors read and approved the final manuscript.
Funding Open access funding provided by University of Sopron. This research received no specific grant from any founding agency in the public, commercial, or not-for-profit sectors.
Data availability Not applicable.
Code availability Not applicable.

Declarations
Ethics approval Not applicable.

Conflict of interest
The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.