These RNA-seq data are available at the NCBI BioProject database with accession number PRJNA627111. The transcriptomes of 11 eyelid skin samples, five samples from 100% pigmented animals, and six samples from 0% pigmented animals were analyzed (Jara et al., 2020). The eyelid transcriptomes were generated using poly-A capture and strand-specific RNA sequencing in an Illumina HiSeq 2500 sequencing system. A total of 542,751,474 (38.5 G) clean reads were obtained. These reads were mapped to the latest cattle (Bos taurus) reference genome ARS-UCD1.2, using the software Hisat2 (v2.1.0) (Kim et al., 2015). The overall mapping rate ranged from 91 to 93%, and only uniquely mapped reads were considered (Supplementary Table S1). Transcriptomes were assembled for each sample using Cufflinks, and then, all the assemblies were merged into one using Cuffmerge.

Potential lncRNAs were identified using successive filters starting from all obtained transcripts. First, transcripts that presented the class code “ = “, “e”, “p”, and “c” from the output generated by Cuffcompare were filtered out (You et al., 2019). Transcripts with less than two exons with low expression levels (ten mapped reads per sample in at least five biological replicates were defined as the expression threshold) and shorter than 200 bp were also filtered out. DNA sequences of the transcripts were extracted using gffread (Trapnell et al., 2010). The sequences of known lncRNAs were downloaded from two multispecies lncRNA databases, ALDB (Domestic-Animal LncRNA Database) (Li A. et al., 2015) and NONCODE (Zhao Y. et al., 2016), which contain 8,250 and 23,515 cattle lncRNAs, respectively. BLASTn (version 2.2.31) (Altschul et al., 1990) was used to align the unannotated transcripts to lncRNAs in these two databases using stringent parameters (e-value ≤ 1 × 10⁶, coverage and identity ≥90%). Accurately aligned transcripts against either ALDB or NONCODE (5.0) were regarded as known lncRNAs. On the other hand, the transcripts that did not align with sequences in either ALDB or NONCODE were considered putative novel lncRNAs. This pipeline was adapted from Konsinska-Selbi et al. (2020).

We calculated the coding potential of each transcript using three complementary tools: Coding Potential Calculator 2 (CPC2, version 0.1) (Kang et al., 2017), Coding-Potential Assessment Tool (CPAT, version 2.0.0) (Wang et al., 2013), and Predictor of long noncoding RNAs and messenger RNAs based on an improved k-mer scheme (PLEK, version 1.2) (Li et al., 2014). PLEK and CPC2 are based on the same support vector machine-based (SVM) classification model, while CPAT is based on a logistic regression (Wang et al., 2013; Antonov et al., 2019). All these programs are alignment-free tools and have been proven to be highly effective in discriminating lncRNAs (Han et al., 2016).

The CPAT program estimates coding probability scores. The optimum cutoff value for protein-coding probability is species-specific, and hence, the CPAT was trained using a set of 10,000 known bovine protein-encoding transcripts (Bos_taurus.ARS-UCD1.2.cds.all.fa, version 95) and a set of 10,000 bovine noncoding sequences (larger than 200 bases). The final reference dataset of noncoding RNAs comprised 37,695 sequences (5,930 from Bos_taurus.ARS-UCD1.2.ncrna.fa (version 95), 23,515 from NONCODEv5_cow.fa, and 8,250 from ALDB.cow.lincRNAs.v1.0.fa). Bovine protein-coding transcripts and noncoding sequences were extracted randomly from each annotation, following previously published studies (Billerey et al., 2014; Gupta et al., 2019). In brief, the two training sets were randomly split into ten different parts to perform a 10-fold cross-validation analysis. The cutoff value was selected to maximize specificity and sensitivity. On the other hand, PLEK uses a sliding window-based approach to analyze the transcripts based on a k-mer frequency distribution. PLEK was trained using the same set of sequences that were used for training the CPAT program.

Transcripts that displayed a CPC2 score lower than 0.5, a CPAT score lower than or equal to 0.36, and a PLEK score lower than 0 were considered noncoding genes and were used in subsequent analyses.

Differentially expressed genes were identified using the R package DESeq2 (version 1.18.1) with default parameters (Love et al., 2014). Putative novel lncRNAs, known lncRNAs, and annotated protein-coding genes were all included in this statistical analysis. Only genes with at least ten reads per sample in at least five biological replicates were considered. Genes with an adjusted p-value ≤ 0.05 (Benjamini and Hochberg, 1995) and a |log₂FC| ≥ 1.5 were considered differentially expressed (DEs) between pigmented and unpigmented samples.

The analyses of the different sequence characteristics were performed on 14,361 coding genes (34,447 transcripts), 4,937 putative novel lncRNAs, and 218 known lncRNAs. The GC content and the length of transcripts were obtained using the infoseq function of the EMBOSS package version: 6.6.0.0. The exon number was estimated using custom scripts in bash. The abundance was estimated using Cuffnorm. The minimum free energy (ME) was calculated using the RNAfold program included in the ViennaRNA package version 2.4.9 (Zuker and Stiegler, 1981; Lorenz et al., 2011). To make the ME value of the different RNA sequences comparable, we normalized the ME by the sequences’ length, yielding the MEN. Thus, MEN was calculated by dividing the ME value by the transcript’s size and multiplying it by 100. Thus, the MEN value relates the ME estimation to a segment of 100 nucleotides: MEN = (ME/sequence length)*100 (Zhang et al., 2006).

Two strategies were used to study the association between lncRNAs and target genes acting in cis and trans. In the first case, all “neighbor” protein-coding genes showing differential expression (DEGs, p-value ≤ 0.01, and |log₂FC| ≥ 1) were identified in a range of 300 kb upstream and downstream of differentially expressed lncRNAs. Genes that showed significant correlation coefficients at the expression level with neighboring lncRNAs were considered cis target genes (Pearson, r ≥ |0.60|, p-value ≤ 0.05).

For potential associations in trans with differentially expressed protein-coding genes (DEGs, p-value ≤ 0.01 and |log₂FC| ≥ 1), sequence complementarity with differentially expressed lncRNAs (p adjust ≤0.05 and |log₂FC| ≥ 1.5) was analyzed using LncTar (Li J. et al., 2015). Note that a more stringent threshold was used to detect DE lncRNAs than to detect DE protein-coding genes to obtain a more comprehensive and diverse sample of genes and their functions. The program LncTar identifies potential lncRNA targets by finding the minimum free energy of lncRNA and mRNA pair joint structures. This program was run with a threshold value of ndG −0.08. All predicted interactions of lncRNAs with target mRNAs that showed a significant coexpression correlation (Pearson, r ≥ |0.60|, p-value ≤ 0.05) were kept for further analysis.

The functional roles of the target genes were evaluated using the R package EnrichKit (https://github.com/liulihe954/EnrichKit). Different gene set databases, including GO, MeSH, Reactome, InterPro, and MsigDB, were interrogated in the enrichment analysis. Terms significantly enriched within target genes were detected using Fisher’s exact test, a test of proportions based on the hypergeometric distribution.

The results of the RNA-seq analysis were validated using real-time polymerase chain reaction (qRT-PCR). We selected three differentially expressed lncRNAs and three differentially expressed protein-coding genes for validation. Primer sequences and expected product lengths are listed in Supplementary Table S2. Real-time PCRs were performed using 7.5 µL of SYBR^® Green master mix (Maxima SYBR Green qPCR Master Mix (2X), with separate ROX vials, Thermo Scientific™, United States of America), equimolar amounts of forward and reverse primers (200 nM, Operon Biotechnologies GmbH, Cologne, Germany), and 20 ng diluted cDNA (1:7.5 in RNase/DNase free water) in a final volume of 15 µL. Samples were analyzed in duplicate in a 72-disk Rotor-GeneTM 6000 (Corbett Life Sciences, Sydney, Australia). Standard amplification conditions were 10 min at 95°C and 40 cycles of 15 s at 95°C, 30 s at 60°C, and 15 s at 72°C. At the end of each run, dissociation curves were analyzed to ensure that the desired amplicon was detected, discarding contaminating DNA or primer dimers. Gene expression was normalized using ACTG1 as a housekeeping gene. Normalized gene expression values (ΔCt) were analyzed using a linear model including the pigmentation group as an independent variable. The association between the normalized gene expression and the pigmentation group was tested using a t-test. The mean and the range of the log2-fold change for each gene were calculated as log₂ (2^−ΔΔCt) using the estimated ΔΔCt value ±standard error.

A total of 111,791 transcripts, including 27,806 unannotated transcripts, were obtained after mapping the sequencing reads to the latest bovine genome reference. Among the novel transcripts, 218 were detected in the ALDB and NONCODE databases and were therefore considered known lncRNAs (Supplementary Data Sheet S1). After training, a coding potential cutoff of 0.36 was selected for the CPAT. This value maximizes specificity and sensitivity (98.3%) (see Supplementary Figure S1), similar to other studies (Billerey et al., 2014; Gupta et al., 2019). All transcripts with a score below 0.36 were retained as potential long noncoding RNAs. After integrating the results from CPC2, CPAT, and PLEK, we identified a total of 4,937 unannotated transcripts as putative lncRNAs (Figure 1, Supplementary Figure S2).

Based on RNA-seq analysis of 18 different tissues, Koufariotis et al. (2015) identified 9,778 lncRNAs in cattle. Of special interest, seven lncRNAs were found to be differentially expressed between white and black skin samples (Koufariotis et al., 2015). A recent paper reported a total of 1,535 expressed lncRNAs encoded by 1,183 putative noncoding genes in bovine oocytes (Wang et al., 2020). More specifically, Weikard et al. (2013) identified 4,848 lncRNAs that are associated with pigmentation in cattle, supporting the idea that lncRNAs play an important role in skin pigmentation. Notably, this number of lncRNAs is similar to the 4,937 transcripts that we reported in the present study, although only 314 were in common. The difference between studies could be explained in terms of the analytical pipeline used to identify putative lncRNAs in each case. Note that at the time of writing this manuscript, there is no widely accepted pipeline to identify lncRNAs.

We identified a total of 4,937 novels and 218 known lncRNAs in our RNA-seq dataset. The lncRNAs were characterized as novel isoforms (51.5%), sense exon overlap (14.5%), antisense intron overlap (12.3%), intergenic (16.5%), and antisense exon overlap (5.2%) forms (Supplementary Figure S3, Supplementary Data Sheet S1). Of the 51.5% characterized as novel isoforms, 94% were generated from regions that harbor protein-coding genes while 6% were novel isoforms from known lncRNAs. Note that similar results were recently reported by Alexandre et al. (2020) working on lncRNAs associated with feed efficiency in cattle.

Novel and known lncRNAs were evenly distributed across the whole genome, ranging from 1.3% (novel lncRNAs) and 0.45% (known lncRNAs) in BTA20 to 6.6% (novel lncRNAs) and 6.8% (known lncRNAs) in BTA3 (Figure 2).

We analyzed the guanine-cytosine content (GC), normalized minimum free energy (MEN), transcript length, exon number, and expression level of all putative novel lncRNAs and compared these metrics with those from known lncRNAs and protein-coding genes. Both groups of lncRNAs displayed significantly lower GC content (Figure 3B), shorter length (Figure 3C), and higher MEN (Figure 3D) than protein-coding genes (p-value ≤ 0.05, Wilcoxon test, Table 1). We found no significant differences in the GC content nor in the MEN between the novel and known lncRNAs (p-value = 0.37 and p-value = 0.55, Wilcoxon test, respectively, Table 1). These results agree with previous studies that showed that lncRNAs have some unique sequence features compared to coding genes (Niazi and Valadkhan, 2012).

Both novel and putative lncRNAs showed significantly higher MEN values than protein-coding genes (p-value ≤ 0.05, Wilcoxon test, Table 1). The sequence itself seems to contribute considerably to MEN since sequences with the same GC content can potentially present different folding energy values based on their secondary structure (Niazi and Valadkhan, 2012). It can be inferred that lncRNAs have a more flexible structure than mRNAs, which may reflect a higher potential to interact with other molecules. Indeed, it has been shown that lncRNAs have GC and MEN values similar to those of 3′ UTR sequences, suggesting that lncRNAs have regulatory functions (Niazi and Valadkhan, 2012).

We found that the novel lncRNAs were significantly shorter than the known lncRNAs (p-value ≤ 0.05, Wilcoxon test, Table 1). Both groups of noncoding transcripts showed lower expression levels (Figure 3A) and a lower number of exons (Figure 3E) than protein-coding genes (p-value ≤ 0.05, Wilcoxon test).

Overall, our results showed that lncRNAs are characterized by a lower number of exons, lower GC content, lower expression level, and higher normalized minimum free energy and tend to be shorter than protein-coding sequences (Derrien et al., 2012; Harrow et al., 2012; Billerey et al., 2014). Note that these sequence feature analyses support the reliability of the putative novel lncRNAs identified. We believe that the marginal differences observed in the length, expression level, and lower number of exons between the novel and known lncRNAs could be explained by the limited catalog of lncRNAs in cattle (Kosinska-Selbi et al., 2020). However, we could not completely discard the possibility that some of our novel lncRNAs are either pseudogenes or misidentified coding sequences.

We found 65 differentially expressed protein-coding genes (p adjusted ≤0.05 and |log₂ FC| ≥1.5) between the pigmented and unpigmented samples (Appendix S2). Among them, MC1R, TYR, PMEL, DCT, MLANA, and KIT showed upregulated expression in pigmented eyelid samples (Jara et al., 2020). These genes are key for the generation, storage, and distribution of melanin (D’Mello et al., 2016). Gene KIT encodes a receptor of tyrosine kinase and is considered a proto-oncogene, involved in the development of melanocytes, including migration, survival, proliferation, and differentiation (Grichnik, 2006). Previous studies showed the important role of KIT in UVB-induced epidermal melanogenesis in humans (Yamada et al., 2013) and in the survival of melanocytes (Mort et al., 2015). Here, KIT showed higher expression in pigmented eyelid samples, suggesting an important role in pigmentation. Interestingly, our analyses suggest that lncRNAs do not interact directly with KIT, and hence, it seems this gene is not directly regulated by these non-coding elements. Moreover, we identified a total of 27 differentially expressed lncRNAs (p adjusted ≤0.05 and |log₂FC| ≥ |1.5|); twenty-four were putative novel lncRNAs, and three were classified as known lncRNAs. Eight lncRNAs showed upregulated and 19 lncRNAs showed downregulated expression in the pigmented eyelid samples. The top lncRNA with significantly downregulated expression was TCONS_00073858 (novel) with log₂FoldChange = −10.3, while the top lncRNA with upregulated expression was TCONS_00088900 (novel) with log₂FoldChange = 7.7.

We analyzed potential target mRNAs of lncRNAs acting in cis, that is, protein-coding genes within 300 kb upstream and downstream of the location of significant differentially expressed lncRNAs. We found four putative target genes, namely, CXCL13, FABP4, GPR143, and UGT1A1, in the flanking regions of four differentially expressed lncRNAs, TCONS_00094682 (novel), TCONS_00021379 (novel), TCONS_00109545 (novel), and TCONS_00078693 (novel), respectively (Table 2). The lncRNAs TCONS_00094682 and TCONS_00021379 showed downregulated expression in the pigmented samples and acted at the cis level with the CXCL13 and FABP4 genes, respectively. The CXCL13 gene encodes a homeostatic chemokine that traffics B cells and is involved in the immune response (Müller et al., 2002). Fatty acid transporters, such as FABP4, were recently shown to increase their expression during metastatic melanoma development, suggesting a higher uptake of fatty acids and metabolic reprogramming that serves as a signature for this condition (Lee et al., 2020). In contrast, the lncRNAs TCONS_00109545 and TCONS_00078693 showed upregulated expression in the pigmented samples. The lncRNA TCONS_00109545 interacts at the cis level with the GPR143 gene, which plays a critical role in retinal health and development (McKay, 2019). Mutations in GPR143 have been associated with the ocular albinism type 1 (OA1) phenotype (Gao et al., 2020), a genetic disorder characterized by reduced ocular pigmentation. Finally, lncRNA TCONS_00078693 was correlated with the expression of UGT1A1, a gene implicated in retinoic acid binding. Retinoid signaling is affected in early carcinogenesis (Tang and Gudas, 2011), and retinoic acid is often used to prevent photoaging of human skin, preventing melanocytic and keratinocytic atypia (Cho et al., 2005). Note that all these lncRNAs showed a positive correlation with their target genes, and hence, we hypothesize that these lncRNAs interact at the cis level, promoting the expression of these genes through the recruitment of proteins that enable transcription loops (Long et al., 2017; Li et al., 2019; Gil and Ulitsky, 2020).

We analyzed the potential target mRNAs of lncRNAs acting in trans from a total of 273 differentially expressed genes that showed sequence complementarity with at least one of the twenty-seven differentially expressed lncRNAs. Of these 273 genes, 193 showed sizable expression correlations with significant lncRNAs (Pearson, r > |0.60|, p-value < 0.05) (Supplementary Data Sheet S3); hence, they were classified as potential trans target genes. We identified target genes involved in melanogenesis, melanosome development (Lin and Fisher, 2007; Kubic et al., 2008; Schmutz and Dreger, 2013), pigmentation processes (Lalueza-Fox et al., 2007; Enriqué Steinberg et al., 2013), tumor pathways (Hoashi et al., 2005), innate immunity and inflammatory signaling (Ito et al., 2015), such as MC1R, PMEL, MLANA, PAX3, IGFBP2, FGF23, and TREM-2. Interestingly, many of these potential trans target genes were reported previously as differentially expressed in pigmented versus non-pigmented samples (Jara et al., 2020) (Supplementary Figure S4, Supplementary Data Sheet S3). The genes MC1R, PMEL, MLANA, and PAX3 showed upregulated expression in the pigmented samples and interacted with lncRNAs with up- and down-regulated expression (Supplementary Data Sheet S2, S3). MC1R is a key pigmentation gene, and its activation in melanocytes stimulates melanogenesis, particularly eumelanogenesis (Beaumont et al., 2007; Cheli et al., 2010; Enriqué Steinberg et al., 2013). The PMEL gene is a key component of mammalian melanosome biogenesis (Clark et al., 2006), and it is required for the generation of cylindrical melanosomes in zebrafish (Watt et al., 2013; Burgoyne et al., 2015). Mutations in PMEL have been shown to regulate hypopigmented phenotypes in vertebrates (Kwon et al., 1995; Kerje et al., 2004). The transcription factor PAX3 interacts with lncRNA TCONS_00088900 and regulates the expression of MITF (Lin and Fisher, 2007), a gene associated with ambilateral circumocular pigmentation in cattle (Pausch et al., 2012). MC1R and MLANA expression levels showed significant positive correlations with the novel lncRNAs TCONS_00088900 and TCONS_00091890, while PMEL also showed a significant positive correlation with TCONS_00091890. MLANA encodes a protein (MART-1) that is localized in melanosomes (De Maziere et al., 2002). Interestingly, MLANA interacts with PMEL and regulates its expression, stability, trafficking, and processing (Hoashi et al., 2005). MLANA showed negative correlations with the lncRNAs ALDBBTAT0000004273 and TCONS_00106295. LncRNAs have different ways of affecting the transcription of target genes in trans, for example, by stabilizing the mRNA (Cao et al., 2017; Long et al., 2017; Li et al., 2019). The participation of lncRNAs in mRNA stabilization and in increasing the expression of specific genes has been recently discussed (Li et al., 2019). From our results, we can infer that ALDBBTAT0000004273 and TCONS_00106295 could be involved in reducing the stability of the mRNA encoded by the MLANA gene and consequently contributing to reduced expression in unpigmented samples. Additionally, the lncRNAs could be involved in the stabilization of the mRNAs of certain genes responsible for pigmentation, such as PEML and MC1R, and thus enhance their expression levels in pigmented eyelids. The genes IGFBP2, FGF23, and TREM-2 showed downregulated expression in the pigmented samples and interacted with lncRNAs with up- and down-regulated expression (Supplementary Data Sheet S2, S3).The FGF23 gene may act as a proinflammatory cytokine (Ito et al., 2015), suggesting a connection between FGF23 and inflammatory processes (Courbebaisse and Lanske, 2018). In fact, high expression of FGF23 is associated with an increased risk of mortality, likely because of its contribution to decreased host defense response to infection, inflammation, and anemia (Courbebaisse and Lanske, 2018). FGF23 showed positive correlations with the ALDBBTAT0000004273 and TCONS_00106295 lncRNAs. The IGFBP2 gene is involved in tumor pathways (Pickard and McCance, 2015) and showed positive correlations with the lncRNAs ALDBBTAT0000001157 and TCONS_00106295. The fatty acid-binding protein 7 (FABP7) gene was also described as a key regulator of cancer metastasis (Cordero et al., 2019), and its expression was upregulated in pigmented samples (Jara et al., 2020). This gene showed significant positive and negative correlations with several lncRNAs. The TREM-2 gene showed a positive correlation with lncRNA TCONS_00073858, which is implicated in innate immunity and inflammatory signaling (Sharif and Knapp, 2008). TREM (TREM-1/TREM-2) gene expression is lower in cutaneous melanoma versus control samples (Nguyen et al., 2015), and these genes could have prognostic and therapeutic value in the treatment of melanoma (Nguyen et al., 2015; Nguyen et al., 2016). From our results, it can be speculated that the lncRNAs described here are involved in stabilizing the mRNAs of genes involved in the immune system and cancer, such as FGF23 (Courbebaisse and Lanske, 2018), IGFBP2 (Pickard and McCance, 2015) and TREM-2.

A total of 193 genes were identified as putative targets of DE lncRNAs in pigmented samples. To identify enriched functions among them, we performed a gene set enrichment analysis using the annotation retrieved from different databases, including GO, KEGG, MeSH, InterPro, Reactome, and MSigDB. Figure 4 shows the most relevant biological terms and pathways associated with eyelid pigmentation. The most significant functional terms were related to skin pigmentation, such as melanosome (GO:0042470), melanin biosynthetic process (GO:0042438), melanosome organization (GO:0032438), melanocyte differentiation (GO:0030318), melanogenesis (bta04916), skin pigmentation (D012880), pigmentation (D010858) and melanogenesis (M7761) (Figure 4, Supplementary Data Sheet S4). Notably, we found several significant terms related to the inflammatory response and infectious and tumoral pathways: immune response (GO:0006955), defense response to bacterium (GO:0042742), leukocyte chemotaxis (GO:0030595), cytokine-cytokine receptor interaction (bta04060), T cell activation (GO:0042110), cytokine-cytokine receptor interaction (bta04060), immune response (M12401), and immune (humoral) and inflammatory response (M8838) (Figure 4, Supplementary Data Sheet S4). Our findings indicate that target genes of lncRNAs with up- and down-regulated expression in eyelid skin are associated not only with pigmentation or melanogenesis but also with the immune response.

We validated the findings of the RNA-seq experiment using qRT-PCR. Three lncRNAs and three protein-coding genes were evaluated. Figure 5 shows the log2-fold differences in gene expression measured by both RNA-Seq and qRT-PCR, confirming that the six genes showed very similar patterns of abundance with both methods.