Anti-angiogenic VEGFAxxxb transcripts are not expressed in the medio-basal hypothalamus of the seasonal sheep

This study investigated Vegfa expression in the pars tuberalis (PT) of the pituitary and medio-basal hypothalamus (MBH) of sheep, across seasons and reproductive states. It has recently been proposed that season impacts alternative splicing of Vegfa mRNA in the PT, which shifts the balance between angiogenic VEGFAxxx and anti-angiogenic VEGFAxxxb isoforms (with xxx the number of amino acids of the mature VEGFA proteins) to modulate seasonal breeding. Here, we used various RT-PCR methodologies and analysis of RNAseq datasets to investigate seasonal variation in expression and splicing of the ovine Vegfa gene. Collectively, we identify 5 different transcripts for Vegfa within the ewe PT/MBH, which correspond to splicing events previously described in mouse and human. All identified transcripts encode angiogenic VEGFAxxx isoforms, with no evidence for alternative splicing within exon 8. These findings led us to investigate in detail how “Vegfaxxxb-like” PCR products could be generated by RT-PCR and misidentified as endogenous transcripts, in sheep and human HEK293 cells. In conclusion, our findings do not support the existence of anti-angiogenic VEGFAxxxb isoforms in the ovine PT/MBH and shed new light on the interpretation of prior studies, which claimed to identify Vegfaxxxb isoforms by RT-PCR.


Introduction
Our understanding of the molecular and neuroendocrine mechanisms that govern seasonal breeding has increased considerably over the last 20 years. In the current model, melatonin acts within the pars tuberalis (PT) of the pituitary to increase Tshb transcription under long days [1][2][3][4]. This PT-derived specific TSH (Thyrotropin-Stimulating Hormone) acts locally on specialized glial cells, known as tanycytes, to trigger the expression of thyroid hormone deiodinase type II (Dio2), which increases T3 levels within the medio-basal hypothalamus (MBH) during long days. This phylogenetically conserved mechanism underlies the critical role of hypothalamic thyroid hormone signaling in the control of seasonal cycles [5][6][7].
Beyond this TSH/DIO2/T3 signaling pathway, it has become clear that season impacts the expression of hundreds of genes in the PT/MBH region. The combination of transcriptomics and gene ontology analyses in sheep has pointed to several key pathways that include cellular plasticity, extracellular matrix remodeling, cell proliferation, epigenetics and angiogenesis [8,9,10]. In an unexpected twist, Castle-Miller et al. [11] recently proposed a role for PTderived VEGFA in the seasonal regulation of pituitary angiogenesis in sheep [11][12][13]. Such changes are proposed to play a role in the control of seasonal breeding through the modulation of pituitary blood flow hence the access of hypothalamic releasing factors such as GnRH to the pars distalis (PD) of the pituitary. Specifically, the authors show that the Vegfa transcript undergoes seasonally-controlled alternative splicing of exon 8, which gives rise either to classical angiogenic Vegfaxxx isoforms (with xxx the number of amino acids of the mature VEGFA proteins) or to anti-angiogenic Vegfaxxxb isoforms [11]. This alternative splicing event of the exon 8 of Vegfa was initially described by Bates et al. for human tissues and would therefore likely be conserved across mammals [14]. Castle-Miller et al. proposed that the balance between both isoforms-angiogenic VEGFAxxx and anti-angiogenic VEGFAxxxb-is integral to the seasonal control of pituitary function [11]. Based on these findings, we investigated the expression profile and splicing of Vegfa transcripts, across seasons and reproductive states in sheep. Here, we report that both Vegfa and VegfR2 display modest seasonal variation in their overall expression level in the PT/MBH, suggestive of enhanced VEGFA signaling during long days. We identify five different pro-angiogenic Vegfaxxx splice variants in the ovine MBH. The main splicing event corresponds to the skipping of exon 6, thereby leading to VEGFA164. We also identify transcripts in which both exons 6 and 7 are skipped (leading to VEGFA120) and transcripts in which cryptic splice sites within exon 6 are used. However, neither our extensive RT-PCR analysis nor our investigation of recent ovine RNAseq datasets [9,10] provided experimental support for the existence of endogenous anti-angiogenic Vegfaxxxb transcripts in the ovine PT/MBH, regardless of season or reproductive status. This suggests that such Vegfaxxxb transcripts are either extremely rare or that they do not exist at all, which makes it very unlikely that such anti-angiogenic isoforms play a significant role in seasonal physiology. Furthermore, using a step-by-step empirical approach to PCR primer design we provide evidence that detection of these isoforms in prior studies, including studies in human cell lines and tissues [14] might be artefactual, as previously suggested by Harris et al. [15].
Collectively, the results of our investigation into Vegfa splicing in the MBH of the sheep concur with prior findings in mouse and humans and suggest that mechanisms are phylogenetically conserved and lead only to the production of pro-angiogenic VEGFAxxx isoforms [12,13,16].

Materials & methods
Apart from experiment 4 (see Experimental animals & procedures), the samples being analysed in the current study were collected in multiple cohorts, which have been used in previous experiments aimed at investigating the roles of photoperiod and thyroid hormone in seasonal breeding. For detailed accounts of these cohorts and experimental procedures, the reader is referred to a recently published study [10].
All surgeries were performed after sodium thiopental anesthesia (Nesdonal1, 1g/80kg), under constant isoflurane administration (Vetflurane1) and all efforts were made to minimize suffering. Following surgery, animals received an injection of antibiotics (oxytetracycline, Terramycine LA1, 1ml/10kg) and an injection of a non-steroidal anti-inflammatory drug (Finadyne1, flumixin megumine, 2ml/50kgs). Animals were followed daily throughout the experiment. Blood samples were collected twice weekly by jugular venipuncture in heparinized tubes. After centrifugation, serum was collected and frozen until assayed for hormones. All animals were killed by decapitation under deep barbiturate anesthesia (Nesdonal1, 5mL). To minimize potential issues linked to time-of-day fluctuation in gene expression (i.e. circadian rhythms) all animals were killed in the early morning (ZT1-4, with ZT0 being the time of lights on, or sunrise).

Experimental animals & procedures
Experiments were conducted in adult Ile-de-France ewes (3-5 years old; weight 60-80 kg) kept under normal husbandry conditions at the research station of the Institut National de la Recherche Agronomique (Nouzilly, Unité Expérimentale PAO n˚1297 (EU0028)). All ewes were ovariectomized (OVX) and estradiol-implanted (E2; 1cm silastic implant) at the beginning of the experiment (three independent experiments; total of n = 74, detailed below). The OVX+E2 model normalizes the level of circulating E2, which uncovers the well-documented central seasonal shift in the negative feedback action of E2 on gonadotropin secretion [17]. In this model, serum levels of the gonadotropins LH and FSH provide a reliable index of the state of the GnRH pulse generator: LH/FSH are low during the non-breeding season, then gradually increase to reach high levels during the breeding season [10,17]. Two groups of intact ewes were also used for a 4 th experiment (total of n = 12, n = 6 / group). In total, four independent experiments were performed. For the 1 st and 2 nd experiments, RNA was extracted from MBH and from the caudal part of the pituitary pars distalis (PD; see S1 Fig for details of the dissection procedure). The MBH blocks not only comprise the PT but also include the median eminence, arcuate nucleus, dorsomedian and ventromedian hypothalamic nuclei. The RNA samples from these 2 experiments were used to carry out RT-PCR and qRT-PCR in the current paper and have also been used previously for RNA-seq analysis as detailed further in this M&M (also see [10]). Importantly, many genes with low level of expression (revealed by the number of reads, RPKM), were identified as being differentially expressed by the RNA-seq, which was confirmed by ISH [10]. Furthermore, ISH showed that each candidate gene is specifically expressed in a distinct compartment of the MBH: PT, tanycytes, arcuate nuclei or dorsomedian hypothalamus [10]. Therefore, our approach effectively allows detection of small variations in expression, even for genes expressed to low levels, and by only a specific subset of cells within the MBH.
For the 3 rd and 4 th experiments ISH was performed on brain sections (see next section and [10]). In the 1 st experiment, 18 ewes were kept in an open barn and killed in May (non-breeding season, n = 6), August (when breeding resumes, n = 6) and November (breeding season, n = 6). In the 2 nd experiment, ewes (total of n = 24; 3 groups with n = 8/group) were kept indoors under a short photoperiod (SP group, 8.5h of light / day) from December onwards and killed in February, or exposed to a long photoperiod (LP group, 15.5h of light / day) for 3 weeks before culling. A third group had been thyroidectomized (THX) months before the LP exposure (LP-THX group). In animals kept either under prolonged SP or exposed to LP, the hypothalamus drives reproductive arrest, as assessed by LH & FSH levels. In contrast, THX ewes maintain high LH & FSH levels (see [10]). In the 3 rd experiment, ISH was performed on brains, which were sampled from ewes maintained outdoors and killed in May and November (total of n = 32; n = 16/group). Half of each group had undergone THX~5 months prior to the end of the experiment (n = 8/group) while the other half had undergone a sham-operation (n = 8/group). THX specifically blocks the spring transition to non-breeding (see [10]). Finally, in the 4 th experiment two groups of intact ewes (total of n = 12; n = 6 / group) were kept outdoors and killed in May and November.

RNA extraction, standard RT-PCR and qRT-PCR analysis
RNA extraction was performed on each MBH block or PD (see S1 Fig) using TriReagent (Sigma). Concentration and purity of individual samples were determined with Nanodrop 2000 (ThermoScientific) and integrity was checked by standard agarose gel electrophoresis. The same procedure was applied for extraction of RNA from HEK293 cells. For both standard RT-PCR and qRT-PCR, cDNA was synthesised using Omniscript RT kit (Qiagen) and Oligo-dT primers (Eurofins, Germany). Standard PCR was performed with Taq DNA polymerase (Qiagen) using an Applied Biosystems 2720 Thermocycler (Thermofisher). Conditions for PCR were as follows: 95˚C for 5 min then 30 cycles: 95˚C for 30s / 58˚C for 30s / 72˚C for 1 min, and a final extension at 72˚C for 10 min. PCR reactions were loaded on a 1.5% agarose gel and migrated at 80-110V in Tris-EDTA-Acetic acid (TEA) buffer for~30min.
For qRT-PCR, the optimal cDNA dilution and calibration curves for Vegfa and VegfR2 primer pairs (Eurofins) were established using cDNA synthesized from an equimolar mix of individual RNA from each experiment. As a negative control, the same mix with water instead of RT was prepared. Quantitative PCR was performed using CFX-96 Real-Time PCR Detection System (Bio-Rad) and Sso Adv Universal SYBR Green Supermix (Bio-Rad). All samples (unknown, standard curves) were assessed in triplicate and Rplp0 (ribosomal subunit P0, a.k.a 36B4) was used as a housekeeping gene; Rlpl0 expression did not display any statistical differences between groups (data not shown). The quantification of mRNA level was obtained by the 2 -DCT method and data are presented as fold-increase compared to the condition with the lowest expression level. For both standard PCR and qPCR, the identity of amplicons was validated by cloning (pGEMT easy vector, Promega) and Sanger sequencing (Eurofins Sequencing services, Germany).

Analysis of Illumina RNA-seq datasets
We analyzed two independent RNA-seq datasets generated by Illumina sequencing. The first dataset was generated using the same MBH RNA samples which are used here to perform qRT-PCR as well as RT-PCR. All details regarding the procedures, thorough validation and outcomes of these RNA-seq experiments (comparisons between May/August/November = Expt 1 under M&M and comparisons between SP/LP/LP-THX = Expt 2 under M&M) have recently been published [10]. The second dataset was independently generated using RNA extracted from the PT of castrated rams (see details in [9]). First, the data were re-mapped to Oar_V4.0 (ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/vertebrate_mammalian/Ovis_aries/ latest_assembly_versions/GCA_000298735.2_Oar_v4.0) using STAR RNA-seq aligner and standard mapping options [18]. STAR allows the identification of de novo and annotated splice junctions and the algorithm uses quality scores to increase its efficacy. Annotation from ftp:// ftp.ncbi.nih.gov/Ovis_aries/GFF/ref_Oar_v4.0_top_level.gff3.gz were used to define known junctions. Novel junctions are identified by STAR based on uniquely mapping reads. Using a custom perl script the junctions were filtered for the Vegfa gene. The data represent all the junctions identified with >10 reads.

Data analysis
Data were analysed using GraphPad Prism and are reported as mean ± sem. qRT-PCR data were analyzed by 1-way ANOVA using time or treatment as variables. The post-hoc Tukey test was used for mulitple comparisons. ISH data were analysed by either t-test or 2-way ANOVA using time and treatment as variables. p<0.05 was considered significant.

Results
A complete list of the primers used in this study is provided in Table 1.

Seasonal expression of Vegfa and VegfR2 in the ovine MBH and pars distalis (PD)
Quantitative RT-PCR or in situ hybridization (ISH) were used to determine level and profile of expression of Vegfa and its cognate receptor VegfR2 in the MBH and PD in 4 experiments.
In the 1 st experiment qRT-PCR was used to evaluate expression of Vegfa and VegfR2 in MBH tissue blocks collected in May (non-breeding season), August (breeding season resumes) and November (breeding season), with n = 6 for each group. Such a block comprises the PT, median eminence, arcuate nucleus, dorsomedian and ventromedian hypothalamic nuclei. Details of the tissue sampling are provided in S1 Fig. Both Vegfa and VegfR2 showed a~2-fold decrease in expression between May and August and remained low through to November ( Fig  1A).
In the 2 nd experiment, we assessed Vegfa and VegfR2 expression in the MBH block and in the PD of the pituitary of 3 groups of ewes (n = 8 / group) kept under SP or exposed to an acute LP stimulus (both sham-operated and THX): neither of these genes showed any difference in expression levels amongst tissues or groups ( Fig 1B).
In the 3 rd experiment, ISH was performed on hypothalamic sections from ewes also sampled in May and November (n = 16/group). In addition, half of each group had undergone THX~5 months prior to the end of experiment. The use of ISH on coronal sections showed strong Vegfa expression in the PT with no difference in relative levels amongst groups ( Fig 1C).
In the 4 th experiment intact ewes were killed either in May or in November. Here again, ISH on coronal sections showed strong Vegfa expression in the PT, which was not different between May and November as revealed by quantification on autoradiograms ( Fig 1D).
In conclusion, Vegfa and VegfR2 are expressed in the MBH and PD of the ewe. Both Vegfa and VegfR2 display moderate seasonal modulation in the MBH, with higher levels during the non-breeding season. Within the MBH, the PT appears to be the main site of Vegfa expression, with no obvious correlation between relative levels and the reproductive status. In these experiments, the PCR primers for qRT-PCR and the probe used for ISH were designed to detect all Vegfa transcripts; i.e. they would not allow discrimination between Vegfaxxx and Vegfaxxxb isoforms. We therefore designed approaches to specifically investigate expression of Vegfaxxx / Vegfaxxxb splice variants across seasons and physiological states.

Primer (strand) 5' -> 3' sequence Function
Vegfa primers  Fig 2B); (ii) XM_012100431, which contains all 8 exons but uses an alternatively spliced form of exon 6, which keeps the ORF and yields VEGFA182; (iii) XM_012100432, in which exon 6 skipping leads to VEGFA164 and (iv) NM_001025110, in which both exons 6 and 7 are skipped, thereby leading to VEGFA120. The existence of one or several additional transcripts, Vegfaxxxb, has been proposed in sheep [11]. Such transcript(s) would comprise an alternatively spliced exon 8 (exon 8b), due to the use of a cryptic acceptor splice site, located downstream of the known acceptor site (Fig 2B,  2C & 2D). Such transcripts would therefore encode VEGFAxxxb proteins with a CLTRKD and PD (right column) cDNA samples from OVX+E2 ewes, which had been kept indoors under prolonged SP or exposed to 3 weeks of LP. A group of LP-exposed ewes had been THX months before photoperiodic transfer (LP-THX). C/ 3 rd Experiment: ISH for Vegfa on coronal brain sections at the level of the caudal PT / infundibulum region. Shamoperated and THX ewes (OVX+E2 model) were sampled in May and November and quantification was performed in the PT. D/ 4 th Experiment: ISH for Vegfa on coronal brain sections at the level of the caudal PT / infundibulum region. Intact ewes were sampled in May and November and quantification was performed in the PT. Representative autoradiograms are shown (scale bar = 2mm). https://doi.org/10.1371/journal.pone.0197123.g001 motif on their C-Term, rather than the CDKPRR motif (Fig 2B & 2C). We first aimed at the identification of splicing variants of the Vegfa gene using a RT-PCR methodology suitable for amplification of both Vegfaxxx and Vegfaxxxb transcripts.

Identification of ovine Vegfa transcripts by standard PCR reveals 5 Vegfaxxx isoforms
We first performed standard RT-PCR using an equimolar mix of RNA extracted from the 18 MBH blocks from OVX+E2 ewes of the 1 st experiment described above (May/August/ November). Therefore, our cDNA preparation is expected to contain a representative set of seasonal Vegfa transcripts, should there be variation in their relative proportions throughout the year. As can be deduced from Fig 2, Vegfaxxxb PCR products would be expected to migrate approximately 66bp below their Vegfaxxx counterparts.
We performed four different PCR reactions using combinations of two upstream primers, located within exon 3 and exon 4, respectively (O12I & O14I, Fig 2B, Fig 3A, Table 1) and two downstream primers, both located in exon 8, >120bp downstream of the predicted alternative splice site yielding exon 8b (O15I & O16I, Fig 2B, Fig 3A, Table 1). As assessed by agarose gel electrophoresis, all four PCR combinations yielded a distinctive "3-bands pattern" (b1, b2 and b3; Fig 3A), consistent with the existence of at least three different transcripts for Vegfa. The three bands obtained with the O14I/O15I primer pair (Fig 3A) were extracted from the gel, cloned and sequenced (Fig 3B).
A total of 84 clones were sequenced: 15 clones for b1, 44 clones for b2 and 25 clones for b3 (Fig 3B). This revealed that band1 comprises at least three distinct transcripts: b1a (500bp; 10 clones) corresponds to XM_012100430; b1b (482bp; 3 clones) corresponds to XM_012100431; b1c corresponds to a novel transcript (492bp; 2 clones) which uses a cryptic acceptor site at the beginning of exon 6 thereby removing 8bp. As a consequence, this b1c transcript would encode a truncated protein product. The 44 clones for Band 2 represent a single transcript, which corresponds to XM_012100432 (428bp, exon 6 skipped). Finally, the 25 clones for Band 3 also represent a single transcript, which corresponds to NM_001025110 (296bp, both exons 6 and 7 skipped). All clones sequenced and transcripts identified correspond to the usage of the "classical" splice site for exon 8.

Specifically searching for Vegfaxxxb transcripts with isoform-specific PCR primers
We reasoned that failure to identify Vegfaxxxb transcripts might be due to a low level of expression and/or a very tight, season-specific, expression pattern. Therefore, we used cDNA pools from each group of animals of the 1 st experiment (i.e. May, August and November, n = 6/condition; Fig 1A) rather than a cDNA pool generated with a RNA mix from all 18 animals. In a strategy broadly similar to that used by Harris et al. for mouse and human tissues [15], we designed two Vegfaxxxb-specific downstream primers: O21I comprises the last 5 bases of exon 5 and the first 17 bases immediately downstream of the splice site for putative exon 8b while O22I comprises the last 5 bases of exon 7 and the 14 bases immediately downstream of the splice site for putative exon 8b. The last 5 bases of exon 7 (5'-TGCAG-3') are identical to the last 5 bases of exon 4 and 3 (see Table 1).
Therefore, O21I is designed to amplify Vegfaxxxb transcripts that lack both exons 6 and 7 while O22I would amplify transcripts that include exon 7. Based on findings for Vegfaxxx, this primer design encompasses the complete set of putative Vegfaxxxb transcripts. We designed an additional downstream primer (O23I), which bridges the putative splice junction for exon 8b. This primer would amplify all Vegfaxxx isoforms and therefore serves as a positive control. These primers were used in combination with the upstream primer O12I for all PCR reactions. Whether performed with cDNA from animals culled in May, August or November, PCR reactions using Vegfaxxxb-specific primers (O21I & O22I) did not yield any PCR product (Fig 4A). In contrast, PCR with the O12I/O23I pair resulted in a very robust "3-bands pattern" of amplification.
The same PCR reactions were also performed using cDNA synthesized from MBH RNA extracts from the 2 nd experiment (see Fig 1B). Here again, the O12I/O23I pair led to robust PCR amplification of the "3-bands pattern" while the Vegfaxxxb-specific primers O21I & O22I did not yield any PCR product (Fig 4B). Analysis of RNA-seq datasets: Strong support for the existence of Vegfaxxx splice variants but not for Vegfaxxxb isoforms We next analyzed two independent RNA-seq datasets generated by Illumina sequencing. A summary of the results is presented in Table 2 Table 2 and Fig 2). Importantly, 408 of these UMR cover the junction between exon 5 and exon 8 (highlighted in grey) while 2285 cover the junction between exon 7 and exon 8 (also highlighted in grey). Thus, a total of 2693 UMR support the existence of the canonical splice site for exon 8. In contrast, no UMR supports the existence of the alternative splice site yielding the putative exon 8b.
To determine whether prior investigations identified "Vegfaxxxb-like" PCR products due to methodological issues in primer design we aimed to define the minimal primer requirements that lead to amplification of "Vegfaxxxb-like" PCR products.

How many bases complementary to exon 7 to yield "Vegfaxxxb-like" PCR products?
We focused our analysis on transcripts that would include exon 7, since these represent the majority of Vegfa transcripts (b1 and b2, see Fig 3B and Table 2). We first wished to define how many bases specific to exon 7 are required to yield "Vegfaxxxb-like" PCR products. As demonstrated by the use of O22I (Fig 4), 5 bases is not enough. Using an approach similar to that used by Harris et al [15], we designed downstream primers with a 3' end that includes 7, 9 or 11 bases complementary to the 3' end of exon 7 (O61I, O62I & O63I), while adjusting the length of the "exon 8bspecific" 5' part of the primer to maintain similar melting temperatures (Tm) across primers.
PCR with O23I (positive control, see Fig 4) led to the typical "3-band pattern" (Fig 5A), while PCR with O22I (negative control, see Fig 4), O61I and O62I yielded at best very faint PCR amplification. In contrast, PCR with O63I resulted in robust amplification of "Vegfaxxxb-like" PCR products, corresponding to b1 and b2. We therefore surmised that the 11 bases at the 3'end, which hybridize to exon 7, suffice to get robust PCR amplification. Accordingly, we designed O64I, which retains these 11 bases while the other 9 bases of O63I are scrambled. Interestingly, very little PCR amplification was obtained with O64I compared to O63I (Fig 5A). In sum, to yield "Vegfaxxxb-like" PCR products, the "Vegfaxxxb-specific" downstream primer must include a minimum of 11 bases complementary to exon 7. Additionally, the 9 bases of the "exon 8b-specific" 5' part of the primer also contain feature(s) required for efficient PCR amplification.

Critical feature(s) of the 9 bases of the "exon 8b-specific": Ruling out cooperativity
We first tested the hypothesis that the 9 bases of O63I are involved in some form of cooperativity: since these 9 bases are complementary to the Vegfa mRNA sequence (contrary to the non-combinations. The 3 bands, labeled b1, b2 and b3 (arrows), were individually extracted on gel for the primer combination O14I/O15I. B/ Schematics summarizing the results of Sanger sequencing after cloning of PCR products contained in b1, b2 and b3 (location of primers O14I/O15I is provided). Five distinct amplicons were identified: 3 for b1 (labeled b1a, b1b and b1c) and a single product for both b2 and b3. The number of clones sequenced for each amplicon is indicated on the right. An annotated sequence of ovine exon 6 is provided and the alternative splice sites used to generate the transcripts b1b and b1c are identified. Note that 0/84 clones included the putative exon 8b. specific scrambled sequence of O64I) they might hybridize to the Vegfa transcript, thereby bringing the 3'end of the primer closer to its target sequence (i.e. 3' end of exon7), which may indirectly promote PCR. To test this hypothesis, we designed four primers (O65I, O66I, O67I and O68I), which include a stretch of 9 bases complementary to sequences contained within exon 8-8b (Fig 5B, sequences in blue and underlined in bottom panel). We kept the AT/GC content of these 9 bases identical to that of O63I (Fig 5B, sequence in green and bold in bottom panel). As seen in Fig 5B, none of these primers led to any significant PCR amplification compared to O63I. We therefore ruled out some form of "cooperativity" as a potential explanation for the amplification specifically seen when using O63I.

Critical feature(s) of the 9 bases of the "exon 8b-specific": Exon 8 and putative exon 8b share a common 4 bases motif at their 5' end
A closer look at exon 8 sequence revealed that the four bases immediately following the splice sites of exon 8 and putative exon 8b are identical: spliced transcripts would both start with the Searching for Vegfaxxxb transcripts with isoform-specific PCR primers. A/ Standard PCR / agarose gel electrophoresis was performed using MBH extracts from ewes culled across seasons (same samples used for qPCR data in Fig 1). O21I and O22I are Vegfaxxxb-specific primers while O23I spans the putative splice site for exon 8b (and is therefore a Vegfaxxx-specific primer). B/ Standard PCR / agarose gel electrophoresis was performed using MBH extracts from ewes culled under SP or after a photoperiodic transfer to LP (intact and THX ewes; same samples used for qPCR data in Fig 1).
https://doi.org/10.1371/journal.pone.0197123.g004 same ATGT motif (bold green, Fig 6A). As a consequence, exon O63I not only comprises 11 bases complementary to exon 7 but also includes 4 bases, which would correspond to either exon 8 or putative exon 8b. Consequently, the 15 bases at the 3' end of this primer would be complementary to either Vegfaxxx or Vegfaxxxb transcripts.
We reasoned that these 4 bases might explain why O63I leads to PCR amplification while other primers do not. As a first step, we designed primers O73I, O74I, O75I and O76I, whose sequence corresponds to the last 15, 14, 13 and 12 bases of the 3'end, respectively. As seen in Fig 6B, no PCR amplification was obtained with these primers. This was not unexpected since these PCR primers are very short and therefore have quite a low theoretical Tm: 44˚C for O73I, down to 36˚C for O76I (using G/C: 4˚C, A/T: 2˚C). We hypothesized that the 5 extrabases missing from O73I are not required for specificity (i.e. primer hybridization) but that their mere presence increases the Tm of the primer, which allows PCR to proceed. If true, addition of virtually any 5 bases upstream of the 15 bases stretch might support PCR amplification. To test this, we designed three primers derived from O73I, with variable stretches of 5 bases at their 5' end: O77I, with a 5'-CACTC-motif (mirror of the 5-GTGAG-motif of O63I), O78I with a 5'-GGGGG-stretch and O79I with a 5'-AAAAA-stretch. All three primers were as efficient as O63I (Fig 6B). This unambiguously validates our hypothesis that the extra bases on the 5'end do not confer specificity but are merely required to increase Tm.
Next, we wondered if, in the context of a 20 bases primer, all 4 bases of the ATGT motif (5'-ACAT for the downstream primer) are actually required for PCR amplification. To test this, we designed primers with point mutations. As seen before, robust PCR amplification was observed using O15I and O63I as positive controls (Fig 6C). Interestingly, PCR also worked with O80I, in which the last base of the motif is mutated to yield a 5'-TCAT motif. In contrast, mutating the third (5'-AGAT, O81I) or both the third and fourth bases of the motif (5'-TGAT, O82I) severely impaired PCR. Interestingly, a primer with a 5'-CACTCGGAT motif (O90I), which spares only the first two bases of exon 8 led to substantial PCR amplification.
We conclude that amplification of "Vegfaxxxb-like" PCR products may be obtained through primer design. We identified the minimal features of an error-prone PCR primer ( Fig  7A): a 3'end with 11 bases complementary to exon 7, preceded by a minimum of 2-3 bases complementary to either exon 8 or "putative exon 8b" and the addition of at least 5-6 bases Table 2. RNA-seq identifies intron-spanning, uniquely mapping reads at the Vegfa gene.

Intron position Uniquely mapping reads (UMR)
Ovis aries Vegfa gene NC_019477.2 1 st base Last base Exon junction # band (Fig 2) Lomet et al. [10] Wood et al. [9]  (any bases) at the 5'end. The 5' end appears to be solely required to increase Tm, which allows PCR to proceed.

Vegfaxxxb transcripts as PCR artefacts: Extending findings to human
Interestingly, the sequence of the error-prone downstrem PCR primer (exemplified by O63I) identified through our empirical approach is virtuallly identical to the one proposed by Bates Fig 5. Defining features of an error-prone primer yielding "Vegfaxxxb-like" products. Standard PCR / agarose gel electrophoresis was performed using MBH extracts from ewes culled across seasons. A/ 20 bases downstream primers including a variable number of bases complementary to exon 7 were tested. O23I was used a positive control. A minimum of 11 bases complementary to exon 7 was required to yield "Vegfaxxxb-like" PCR products. The sequence of the 9 bases complementary to putative exon 8b also matters (9 bases scrambled primer O64I, see text). Neg: cDNA was omitted from the PCR reaction. B/ PCR was performed with various 20 bases downstream primers, which all contain the 3'end with the 11 bases complementary to exon 7 and a 5'end with 9 bases stretches complementary to variable sequences within exon 8. Bottom panel provides the sequence of exon 8. Locations of the 9 bases stretches used for designing 5'ends of O65I-O68I are underlined in blue. A PCR product was obtained only with the "exon 8b-specific" O63I primer; the 9 bases stretch of this primer is indicated in bold green. O23I was used a positive control.
https://doi.org/10.1371/journal.pone.0197123.g005 et al. [19] for specific amplification of Vegfaxxxb transcripts in human (O81I; Fig 7B). Indeed, Vegfa sequence is broadly conserved amongst vertebrates and Exon 8 and Exon 8b in hVegfa would both start with an AT motif (Fig 7B). We therefore tested the hypothesis that amplification of "Vegfaxxxb-like" PCR products may be a PCR artefact using human cDNA from HEK293 cells. As seen in Fig 7C, only faint PCR amplification was observed with the "exon 8b-specific" O81I primer, which would be entirely specific (20 bases / 20 bases) to "Vegfaxxxblike" transcripts, including the 9 bases complementary to the 5'end of putative exon 8b ( Fig   Fig 6. Exon 8 and putative exon 8b share a common 4 bases motif at their 5' end. A/ Nucleotide sequence of ovine exon 8, splice sites are indicated. Both exon 8 and exon 8b would start with the same 4 bases motif ATGT, in green and bold. B/ Standard PCR / agarose gel electrophoresis was performed using MBH extracts from ewes culled across seasons. The use of downstream primers, which all contain the 3'end with the 11 bases complementary to exon 7 and a 5'end that includes these 4 bases or less (O73I-O76I), did not yield any PCR products. However downstream primers, which all contain the 3'end with the 11 bases complementary to exon 7, along with a 5'end that includes these 4 bases preceeded by a random stretch of 5 bases (O77I-O79I), all yielded robust amplification. C/ Single or double point mutation were made to the ATGT motif in the backbone of "exon 8b-specific" O63I primer. Note that a single point mutation (ATCT, O81I) is enough to blunt PCR (see text for details). 7C, sequence provided above the gel). In contrast, very robust amplification was obtained with O63I and O80I, which both carry a single mismatch to putative "Vegfaxxxb-like" transcripts. Interestingly, such point mutations result in primers that respectively possess 4 and 3 bases complementary to the proximal exon 8 motif (indicated in bold green in Fig 7B, Fig 7C), while O81I retains only 2 such bases (AT). Furthermore, the non-specific O90I primer (with a 5'-CACTCGGAT motif), which solely retains the AT motif consistently gave a stronger PCR amplification than the "exon 8b-specific" O81I primer.
In conclusion, our data are fully consistent with the notion that the weak PCR amplification obtained with the "Vegfaxxxb-specific" primer O81I arises from hybridization to the 5' of exon 8 of angiogenic Vegfaxxx transcripts, partly due to the common AT motif. The stronger PCR amplification observed with O80I and O63I, which both bear mismatches to Vegfaxxxb sequence, yet support robust PCR amplification is in line with this hypothesis (i.e. amplification of Vegfaxxx transcripts), as is the consistent amplification obtained with O90I.

Discussion
Our analysis of Vegfa expression in the PT/MBH of sheep, across seasons and reproductive states, using RT-PCR and analysis of RNAseq datasets, led to the identification of five Vegfa splice variants. These different splice variants have been previously identified in mouse and human tissues, which suggests phylogenetical conservation of Vegfa splicing events. All the sequences identified in this work correspond to Vegfaxxx transcripts, with no evidence for usage of an alternative splice site in exon8. This suggested that identification of anti-angiogenic Vegfaxxxb transcripts in prior studies may have been caused by methodological issues with PCR. We therefore devised an empirical approach, which led to the identification of the minimal sequence requirements of an error-prone "exon 8b-specific primer".
The VEGFA165b isoform was initially identified by Bates et al. [14] in human renal carcinomas. In striking contrast to other VEGFA isoforms, which promote angiogenesis, functional tests showed that VEGFA165b inhibits angiogenesis. Since then, additional VEGFAxxxb isoforms have been identified, which are all splice variants that comprise or lack exons 6 and/or 7 [19]. In human, isoforms that include exon 8a harbor a CDKPRR sequence at their C-term while isoforms that include exon 8b bear the SLTRKD sequence. It has been proposed that this difference confers either angiogenic or anti-angiogenic properties (VEGFAxxx and VEG-FAxxxb isoforms, respectively). Consequently, it has been suggested that the balance between these two families of VEGFA sets the angiogenic outcome [20,21]. This has considerable therapeutic implications as current anti-VEGFA treatments against multiple froms of cancer (Bevacizumab, Avastin1) or diabetic retinopathies and macular degeneration (Ranibizumab, Lucentis1) potentially target both pro-and anti-angiogenic forms of VEGFA [13,22]. Using an anti-VEGFAxxxb antibody developed against human isoforms and RT-PCR, Castle-Miller et al. [11] recently showed that anti-angiogenic VEGFAxxxb isoforms are also present within the ovine PT and PD. They reported that Vegfa transcripts undergo a seasonal switch in splicing, such that production of VEGFAxxxb is increased during the breeding season and leads to reduced angiogenesis. According to Castle-Miller et al., these anti-angiogenic Vegfaxxxb transcripts represent a substantial part of the total pool of Vegfa mRNA, as they are up to 4 times more abundant than Vegfaxxx mRNA in the ovine PT/ME [11]. This finding in sheep is in line with prior reports in various normal human tissues, wherein Vegfaxxxb transcripts may account for >70% of total Vegfa [20,21,[23][24][25].
To investigate this splicing mechanism we took advantage of RNA samples and brain sections, which were previously used to establish the ovine MBH transcriptome across seasons and reproductive states (Fig 1, [10]). Our qPCR data in the MBH, using non-isoform specific PCR primers which allow amplification of Vegfaxxx and Vegfaxxxb isoforms, show that both Vegfa and VegfR2 transcripts undergo a~2-fold decrease from the non-breeding season (May) to the breeding season (August/November). This is suggestive of enhanced VEGFA signaling during long days. However, neither Vegfa nor VegfR2 are acutely regulated by long days or by THX, which blocks the transition to anestrus at the end of winter. Consequently, there is no obvious correlation between levels of Vegfa or VegfR2 and the reproductive state. Our ISH data refined the pattern of Vegfa expression, showing it is mostly expressed in the PT, and its expression is not affected by photoperiod or THX. The apparent inconsistency between qPCR and ISH data might be due to the relatively small transcriptional changes (~2-fold) and to the semi-quantitative nature and lower sensitivity of ISH compared to qPCR. Overall, our data are consistent with prior findings, by both Jabbour et al. [26] and Castle-Miller et al. [11] which show that total Vegfa mRNA and VEGFA protein levels do not display large seasonal variation in the PT.
We then used non-isoform-specific primers to allow RT-PCR amplification and characterization of both Vegfaxxx and putative Vegfaxxxb transcripts (Fig 3). The sequences identified are consistent with predicted sequences in databases and findings in other species, including human and mouse, which suggest phylogenetical conservation of Vegfa splicing events [13,16]. As assessed by RT-PCR and agarose gel electrophoresis, the transcripts with exon 6 skipping represent the major isoform (b2 band). We also identify four minor isoforms. In the shortest (b3 band), both exons 6 and 7 are skipped. For the three longer isoforms, exon 7 is retained while exon 6 may be retained, or undergo usage of cryptic donor or acceptor sites. This latter isoform (denoted b1c), which was not predicted in sheep, has recently been identified in human tissues [16]. Therefore, our approach successfully identified known, novel, and, rare Vegfa splice variants, all of which correspond to angiogenic Vegfaxxx isoforms as none of the 84 clones included exon 8b. These findings are not congruent with findings by Castle-Miller et al. that Vegfaxxxb transcripts are very abundant in the MBH [11].
In contrast, our data are fully consistent with findings by Eswarappa et al. [27], who used a similar RT-PCR approach in the cow, for which the proximal part of exon 8 is virtually identical to that of sheep. Echoing our findings, out of 74 cDNA clones sequenced, none included the putative exon 8b [27]. These authors proposed that anti-angiogenic isoforms may nevertheless arise as a consequence of programmed translational read-through leading to another isoform, dubbed VEGFA-x [27]. These isoforms may help explain the immunostaining obtained with the anti-VEGFAxxxb antibody, in spite of the absence or scarcity of specific Vegfaxxxb transcripts. However, the existence of this VEGFA-x isoform is disputed and its anti-angiogenic properties could not be replicated [28]. Furthermore, VEGFAxxxb isoforms have been found to be weakly angiogenic compared to VEGFAxxx isoforms, rather than antiangiogenic [29,30]. Therefore, the existence of anti-angiogenic VEGFA isoforms remains controversial.
Interestingly, the anti-VEGFAxxxb antibody used by Castle-Miller et al [11] in sheep was initially developed against an epitope of 6 amino acids, SLTRKD, which corresponds to the C-Term shared by all putative human VEGFAxxxb isoforms. The C-Term sequence of presumptive VEGFAxxxb isoforms in sheep is slightly divergent, as these isoforms would harbor a CLTRKD sequence. Therefore, we infer that the epitope recognized by the anti-VEGFAxxxb antibody might be restricted to the sequence LTRKD. Considering that Vegfaxxxb transcripts and isoforms are non-existent or scarce, an alternative explanation for the VEGFAxxxb-like immunostaining might be cross-reactivity of the antibody with a protein bearing a similar epitope. To substantiate this hypothesis, we searched human protein databases using blastp (protein-protein Blast, see [31]), which led to the unambiguous identification of 10 proteins harboring an LTRKD motif (S1 Table). Interestingly, this anti-VEGFAxxxb antibody was also used to detect VEGFAxxxb isoforms in tissues from mouse [32] and rat [33,34], while these isoforms would harbor a distinct PLTGKTD motif at their C-term, as already pointed by others [15,35]. Taking these considerations into account, it is not possible to rule out that the anti-VEGFAxxxb antibody, which has been extensively used in multiple species to define the roles of VEGFAxxxb isoforms, may cross-react with a substantial number of proteins (>10). This calls for the development of more specific tools.
It has been proposed that the RT-PCR approach described above, making use of non-isoform-specific primers, may be biased towards amplification of certain Vegfa transcripts with respect to their abundance or potential issues with secondary structures [20]. We therefore designed isoform-specific primers, to specifically target and amplify Vegfaxxxb transcripts (Fig  4). We tested these primers on cDNA samples obtained across seasons and reproductive states. No PCR product could be obtained. Taken together, these data reinforce our conclusion that Vegfaxxxb transcripts are absent from the ewe MBH, or expressed to very low levels, i.e. below the sensitivity threshold of our RT-PCR methodology. As alluded to above, our findings in sheep echo multiple failed attempts by different research teams to reveal Vegfaxxxb mRNA in various animal models (cow, mouse and human) and using multiple methodologies [15,16,27,[35][36][37]. In the most recent study, Bridgett et al. [16] used publicly available RNA-seq datasets to search for Vegfaxxxb transcripts in human tissues. Their extensive analysis identified >40000 Vegfa transcripts that cover the junction between exon 7 and exon 8 in >10 different human tissues. None of the reads supported the existence of an exon 8b splice site, hence the existence of Vegfaxxxb mRNA [16]. Our current analysis of independent MBH/PT RNAseq datasets, which were generated in ewes [10] and rams [9], fully concurs with findings by Bridgett et al. in human tissues (Table 2). This provides further support to the notion that previous detection of Vegfaxxxb mRNA may be accounted for by methodological issues. Indeed, Harris et al [15] had already dedicated a thorough study to the identification of Vegfaxxxb in mouse and human tissues and cells and concluded that detection of Vegfaxxxb mRNA might arise as a PCR artefact through "5' tailing".
In reply to the study by Harris et al. [15], Bates et al. [19] argued that inability to detect Vegfaxxxb mRNA was due to inadequate primer design and lack of appropriate controls. According to Bates et al., the 3'end of the adequate Vegfaxxxb-specific primer includes 11 bases complementary to exon 7, while the 9 bases at the 5'end are complementary to the putative exon 8b [19]. The authors validated their primer specificity using a construct that bears the putative Vegfaxxxb sequence. However, the authors did not valide their primers on tissue cDNA, nor did they assess whether the 9 bases on the 5'end provide specificity.
With these considerations in mind, we used a step-by-step approach to identify an "exon 8b-specific" downstream primer able to yield a "Vegfaxxxb-like" PCR product (Figs 5 and 6). Collectively, our experiments define the minimal requirements for such a primer: a 3'end with 11 bases complementary to exon 7 followed by the first 2-3 bases of exon 8, which together provide specificity, then a stretch of random nucleotides, which increases the melting temperature of the oligonucleotide (Fig 7). This error-prone primer appears identical to the one proposed by Bates et al. [19] to specifically detect Vegfaxxxb RNA in human tissues, which questions interpretation of prior findings. Overall, the RT-PCR approaches we developed using sheep MBH are quite similar to those developed for mouse and human tissues and cells (including HEK293) by Harris et al. [15] and our conclusions concur with theirs: "VEGFxxxb products were consistently amplified only when the reverse primer contained more bases of complementary sequence to exon 7 or exon 5 as opposed to exon 8b (i.e. 5' tailing)." In this study [15], amplification of Vegfaxxxb-like transcripts in mouse and human could be obtained when the number of bases complementary to exon 7 or exon 5 was equal or >13, which is in line with our current findings in sheep.

Conclusions
Our results (i) show that the ovine PT may use enhanced VEGFA signaling during long days, the consequences of which in terms of breeding control-if any-remain to be established, (ii) reveal the existence of at least five Vegfa splice variants in the ovine MBH, all of which belong to the pro-angiogenic Vegfaxxx family and (iii) unveil issues in Vegfaxxxb detection methodology by RT-PCR. In conclusion, our findings do not support the existence of anti-angiogenic Vegfaxxxb isoforms in the ovine PT/MBH and shed new light on the interpretation of prior studies, which claimed to identify Vegfaxxxb isoforms by RT-PCR.  Table 2 for details). The Ovis aries v4.0 genome and annotation was uploaded into the Broad Institute IGV browser (http://software.broadinstitute. org/software/igv/). STAR junction tables were converted to bed format using a custom bash script and the resulting bed file was uploaded along with the BAM and BAM index files for a representative sample. (PDF) S1