Development of computational pipelines for transcriptome and miRNome characterization from RNA-seq data applied to swine adipose tissue

High throughput technologies for DNA sequencing are used more and more frequently for gene expression profiling studies (RNA-seq). With respect to other techniques such as microarrays, RNA-seq has higher sensitivity in retrieving the expressed molecules and presents the advantageous feature of allowing the detection of unknown or uncharacterized transcripts. RNA-seq data processing involves several computational steps (input preprocessing for quality evaluation and cleaning; read alignment to reference genome; transcript identification, quantification, and annotation; differential expression assessment) that have to be performed in sequential order, thus resulting in a computational pipeline. Each single RNA-seq experiment can produce large amounts of data that require the use of efficient computational methods to obtain transcriptome qualitative and quantitative characterization. There are different methods that implement each conceptual pipeline step, and new ones are continuously proposed. However, because of the variety of biological questions and study designs to which RNA-seq experiments can be applied to, there is not a commonly adopted implementation of the processing workflow. In this thesis, we developed a computational pipeline for the analysis of RNA-seq data focused on the linear transcriptome, extended an existing pipeline that analyzes RNA-seq data of microRNAs (miRNAs) and miRNA-like small RNAs, and started to develop a computational pipeline for the detection and quantification of circular RNAs. The main objectives of the first two pipelines were the profiling of the set of the transcripts (transcriptome) and small RNAs (miRNome) expressed in the considered samples, by the identification of known and new RNAs. They allowed as well to investigate RNA sequence variations (such as miRNA isomiRs), transcripts and small RNAs expression levels, and to compare expression profiles between different sample groups. The pig (Sus scrofa) is a model organism for human diseases, and very important per se for the meat industry. Fat and backfat tissues are subject of very active research since fat attributes and deposition traits are in strong connection with technological aspects and quality of pig products. However, the global framework of the biological and molecular processes regulating backfat deposition in pig is still incomplete. We applied our pipelines to RNA-seq data of polyadenylated and of small RNAs from pig subcutaneous adipose tissue samples from 20 Italian Large White (ILW) individuals. Selected animals were reared under very standard conditions but presented, for fat traits, extreme and divergent phenotypes (FAT and LEAN pigs) and genetic merits. The backfat transcription profile was characterized by the expression of 23,483 genes, of which only 54.1% were represented by known genes. Of 63,418 expressed transcripts, about 80% were non-previously annotated isoforms. By comparing the expression level of FAT vs. LEAN pigs, we detected 86 robust differentially expressed transcripts, 72 more expressed in fat pigs (including ACP5, BCL2A1, CCR1, CD163, CD1A, EGR2, ENPP1, GPNMB, INHBB, LYZ, MSR1, OLR1, PIK3AP1, PLIN2, SPP1, SLC11A1, STC1) and 14 less expressed (including ADSSL1, CDO1, DNAJB1, HSPA1A, HSPA1B, HSPA2, HSPB8, IGFBP5, OLFML3). Overexpressed genes were implied particularly in immune system processes, response to stimulus, cell activation and skeletal system development. Underexpressed genes included five heat shock proteins and were involved in unfolded protein binding and stress response functional categories. Adipose tissue alterations and impaired stress response are linked to inflammation and, in turn, to adipose tissue secretory activity, similar to what is observed in human obesity. MiRNAs play important roles in cell differentiation and physiology acting as post-transcriptional regulators of gene expression by silencing targeted transcripts. The pig backfat miRNome showed the expression of hundreds of small RNAs, including putative new miRNAs, new miRNA isoforms (isomiRs), and new moRNAs, likely produced from the terminal regions of non-canonically processed hairpin precursors. From a first study on two samples, we detected 222 known miRNAs, 68 new miRNAs and 17 moRNAs expressed from known hairpins, and 312 new miRNAs expressed from 253 new hairpins. The expression of five small RNAs, including moRNA ssc-moR-21-5p and a miRNA from a new hairpin, was validated by a qRT-PCR assay, thus confirming the robustness of our results. A second study on 18 samples identified a largely overlapping miRNome in terms of expressed elements and variations, and was important to identify differentially expressed miRNAs and moRNAs in FAT and LEAN subjects. We predicted putative regulatory interactions between small RNAs and transcripts by sequence analysis, using custom target predictions on reconstructed transcript sequences and miRNA isomiRs. We then integrated target prediction results with combined analysis of miRNA and transcript expression data, to eventually select miRNA-transcript relations most supported by negative correlation of expression profiles. Further, the predicted network of miRNA-transcript interactions was enriched by information on transcript differential expression, functional annotations and coding potential predictions, and transcript overlap with pig QTL genomic regions. In this way we were able to focus on a restricted and possibly most significant number of interactions that need to be experimentally investigated. Additional considerations are coming from the study of the possible impact of specific differentially expressed miRNAs to genes belonging to the pathways most germane to adipose tissue features. The applicative results of these studies enlarged the knowledge of transcripts and small RNAs expressed in the pig adipose tissue, as well as small RNA-transcripts regulatory interactions, providing information helpful for a better understanding of ILW pig backfat and future studies on gene expression regulation in this tissue. Moreover, the methods presented here are currently undergoing further development and extension, and have applications well over and above those presented in this thesis.

Le tecnologie per il sequenziamento massivo del DNA sono spesso usate per studiare il trascrittoma e ottenre profili d’espressione genica su scala genomica (RNA-seq). Rispetto ad altre tecnologie come i microarray, l’RNA-seq ha una maggiore sensibilità nel campionare e quantificare le molecole espresse e permette inoltre l’identificazione di trascritti sconosciuti o non caratterizzati. Il processamento di dati RNA-seq prevede molteplici passaggi di analisi (preprocessamento degli input per la valutazione della qualità e pulizia, allineamento delle read al genoma di riferimento, identificazione, quantificazione e annotazione dei trascritti, stima di espressione differenziale) che devono essere eseguiti in ordine sequenziale, mediante pipeline computazionali. Ogni singolo esperimento di RNA-seq può produrre grandi quantità di dati che richiedono l’impiego di metodi efficienti per ottenere la caratterizzazione qualitativa e quantitativa del trascrittoma. Esistono diversi metodi che implementano ogni passaggio concettuale di analisi e nuovi ne vengono continuamente proposti. Questo e’ anche dovuto alla varietà dei quesiti biologici e disegni sperimentali a cui gli esperimenti di RNA-seq possono essere applicati. Di converso, non esiste un’implementazione comunemente adottata dello schema di processamento. In questa tesi, abbiamo sviluppato una pipeline computazionale per l’analisi di dati RNA-seq focalizzata sul trascrittoma lineare; abbiamo esteso una pipeline esistente che analizza dati di RNA-seq di microRNA (miRNA) e piccoli RNA simili ai miRNA ed abbiamo iniziato a sviluppare una pipeline computazionale per l’identificazione e la quantificazione di RNA circolari. Gli obiettivi principali delle prime due pipeline sono il profiling dell’insieme dei trascritti (trascrittoma) e piccoli RNA (miRNoma) espressi, con l’identificazione di RNA noti e nuovi. Inoltre, è stato possibile studiare le variazioni di sequenza degli RNA (come gli isomiR dei miRNA), dei livelli di espressione di trascritti e piccoli RNA, e confrontare i profili di espressione tra diversi gruppi di campioni biologici. Il maiale (Sus scrofa) è un organismo modello per numerose malattie o condizioni umane, ma anche molto importante di per sé per l’industria di carne e derivati di alto pregio economicamente importanti. Il tessuto adiposo e il lardo dorsale sono oggetto di attiva ricerca, poichè alcune caratteristiche qualitative e quantitative del grasso e i meccanismi e tassi di deposito e accumulazione del grasso sono in stretta connessione con aspetti tecnologici e risultati qualitativi dei prodotti finali, come il prosciutto crudo. Tuttavia, il quadro complessivo dei processi biologici e molecolari che regolano il deposito del lardo dorsale nei maiali è ancora incompleto. In questa tesi, abbiamo applicato i metodi di analisi sviluppati a dati RNA-seq di RNA poliadenilati e piccoli RNA da campioni di tessuto adiposo sottocutaneo di 20 soggetti di razza Italian Large White (ILW). Gli animali selezionati sono stati allevati in condizioni molto standardizzate, ma presentano, riguardo i tratti del grasso, fenotipi e corrispondenti meriti genetici estremi e divergenti (maiali FAT e LEAN). L’analisi del profilo trascrizionale del lardo dorsale ha identificato l’espressione di 23.483 geni, dei quali solo il 54,1% rappresentato da geni noti. Dei 63.418 trascritti espressi, circa l’80% erano isoforme non precedentemente annotate. Confrontando i livelli di espressione dei maiali FAT contro i maiali LEAN, abbiamo poi identificato, con criteri molto stringenti, 86 trascritti differenzialmente espressi: 72 espressi a livelli più alti nei maiali obesi (tra cui ACP5, BCL2A1, CCR1, CD163, CD1A, EGR2, ENPP1, GPNMB, INHBB, LYZ, MSR1, OLR1, PIK3AP1, PLIN2, SPP1, SLC11A1, STC1) e 14 meno espressi (inclusi ADSSL1, CDO1, DNAJB1, HSPA1A, HSPA1B, HSPA2, HSPB8, IGFBP5, OLFML3). I geni sovraespressi sono implicati in processi del sistema immunitario, di risposta allo stimolo, attivazione cellulare e sviluppo dell’apparato scheletrico. I geni sottoespressi includono cinque proteine heat shock e sono associati a categorie funzionali quali il legame di proteine mal ripiegate, e la risposta allo stress. Nel tessuto adiposo un’eccessiva adiposità combinata a carenze nei meccanismi di risposta allo stress sono collegate ad uno stato infiammatorio del tessuto e, di conseguenza, ad alterazioni dell’attività secretoria del tessuto adiposo, similmente a quanto è stato osservato nell’obesità umana. I miRNA sono importanti regolatori dell’espressione genica nel differenziamento, nell sviluppo e nella fisiologia cellulare dei diversi tessuti. Essi agiscono come regolatori post-trascrizionali dell’espressione genica, silenziando i trascritti bersaglio. Lo studio del miRNoma del lardo dorsale di maiale ha identificato l’espressione di centinaia di piccoli RNA, includendo potenziali nuovi miRNA, nuove isoforme di miRNA (isomiR) e nuovi microRNA-offeset RNA (moRNA), probabilmente prodotti dalle regioni terminali di precursori a forcina processate in modo non canonico. Da uno studio preliminare condotto su due campioni abbiamo rilevato 222 miRNA noti, 68 nuovi miRNA e 17 moRNA espressi da forcine note, e 312 nuovi miRNA espressi da 253 nuove forcine. L’espressione di cinque piccoli RNA, inclusi il moRNA ssc-moR-21-5p e un miRNA prodotto da un precursore da noi predetto, è stata validata mediante qRT-PCR, confermando l’affidabilità dei nostri risultati. In accodo con questi dati, un secondo studio condotto su 18 campioni ha identificato un miRNoma molto simile in termini di elementi espressi e varianti. Questo ha inoltre permesso di identificare miRNA e moRNA differenzialmente espressi tra soggetti FAT e LEAN, potenziali regolatori di trascritti la cui modulazione dell’espressione potrebbe essere implicata nelle variazioni fenotipiche dei soggetti considerati. Abbiamo predetto i potenziali bersagli dei miRNA e dei moRNA (nell ipotesi che i moRNA possano funzionare come miRNA) modulati prendendo in considerazione, per analisi ad hoc le sequenze dei trascritti ricostruite in precedenza e gli isomiR dei miRNA risultati maggiormente espressi e quindi rilevanti. Abbiamo integrato i risultati di queste predizioni con l’analisi combinata dei profili d’espressione di miRNA e trascritti, per selezionare le relazioni miRNA-trascritto maggiormente supportate dai dati d’espressione. La rete di interazioni miRNA-trascritti ottenuta in questo modo è stata arricchita dall’informazione su espressione differenziale, annotazione funzionale e predizioni del potenziale codificante e sovrapposizione dei trascritti con regioni genomiche di QTL di maiale. In questo modo siamo stati in grado di identificare un numero ristretto di interazioni potenzialmente molto significative che necessitano di essere investigate sperimentalmente. Ulteriori considerazioni stanno emergendo dallo studio del potenziale impatto di specifici miRNA differenzialmente espressi su geni appartenenti a pathway molto attinenti alla biologia del tessuto adiposo. I risultati applicativi di questi studi hanno allargato la conoscenza dei trascritti e dei piccoli RNA espressi nel tessuto adiposo di maiale, e anche delle interazioni regolative tra piccoli RNA e trascritti, fornendo utili informazioni per una miglior comprensione del lardo dorsale di maiali ILW e nuove ipotesi per studi futuri sulla regolazione dell’espressione genica in questo tessuto. In aggiunta, stiamo attualmente sviluppando ed estendendo ulteriormente i metodi qui presentati, con applicazioni e obiettivi ulteriori rispetto a quelli descritti in questa tesi.