biogenesis of most eukaryotic mrnas involves splicing and cleavage and polyadenylation ( 3 end processing ) ( derti et  al . , 2012 , tian et  al . , 2005 ) .
both mechanisms are required to produce functional mrnas and are also important to regulate gene expression by producing alternative mrna isoforms and for efficient transcription termination ( di giammartino et  al . , 2011 ,
the alternative isoforms are produced by alternative splicing or by use of alternative polyadenylation ( apa ) sites .
both alternative splicing and apa are regulated by rna - binding proteins ( rbps ) ( di giammartino et  al . , 2011 ,
, 2013 , fu and ares , 2014 , shi , 2012 , witten and ule , 2011 ) .
however , few tools are available to study both processes in an integrated manner , and the overlap between their regulatory programs is poorly understood . the regulatory function of many rbps depends on their binding position in respect to the regulated exon ( witten and ule , 2011 ) .
such position - dependent regulatory principles have been visualized at high - resolution in the form of rna maps ( ule et  al . ,
2006 ) , and were exploited to derive codes that can predict tissue - specific splicing patterns ( alipanahi et  al .
it is clear that the rbp binding on nascent rna can also affect apa in a position - dependent manner ( batra et  al . , 2014 ,
three studies have used crosslinking and immunoprecipitation ( clip ) to define the rna map of apa ( licatalosi et  al . , 2008 , masuda et  al . , 2015 ,
batra et  al . , 2014 ) , but these studies have plotted the position of full clip reads around the regulated poly(a ) sites , rather than the position of crosslink sites , and have not evaluated the statistical significance of identified enrichments . therefore , the importance in binding position for guiding the repressing or enhancing effects of rbps remains unclear .
the nucleotide resolution that is obtained by the truncated cdnas in iclip , and the quantitative nature gained by the analysis of unique molecular identifiers ( umis ) , allowed us to examine the rna maps of apa regulation in much greater detail .
our study examines how the tar ( transactive response ) dna - binding protein 43 ( tdp-43 , also referred to as tardbp ) regulates splicing and apa via position - dependent principles .
tdp-43 is an rbp involved in several neurodegenerative diseases , including frontotemporal lobar degeneration ( ftld - tdp ) and amyotrophic lateral sclerosis ( als ) ( ratti and buratti , 2016 ) .
tdp-43 regulates alternative splicing in a position - dependent manner by binding to intronic dinucleotide of uridine and guanosine ( ug)-rich motifs , such that binding close to the splice sites represses splicing , whereas binding further downstream of the exon enhances splicing ( tollervey et  al . , 2011 ) .
the binding of tdp-43 is also enriched in 3 utrs , and it can regulate apa of its own transcript ( ratti and buratti , 2016 , tollervey et  al . , 2011 ) , indicating a possible role for regulating apa of other transcripts .
we identified poly(a ) sites regulated by tdp-43 in hek293 cells  by 3 mrna sequencing and defined tdp-43 binding sites  with individual nucleotide crosslinking and immunoprecipitation ( iclip ) . to perform 3 end data analysis
, we developed  expressrna , a modular bioinformatics research platform that manages data and inter - connects recently developed ( expressrna.apa module , rnamotifs2 ) and existing ( star , dexseq ) analysis software .
the expressrna.apa module identifies the sites where cleavage and polyadenylation takes place ( poly(a ) sites ) , marks and annotates the differentially regulated poly(a ) sites , and visualizes binding patterns of rbps around these sites .
the resulting rna maps defined the binding positions where tdp-43 can enhance or repress multiple types of poly(a ) sites , and we validated its direct action with a minigene reporter .
we demonstrated that the role of tdp-43 can be replaced in this minigene by tia1 , another known splicing regulator , if ug - rich motifs are replaced by ua - rich motifs . to identify enriched clusters of diverse regulatory motifs around the regulated exons or poly(a ) sites
, we extended and upgraded rna motifs ( cereda et  al . , 2014 )
this identified enriched clusters composed of multiple types of motifs around tissue - specific exons or poly(a ) sites , which can indicate a potential for combinatorial regulation of splicing and apa .
to study how tdp-43 regulates apa , rna  was isolated from six independent tdp-43 knockdown ( kd ) and control hek293 cells .
the 3 ends of mrnas were amplified with the quantseq rev 3 mrna sequencing ( mrna - seq ) library prep kit ( lexogen ) , which uses a poly(t ) primer to reverse transcribe the mrnas .
the library was sequenced with hiseq , producing 60-nt single - end reads and 10-nt index reads . for data analysis , we developed the expressrna web application and analysis platform ( figures 1 and s1b ) .
the first step of expressrna is to process and map the 3 mrna - seq data to the genome , classify the sites where cleavage and polyadenylation takes place ( poly(a ) sites ) , and identify the differentially regulated poly(a ) sites .
we first aligned the reads to the hg19 reference genome with star ( dobin et  al . , 2013 ) , allowing soft clipping from both 5 and 3 end .
of the total 52.96 million reads , 54.83% aligned uniquely to the human genome ( table s1 ) . to construct the database of identified poly(a ) sites , we considered sequence data from all experiments as one dataset .
some of the aligned reads were soft - clipped due to the sequencing running into the poly(a ) tail or imperfect primer annealing ( figure  s2a).figure  1expressrna research platformthe expressrna platform performs analyses of alternative polyadenylation datasets and can include external alternative splicing datasets .
the identified or provided regulated features ( poly(a ) sites , alternative exons ) are then combined with rna protein binding information ( iclip ) .
expressrna research platform the expressrna platform performs analyses of alternative polyadenylation datasets and can include external alternative splicing datasets . the identified or provided regulated features ( poly(a ) sites , alternative exons )
since the longest 3 utr isoforms are in some cases not fully annotated , we added 5 kb of the intergenic region downstream of each gene . if two genes are closer than 10 kb , only the region up to the middle was added .
we found that alternative positions of cleavage and polyadenylation can be clustered in close proximity of a dominant poly(a ) site , with most variation occurring within 5 nt of the dominant site , indicating that the cleavage position is not always precise and can vary by a few nucleotides ( figure  s2b ) .
such cleavage positions occur downstream of a single polyadenylation signal ( pas ) , and therefore their variation likely reflects lack of cleavage precision by the cleavage and polyadenylation machinery . for quantification of poly(a ) sites
, we therefore summed up the counts of reads that identified cleavage up to 5  nt away from each dominant poly(a ) site ( figure  s2c ) .
published studies also indicate that auxiliary rna motifs tend to be enriched in the region approximately up to 75 nt upstream ( ure , upstream regulatory elements ) and 50 nt downstream ( dre , downstream regulatory elements ) of each cleavage site ( beaudoing et  al . ,
we wished to focus our study on fully independent cleavage sites that contain their own pas and auxiliary motifs .
therefore , we identified the dominant poly(a ) sites based on read count , such that all resulting sites were at least 125 nt apart .
quantseq relies on annealing a poly(t ) primer to the poly(a ) tail of mrnas to identify the 3 end of the mrnas .
however , as has been shown previously , poly(t ) primers often also anneal to internal a - rich sites in mrnas ( derti et  al . ,
therefore , poly(a ) sites with a - tracts in the vicinity [ 10 .. 10 ] were first filtered out .
the poly(a ) sites were further classified into three classes by using the software poly(a)r ( akhtar et  al . , 2010 ) , which evaluates the presence of preceding pas .
this identified 17,102 pas - strong , 2,599 pas - weak , and 10,513 pas - less sites ( table s4 ) .
we examined the nucleotide composition , the efficiency of cleavage , and the level of overlap with two past studies ( derti et  al .
, 2012 , gruber et  al . , 2012 ) , which confirms that the pas - strong and pas - weak sites are the most reliable and efficiently used sites ( figures s2d and s2e  s2h ) . classification and broad identification of poly(a ) sites in hek293 cells has been achieved earlier by a comprehensive analysis of 3 end targeted sequencing datasets ( gruber et  al . , 2016 ) .
in contrast , the aim of our study is not to further characterize newly identified sites , but rather to use a stringent filtering approach to examine the regulatory patterns at the functional poly(a ) sites .
we therefore focused only on the pas - strong and pas - weak sites , which identified a total of 19,701 poly(a ) sites , 16,599 sites were annotated to genes ( ensembl v.74 ) . to ensure robust analysis , poly(a ) sites with less than ten read counts in either control or kd experiments
we then identified the poly(a ) site in each gene that contained the highest read count , which we consider as the  major poly(a ) site .  to avoid poorly used sites , we then discarded those poly(a ) sites that have less than 5% of reads compared to the number of reads present at the major poly(a ) site in the same gene .
15,113 sites in 10,069 genes remained after these filtering steps ( figure  s2d ) .
comparisons of our 15,113 sites to the published poly(a ) database in human tissues ( derti et  al . , 2012 ) and hek293 cells ( gruber et  al . , 2012
) demonstrates that the majority of sites overlap at nucleotide precision ( figure  2d ) . allowing for up to a 10-nt spacing
, approximately 70% of our sites overlap with the sites in derti et  al .
( 2012 ) and 65% with the previously defined sites in  hek293 cells from gruber et  al .
the overlapping sites have the characteristic nucleotide signature of a - rich pas upstream of the cleavage site , which is followed by t - rich sequences , in agreement with past studies ( sheets et  al . , 1990 ) ( figure  s2f ) . a larger number of poly(a ) sites were identified in hek293 cells by a previous study due to the use of milder filtering criteria ( gruber et  al . , 2012 ) , but the sites that are not shared with our study have weaker nucleotide signature compared to the overlapping sites ( figure  s2 g )
. however , the sites identified in our study have similar nucleotide signature to the ones that overlap with the published sites , indicating that the stringent filtering criteria applied by our study are sufficient to ensure that most of the identified poly(a ) sites are valid ( figure  s2 g ) .
moreover , de novo search of the most enriched hexamers in the region where the pas is normally located ( [ 30 .. 18 ] upstream of the poly(a ) site ) found the expected consensus sequences ( beaudoing et  al . , 2000 )
( figure  s2h ; supplemental experimental procedures).figure  2mapping reads and evaluating poly(a ) site loci and expression(a ) read alignment to the hg19 reference genome .
reads are soft clipped ( red ) due to poly(a ) tail sequence at the 5 end or imperfect primer annealing at the 3  end.(b ) classification of proximal and distal poly(a ) site pairs ( 5 s1  = 5 splice site 1 ) .
the same - exon pairs are not limited to the last exons in the gene , 9% ( 201 of 2,202 ) of the same - exon pairs are annotated to a non - terminal exon.(c ) dexseq computes direction ( log2 fold change ) and significance ( adjusted p value ) of poly(a ) sites in control versus tdp-43 kd . the proximal - distal site pair with highest fold change is selected among significantly changed poly(a ) sites ( adjusted p value < 0.05 ) .
control distributions are drawn as black lines in all rna maps.(d ) overlap with derti et  al .
, 2012 , gruber et  al . , 2012 , and nam et  al .
( 2014 ) poly(a ) sites around the poly(a ) sites defined by the present study .
mapping reads and evaluating poly(a ) site loci and expression ( a ) read alignment to the hg19 reference genome .
reads are soft clipped ( red ) due to poly(a ) tail sequence at the 5 end or imperfect primer annealing at the 3  end .
( b ) classification of proximal and distal poly(a ) site pairs ( 5 s1  = 5 splice site 1 ) .
the same - exon pairs are not limited to the last exons in the gene , 9% ( 201 of 2,202 ) of the same - exon pairs are annotated to a non - terminal exon .
( c ) dexseq computes direction ( log2 fold change ) and significance ( adjusted p value ) of poly(a ) sites in control versus tdp-43 kd . the proximal - distal site pair with highest fold change is selected among significantly changed poly(a ) sites ( adjusted p value < 0.05 ) .
, 2012 , gruber et  al . , 2012 , and nam et  al .
( 2014 ) poly(a ) sites around the poly(a ) sites defined by the present study .
we further tested the validity of our detected poly(a ) sites by  plotting the coverage of rna sequencing ( rna - seq ) experiments from published studies performed in hek293 cells ( pham et  al .
, 2016 , trakman et  al . , 2016 ) ( figures s2j  s2l ) .
we computed the read coverage in the region 500 .. 200 around proximal and distal poly(a ) sites , defined by the present study , which were separated into the 85% that are shared with the previous study by nam et  al .
we find a peak of increased rna - seq average coverage upstream of both types of poly(a ) sites , and a drop to negligible levels of rna - seq coverage downstream of the distal poly(a ) sites ( figure  s2j ) .
this significant difference in read coverage just upstream and downstream of both the annotated and poly(a ) sites that we have annotated additionally supports the validity of poly(a ) sites uncovered here . for further comparison
, we also plotted heatmaps of rna - seq read coverage around the 50 known and poly(a ) sites uncovered by this analysis that contain most quantseq reads by displaying the quantseq 3 end targeted sequencing coverage ( data from our study ) alongside the rna - seq coverage data ( figures s2k and s2l ) .
this shows that the rna - seq read distribution at individual poly(a ) sites is similar for known and sites uncovered by this analysis , since most sites contain enrichment of rna - seq reads upstream , and lack the reads downstream of the poly(a ) sites .
taken together , this indicates that the poly(a ) sites are accurately identified by the expressrna analysis of quantseq data in hek293 cells . to define the regulatory principles with high fidelity
, we focused our analyses on 3,291 genes where we can robustly annotate multiple poly(a ) sites .
this represents about 33% of detected genes , which is slightly lower compared to most past studies , which reported approximately 40% of protein coding genes with > 1 poly(a ) site ( ni et  al . , 2013 ) .
it is much lower than the report of 70% with multiple poly(a ) sites that were detected with the cross - tissue analysis ( derti et  al . , 2012 ) , but that is expected given that a single cell type has less rna diversity compared to all the tissues .
the lower proportion of genes with multiple poly(a ) sites likely reflects also our stringent requirement for that each poly(a ) site contains at least 5% of the reads that map to each gene .
next , we wished to identify those genes with multiple poly(a ) sites where the use of the poly(a ) sites changes upon tdp-43 kd .
statistically significant changes in poly(a ) site use between control and kd conditions were identified using dexseq ( anders et  al . , 2012 )
if more than two poly(a ) sites were identified in a gene , then the two poly(a ) sites most significantly changed ( adjusted p value < 0.05 ) were considered for further analyses of regulated sites . if only one site had an adjusted p  value < 0.05 , then the second site was selected based on highest read count . if no site had p value < 0.05 , then both sites were selected based on highest read count .
our approach identified 3,291 genes ( poly(a ) site pairs ) that we then classified based on the position of the two poly(a ) sites in the gene ( figure  2b ) . if both poly(a ) sites are in the same exon , we classified them as same exon ( the most common apa type , also referred as tandem in previous studies ) ( elkon et  al . , 2013 ) . if the proximal poly(a ) site is part of a composite exon that contains an internal 5 splice site , and the two poly(a ) sites were generated by alternative 5 splice site use , we classified them as a composite exon . finally , if the proximal poly(a ) site is in an exon that is fully skipped when the distal poly(a ) site is used , we classified them as skipped exon ( figure  2b ) .
the final filtered poly(a ) data included 3,291 poly(a ) site pairs , of which 2,202 belonged to the same - exon , 662 to the skipped - exon , and 427 to the composite - exon class . to identify changes in apa between control and tdp-43 kd
, we examined the changes in relative read counts at the proximal and distal poly(a ) sites in each gene , which is referred to as fold  change ( reported by dexseq ) .
the fold change is used to determine the direction of change as repressed ( fold change < 0 ) or enhanced ( fold change > 0 ) . for genes with no significantly regulated poly(a ) sites , two sites with highest read counts across both control and kd conditions were selected and classified into control enhanced and control repressed ( figure  2c ) . to identify tdp-43 binding sites in hek293 cells , we performed iclip of tdp-43 in two replicate experiments , which together generated 12,622,661 uniquely mapped iclip cdnas .
we report the number of iclip cdnas around each regulated poly(a ) site ( table s2 ) .
we then identified 415,238 significant crosslink clusters , 9.4% of which map to annotated 3 utrs ( figure  s2i ) .
we visualized the positions of these crosslink clusters around the same - exon class of regulated poly(a ) sites , including both the proximal ( figures 3a3c ) and the distal poly(a ) sites ( figure  3d ) .
we also examined the ug - rich binding motifs ( see defining tdp-43 binding positions with iclip ) , which are known to bind tdp-43 . reassuringly , crosslink clusters and enriched ug - rich motifs are abundant at similar positions ( figures 3a and s3a  s3d ) , further indicating that these two independent approaches correctly define the binding sites of tdp-43.figure  3rna maps of tdp-43 around proximal and distal same - exon poly(a ) sites with iclip clusters(a ) proximal poly(a ) site rna map of tdp-43 iclip at regulated genes .
tdp-43 is bound around repressed poly(a ) sites and to a lesser extent ( sparsely ) further upstream and downstream of enhanced poly(a ) sites.(b ) top 20 iclip mrna targets contributing to proximal poly(a ) site repression.(c ) top 20 iclip mrna targets contributing to proximal poly(a ) site enhancement.(d ) distal site rna map of tdp-43 binding around regulated genes .
the level of binding is lower compared to the proximal sites ( table s6 ) .
rna maps of tdp-43 around proximal and distal same - exon poly(a ) sites with iclip clusters ( a ) proximal poly(a ) site rna map of tdp-43 iclip at regulated genes .
tdp-43 is bound around repressed poly(a ) sites and to a lesser extent ( sparsely ) further upstream and downstream of enhanced poly(a ) sites .
( b ) top 20 iclip mrna targets contributing to proximal poly(a ) site repression . ( c )
the level of binding is lower compared to the proximal sites ( table s6 ) . to further assess the specificity of the rna map
, we redrew it by using just the 13,034 poly(a ) sites that overlap with the ones identified also with the 3p - seq method in the previous study ( figure  2d ) .
the resulting same - exon proximal rna map is almost identical to the map that used all the sites ( figures 3a and s3e ) , which demonstrates that the map is robust and largely unaffected by the identification of additional poly(a ) sites .
moreover , we compared iclip binding of tdp-43 , tial1 , and stau1 around the same - exon proximal poly(a ) sites regulated by tdp-43 ( figure  s3f ) .
we chose tia1 and stau1 as controls , since both of these rbps also have enriched crosslinking to 3  utrs .
we plotted the enrichment of crosslink clusters for each protein by comparing regulated versus control poly(a ) sites , which demonstrated much stronger and position - specific enrichment for tdp-43 compared to the other control rbps .
since tdp-43 is a known regulator of splicing , we also examined whether use of the poly(a ) sites might be indirectly affected due to co - regulation ( binding ) at nearby splice sites .
this is possible in the case of composite - exon and skipped - exon poly(a ) sites , which contain at least one splice between the regulated poly(a ) sites ( s1 and s2 , as marked in figure  2b ) .
we found position - dependent tdp-43 crosslinking at the regulated composite - exon and skipped - exon poly(a ) sites ( figures s4a  s4d ) and much less at splice sites flanking these poly(a ) sites ( figures s5a  s5f ; table s6 ) .
this indicates that tdp-43 rarely regulates apa indirectly via splicing but rather directly regulates both types of competing poly(a ) sites . to demonstrate in a simple way that tdp-43 binds at different positions to repress or enhance the poly(a ) sites
, we defined the 40-nt window around each type of regulated poly(a ) sites that had the strongest enrichment of crosslink clusters compared to controls ( repressed / enhanced versus control genes ) ( table s6 ) . to determine the number of genes that contain a specific class of alternative poly(a ) site that is directly regulated by tdp-43
, we counted the number of iclip crosslink clusters at this 40-nt window , which allowed us to calculate the bound regulated genes score ( brg ) ( see experimental procedures and table s6 ) .
this shows that tdp-43 most often binds next to the proximal same - exon poly(a ) sites ( 97 .. 57 relative to proximal poly(a ) site , brg  = 33 , p value 1e-6 , figure  3 ) , while enriched binding is also seen further downstream of the enhanced sites ( 72 .. 112 relative to proximal poly(a ) site , brg  = 19 , p value 3e-5 , figure  3 ) .
this pattern is reminiscent of the rna maps of splicing regulation , where tdp-43 binds directly upstream or within the exon to repress and further downstream of the exon to enhance splicing ( tollervey et  al . , 2011 ) . to understand whether the rna map of apa regulation by tdp-43 can be discovered without any knowledge of its binding specificity
so that it could assess regulated poly(a ) sites , in addition to alternative exons .
moreover , while the original rnamotifs could only identify clusters of highly similar motifs around alternative exons , the version ( named rnamotifs2 ) can identify clusters of diverse motifs enriched around groups of exons or poly(a ) sites .
this can assess the potential for multiple rbps to combinatorially regulate pre - mrna processing .
we applied rnamotifs2 to  analysis of the same - exon class of poly(a ) sites regulated by  tdp-43 .
the enrichment of the detected significant motif clusters ( fisher s exact test , p  < 0.01 ) is plotted in blue at the  repressed and in red at enhanced poly(a ) sites ( figure  4 ) .
ug - rich motif clusters were enriched mainly around the proximal regulated sites at similar positions as the iclip crosslink sites ( figures 3a and s3a ; table s6 ) .
interestingly , the distal sites mainly had enrichment of u - rich and ya - rich motifs ( figure  4b ) , indicating a potential for regulation of competing poly(a ) sites by different rbps.figure  4motif analysis around tdp-43-regulated same - exon poly(a ) sitesup to two of the most significant motif clusters are shown for each search region ( r1 , r2 , r3).(a ) motif analysis around proximal poly(a ) sites show significant ug - rich clusters in r1 [ 100 .. 40 ] and r2 [ 40 .. 20 ] around repressed proximal poly(a ) sites .
more distal binding in r3 [ 20 .. 80 ] results in enhancement.(b ) motif analysis around distal poly(a ) sites reveals less pronounced regulatory effects , mostly enhancement guided by uc and ua - rich clusters .
motif analysis around tdp-43-regulated same - exon poly(a ) sites up to two of the most significant motif clusters are shown for each search region ( r1 , r2 , r3 ) .
( a ) motif analysis around proximal poly(a ) sites show significant ug - rich clusters in r1 [ 100 .. 40 ] and r2 [ 40 .. 20 ] around repressed proximal poly(a ) sites .
( b ) motif analysis around distal poly(a ) sites reveals less pronounced regulatory effects , mostly enhancement guided by uc and ua - rich clusters . to study whether the binding of tdp-43 alone is sufficient to  regulate apa
, we produced a minigene reporter with the 3  utr of structural maintenance of chromosomes 1a ( smc1a ) gene , which contains regulated poly(a ) sites from the same - exon class ( figure  5a ) .
the 3 utr of smc1a gene was cloned downstream of the firefly luciferase orf in a modified pcdna3 plasmid , which does not contain any additional poly(a ) sites .
qrt - pcr was then used to quantify the use of distal poly(a ) sites , which was normalized by the total amount of mrna that was produced . to distinguish the minigene expression from the endogenous smc1a gene
, we quantified the total mrna with a forward primer that recognizes part of the luciferase coding sequence . to monitor the distal poly(a ) site use ,
the forward sequence was designed across the artificial junction that was introduced with production of the shortened smc1a 3 utr .
this confirmed the significant increase in the use of the proximal poly(a ) site upon tdp-43 kd ( figure  5b).figure  5smc1a mini - gene(a ) iclip tdp-43 binding along smc1a 3-utr and zoomed in binding upstream of proximal poly(a ) site . the intact smca1 sequence around the proximal polya site is shown in the lower part ( wild - type [ wt ] smca1 ) , which includes the region upstream of the polya signal ( in bold ) that contains a ug - rich region ( in bold ) that crosslinks to tdp-43 ( shown above the sequence with blue bars ) .
we introduced mutations ( mut ) or deletion of tg dimers ( del ) to prevent tdp-43 binding to the ug - rich region in the rna .
the final minigene was designed with mutations that replaced g in tg dimers into t or a to convert the sequence into a binding site for tia proteins ( tia ) .
the mutant nucleotides are marked in red.(b ) ratio in the use of distal vs. proximal polya site in control cells , tdp-43 kd or the double kd of tia1 and tial1 ( tia kd ) .
( a ) iclip tdp-43 binding along smc1a 3-utr and zoomed in binding upstream of proximal poly(a ) site . the intact smca1 sequence around the proximal polya site is shown in the lower part ( wild - type [ wt ] smca1 ) , which includes the region upstream of the polya signal ( in bold ) that contains a ug - rich region ( in bold ) that crosslinks to tdp-43 ( shown above the sequence with blue bars ) .
we introduced mutations ( mut ) or deletion of tg dimers ( del ) to prevent tdp-43 binding to the ug - rich region in the rna .
the final minigene was designed with mutations that replaced g in tg dimers into t or a to convert the sequence into a binding site for tia proteins ( tia ) .
( b ) ratio in the use of distal vs. proximal polya site in control cells , tdp-43 kd or the double kd of tia1 and tial1 ( tia kd ) .
and thus it is theoretically possible that it binds to the longer mrna isoform and stabilizes this isoform , which could explain why the shorter mrna isoform is increased upon tdp-43 kd . to rule out this possibility , we mutated or deleted a tdp-43 binding site that is located upstream of the proximal poly(a ) site ( figure  5a ) .
this binding site is present in both short and long isoforms , and it is unlikely that this site could lead to differential stability of the two isoforms . in contrast
, the site is positioned upstream of the repressed poly(a ) site , thus representing the pattern where the rna map detects most binding enrichment ( figure  3a ) .
tdp-43 binding in this region is ideally positioned to repress the nearby poly(a ) site by blocking the recruitment of cleavage and polyadenylation factors .
indeed , disruption and deletion of the tdp-43 binding site in the minigene caused a strong increase in proximal poly(a ) site use under control conditions , which is comparable with the derepression of the site that is seen upon tdp-43 kd in the wild - type minigene ( figure  5b ) .
tdp-43 kd caused no further effect on the mutated or deleted minigene , confirming that these mutations abolished the capacity of tdp-43 to regulate the proximal poly(a ) site .
another protein , tia1/tial1 , can also bind to u - rich sequences and has been shown to bind to the 3 utr ( wang et  al . , 2010 ) .
to test whether the position - dependent activity of tdp-43 is shared by other rbps , we therefore replaced the ug - rich motifs with a ua - rich sequence that is designed to promote binding of cytotoxic granule - associated rna binding ( tia ) proteins based on its similarity to the tia - binding sites ( figure  5a ) identified previously by iclip ( wang et  al . , 2010 ) .
addition of this site restored repression of the proximal poly(a ) site , and kd of tia1/tial1 proteins confirmed that this repression is caused by binding of tia proteins rather than tdp-43 ( figure  5b ) .
this demonstrates that diverse rbps can bind in the vicinity of the poly(a ) site to block its use .
both tdp-43 and tia proteins are also important regulators of splicing , thus indicating that splicing and 3 end processing could often be regulated by the same proteins . to provide insight into shared mechanisms of splicing and 3 end processing , we used rnamotifs2 to further examine the poly(a ) sites and exons that are differentially used in the brain compared to other tissues .
we first examined the poly(a ) sites that are differentially used between brain and universal human reference  ( derti et  al . , 2012 ) .
this detected several types of motifs , including enrichment of tcat immediately downstream of the proximal poly(a ) sites that are repressed in the brain ( figure  s4e ) , which agrees with the previous finding that nova ( neuro - oncological ventral antigen ) proteins generally repress poly(a ) sites when binding in close vicinity ( licatalosi et  al .
, 2008 ) and indicates that nova proteins promote formation of mrna isoforms with longer 3 utrs , which are known to be more common in neurons ( derti et  al . , 2012 , miura et  al . , 2013 ) . finally , we analyzed motifs at alternative exons that are differentially spliced between brain and heart as defined by splicing microarray ( arrayexpress emtab-1911 ) .
as in our previous study , this detected tc - rich motifs and tcat - related motifs , which correspond to the tc - rich and tcat - related preferences of polypyrimidine tract - binding protein ( ptbp ) and nova proteins , respectively ( cereda et  al .
the positional enrichment agreed with the known rna maps of the two proteins , since ptbp proteins mainly bind upstream of brain - specific exons to repress them , while nova proteins bind upstream and within the exons that are repressed , and downstream of exons that are enhanced in the brain ( ule et  al .
, 2006 , witten and ule , 2011 ) . by detecting the nova - binding motifs both at the poly(a ) sites and exons that are differentially regulated in the brain ,
rnamotifs2 confirmed the important role of nova proteins in both splicing and 3 end processing.figure  6motif clusters regulating alternative splicing in the human brainsignificant motif clusters involved in alternative splicing ( comparing brain and heart tissue ) obtained with rnamotifs2 .
alternative exons flanked by upstream tc and t - rich clusters ( ptbp targets ) are enhanced in the brain , where ptbp expression is low compared to other tissue .
contrary , nova proteins are highly expressed in the brain , and exons flanked by upstream tcat - rich clusters ( nova targets ) are therefore repressed in this tissue . for previous analysis of data ,
motif clusters regulating alternative splicing in the human brain significant motif clusters involved in alternative splicing ( comparing brain and heart tissue ) obtained with rnamotifs2 .
alternative exons flanked by upstream tc and t - rich clusters ( ptbp targets ) are enhanced in the brain , where ptbp expression is low compared to other tissue .
contrary , nova proteins are highly expressed in the brain , and exons flanked by upstream tcat - rich clusters ( nova targets ) are therefore repressed in this tissue . for previous analysis of data ,
our upgraded rnamotifs2 can discover clusters of diverse motifs , and therefore it can provide insights into potential combinatorial regulation of poly(a ) sites or exons .
indeed , it showed that ug - rich motifs tend to cluster close to both uc - rich and ucau - rich motifs around the regulated exons ( figure  6 ) and regulated poly(a ) sites ( figure  s4e ) in the brain . to further examine this clustering
, we plotted the co - occurrence of ug and uc - rich motif clusters at 3 splice sites of silenced alternative exons ( figure  s4f ) in the brain .
moreover , analysis of the alternative exons regulated by ptb , nova , and tdp-43 proteins discloses that ug and ucau - rich motifs tend to cluster together around the exons regulated by nova proteins ( figure  s6 ) . these results indicate that proteins such as tdp-43 might cooperate with tissue - specific rbps to regulate pre - mrna processing ; however , this hypothesis will need further experimental validation .
in this study , we developed several tools to examine how cis - acting elements recruit rbps to regulate splicing and apa in  a transcriptome - wide manner ( figure  s7 ) .
we developed expressrna , a platform that processes 3 mrna sequencing data , classifies the sites where cleavage and polyadenylation takes place ( poly(a ) sites ) , and identifies the differentially regulated poly(a ) sites .
we also integrated the expressrna platform with an upgraded rnamotifs2 software to identify enrichment of combinatorial motif clusters around regulated exons or poly(a ) sites .
we used the motifs and the iclip data to define high - resolution rna maps of apa regulation by tdp-43 .
this showed that tdp-43 can directly regulate both proximal and distal poly(a ) sites of same - exon or skipped classes ( figure  7 , table s6 , p value < 0.01 , brg > 10 ) . when binding close to the poly(a ) signal or the poly(a ) site , tdp-43 represses use of the site , whereas it can enhance the site when binding further downstream .
similarly , tdp-43 generally represses splicing when binding close to 3 splice site or the exon , while enhancing when binding further downstream of the regulated alternative exon ( tollervey et  al . , 2011 ) .
this indicates common position - dependent regulatory principles of both mechanisms.figure  7summary of the main position - dependent modes of apa regulation by tdp-43(a ) the binding patterns that are most enriched for the same - exon type of apa , as shown in figures 3a3d and s3b - e.(b ) the binding patterns that are most enriched for the skipped - exon type of apa , as shown in figures s4c and s4d ; the arrow marks the position of the  regulated poly(a ) site , and the circle marks the  main position of tdp-43 binding .
the blue color  denotes the repressive , and the red the enhancing effect of tdp-43 , and the positions of the regulatory patterns can be found in the table s6 ( p  value  < 0.01 , brg > 10 ) . summary of the main position - dependent modes of apa regulation by tdp-43 ( a ) the binding patterns that are most enriched for the same - exon type of apa , as shown in figures 3a3d and s3b - e .
( b ) the binding patterns that are most enriched for the skipped - exon type of apa , as shown in figures s4c and s4d ; the arrow marks the position of the  regulated poly(a ) site , and the circle marks the  main position of tdp-43 binding .
the blue color  denotes the repressive , and the red the enhancing effect of tdp-43 , and the positions of the regulatory patterns can be found in the table s6 ( p  value  < 0.01 , brg > 10 ) . our findings are consistent with previous studies of rna maps for apa regulation by several other rbps ( batra et  al . , 2014 , licatalosi et  al . , 2008 , liu et  al . , 2013 , masuda et  al . , 2015 ) .
we extend these findings by assessing the regulatory principles at high resolution and in a quantitative manner , which shows the strongest direct role of tdp-43 is in repression of proximal poly(a ) sites , especially by binding in close proximity upstream and downstream of the poly(a ) site .
while pas is the primary element recruiting the cleavage and polyadenylation specificity factor ( cpsf ) , other auxiliary sequences tend to be ug or u - rich in mammalian transcripts ( yang and doubli , 2011 ) .
these include the tgta motif upstream of the pas that recruits the cleavage factor i(m ) ( cfim ) and the downstream ug - rich motifs that recruit the cleavage stimulation factor ( cstf ) ( hoque et  al . , 2013 , martin et  al . , 2012 , yao et  al . , 2012 )
tdp-43 thus appears to act as a competitor to displace cfi and cstf from the ug - rich sites on pre - mrnas .
interestingly , we also find that tdp-43 binding is pronounced further downstream of the enhanced poly(a ) sites , in which cases it might be able to stabilize the binding of processing factors at the nearby poly(a ) site . the rna maps of splicing suggest that tdp-43 might compete with binding of u2af and other splicing factors when binding close to 3 splice sites or recruit u1 small nuclear ribonucleoprotein ( snrnp ) when binding downstream of 5 splice sites .
it remains to be seen how tdp-43 manages to block or enhance pre - mrna binding of such a variety of factors .
the use of rnamotifs2 confirms that apa is , like alternative splicing , often regulated by rbps that bind clusters of closely spaced short motifs on pre - mrnas .
it is clear that the precise position of ug - rich motifs defines the function of tdp-43 in repressing or enhancing either alternative exons or poly(a ) sites ( tollervey et  al . , 2011 ) .
mutation of the ug to ua - rich motifs replaced the function of tdp-43 with tia proteins , which are also otherwise major regulators of splicing .
this indicates that many rbps may have pleiotropic functions in splicing and apa , with the position of their rna - binding site being the primary determinant of their function .
our preliminary analyses indicate that an integrative analysis of tissue - specific splicing and 3 end processing could uncover additional rbps that can regulate both processes .
for example , we detected the known function of the cis - acting element tcat in the brain - specific splicing and apa .
moreover , we find that ug - rich motifs are present within clusters of uc motifs at the tissue - specific exons , indicating a potential for crosstalk between proteins binding these motifs , such as tdp-43 , celf , and ptbp proteins .
enrichment of ug - rich motifs was recently also identified next to binding sites of rbfox proteins ( damianov et  al . , 2016 ) . since our software can detect clusters of diverse motifs , it is well suited for identification of such sites of overlapping binding motifs for multiple rbps , which might allow cooperative or competitive regulation of alternative splicing or alternative polyadenylation across various tissues and conditions .
the platform can examine published poly(a ) sites and exons as well as discover additional ones .
we apply state - of - the - art statistical methodology to identify differentially polyadenylated genes ( dexseq ) . in this context ,
the ability of successfully identifying the regulatory rna maps could also be an estimator ( benchmark ) of the validity and success of the presented approaches .
therefore , expressrna provides a flexible ( modular ) data integrative research platform , making computational analysis highly reproducible and allowing user - friendly visualization with sharing of data and results .
hek293 flpin t - rex cells were maintained in dmem with 10% fetal bovine serum ( fbs ) , supplemented with 3  g / ml blasticidine and 50  g / ml zeocin .
for the small interfering rna ( sirna)-induced tdp-43 kd , 20  nm of tardbp stealth sirna ( invitrogen , a-012394 - 14 , 5-ggcucaagcauggauucua-3 ) was mixed with 10  l of rnaimax following the manufacturer s reverse transfection protocol and added to a 10-cm dish of hek293  flpin cells . for the non - targeting sirna (
stealth rnai sirna negative control med gc content , invitrogen 12935 - 300 ) , 20  nm was also used to distinguish off - target effects from biologically relevant ones .
after the first 24  hr of transfection , the medium was replaced with dmem with 10% fbs , and , after an additional 24  hr , the cells were collected for analysis .
rna was extracted using the direct - zol rna kit , and an in - column dnase digestion step was performed at room temperature for 15  min .
the poly(a)seq libraries for samples that were transfected with either non - targeting sirna or tardbp sirna were generated using the reverse quantseq 3 mrna - seq kit .
one fragment per transcript was generated , which resulted in extremely accurate gene expression values .
for the initial step of this kit , oligodt priming , including illumina - compatible linker sequences , was carried out .
single - end sequencing ( 60 nt ) was performed on a illumina ga-2 with a rapid run flow - cell .
this kit makes use of a custom sequencing primer that anneals to a linker sequence previously introduced in the oligo(dt ) priming step for reverse transcription .
we developed two independent research platforms , expressrna for computational analysis of post - transcriptional modifications including analysis of 3  end sequence data ( expressrna.apa ) and rnamotifs2 for clustered motif analysis ( see rnamotifs2 ) .
the web application part of expressrna ( jquery and javascript ) allows the exploration of combined analysis results over an interactive web interface ( figure  1 ) .
expressrna supports several open - source and commercial 3 end sequencing protocols ( lexogen forward / reverse , 3p - seq , pa - seq , polyaseq ) .
the platform is modular and scalable and runs on desktop computers for  small sequence samples and on multi - core machines for large datasets .
it provides a complete analytical framework incorporating several tools for reads alignment , genome annotation , calling of differentially polyadenylated genes , and integration of rna - protein binding ( iclip ) datasets ( figure  s1 ) .
we established an online server ( http://expressrna.org/paper ) to provide interactive browsing of results presented in this publication .
we processed each of the 12 experimental datasets ( table s1 ) by aligning the reads to the reference human genome ( hg19 ) using star aligner ( dobin et  al . , 2013 ) with default parameters . tagging only one position per alignment (
the first 5 aligned nucleotide ) , we constructed the database of genome - wide polyadenylation events ( figure  2a ) .
since internal priming ( annealing to the genomic sequence instead of the poly(a ) tail ) is a major problem in 3 end sequencing protocols , we checked the genomic sequence in the region [ 10 .. 10 ] surrounding the polyadenylation events and filtered out alignments containing stretches of six consecutive a or with 70% a coverage in any 10-nt sub - window in this region . published studies indicate that auxiliary rna motifs tend to be enriched in the region approximately up to 75 nt upstream ( ure ) and 50 nt downstream ( dre ) of each cleavage site ( beaudoing et  al .
we wished to focus our study on fully independent cleavage sites that contain their own pas and auxiliary motifs .
therefore , we identified the dominant poly(a ) sites based on read count , such that all resulting sites were at least 125 nt apart .
for this purpose , we ranked the polyadenylation events by read count in descending order and considered only the high - ranking events that are more than 125 nt apart .
it has been shown that cleavage is not an exact process ( pauws et  al . ,
2001 ) , and we find that cleavage can occur within a small window of positions around the dominant cleavage site ( figure  s2b ) .
to allow for the variation in cleavage precision ( as seen in figure  s2b ) , we computed the per - experiment expression of each poly(a ) site by summing the read counts that identify any position in the region [ 5 .. 5 ] . to identify regulated poly(a ) sites in genes ,
we input the read counts of all poly(a ) sites remaining after filtering into dexseq .
the analogy to alternative exons counts usually input to dexseq are read counts at poly(a ) sites .
we input count values for all replicates for control and tdp-43 kd 3 end sequencing experiments .
dexseq returns fold change ( log2 ) and adjusted p  value for each site .
the genes where no poly(a ) site reaches significance ( p  < 0.05 ) are classified as controls . in genes with more than one poly(a ) site , only two poly(a ) sites are selected for further analysis . in control genes , the two poly(a ) sites with highest read count
two significantly changed sites ( adjusted p value < 0.05 ) with highest difference in fold change are selected for each gene , additionally requiring that fold changes are of opposite signs .
if the gene was marked as control , the proximal and distal control poly(a ) sites are further labeled as control - down and control - up ( dependent on their fold change ) .
for regulated sites , if a proximal site has fold change < 0 , the site is marked as repressed , and if fold change > 0 , marked as enhanced ( the reverse holds for distal sites ) .
gene site pairs are then used for further analyses , including rna maps , rnamotifs2 , and gene ontology ( go ) term analysis .
after we select the proximal and distal site for each gene , we classify the pair ( and consequently the gene ) into same exon , composite exon , and skipped exon .
for this classification , we use the gene level annotation that is computed by linearizing the ensembl gene annotation by merging the transcript annotation .
when no annotated splice site is present between the selected poly(a ) sites , then the gene is assigned to the same - exon category ( figure  2b ) .
if there is a splice site present between the selected poly(a ) sites , further classification depends on the type of splice site preceding the proximal poly(a ) site .
if there is a 5 splice site immediately upstream of the proximal site , the gene is assigned to the composite - exon category , and , if there is a 3 splice site immediately upstream of the proximal site , the gene is assigned to the skipped - exon category ( figure  2b ) .
rbps play a significant role not only in the regulation of alternative splicing but also in the regulation of apa .
we designed experiments with intact cells as controls and kd samples as test experiments . after identifying alternatively polyadenylated genes , we defined three sets of genes ( repressed , enhanced , and controls ) and cumulatively plotted the iclip data around poly(a ) sites of these three categories .
the position of the polya site ( also referred to as the cleavage site ) is marked with a red line at the center of the rna maps ( figures 3 and s3s5 ) .
the approach for visualizing an rna map of apa is analogous to any other region of interest , such as , for instance , the intron - exon boundaries in the context of alternative splicing ( ule et  al . , 2006 ) .
to display the variability of positional binding enrichment detected by rna maps , we performed 10,000 bootstraps ( sampling with replacement ) of each gene class to obtain the sd ( gray ) .
the colored areas ( red , blue ) on rna maps are drawn between the signal ( repressed , enhanced ) and control regions , excluding sds , which conveys a picture of the most reliable enrichment of candidate regulatory binding sites ( figure  3a ) . to determine the number of regulated genes that contain tdp-43 binding within a specific set of positions , we calculated the brg , where the binding sites are defined either by the significant iclip crosslink clusters or rna motif clusters ( table s6 )
. brg corresponds to the number of genes that contain at least one binding site in a specific window within the rna map at a specific class of poly(a ) sites .
the brg allows a comparative analysis of rna maps , for example , showing the brg in the entire map ( reported brg400 , since the entire rna map window is 400 nt ) .
furthermore , to determine positions on the rna map that have the strongest regulatory significance ( fold change between repressed / enhanced and control genes ) , we performed 1 m bootstraps ( sampling with replacement ) on the entire dataset .
we then used a 40-nt sliding window and computed the p value by considering data versus bootstrap fold changes . for same - exon proximal poly(a ) sites , we found that the most significant 40-nt window is between the 97 and 57 nucleotide relative to the repressed proximal poly(a ) sites ( p value  = 1e-6 , log2 fold change  = 2.04 , brg  = 33 ) .
the positions of the 40-nt most significant windows with respective brgs and fold changes are available for each type of poly(a ) site ( table s6 ) .
we additionally plotted heatmaps of the repressed and enhanced top 20 bound genes , which demonstrate positional contributions of individual sites ( figures 3b , 3c , s3b , and s3c ; table s2 ) .
iclip data were processed as described previously by icount web server ( http://icount.biolab.si ) ( wang et  al . , 2010 ) .
the crosslink clusters were identified considering all crosslink sites that were significant with a false discovery rate ( fdr ) < 0.05 at a maximum spacing of 20 nt between crosslink sites .
each binding site position was then defined by the exact position of the  bound gu - rich motifs 
( gtgtg , tgtgt , tgcgt , tgtgc , cgtgt , atgtg , gtatg , gtgta , gcgtg , gtgcg , tgtga , tgtat , gtgtt , ctgtg , tatgt , ttgtg , tgagt , ggtgt , gagtg , gtgtc , tgtgg , agtgt , gtttg ) that were present in iclip crosslink clusters .
rnamotifs2 is an upgraded version of the previously published rnamotifs software ( cereda et  al . , 2014 ) .
we extended the motif search from alternative splicing regions to three regions surrounding poly(a ) sites and modified the search algorithm with an iterative clustering method .
the approach finds the motif with the strongest signal that differentiates up / downregulated sequences to control sequences first .
next , it removes sequences where the motif signal is strongest and compensates the loss of signal by clustering the strongest motif with additional motifs , iteratively constructing a cluster consisting of up to four different motifs .
the input of rnamotifs2 is a list of genomic positions , either alternative exons or alternative poly(a ) sites .
each input needs to be assigned to one of the three classes : control , enhanced , or repressed . after computing cluster motif analysis ,
the region used for the smc1a minigene is composed of two sequences surrounding both poly(a ) sites .
the first part of smc1a minigene contains  an 812-bp - long sequence surrounding the proximal poly(a ) site ( positions 53406176 .. 53406987 ) , which is coupled to the 537-nt - long 3 utr sequence close to the distal poly(a ) site ( positions 53400970 .. 53401506 ) .
the intervening part of the 3 utr was not incorporated in the minigene reporter due to the limiting size of the plasmid .
the minigene ( 0.3  g ) was co - transfected into hela cells together with sirna against tdp-43 or control or tia , and the cells were harvested after 48  hr .
reverse transcription was performed using revertaid ( fermentas ) , and qpcr was performed to assess the level of different poly(a ) site usage using sybrgreen according to manufacturer s protocol .
the expression of 3 utr of minigene was assessed using forward primer in the luciferase gene and reverse primer in the short 3 utr .
the use of distal poly(a ) site was assessed by using forward primer across the junction between proximal poly(a ) site and the downstream sequence surrounding the distal poly(a ) site , and reverse primer around the distal poly(a ) site .
ratio between the distal and proximal poly(a ) site use in each condition was normalized to the ratio in the control cells .
 summarymany rna - binding proteins ( rbps ) regulate both alternative exons and poly(a ) site selection . to understand their regulatory principles 
 , we developed  expressrna , a web platform encompassing computational tools for integration of iclip and rna motif analyses with rna - seq and 3 mrna sequencing . 
 this reveals at nucleotide resolution the  rna maps  describing how the rna binding positions of rbps relate to their regulatory functions . 
 we use this approach to examine how tdp-43 , an rbp involved in several neurodegenerative diseases , binds around its regulated poly(a ) sites . binding close to the poly(a ) 
 site generally represses , whereas binding further downstream enhances use of the site , which is similar to tdp-43 binding around regulated exons . 
 our rnamotifs2 software also identifies sequence motifs that cluster together with the binding motifs of tdp-43 . 
 we conclude that tdp-43 directly regulates diverse types of pre - mrna processing according to common position - dependent principles . 
Introduction
Results
Discussion
Experimental Procedures
Author Contributions