Total RNA was extracted from fresh tissue samples from 1-year-old P. tomentosa clone LM50, including juvenile leaf, mature leaf, apical shoot meristem, phloem, cambium, developing xylem, mature xylem, and root tissue as described (Du et al., 2013a). The RNA was reverse-transcribed into cDNA using the SuperScript First-Strand Synthesis system with the supplied polythymine primers (Invitrogen, Life Technologies, New York, NY, United States) (Zhang et al., 2010a). All cDNA samples were used to investigate the tissue-specific expression of PtoNERD and NERDL through reverse-transcription quantitative PCR (RT-qPCR). The qPCR was performed on the 7500 Fast Real-Time PCR System with SYBR Premix EX Taq (TaKaRa, Dalian, China) with gene-specific primers (Supplementary Table S1) designed with Primer Express 3.0 (Applied Biosystems). The PCR conditions were as follows: 94°C for 5 min; 40 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 30 s; and a final cycle at 70–95°C for melting curve analysis to confirm the specificity of the amplification. The qPCR was conducted as described by Du et al. (2013a). All reactions were conducted with three technical and three biological replicates. The results for the various tissues were standardized to the internal control (Actin, Accession number: EF145577). Data analysis was performed with Opticon Monitor Software 3.1 (MJ Research, Bio-Rad, Hercules, CA, United States).

Total genomic DNA was extracted from fresh young leaves with a Plant DNeasy Mini kit (Qiagen China, Shanghai) according to the manufacturer’s instructions. The primer sets used for amplification and sequencing were specifically designed based on the sequence of the lncRNA NERDL and the cDNA sequence of PtoNERD. The 16,119 bp genomic DNA sequence of PtoNERD, including the flanking region (2000 bp), was obtained by direct sequencing of the P. tomentosa (LM50 clone) using conserved primers and the BigDye Terminator Cycle Sequencing kit version 3.1 (Applied Biosystems, Beijing, China) on a Li-Cor 4300 genetic analyzer (Li-Cor Biosciences, Lincoln, NE, United States). The Prosite database² was used to identify protein domains, families, and functional sites, and SMART³ was used to analyze the protein sequence. To identify SNPs, NERDL, and PtoNERD were sequenced and analyzed in 44 unrelated individuals from the association population with MEGA 4.0 and DnaSP v4.90.1 (Rozas et al., 2003). Insertions/deletions (InDels) were disregarded in the analysis. All 88 sequences have been deposited in the GenBank database (GenBank Accession No. KY344444-KY344487 for NERDL and KY344488-KY344531 for PtoNERD). Sequence alignment was performed with BioEdit, and manual editing was performed to remove primer sequences and to confirm sequence quality. All common SNPs (minor allele frequencies > 0.05) detected in NERDL and PtoNERD were genotyped in the association population using the Beckman Coulter sequencing system (Franklin Lakes, NJ, United States).

The overlapping significant SNPs were used to investigate the roles of different genotypes on the expression of NERDL and PtoNERD. For each SNP genotypic class, 10 individuals were selected and used to measure the expression of the lncRNA NERDL and PtoNERD via RT-qPCR. The differential expression across different genotypes was calculated with ANOVA. The effects of different genotypes of the shared locus on phenotypic variation were also investigated.

To identify SNPs within PtoNERD, we cloned a cDNA from PtoNERD. This cDNA is 9527 bp long, including a 5163 bp open reading frame and 364 bp of 3′UTR. Based on the alignment of the full-length cDNA sequence to the genomic sequence, PtoNERD (16119 bp) has 10 introns and 11 exons flanked by 2000 bp of promoter region, 364 bp of 3′UTR, with a 2000 bp flanking region (Figure 1A); the length of each region is listed in Table 1. PtoNERD encodes a putative polypeptide of 1720 amino acids, with an estimated molecular mass of 188.4 kD and a pI of 4.82. Characterization of the amino acid sequence of PtoNERD revealed that it contains the essential domain PHD, SWIB, Plus-3, and GYF (Figure 1A). Phylogenetic analysis indicated that PtoNERD is an ortholog of AtNERD and PtrNERD (Figure 1B). Phylogenetic analysis using the neighbor-joining method showed that PtoNERD in P. tomentosa shares a close relationship with the homologous sequences in P. trichocarpa and P. euphratica, indicating that PtoNERD is highly conserved in Populus, whereas the conservation of PtoNERD in distantly related species is weak.

To further explore the roles of PtoNERD in plant development, using genome annotation and a genome browser, we identified NERDL, an lncRNA located in the antisense strand and overlapping with the promoter region of PtoNERD, which acts as a cis-regulator of PtoNERD. Sequencing analysis revealed the 2286-nt sequence of NERDL, which shares an overlapping region of 961 bp with PtoNERD (Figure 1A). Therefore, PtoNERD was regarded as the candidate cis target gene of the overlapping lncRNA NERDL (lncRNA related to PtoNERD). Several functional motifs are found in the overlapping region, such as light responsive element 3-AF1 binding site, Box I, and GT1-motif, as well as as-2-box involved in shoot-specific expression and light responsiveness.

We explored the tissue-specific expression patterns of the lncRNA NERDL and PtoNERD in eight different plant tissues via RT-qPCR (Figure 1C). NERDL was highly abundant in phloem (1.82 ± 0.09) (arbitrary units normalized to control), followed by mature xylem (1.33 ± 0.07), with the lowest abundance detected in cambium (0.56 ± 0.03), suggesting that NERDL plays a regulatory role in wood formation. PtoNERD mRNA was enriched in juvenile leaves (4.23 ± 0.20), followed by mature xylem (3.10 ± 0.14), with the lowest abundance in cambium (0.41 ± 0.19). Overall, PtoNERD was more highly expressed than NERDL. The Pearson’s product moment correlation was 0.62 (P < 0.05), indicating a strong correlation in transcript levels between NERDL and PtoNERD. The reads numbers of NERDL and PtoNERD obtained by RNA-seq showed a similar trend in cambium, developing xylem, and mature xylem tissue (Supplementary Figure S3).

To characterize the intraspecific polymorphisms for use in association mapping, we obtained the full-length genomic sequences of NERDL and PtoNERD from 44 unrelated individuals of P. tomentosa. In total, 26 SNPs were detected in NERDL at a density of 1/88, and 142 SNPs in PtoNERD were detected at a frequency of approximately one SNP every 114 bp. The nucleotide diversity (excluding InDels) was higher in NERDL (π = 0.00257 and 𝜃w = 0.00261) than in PtoNERD (π = 0.00174 and 𝜃w = 0.00203). As expected, the average nucleotide diversity in PtoNERD was higher in introns (π = 0.00096) than in exons (π = 0.00030), since under selection pressure, the coding region of a gene shows more conservation than the non-coding regions (Table 1). Within the coding region of PtoNERD, the average nucleotide diversity of synonymous polymorphisms (d_S, π = 0.00099) was higher than that of non-synonymous polymorphisms (d_N, π = 0.00078) at a ratio of d_N/d_S (0.788 < 1), indicating that non-synonymous sites in the exons are under purifying selection.

To perform LD analysis, 24 and 135 common SNPs (minor allele frequency > 0.05) from NERDL and PtoNERD, respectively, were successfully genotyped across 435 individuals in the association population. Four common SNPs were located in the overlapping region and were shared by NERDL and PtoNERD (SNP21–24 in NERDL and SNP1–4 in PtoNERD). A plot of the r² combined with physical distance between SNPs within NERDL and the PtoNERD genes is shown in Supplementary Figure S2. In this population, the rate of LD declined rapidly, and the r² values decreased to 0.1 within 1200 and 8000 bp in NERDL (Supplementary Figure S2A) and PtoNERD (Supplementary Figure S2B), respectively, indicating that LD does not extend across the entire gene region.

Unlike microRNAs, whose roles in post-transcriptional gene regulation are well-characterized (Bartel, 2009), the biological functions of lncRNAs are still largely unknown, although studies have increasingly focusing on the functional mechanisms of lncRNAs. LncRNAs function in a variety of process in plants, such as flowering, RdDM, and responses to biotic and abiotic stress (Quan et al., 2015). LncRNAs function through various mechanisms, including mutual interactions with RNAs and proteins (Guhaniyogi and Brewer, 2001). Base pairing-dependent lncRNA–mRNA interactions are a common mode of action for such interactions. In this study, we found that the lncRNA NERDL partly overlaps with the promoter region of PtoNERD and shares a 961 bp fragment containing four common SNPs and two haplotypes (Figure 1A and Tables 2, 3), suggesting a potential correlation between this lncRNA and PtoNERD. The overlap of lncRNAs with the promoter regions of protein-coding genes can affect the expression of these genes (Hung et al., 2011). Transcript expression analysis revealed a Pearson Correlation Coefficient of 0.62 (P < 0.05), providing evidence for a strong correlation between NERDL and its putative target gene, PtoNERD (Figure 1C). LncRNAs are usually shorter and are expressed at lower levels than protein-coding transcripts in P. tomentosa (Chen et al., 2015) and in mammals (Pauli et al., 2012). Here, RT-qPCR analysis and RNA-seq showed that the expression level of PtoNERD was higher than that of NERDL and that they both exhibit tissue-specific expression in P. tomentosa (Figure 1C and Supplementary Figure S3). NERDL was expressed at the highest level in phloem, suggesting that it may play a regulatory role in secondary wood formation. Furthermore, NERDL and PtoNERD were highly expressed in xylem (Figure 1C), suggesting that they play important roles in wood formation.

To understand the forces of evolutionary change and to evaluate the precision and power of association mapping, thorough knowledge of the levels of nucleotide diversity and the extent of LD in natural populations is required (Ingvarsson, 2008; Zhang et al., 2010b). In general, we found that the levels of average nucleotide diversity were lower in coding regions than in the non-coding regions (Table 1), implying that the coding regions are conserved relative to the other regions under natural pressure, with a dN/dS ratio (0.7879) less than 1 for PtoNERD (Neale and Kremer, 2011). Indeed, synonymous mutations occurring during evolution are more likely to be fixed than neutral ones owing to purifying selection (Zhang et al., 2010b). Also, the nucleotide diversity of NERDL (π = 0.00257) is slightly higher than that of PtoNERD (π = 0.00174), further confirming the conservative nature of the protein-coding gene. SNP diversity in genes has primarily been explored in exons, introns, and UTRs, while less attention has been paid to promoter regions and lncRNAs. It would be quite useful to analyze the levels of nucleotide diversity within the promoters of numerous genes compared with other gene regions. In the current study, we found that the promoter region possessed a significantly higher frequency of SNPs than the other regions in PtoNERD (π = 0.0032; Table 1). This finding suggests that this region may be an unstable regulatory region that plays important roles during the evolutionary process (Weickert et al., 2012).

We identified 4 significant SNPs in NERDL and 26 in PtoNERD associated with D, FW, AC, HEC, and/or LC, and 15 haplotype-based associations that were strongly supported by single-marker associations (Table 3), leading to the suggestion that NERDL and its potential target might play a role in tree growth and wood formation. We identified 35 associations with both additive and dominant effects associated with five traits (D, FW, AC, HEC, and LC), representing 30 significant SNPs from NERDL and PtoNERD (Figure 2A and Table 2), suggesting that the lncRNA–mRNA pair might have a common effect on certain traits. Eight SNPs from NERDL and PtoNERD were associated with HEC, with additive and dominant effects (r², 4.10–6.54%), an association considered to indicate generalized pleiotropy (Mackay et al., 2009). Potential pleiotropy within trait categories and within different phenotypic traits has been reported in P. trichocarpa (Porth et al., 2014). Together, these findings suggest that NERDL and its putative target PtoNERD might affect HEC through a shared pathway. A similar trend was found for D, which is associated with nine SNPs from NERDL and PtoNERD (r², 4.21–6.10%). Considering the associations based on additive and dominant effects, our findings suggest that NERDL and PtoNERD co-regulate HEC and D and that the role of NERDL might be similar to that of its target, PtoNERD. Furthermore, under the additive and dominant models, several SNPs are associated with multiple traits but with different contributions, such as SNP9 in NERDL, which is associated with two traits, with r² varying from 11.28 to 13.08%, revealing the pleiotropy of specific SNPs, as observed in a previous study (Yang et al., 2015).

Three non-synonymous markers (SNP33, SNP49, and SNP75), one synonymous marker (SNP105), and 10 non-coding markers in promoter and introns were associated with FW simultaneously (Supplementary Table S2), confirming the necessity of using an integrated approach to characterize the genetic basis of wood traits. In addition, the existence of seven markers in exons suggests that these differences in PtoNERD might affect phenotypic variation by causing amino acid changes or influencing gene expression levels, whereas the five SNPs in introns might affect the regulatory roles of introns or exon splicing, thereby leading to phenotypic variation. One common SNP (L23/N3) overlapping in NERDL and PtoNERD was significantly associated with D (r² = 4.7%). Notably, genotypes of the overlapping SNP L23/N3 led to significant differences in the transcript levels of NERDL and PtoNERD (Figure 2B). The transcript levels of the AA and AG genotypes showed similar trends in both NERDL and its putative target PtoNERD, indicating that this overlapping locus can influence the expression of different transcripts. Overall, the average D was higher in the AA genotype (22.74 cm) than in the AG genotype (17.28 cm) (Figure 2C), indicating that variations in genotype may lead to changes in phenotype. Moreover, two shared common haplotypes from NERDL_SNP22-24/PtoNERD_SNP2-4 associated with D surround L23/N3. The locus might serve as a functional marker involved in the regulation of D, which should be explored in the future.

Epistasis analysis is a crucial method for deciphering the effects of genetic variation and interactions (Mackay, 2014). Epistasis analysis can address the limitations of single-locus analysis and is a convenient method for further investigating the interactions of lncRNAs and mRNA and their effects on quantitative traits (Quan et al., 2016). Complex gene interactions have a stronger effect on disease susceptibility than single genes (Chou et al., 2011). Here, we characterized nine significant lncRNA–mRNA associations involved in two traits (D and FW) and 11 unique SNPs from NERDL and PtoNERD (Figures 4D,E), reflecting the possible genetic interaction of the lncRNA NERDL with PtoNERD. Furthermore, 55.56% of these SNP–SNP pairs had negative IG values, implying that redundant information is carried in the pairs for each trait and that their similar functional modes contribute to phenotypic variation. In associations with positive IG values, the contribution of pairwise interactions to phenotypic variation is more than the sum of the contributions of two SNP individuals. The genotypic combinations with the highest and lowest IG values for FW were found for SNP pair L3-N11 (Figures 4A–C), where the value for the genotypic combinations was significantly higher than that of the single-locus genotype. The identification of positive and negative IGs contributes to our understanding of the close interactions of NERDL with its target.

Of the SNPs showing lncRNA–mRNA interactions identified under the epistasis model, L3 in NERDL formed five SNP pairs with N11, N12, N20, N98, and N135 in PtoNERD under the epistasis model associated with FW with diverse IGs, revealing the complexity of the interactions between NERDL and PtoNERD (Figure 2B). L20 also showed similar associations with four SNPs associating with D from PtoNERD (N72, N80, N81, and N135) (Figure 2C). These observations further demonstrate the important functions of L3 and L20 in NERDL and the importance of the SNP pairs under the epistasis model for phenotypic variation, and they help confirm the notion that the interaction of NERDL with PtoNERD contributes to wood formation. In addition, 3 SNPs (N11, N12, and N20) out of 11 also showed significant effects in the additive and dominant models, while the remaining epistatic effects involved pairs of SNPs that showed deficiency for both additive and dominant effects. These findings are consistent with the suggestion by Xu and Jia (2007) that epistatic variance affects genetic variance. Moreover, whether or not two SNP loci have additive and dominant effects does not influence whether they have epistatic effects. Interestingly, we detected an intragenic SNP pair, N4-N131, in which N4 had no effect on HEMC, whereas it did have pairwise effects (4.03%) on HEMC in combination with N131 (Supplementary Table S3). Epistasis analysis indicated that the lncRNA NERDL might interact with PtoNERD via similar pathways and that their complex interactions play a role in regulating tree growth and determining wood properties.

Our RT-qPCR analysis uncovered a positive correlation in expression of NERDL and its putative target (r = 0.62; Figure 1C). This finding is consistent with the results of association studies of additive and dominance effects, strongly suggesting that NERDL and its potential target affect phenotypic variation in a similar manner. NERDL and PtoNERD were more highly expressed in secondary tissues (xylem) than in primary tissues (cambium), indicating that they may be involved in secondary wood formation. The high abundance of the two transcripts in phloem and xylem and the 35 associations with growth and wood formation traits suggest that NERDL and PtoNERD together perform a regulatory function in secondary xylem formation in trees. The overlapping SNP L23/N3 also had significant effects on the transcript levels of NERDL and PtoNERD (Figure 2B). Analysis of genotype data indicated that the 15 haplotypes explained 1.06–13.16% of phenotypic variation, which was confirmed by the finding that significant SNPs are located in blocks of haplotypes associated with the same traits.