Yon Ho Jee

Florian Thibord

Alicia Dominguez

Corriene Sept

Kristin Boulier

Vidhya Venkateswaran

Yi Ding

Tess Cherlin

Setia Shefali

Valeria Lo Verma

Faro

Verma

Valeria Lo Faro

Traci M Bartz

Anne Boland

Jennifer A Brody

Jean-Francois Deleuze

Joseph Emmerich

Marine Germain

Andrew D Johnson

Charles Kooperberg

Pierre- Emmanuel Morange

Nathan Pankratz

Bruce M Psaty

Alexander P Reiner

David M Smadja

Colleen M Sitlani

Pierre Suchon

Weihong Tang

David-Alexandre Trégouët

Sebastian Zöllner

Bogdan Pasaniuc

Scott M Damrauer

Serena Sanna

Harold Snieder

Cohort Lifelines

Christopher Study

Kabrhel

Nicholas L Smith

Peter Kraft

Multi-ancestry polygenic risk scores for venous thromboembolism

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Venous thromboembolism (VTE) is among the top five most common vascular diseases in most countries [START_REF] Wendelboe | Global Burden of Thrombosis: Epidemiologic Aspects[END_REF]. The estimated lifetime risk of VTE is 8% among US adults [START_REF] Bell | Lifetime Risk of Venous Thromboembolism in Two Cohort Studies[END_REF]. Approximately 20% of individuals die within 1 year of a VTE diagnosis often from the provoking conditions, and complications are common among survivors [START_REF] Søgaard | 30year mortality after venous thromboembolism: a population-based cohort study[END_REF]. Thus, the development of tools that stratify people according to their risk of developing VTE is helpful, which could inform risk-stratified prevention strategies that contribute to reducing the burden of VTE.

Polygenic risk scores (PRS) are useful tools for approximating the cumulative genetic susceptibility to complex traits and diseases. PRSs based on the independent genome-wide significant variants discovered in genome-wide association studies (GWAS) European-ancestry samples [START_REF] Heit | A genome-wide association study of venous thromboembolism identifies risk variants in chromosomes 1q24.2 and 9q[END_REF][START_REF] Tang | A Genome-Wide Association Study for Venous Thromboembolism: The Extended Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium[END_REF][START_REF] Germain | Meta-analysis of 65,734 individuals identifies TSPAN15 and SLC44A2 as two susceptibility loci for venous thromboembolism[END_REF][START_REF] Klarin | Genome-wide association analysis of venous thromboembolism identifies new risk loci and genetic overlap with arterial vascular disease[END_REF][START_REF] Lindström | Genomic and transcriptomic association studies identify 16 novel susceptibility loci for venous thromboembolism[END_REF][START_REF] Thibord | Cross-Ancestry Investigation of Venous Thromboembolism Genomic Predictors[END_REF] have been demonstrated to identify individuals at high risk of VTE [START_REF]Genome-wide association analysis of venous thromboembolism identifies new risk loci and genetic overlap with arterial vascular disease[END_REF][START_REF] Kolin | Prediction of primary venous thromboembolism based on clinical and genetic factors within the U[END_REF]. However, there is limited evidence on whether high-dimensional PRS that are not restricted to genome-wide significant variants can enhance the predictive ability.

In the USA, the incidence of VTE is approximately 65% higher in those who identify as Black Americans than White Americans [START_REF] Silverstein | Trends in the incidence of deep vein thrombosis and pulmonary embolism: a 25-year population-based study[END_REF][START_REF] Zakai | Racial differences in venous thromboembolism[END_REF]. Polygenic risk prediction models for VTE could be particularly important among Black Americans, as a clinical tool to reduce this disparity in VTE risk. (This does not preclude research into structural inequities and social determinants of health, which might inform policy interventions to reduce disparities between Black and White Americans.) However, previously developed VTE PRS have been optimized for Europeanancestry populations, and their utility in other populations is unknown. In particular, we are unaware of any efforts to develop VTE PRS specifically for Black Americans.

We developed ancestry-specific and multi-ancestry PRSs for VTE leveraging large GWAS metaanalyses in European-and African-ancestry samples. We validated these PRSs by estimating relative VTE risks across PRS quintiles in five independent U.S.-based studies. We focus on PRS including common variants (minor allele frequencies above 1%) due to difficulties measuring or imputing low frequency or rare variants from GWAS data or imprecision of estimating rare variant associations. Thus our PRSs complement known low frequency variants (such as rs6205 in F5) or known clinical and behavioral risk factors. Here we concentrate on developing PRSs that perform well in diverse populations. Future work will be needed to (a) develop and evaluate models that combine these PRSs with low-frequency and rare variants and other risk factors and (b) assess the clinical utility of VTE risk models for targeted prevention, screening, or treatment [START_REF] Privé | LDpred2: better, faster, stronger[END_REF][START_REF] Ruan | Improving Polygenic Prediction in Ancestrally Diverse Populations[END_REF].

Results

Study sample

The overall study design is illustrated in Figure 1. Our PRS development consisted of two steps: training ancestry-specific PRS and tuning multi-ancestry PRS. We trained ancestry-specific PRSs using European-and African ancestry GWAS summary statistics from the INVENT consortium and two Bayesian methods (LDPRED2 [START_REF] Privé | LDpred2: better, faster, stronger[END_REF] and PRSCSx [START_REF] Ruan | Improving Polygenic Prediction in Ancestrally Diverse Populations[END_REF]). We then tuned the constructed multi-ancestry PRSs by regressing VTE case-control status on a linear combination of the two ancestry-specific PRSs in two separate tuning samples: one European-ancestry tuning sample (1,329 cases and 1,324 controls) and one African-ancestry tuning sample (238 cases and 3,589 6 controls). The testing data set comprised 6,781 cases and 103,016 controls of European ancestry and 1,385 cases and 12,569 controls of African ancestry from five independent studies. Table S1 presents a brief summary of participating studies and biobanks, including basic information about each study or biobank (location, institute, cohort size, and sample recruiting approach), participants (ancestry and age), and genotypes (genotyping platforms and imputation reference).

Comparing PRS distributions across populations

Four single-ancestry PRSs and four multi-ancestry PRSs for VTE were constructed using LDpred2 and PRSCSx and validated in independent European ancestry and African ancestry individuals:

(i) LDpred2 trained using European-ancestry GWAS summary statistics (LDpred2 EUR ); (ii) LDpred2 trained using African-ancestry summary statistics (LDpred2 AFR ); (iii) PRS-CS trained using European-ancestry summary statistics (PRSCSX EUR ); (iv) PRSCS trained using African ancestry summary statistics (PRSCSX AFR ); and (v) LDpred2 EUR + LDpred2 AFR with weights tuned in an independent European-ancestry tuning sample; (vi) LDpred2 EUR + LDpred2 AFR with weights tuned in and independent African-ancestry tuning sample (LDpred2_combined AFR ); (vii) PRSCSX EUR + PRSCSX AFR with weights tuned in the European-ancestry tuning sample (PRSCSX_combined EUR ); (viii) PRSCSX EUR + PRSCSX AFR with weights tuned in the African-ancestry tuning sample (PRSCSX_combined AFR ). All PRSs had higher means in cases than controls in the test data sets (Table 1). Among the European-ancestry VTE cases, the mean PRS was higher for the PRS tuned in European-ancestry samples than for the PRS tuned in African-ancestry samples. The difference was higher for the ancestry-specific PRS (LDpred2 EUR : 0.39 vs LDpred2 AFR : 0.07, PRSCSX EUR : 0.42 vs PRSCSX AFR : 0.31) than for the multi-ancestry PRS (LDpred2_combined EUR : 7 0.39 vs Dpred2_combined AFR : 0.38, PRSCSX_combined EUR : 0.44 vs PRSCSX_combined AFR : 0.41).

Similarly, among the African-ancestry VTE cases, the mean PRS was higher for the Africanancestry-tuned PRS than for the European-ancestry-tuned PRS, with larger difference for the population-specific PRS (LDpred2 EUR : 0.18 vs Dpred2 AFR : 0.19, PRSCSX EUR : 0.22 vs PRSCSX AFR : 0.28) than the multi-ancestry PRS (LDpred2_combined EUR : 0.19 vs Dpred2_combined AFR : 0.23, PRSCSX_combined EUR : 0.26 vs PRSCSX_combined AFR : 0.30).

Evaluation of PRS and VTE risk across populations

Table 2 shows the estimated OR per SD increase of PRS and AUC for VTE in the test set individuals of European-and African ancestry. For the ancestry-specific PRS, LDpred2 EUR and LDpred2 AFR were constructed using 604,741 SNPs and 1,184,805 SNPs, respectively, and PRSCSX EUR and PRSCSX AFR were constructed using 591,788 SNPs and 586,660 SNPs, respectively.

Multi-ancestry PRS were developed as a linear combination of the ancestry-specific PRS, resulting in 1,212,566 SNPs for LDpred2 and 598,977 SNPs for PRSCSX. The multi-ancestry PRSs outperformed ancestry-specific PRSs in both European-and African-Ancestry test samples and across training methods (LDpred2, PRSCSx) (Figure 2,S.Figure 1). In the European-ancestry test set, multi-ancestry PRS in which the weights were tuned in European ancestry samples performed the best (PRSCSX_combined EUR : AUC=0.61 (0.6, 0.62), OR=1.48 (1.45, 1.52), LDpred2_combined EUR : AUC=0.60 (0.59, 0.61), OR=1.42 (1.39, 1.46)). Similarly, in the Africanancestry test set, a multi-ancestry PRS in which the weights were tuned in African-Ancestry samples performed the best (PRSCSX_combined AFR : AUC=0.59 (0.57, 0.60), OR=1.38 (1.30, 1.45); LDpred2_combined AFR : AUC=0.57 (0.55, 0.58), OR=1.26 (1.20, 1.33)).

8

The association between the PRSs and VTE risk by PRS percentile are shown in Figure 3. The association between the highest fifth percentile of PRSCSX EUR (RR=1.89) and LDpred2 EUR (RR=1.79) and VTE risk was greater than that of genome-wide significant PRS (RR=1.78). The highest fifth percentile of the best-performing PRS (PRSCSX_combined EUR ) was associated with 1.9-fold increased risk for VTE among European ancestry subjects compared to the middle stratum (40-50%). Among the African-ancestry samples, the corresponding risk was about 1.68fold (PRSCSX_combined AFR ), which is smaller than that in European ancestry samples.

Inclusion of known low frequency alleles

When we reconstructed PRS including the five genome-wide significant variants, the new PRS performed worse than our original PRS without the five SNPs in European-(PRSCSX_combined EUR : AUC= 0.57 (0.56, 0.59), LDpred2_combined EUR : AUC= 0.52 (0.50, 0.53)) and in Africanancestry test samples (PRSCSX_combined AFR : AUC= 0.59 (0.58, 0.60), LDpred2_combined AFR : AUC= 0.56 (0.55, 0.57)) (S.Figure 2). This is likely because the five SNPs are rare in one or both populations (average MAF in European ancestry=0.1, African ancestry=0), and our tuning samples are small, resulting in noisy weights. Future studies with larger and more diverse training samples and further tuning steps are needed to learn better multi-ancestry PRS weights.

Discussion

Multi-ancestry PRSs outperformed population specific PRSs in U.S. European-and Africanancestry samples, with a greater improvement in African-ancestry samples. The highest fifth percentile of the best performing multi-ancestry PRS in the European ancestry test samples was associated with an approximately 2-fold increased risk for VTE relative to the middle stratum among European-ancestry subjects. The corresponding risk was smaller (1.7-fold) among the African-ancestry subjects, but still non-negligible and higher than any single-ancestry PRS, highlighting that multi-ancestry PRS may be used to identify individuals at highest risk for VTE events. These data may also be useful in guiding primary prevention and treatment strategies across populations, although we stress that demonstrating PRS discrimination is not sufficient to establish clinical utility, which requires consideration of risks and benefits of specific proposed interventions [START_REF] Privé | LDpred2: better, faster, stronger[END_REF][START_REF] Ruan | Improving Polygenic Prediction in Ancestrally Diverse Populations[END_REF].

To our knowledge, this is the first attempt to develop PRS of VTE specific to African-ancestry populations. Clinical evaluation of PRS is needed in African-ancestry populations, where the burden of VTE is growing due to its increase in VTE incidence. Our PRS, developed and validated in African-ancestry samples, could be a step towards risk-based clinical management of VTE among Black Americans, as a complement to primary prevention efforts. Black Americans and other population groups suffer social disadvantage and lifestyle risk factors that could be a strong contributors to the disparities in VTE [START_REF] Folsom | Reasons for Differences in the Incidence of Venous Thromboembolism in Black Versus White Americans[END_REF]. Encouragingly, healthy lifestyle factors were associated with a lower incidence of VTE among people at high genetic risk for VTE [START_REF] Evans | Lifestyle Moderates Genetic Risk of Venous Thromboembolism: The ARIC Study[END_REF]. Hence, as with most diseases, primary prevention efforts directed at lifestyle interventions to reduce weight or increase activity would have the great potential to reduce the societal burden of VTE.

Further research should determine best approaches to VTE prevention that improve health equity.

10

A recent GWAS meta-analysis demonstrated that European-ancestry individuals at or above the top fifth percentile of a PRS comprised of 37 genome-wide significant variants had a 3.2-fold greater risk for VTE (OR: 3.19; 95% CI: 2.89-3.52) relative to half of the population in the middle of the range [START_REF] Lindström | Genomic and transcriptomic association studies identify 16 novel susceptibility loci for venous thromboembolism[END_REF]. More recently, a PRS using the 100 lead variants identified in a larger European ancestry meta-analysis showed AUC=0.620 (95% CI, 0.616-0.625) [START_REF] Thibord | Cross-Ancestry Investigation of Venous Thromboembolism Genomic Predictors[END_REF]. Since these previous PRS include low MAF variants with large effect sizes (e.g., rs6025: transancestry OR=2.39 (8) on F5 gene), the performance of these previous PRSs and our PRSs is not directly comparable. It is worth noting that our PRS was built using genome-wide common variants and was designed to be transportable between European-and African-ancestry individuals, which can be useful for settings with diverse genetic background. The PRSs presented here complement the lowfrequency, large-effect variants and clinical and behavioral risk factors; future work should develop and evaluate comprehensive risk models combining multi-ancestry PRS, low-frequency variants and other risk factors.

The major strength of the study is that it is the first attempt to develop and validate multiancestry PRS for VTE, providing potential utility of PRS in VTE prevention among Africanancestry populations, where the VTE burden is high. In addition, we validated the PRS in the five independent biobanks from GBMI using harmonized analysis framework (e.g. phenotype definitions, ancestry assignments, and PRS construction).

There are several limitations in our study. First, we have focused on common SNPs, specifically HapMap3 SNPs for VTE PRS construction. As a result, information from rarer variants missing in the LD reference panel may not be captured in other non-European ancestries. Second, the lower predictive ability of VTE PRS in African-ancestry samples can be explained by smaller sample size of African-ancestry VTE meta-analysis GWAS, which is 10 times smaller than European GWAS. Third, there remains a multitude of factors that may contribute to crossbiobank heterogeneity including phenotype precision, cohort-level disease prevalence, and environmental factors. We have provided analysis results by cohort (Supplementary Figure 1).

Conclusions

We found that multi-ancestry PRS for VTE outperformed population-specific PRS, especially in African ancestry populations with relatively small GWAS sample sizes. These findings suggest that the multi-ancestry PRS may be used to identify individuals at highest risk for VTE event and provide guidance for the most effective treatment strategy across populations.

Materials and Methods

Study populations

We trained the PRS using summary statistics from the International Network against Venous Thrombosis (INVENT) consortium cross-ancestry GWAS meta-analyses of European-(71,771 VTE cases and 1,059,740 controls) and African-ancestry samples (7,482 VTE cases and 129,975 controls) [START_REF] Thibord | Cross-Ancestry Investigation of Venous Thromboembolism Genomic Predictors[END_REF]. The meta-analysis is based on prospective cohorts and case-control data from 30 studies.

Tuning [START_REF] Wendelboe | Global Burden of Thrombosis: Epidemiologic Aspects[END_REF]329 European-ancestry populations included (Figure 1). These tuning and validation data were not included in the GWAS used in the training step. The definitions of African-and Europeanancestry populations in each study are provided in the Supplementary Materials; these definitions typically involve both self-reported race and ethnicity and genetic similarity to a set of (study-specific) labeled reference samples.

Supplementary Table 1 summarizes the study design, genotyping arrays, and the sample size in each study. All studies were approved by the relevant institutional ethics committees and review boards, and all participants provided written informed consent.

Statistical methods

PRS training and tuning

PRS training and tuning using LDpred2. We ran LDpred2-auto [START_REF] Privé | LDpred2: better, faster, stronger[END_REF] to construct PRS on HapMap3 variants using the INVENT GWAS meta-analysis summary statistics corresponding to each population. We constructed linkage disequilibrium (LD) reference panels for the development of the European-ancestry PRS (LDpred2 EUR ) and African-Ancestry PRS (LDpred2 AFR ) using the EUR and AFR supersamples from the 1000 Genomes Project (Phase 3), respectively.(18) These population-specific PRSs were then linearly combined to construct multi-ancestry PRS 13 (LDpred2 EUR + LDpred2 AFR ) in which the relative contribution of each PRS was estimated by logistic regression in the tuning dataset of European-ancestry samples (LDpred2_combined EUR ) and African-ancestry samples (LDpred2_combined AFR ). Analyses were run using R; code is available at https://github.com/yonhojee/VTE_PRS.

PRS training and tuning using PRSCSx. We separately applied PRSCSx [START_REF] Ruan | Improving Polygenic Prediction in Ancestrally Diverse Populations[END_REF] to the summary statistics from the European-and African-ancestry INVENT VTE GWAS, using the EUR and AFR LD reference panels from the 1000 Genomes Project (Phase 3). The global shrinkage parameter was learnt from the data using a fully Bayesian approach. Ancestry-specific PRSs generated using European-(PRSCSx EUR ) and African-specific posterior weights (hereafter denoted as PRSCSx AFR ) were linearly combined to construct multi-ancestry PRS (PRSCSx EUR + PRSCSx AFR ). The regression coefficients for the linear combination were obtained by fitting a logistic regression model in the tuning data set of European ancestry samples (PRSCSx_combined EUR ) and African American samples (PRSCSx_combined AFR ). Analyses were run using Python; code is available at https://github.com/yonhojee/VTE_PRS.

PRS validation

In each test dataset, population-specific PRSs were calculated as ܴܲܵ ாோ ൌ ∑ ߚ ݔ and ܴܲܵ ிோ ൌ ∑ ߙ ݔ , where x ik is the dosage of risk allele (0-2) at genetic variant k for subject i, and ߚ and ߙ are the corresponding weight in European and African PRS, respectively. The estimates of ߚ and ߙ were trained using summary statistics from the INVENT consortium and LDpred2 and PRSCSx as described above.

14

We calculated the multi-ancestry PRSs as the linear combination of European-and Africanancestry specific PRS:

ܾ݀݁݊݅݉ܿ_ܴܵܲ ாோ ൌ ߛ ிோ_ாோ ܴܲܵ ிோ ߛ ாோ_ாோ ܴܲܵ ாோ ܾ݀݁݊݅݉ܿ_ܴܵܲ ிோ ൌ ߜ ிோ_ிோ ܴܲܵ ிோ ߜ ாோ_ிோ ܴܲܵ ாோ
where PRS AFR and PRS EUR are the PRSs trained in single-ancestry GWAS and the ߛ and ߜ are "meta-weights" tuned in European-and African-ancestry samples, respectively. SNPs with imputation R 2 > 0.9 in training dataset were retained for subsequent analyses. The lists of SNPs and the weights for the PRS computation are available at https://github.com/yonhojee/VTE_PRS. PRSs were standardized within each validation sample to have unit SD in the control subjects.

Logistic regression, adjusting for ten principal components and sex, was used to estimate odds ratios (ORs) for association between the standardized PRSs and VTE risk in each testing set. The discrimination of PRS was assessed using area under the receiver operating curve (AUC). The OR per SD and AUC were obtained individually for each study and combined separately for European-and African-ancestry samples using a fixed-effect meta-analysis.

All statistical analyses were conducted using R v.4.3.0. Logistic regression and AUC were done using glm() and roc() in R.

The distribution of relative risk of VTE by PRS across populations.

We simulated 100,000 individuals with PRS distribution of N(0,1) multiplied by log OR per SD estimates for each PRS. The simulated PRS was then exponentiated to estimate relative risk estimates and split into the percentile categories: [0-1%] (1-5%], (5-10%], (10-20%], (20-30%],

(30-40%], (40-50%] (reference group), (50-60%], (60-70%], (70-80%], (80-90%], (90-95%],

(95-99%] and (99-100%].

Sensitivity analysis of including known low frequency alleles

Out of the 37 genome-wide significant variants, our current PRSs do not include five variants (rs6025, rs145470028, rs1799963, rs6048, and rs143478537), which would have been filtered out of our analyses for one reason or another (e.g., on the X chromosome, low minor allele frequency [MAF]). These variants are important to be considered in VTE PRS given their large effect sizes (e.g., rs6025: transancestry OR=2.39(8) on F5 gene). As a sensitivity analysis, we constructed new PRSs, which additionally include these previously reported variants that are i) not included in our PRS due to the low frequency and ii) not in LD with the variants already included in our PRS. The final PRSs were obtained by the linear combination of the original PRS (constructed using common variants only) and the additional SNPs where the coefficients for the original PRS and the additional SNPs were tuned in the independent ancestry-specific samples (See more details in the Supplementary Materials). Legends to Tables Table 1. Mean and standard deviation of standardized polygenic risk scores with VTE risk in the test set individuals of European and African ancestry. a Combined PRSs were generated using the formula α0 + α1PRS1 + α2PRS2 where α0, α1 and α2 are the weights obtained by fitting a logistic regression model with VTE as outcome, PRS1 and PRS2 as explanatory variables using the validation data set. The weights for the considered combination of PRSs can be found at https://github.com/yonhojee/VTE_PRS.

Legends to Figures

  cases and 1,324 controls of European-ancestry and 238 cases and 3,589 controls of African-ancestry) and validation data (6,781 cases and 103,016 controls of European ancestry and 1,385 cases and 12,569 controls of African ancestry) came from Nurses' Health Study [NHS] and Health Professional Follow-up Study [HPFS] and 4 Global Biobank Meta-analysis Initiative (GBMI) biobanks (Michigan Genomics Initiative [MGI], UCLA Precision Health Biobank [UCLA], Penn Medicine Biobank [PMBB], and Lifelines) with representation across African and

Figure 1 .

 1 Figure 1. Overview of development and validation of population-specific and multi-ancestry PRS for venous thromboembolism.

Figure 2 .

 2 Figure 2. AUC and OR for population-specific and multiancestry PRS across populations.

Figure 3 .

 3 Figure 3. Distribution of relative risk of VTE by PRS across populations.

Figure 1 .Figure 2 .

 12 Figure 1. Overview of development and validation of population-specific and multi-ancestry PRS for venous thromboembolism.

Figure 3 .

 3 Figure 3. Distribution of relative risk of VTE by PRS across populations.

Table 2 .

 2 Association of polygenic risk scores and VTE risk in the test set individuals of European and African ancestry.

	21

Table 1 .

 1 Mean and standard deviation of standardized polygenic risk scores with VTE risk in the test set individuals of European and African ancestry.

	European

Table 2 .

 2 Association of polygenic risk scores and VTE risk in the test set individuals of European and African ancestry.

	PRS testing population
	European	African

Acknowledgements

The views expressed in this manuscript are those of the authors and do not necessarily represent the views of the National Cancer Institute; National Heart, Lung, and Blood Institute; the National Institutes of Health; or the US Department of Health and Human Services. This work was supported by the National Institutes of Health grant 5U01CA261339.

The INVENT Consortium acknowledges all of the participants across the studies who provided their health information to support these analyses.

The authors would like to thank the participants and staff of the NHS and NHSII for their valuable contributions. The NHS, NHS-II and HPFS were supported by grants UM1 CA186107, U01 CA176726, and U01 CA167552 from the National Institutes of Health.

The authors acknowledge the Michigan Genomics Initiative participants, Precision Health at the University of Michigan, the University of Michigan Medical School Central Biorepository, and the University of Michigan Advanced Genomics Core for providing data and specimen storage, management, processing, and distribution services, and the Center for Statistical Genetics in the Department of Biostatistics at the School of Public Health for genotype data curation, imputation, and management in support of the research reported in this publication.

We gratefully acknowledge the resources provided by the Institute for Precision Health (IPH) and participating UCLA ATLAS Community Health Initiative patients. The UCLA ATLAS Community Health Initiative in collaboration with UCLA ATLAS Precision Health Biobank, is a program of IPH, which directs and supports the biobanking and genotyping of biospecimen samples from participating UCLA patients in collaboration with the David Geffen School of Medicine, UCLA CTSI and UCLA Health. Members of the UCLA ATLAS Community Health Initiative include Ruth Johnson, Yi Ding, Vidhya Venkateswaran, Arjun Bhattacharya, Alec Chiu, Tommer Schwarz, Malika Freund, Lingyu Zhan, Kathryn S. Burch, Christa Caggiano, Brian Hill, Nadav Rakocz, Brunilda Balliu, Jae Hoon Sul, Noah Zaitlen, Valerie A. Arboleda, Eran Halperin, Sriram Sankararaman, Manish J. Butte, Clara Lajonchere, Daniel H. Geschwind, and Bogdan Pasaniuc, on behalf of the UCLA Precision Health Data Discovery Repository Working Group and UCLA Precision Health ATLAS Working Group.

We acknowledge the PMBB for providing data and thank the patient-participants of Penn Medicine who consented to participate in this research program. We would also like to thank the Penn Medicine BioBank team and Regeneron Genetics Center for providing genetic variant data for analysis. The PMBB is approved under IRB protocol# 813913 and supported by Perelman School of Medicine at University of Pennsylvania, a gift from the Smilow family, and the National Center for Advancing Translational Sciences of the National Institutes of Health under CTSA award number UL1TR001878. The Lifelines Biobank initiative has been made possible by funding from the Dutch Ministry of Health, Welfare and Sport, the Dutch Ministry of Economic Affairs, the University Medical Center Groningen (UMCG the Netherlands), University of Groningen and the Northern Provinces of the Netherlands. The generation and management of GWAS genotype data for the Lifelines Cohort Study is supported by the UMCG Genetics Lifelines Initiative (UGLI). UGLI is partly supported by a Spinoza Grant from NWO, awarded to Cisca Wijmenga. The authors wish to acknowledge the services of the Lifelines Cohort Study, the contributing research centers delivering data to Lifelines, and all the study participants. UGLI Lifelines Cohort Study group author: Raul Aguirre-Gamboa (1), Patrick Deelen (1), Lude Franke (1), Jan A Kuivenhoven (2), Esteban A Lopera Maya (1), Ilja M Nolte (3), Serena Sanna (1), Harold Snieder (3), Morris A Swertz (1), Peter M. Visscher (3,4), Judith M Vonk (3), Cisca Wijmenga (1), Naomi Wray (4); (1) Department of Genetics, University of Groningen, University Medical Center Groningen, The Netherlands; (2) Department of Pediatrics, University of Groningen, University Medical Center Groningen, The Netherlands; (3) Department of Epidemiology, University of Groningen, University Medical Center Groningen, The Netherlands; (4) Institute for Molecular Bioscience, The University of Queensland, Brisbane, Queensland, Australia.

Conflict of Interest Statement


Philip Barron Abitua

T Blair Gainous

Angela N Kaczmarczyk

Christopher J Winchell

Clare Hudson

Kaori Kamata

Masashi Nakagawa

Motoyuki Tsuda

Takehiro G Kusakabe

Michael Levine

The pre-vertebrate origins of neurogenic placodes

The sudden appearance of the neural crest and neurogenic placodes in early branching vertebrates has puzzled biologists for over a century 1 . These embryonic tissues contribute to the development of the cranium and associated sensory organs, which were crucial for the evolution of the vertebrate "newhead" 2,3 . A previous study suggests that rudimentary neural crest cells existed in ancestral chordates 4 . However, the evolutionary origins of neurogenic placodes have remained obscure owing to a paucity of embryonic data from tunicates, the closest living relatives to those early vertebrates 5 . Here we show that the tunicate Ciona intestinalis exhibits a proto-placodal ectoderm (PPE) that requires inhibition of bone morphogenetic protein (BMP) and expresses the key regulatory determinant Six1/2 and its co-factor Eya, a developmental process conserved across vertebrates. The Ciona PPE is shown to produce ciliated neurons that express genes for gonadotropin-releasing hormone (GnRH), a G-protein-coupled receptor for relaxin-3 (RXFP3) and a functional cyclic nucleotide-gated channel (CNGA), which suggests dual chemosensory and neurosecretory activities. These observations provide evidence that Ciona has a neurogenic proto-placode, which forms neurons that appear to be related to those derived from the olfactory placode and hypothalamic neurons of vertebrates. We discuss the possibility that the PPE-derived

GnRH neurons of Ciona resemble an ancestral cell type, a progenitor to the complex neuronal circuit that integrates sensory information and neuroendocrine functions in vertebrates.

Neurogenic placodes contribute to cranial sensory systems that mediate hearing, smell and taste [START_REF] Steventon | Neural crest and placode interaction during the development of the cranial sensory system[END_REF] . They are considered a vertebrate innovation [START_REF] Shimeld | Vertebrate innovations[END_REF] and the full repertoire, including olfactory, otic, epibranchial and trigeminal placodes, are already present in jawless hagfish and lampreys. It has been argued that the evolution of neurogenic placodes was a crucial event for facilitating the transition from filter-feeding invertebrate chordates to predatory vertebrates [START_REF] Gans | Neural crest and the origin of vertebrates: a new head[END_REF] .

Comparative embryological studies have failed to identify any clear evolutionary precursors of neurogenic placodes in invertebrate chordates. It has been suggested that the pre-oral organ of amphioxus, Hatschek's pit, is related to the adenohypophyseal (anterior pituitary) placode of vertebrates; however, the fate map of this tissue remains unresolved [START_REF] Patthey | The evolutionary history of vertebrate cranial placodes -I: cell type evolution[END_REF] . In addition, the atrial primordia of Ciona have been likened to the otic and lateral-line placodes of fish, but they are not specified until tailbud stages and do not derive from the anterior neural plate border [START_REF] Mazet | Molecular evidence from Ciona intestinalis for the evolutionary origin of vertebrate sensory placodes[END_REF] . Also, previous studies provided evidence for a placodal-like tissue in Ciona, but it was homologized to the adenohypophyseal placode [START_REF] Mazet | Molecular evidence from Ciona intestinalis for the evolutionary origin of vertebrate sensory placodes[END_REF][START_REF] Christiaen | A modular cis-regulatory system controls isoform-specific pitx expression in ascidian stomodaeum[END_REF] , which lacks neurogenic potential. We show that the latter tissue does in fact produce ciliated sensory neurons that express GnRH and RXFP3, signatures of neuroendocrine cell types found in vertebrates [START_REF] Sutton | Distribution of G-protein-coupled receptor (GPCR) 135 binding sites and receptor mRNA in the rat brain suggests a role for relaxin-3 in neuroendocrine and sensory processing[END_REF][START_REF] Schwanzel-Fukuda | Origin of luteinizing hormone-releasing hormone neurons[END_REF][START_REF] Wray | Evidence that cells expressing luteinizing hormone-releasing hormone mRNA in the mouse are derived from progenitor cells in the olfactory placode[END_REF] .

Neurogenic placodes develop from a specific domain anterior to the neural plate of vertebrate embryos [START_REF] Steventon | Neural crest and placode interaction during the development of the cranial sensory system[END_REF] . All placode subtypes arise from this U-shaped domain, termed the preplacodal ectoderm [START_REF] Kwon | Identification of early requirements for preplacodal ectoderm and sensory organ development[END_REF] . This tissue expresses the transcription factor Six1, and mutations in Six1 cause craniofacial defects in humans, including deafness [START_REF] Ruf | SIX1 mutations cause branchio-oto-renal syndrome by disruption of EYA1-SIX1-DNA complexes[END_REF] . We found that the Ciona homologue of Six1, Six1/2, as well as its co-factor Eya, is initially expressed in a row of eight cells located immediately anterior to the neural plate, as marked by ZicL (Fig. 1a and Extended Data Fig. 1). This pattern is recapitulated by a Six1/2 regulatory sequence driving expression of an mCherry reporter gene (Six1/2. mCherry; Fig. 1b and Supplementary Video 1). At tailbud stages the Six1/2 + cells intercalate and become situated above the anterior brain (that is, the sensory vesicle), surrounding the future oral opening (Extended Data Fig. 2).

In vertebrates, the preplacodal ectoderm is specified by a two-step process during which a gradient of BMP signalling establishes preplacodal competence in the non-neural ectoderm during the blastula/early gastrula stage, followed by subsequent inhibition of BMP signalling during gastrulation [START_REF] Kwon | Identification of early requirements for preplacodal ectoderm and sensory organ development[END_REF] . BMP signalling may not be required for establishing proto-placodal competence in Ciona, but inhibition of BMP signalling seems to be crucial for its specification. We observe conserved expression of BMP2/4 and the Chordin (also known as Ci-Chordin) antagonist in ventral and dorsal regions, respectively, at the time when the PPE is specified during gastrulation and neurulation (Extended Data Fig. 3).

Localized Chordin expression at the lateral edges of the Ciona PPE (Fig. 1b) suggests a role in antagonizing BMPs emanating from ventral regions, similar to the situation seen in vertebrates [START_REF] Reichert | A BMP regulatory network controls ectoderm al cell fate decisions at the neural plate border[END_REF] . We therefore explored the possibility that inhibition of BMP signalling might be required for the activation or maintenance of Six1/2 expression. We sought to overcome the presumed inhibitory effect of chordin by the targeted misexpression of BMP ligands or constitutively active forms of BMP receptors. These experiments were performed using a previously described Dmrt regulatory DNA (that is, Dmrt. GFP) [START_REF] Wagner | FGF signaling establishes the anterior border of the Ciona neural tube[END_REF] , which drives expression throughout the anterior neural plate and adjacent anterior ectoderm (Extended Data Fig. 1).

The co-electroporation of Dmrt. BMP2/4 with Dmrt. GFP and Six1/2. mCherry resulted in loss of Six1/2 expression in the anterior neural plate border, consistent with similar experiments in vertebrates [START_REF] Ahrens | Tissues and signals involved in the induction of placodal Six1 expression in Xenopus laevis[END_REF] (Fig. 1c). Conversely, misexpression of BMP5/7 had no noticeable effect on expression of the Six1/2 reporter or embryo morphology (Fig. 1d). These observations suggest that chordin may be specifically required to exclude BMP2/4 signals from the PPE.

Additional misexpression assays were performed with a constitutively active form of the BMP receptor (BMPR1 CA ) for cell-autonomous activation of BMP signalling. In these experiments, activation is restricted to Dmrt + cells, which form the anterior neural plate and adjacent anterior ectoderm. Ectopic BMPR1 activity results in reduced expression of Six1/2 (Fig. 1e), similar to the results obtained upon misexpression of the BMP2/4 ligand. However, the phenotype is less severe, with persistence of Six1/2 expression in the medially derived regions of the PPE (compare Fig. 1d with Fig. 1e). Expression is lost in lateral regions, which are the source of the GnRH neurons.

To demonstrate that disrupting PPE formation is a specific result of ectopic BMP signalling, we misexpressed a constitutively active form of the transforming growth factor-β receptor (TGF-βR CA ) throughout the anterior brain and adjacent ectoderm (Dmrt. TGF-βR CA ). This had very little effect on the expression of Six1/2, despite severe defects in brain morphology (Fig. 1f). We therefore conclude that inhibition of BMP signalling is essential for the specification of the Ciona PPE, as observed in vertebrates [START_REF] Reichert | A BMP regulatory network controls ectoderm al cell fate decisions at the neural plate border[END_REF][START_REF] Ahrens | Tissues and signals involved in the induction of placodal Six1 expression in Xenopus laevis[END_REF] . However, it is possible that this inhibition indirectly regulates Six1/2 expression.

The olfactory placode is the anterior-most neurogenic placode in vertebrates and develops above the oral opening [START_REF] Patthey | The evolutionary history of vertebrate cranial placodes -I: cell type evolution[END_REF] . GnRH neurons arising from this placode migrate along the axonal tracts of olfactory neurons towards the hypothalamus [START_REF] Schwanzel-Fukuda | Origin of luteinizing hormone-releasing hormone neurons[END_REF][START_REF] Wray | Evidence that cells expressing luteinizing hormone-releasing hormone mRNA in the mouse are derived from progenitor cells in the olfactory placode[END_REF] . They subsequently innervate the anterior pituitary and stimulate the release of gonadotropins, hormones critical for growth, sexual development and reproduction [START_REF] Wray | From nose to brain: development of gonadotrophin-releasing hormone-1 neurones[END_REF] . Previous studies have documented the expression of six different GnRH peptides encoded by two duplicated genes in Ciona [START_REF] Kusakabe | A conserved non-reproductive GnRH system in chordates[END_REF] . GnRH2 (hereafter referred to as GnRH) is expressed near the oral opening (Extended Data Fig. 4) and we used the 59 regulatory sequences (GnRH. GFP) to label this lineage during development. During the tailbud stage, the posterior-most Six1/2 + cells acquire a distinctive cone-shaped morphology, permitting their direct visualization during morphogenesis (Fig. 2a and Extended Data Fig. 2). By the larval stage, these cells have differentiated into neurons that co-express GnRH. GFP and Six1/2. mCherry (Fig. 2b and Supplementary Video 2). The production of these neurons is strongly inhibited by BMP signalling (Extended Data Fig. 5).

Ciona PPE-derived GnRH neurons express two genes encoding G-protein-coupled receptors (GPCRs) implicated in neurosecretion in vertebrates, a somatostatin/opioid/galanin/ chemokine-like receptor (SOG/Chemokine-like) and a relaxin-3 receptor (RXFP3) (Fig. 2c,d, Extended Data Fig. 4 and Supplementary Video 2). The relationship of these Ciona GPCRs to human rhodopsin-class GPCRs is shown in a summary tree (Fig. 2e and Extended Data Figs 6 and7). In vertebrates, RXFP3 is expressed in hypothalamic paraventricular nuclei that provide synaptic inputs for the secretion of GnRH within the hypothalamicpituitary-gonadal (HPG) axis [START_REF] Mcgowan | Relaxin-3 stimulates the hypothalamic-pituitary-gonadal axis[END_REF] . Notably, the PPE-derived Ciona GnRH neurons survive through metamorphosis despite extensive cell death and reorganization of most of the larval nervous system (Fig. 2f). Moreover, injection of GnRH into mature Ciona causes the release of gametes and is likely to trigger spawning in wild populations [START_REF] Terakado | Induction of gamete release by gonadotropin-releasing hormone in a protochordate, Ciona intestinalis[END_REF] .

As both GnRH1 and chemosensory neurons arise from the olfactory placode of vertebrates, we were curious whether Ciona PPE-derived neurons might also possess chemosensory properties. Chemosensation by olfactory neurons occurs via odorant receptors in the tips of their cilia. Swimming larvae possess a group of neurons located behind the oral opening, previously described as the anterior apical trunk epidermal neurons (aATENs). These neurons possess cilia that are clearly visualized using an acetylated α-tubulin antibody (Fig. 3a). Electroporation with GnRH. mCherry followed by acetylated α-tubulin staining confirmed that the cilia project from the PPE-derived GnRH neurons (Fig. 3b). Moreover, comprehensive visualization of Six1/2 + cell morphogenesis reveals that the aATEN-GnRH neurons are derived from the a11.205 lineage (Extended Data Fig. 2). Ultrastructural analysis has shown that cilia of these neurons have a 9 + 2 microtubule arrangement [START_REF] Konno | Distribution and structural diversity of cilia in tadpole larvae of the ascidian Ciona intestinalis[END_REF] , a shared characteristic attributed to non-motile cilia of vertebrate olfactory neurons [START_REF] Reese | Olfactory cilia in the frog[END_REF] .

Conventional odorant receptors have not been identified in the Ciona genome [START_REF] Patthey | The evolutionary history of vertebrate cranial placodes -I: cell type evolution[END_REF] , but it is unlikely that a marine chordate would lack a means to detect odorants. Instead of examining divergent GPCRs, we focused on downstream effectors of chemosensory transduction. Cyclic nucleotide-gated (CNG) cation channels are critical for mediating sensory perception in photoreceptors and olfactory neurons. The binding of odorants to receptors causes an intracellular increase in cyclic nucleotides that, in turn, results in an influx of calcium, which ultimately leads to the depolarization of sensory neurons [START_REF] Kaupp | Olfactory signalling in vertebrates and insects: differences and commonalities[END_REF] . In rats, CNGA2 expression is detected in GnRH1 neurons but is not necessary for their activity. This expression might be due to prior activation in the shared progenitor of GnRH1 and olfactory neurons [START_REF] Constantin | Gonadotropin-releasing hormone-1 neuronal activity is independent of cyclic nucleotide-gated channels[END_REF][START_REF] El-Majdoubi | Localization of olfactory cyclic nucleotide-gated channels in rat gonadotropin-releasing hormone neurons[END_REF] .

A homologue of vertebrate CNGA1 is expressed in fully differentiated aATEN-GnRH neurons of Ciona (Fig. 3c). To determine whether Ciona CNGA forms a functional channel, we transfected a CNGA expression construct into human embryonic kidney 293 cells incubated with the intracellular calcium indicator Fura-2 AM to test for a response to cyclic nucleotides. Indeed, application of either 8-bromoguanine-39,59-cyclic monophosphate (8-Br-cGMP) or 8-bromoadenosine-39,59-cyclicmonophosphate (8-Br-cAMP) elicited an increase in intracellular calcium concentration, with a higher sensitivity to cAMP than cGMP (Fig. 3d and Extended Data Fig. 8). The occurrence of a functional CNGA in PPEderived neurons, in addition to the presence of a 9 + 2 cilium, suggests that they are likely to possess chemosensory activities. These observations provide the first evidence, to our knowledge, that putative sensory neurons arise from a placodal-like territory in a nonvertebrate.

Our data suggest that the aATEN-GnRH neurons have dual neurosecretory and chemosensory properties. In vertebrates, neurosecretory and chemosensory cells both arise from the olfactory placode and are intimately linked. GnRH neuroblasts use the axon tracts of chemosensory neurons to guide them to their final destination in the hypothalamus. In adults, pheromones detected by chemosensory neurons cause a release of gonadotropins via hypothalamic GnRH neurons, a function crucial to sexual reproduction. Furthermore, transneuronal tracing experiments have shown that olfactory and GnRH neurons form a coherent neuronal circuit [START_REF] Boehm | Feedback loops link odor and pheromone signaling with reproduction[END_REF][START_REF] Yoon | Olfactory inputs to hypothalamic neurons controlling reproduction and fertility[END_REF] , a connection that likely arose long ago.

It has been proposed that neuronal circuits evolve from the functional segregation of specialized cell types from multifunctional ancestral cells, akin to the well-known mechanism of gene duplication and subfunctionalization. A chemosensory-neurosecretory precursor is likely to have existed in ancestral chordates before the diversification of the HPG axis seen in vertebrates [START_REF] Gorbman | Olfactory origins and evolution of the brain-pituitary endocrine system: facts and speculation[END_REF][START_REF] Arendt | The evolution of cell types in animals: emerging principles from molecular studies[END_REF] . We have presented evidence that the PPE-derived aATEN-GnRH neurons in Ciona provide a particularly vivid example of such an ancestral cell type. Over the course of evolution within the vertebrate lineage, the proposed chemosensory and neurosecretory functions of these neurons may have become segregated into dedicated cell types that work together within a coherent circuit (Fig. 3e). Such cellular subfunctionalization might be a generally important mechanism of neuronal circuit evolution in vertebrates.

Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

METHODS

Embryo preparation and imaging. Ciona intestinalis adults were obtained, fertilized in vitro and electroporated for transient transgenesis as previously described [START_REF] Wagner | FGF signaling establishes the anterior border of the Ciona neural tube[END_REF] . For each electroporation, typically 70 µg of DNA was resuspended in 100 µl buffer. Embryos were fixed at the appropriate developmental stage for 15 min in 4% formaldehyde. The tissue was then cleared in a series of washes of 0.01% Triton-X in PBS. Actin was stained overnight with Alexa-647-conjugated phalloidin (Fisher, A22287) at a dilution of 1:500.

For the acetyl-α-tubulin staining, larvae were fixed for 10 min in 2% para-formaldehyde, then washed three times in 0.01% Triton-Xin PBS for 10 min each, and blocked with 1% normal goat serum for 10 min. Afterwards samples were incubated with a monoclonal antiacetylated tubulin antibody (Sigma-Aldrich, T6793) at a dilution of 1:1,000 overnight at 4°C. Next, the tissue was washed in 0.01% Triton-X in PBS, three times for 30 min each. The samples were then incubated for 2 h at room temperature with an Alexa Fluor 488 donkey anti-mouse IgG (H + L) antibody (Fisher, A21202) at a dilution of 1:500 in 1% normal goat serum, and finally washed in 0.01% Triton-Xin PBS (three times for 30 min each).

Samples were mounted in 50% glycerol in PBS, with 2% 1,4-diazabicyclo[2,2,2]octane (Sigma) for microscopy. Confocal images were acquired on a Zeiss LSM 700 microscope using a Plan-Apochromat 20×, 40× or 63× objective. Confocal stacks contained approximately 50 optical slices at a thickness of 1-2 µm each. Images were rendered using Volocity 6 with the 3D opacity, extended focus or XY plane visualization tools. Time-lapse images were taken on a Zeiss LSM 700 microscope at intervals of 3-4 min. Bright field in situ micrographs of larvae were acquired on a Zeiss Axio Imager microscope using an EC Plan-Neofluar 40× objective with a ProgRes C14 plus camera.

Experimental perturbations were performed in triplicate and two biological replicates were scored. No statistical methods were used to predetermine sample size. The scored replicates were totalled and these are indicated in the appropriate figure legends. Phenotypes were scored in animals exhibiting transgenesis (that is, using a reporter unaffected by a given perturbation such as Dmrt. GFP in Fig. 1). Unfertilized animals or those that failed to gastrulate were excluded from the analysis. No randomization or blinding was performed on samples.

Molecular cloning

The Genome Browser Gateway (University of California, Santa Cruz) facilitated the identification of conserved non-coding sequences between Ciona intestinalis and Ciona savignyi. Putative enhancer sequences were amplified by polymerase chain reaction (PCR) and cloned into a pCESA vector using AscI/NotI restriction sites for Chordin (KH.C6.145; 59-CGGTTACGTTAAGTTCGTGCG-39 and 59-CTTTGTTCCTGCTCGAAGTGAG-39), Six1/2 (KH.C3.553; 59-GTGACAGGAAACGTCTAGG-39 and 59-TGGCTCTGAGGCCGTGATTGTAG-39), GnRH (KH.C9.484; 59-GTGTACGATTGACAAGTTGG-39 and 59-TGTTACGTTATCTCTCTAGAAG-39), SOG/ Chemokine-like (ci0100133186; 59-GGCAGATTTGACCTCAACTTG-39 and 59-GGCGTTTCC GAAAAGCCCTTT-39), RXFP3 (KH.C6.184; 59-CCGTGATTGTAATACGTATGAC-39 and 59-GCACAGTATTGTGTTATATACC-39), and CNGA (KH.C2.249; 59-GTGGAGATGCCGATTTCAACC-39 and 59-CACTTTGTTTGGCATTATTCC-39). The Dmrt and ZicL enhancers have been described previously [START_REF] Wagner | FGF signaling establishes the anterior border of the Ciona neural tube[END_REF] . A similar cloning method was used to create misexpression vectors using NotI/ EcoRI sites for the BMP2/4 (KH.C4.125; 59-ATGGTGGCGCTTACGGATTGGAC-39 and 59-CTATCTACACCCACAAGCTTGC-39) and BMP5/7 (KH.C2.336; 59-ATGACTTGTCATCGCAATAAAG-39 and 59-TCAGTTGCACCCGCATGACAC-39) plasmids. The TGF-βR CA construct was made by site-directed mutagenesis to introduce a Q154D missense mutation in the regulatory glycine-and serine-rich domain of the wild type receptor. A similar strategy was used to create the previously described BMPR1 CA (ref. 31). The coding sequence of SOG/Chemokine-like was amplified from cDNA (59-ATGTTGCTCGGAATCATGAAATC-39 and 59-GTATAAATATTCATACTTGTTTCTG-39) and cloned into a pCESA vector using NotI/EcoRI restriction sites.

For the calcium imaging experiments, cDNA fragments containing the full-length coding regions of C. intestinalis CNGA, CNGB and CNGC (Ghost Gene IDs: KH.C2.249,KH.L42.6, and KH.C7.605, respectively) were amplified from a cDNA pool of mid tailbud embryos (CNGA) or larvae (CNGB and CNGC) by PCR using a thermostable DNA polymerase exhibiting proofreading activity (Takara LATaq; Takara Bio) with gene-specific primers (59-AGCTAACACAGATTTTAGTGTAATG-39 and 59-GTTGCGGGTAAATTACATGTC-39 for CNGA; 59-ATGGCATTGCACATAAATGCAA-39 and 59-AGCAAAGACTTTGGTAAACATCAG-39 for CNGB; 59-TAAGGTTGATACAGGTTTCATTGG-39 and 59-GACGTAATCATGACGAAACTGTG-39 for CNGC). The PCR products were subcloned into pcDNA6 plasmid vectors and sequenced on both strands using the cycle sequencing method with an Applied Biosystems 3100 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA).

In situ hybridization

The fluorescent in situ hybridizations were performed as described previously [START_REF] Wagner | FGF signaling establishes the anterior border of the Ciona neural tube[END_REF] . The colorimetric in situ hybridizations on hatched larvae were performed as previously detailed [START_REF] Rehm | Fixation and dissection of Parhyale hawaiensis embryos[END_REF] . The in situ hybridizations for Extended Data Fig. 3 were performed as specified [START_REF] Hudson | Patterning across the ascidian neural plate by lateral Nodal signalling sources[END_REF] . mRNA probes were synthesized using linearized cDNA clones for Six1/2 (ciad46p13), Eya (cilv29d14), Chordin (ciclo16e09), BMP2/4 (cicl060n01), GnRH2 (cilv25p08) and RXFP3 (cinc033n18) from Nori Satoh's cDNA gene collection (OIST, Okinawa, Japan). Probe templates for CNGA (59-ATGGCATTGCACATAAATGCAA-39 and 59-AGCAAAGACTTTGGTAAACATCAG-39) and SOG/Chemokine-like (59-ATGTTGCTCGGAATCATGAAATC-39 and 59-GTATAAATATTCATACTTGTTTCTG-39) were amplified from a cDNA pool of mid tailbud embryos and subcloned into a pBSKM vector.

Phylogenetic analysis

For the initial survey tree, BLASTP [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] was used to assign the two Ciona GPCR sequences of interest to the rhodopsin class. We used the MAFFT E-INS-i algorithm [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF] to align their seven-transmembrane domains with those of the proteins constituting the Pfam seed alignment for rhodopsin GPCRs [START_REF] Finn | Pfam: the protein families database[END_REF] , as well as most Ciona rhodopsin-receptor-like proteins identified previously [START_REF] Kamesh | The repertoire of G protein-coupled receptors in the sea squirt Ciona intestinalis[END_REF] . After removing ambiguously aligned amino acids with Gblocks [START_REF] Talavera | Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments[END_REF] , the 118-sequence data set (including five outgroup sequences for Ciona Glutamate, Smoothened, Frizzled and cAMP receptor-like proteins) contained 373 positions (including gaps). For optimal tree searches and measurement of nodal support, minimum evolution (ME), maximum likelihood (ML) and Bayesian inference (BI) criteria were used [START_REF] Winchell | Expression of the Lhx genes apterous and lim1 in an errant polychaete: implications for bilaterian appendage evolution, neural development, and muscle diversification[END_REF] .

To classify the two Ciona GPCRs of interest with greater accuracy and support, a second set of analyses was performed using a more focused sequence sampling. Based on the initial survey tree and the identifications made by Kamesh et al. [START_REF] Kamesh | The repertoire of G protein-coupled receptors in the sea squirt Ciona intestinalis[END_REF] , we selected presumed members of the specific clade containing the two Ciona GPCRs of interest, its large 'Chemokine Cluster' sister clade, and its two closest rootward sister clades (that is, members of the 'SOG' rhodopsin subclass). For these selections, GPCRDB Tools 40 was used to build an alignment of human seven-transmembrane domain sequences from the Swiss-Prot database. Using Clustal Omega [START_REF] Sievers | Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega[END_REF] , we integrated seven relevant Ciona sequences described by Kamesh et al. [START_REF] Kamesh | The repertoire of G protein-coupled receptors in the sea squirt Ciona intestinalis[END_REF] into the GPCRDB Tools alignment. After hand-editing in Jalview [START_REF] Waterhouse | Jalview Version 2-a multiple sequence alignment editor and analysis workbench[END_REF] , this 75-sequence data set contained 208 positions. As stated earlier, ME, ML and BI were used to analyse the final alignment.

Calcium imaging

HEK 293 cells, a gift from H. Yagisawa (Graduate School of Life Science, University of Hyogo), were grown at 5% CO 2 and 100% relative humidity in DMEM medium supplemented with 10% (v/v) fetal bovine serum and 0.1 mM non-essential amino acids. HEK 293 cells do not express endogenous CNG channels [START_REF] Dhallan | Primary structure and functional expression of a cyclic nucleotide-activated channel from olfactory neurons[END_REF] . Each Ciona CNG construct was transfected into HEK 293 cells using Targefect F-1 reagent (Targeting Systems, Santee, CA, USA) according to the manufacturer's instructions, and the cells were incubated for 48 h.

To quantitatively visualize cyclic nucleotide-induced intracellular calcium ions, transfected HEK 293 cells were incubated for 30 min in a solution containing 140 mM NaCl, 5.4 mM KCl, 5 mM CaCl 2 , 5.6 mM glucose, 10 mM HEPES (pH 7.4) and 2 µM Fura-2 AM (Dojindo Laboratories, Kumamoto, Japan). After washing with PBS, cells were placed in a solution containing 140 mM NaCl, 5.4 mM KCl, 5 mM CaCl 2 , 5.6 mM glucose, and 10 mM HEPES (pH 7.4). 8-Br-cAMP (Sigma B7880) or 8-Br-cGMP (Sigma B1381) at various concentrations was administrated to the cells in a perfusion system. Fluorescence signals were monitored using an Argus-50/CA system (Hamamatsu photonics, Shizuoka, Japan). Each experiment was performed in at least two biological replicates. Authentication and contamination tests were not performed; however, this should not affect the reliability of the results as the untransfected negative control did not show any response to cyclic nucleotides.

Extended Data

Extended Data Figure 1.

Lineage information for Six1/2 expression in Ciona intestinalis from the gastrula to initial tailbud stage. a, Schematic of the anterior neural plate border at the mid-gastrula stage, including cell lineage nomenclature. Dmrt is initially activated in six blastomeres of 64-cell embryos (a7.9, a7.10 and a7.13). This lineage produces only the anterior neural plate (and the adjacent anterior neural plate border), which forms the PPE. The anterior-most ZicL + counterstained with phalloidin (blue). a, At the initial tailbud stage I (according to Hotta et al. [START_REF] Hotta | A web-based interactive developmental table for the ascidian Ciona intestinalis, including 3Dreal-image embryo reconstructions: I From fertilized egg to hatching larva[END_REF] ), the Six1/2 + cells are arranged in a U shape, anterior to the neural plate. a9-c9, Underlined cell lineages are derived from the left side of the embryo. a9, Schematic indicates that the U shape is composed of 16 cells in total at the 11th generation (that is, a11.154, a11.138, etc.). There are no further divisions of these cells until after the late tailbud I stage Extended Data Figure 6.

Initial phylogenetic analysis of two GPCRs expressed in the PPE-derived neurons of Ciona intestinalis. This broad survey tree, constructed according to maximum likelihood, shows the approximate placement of the two GPCRs of interest (boxed in blue) within the rhodopsinclass G-protein-coupled receptors. The numbered Ciona sequences and their receptor-type identifications in parentheses are from Kamesh et al. [START_REF] Kamesh | The repertoire of G protein-coupled receptors in the sea squirt Ciona intestinalis[END_REF] . All non-Ciona sequences were downloaded in bulk from http://pfam.xfam.org; they comprise the Pfam 'seed' alignment of seven-transmembrane receptors for the GPCRs. The rhodopsin subclass names are given to the right of each coloured group. The asterisk and double asterisk indicate the genes labelled RXFP3 and SOG/Chemokine-like, respectively, in Fig. 2e.We judge robust nodal support as bootstrap percentages. 70 for minimum evolution (ME) and maximum likelihood (ML) and posterior probability percentages. 95 for Bayesian inference (BI).

Extended Data Figure 7.

Refined maximum likelihood phylogeny of two GPCRs expressed in the PPE-derived neurons of Ciona intestinalis. Relative to the initial survey tree (see Extended Data Fig. 6), analysis of this focused set of 'chemokine cluster' sequences further clarifies the phylogenetic affinities of two Ciona GPCRs expressed in the aATEN-GnRH neurons. The latter sequences are highlighted with violet boxes; all other Ciona sequences are from Kamesh et al. [START_REF] Kamesh | The repertoire of G protein-coupled receptors in the sea squirt Ciona intestinalis[END_REF]

  cells of row IV (yellow) mark the boundary of the neural plate, which gives rise to the anterior sensory vesicle in tailbud embryos. The dotted line indicates the oral opening. b, Schematic of the anterior neural plate border during the gastrula-neurula transition, indicating the lineage-specific expression of Six1/2 (magenta). Six1/2 is initially expressed in eight cells comprising the posterior cells (row V) and the posterior lateral cells (row VI). c, Dorsal view of an initial tailbud embryo co-electroporated with Six1/2. mCherry, Dmrt. GFP and ZicL. CFP. At this stage, the cells initially expressing Six1/2 have divided once, giving rise to 16 cells in total. Brackets indicate the derivatives of the annotated lineages shown in b. Anterior is to the left. Extended Data Figure 2. Six1/2 + cell morphogenesis in Ciona intestinalis from the initial tailbud stage to late tailbud I. a-c, Dorsal views of tailbud embryos electroporated with Six1/2. mCherry (magenta) and Abitua et al. Page 10 Nature. Author manuscript; available in PMC 2016 September 01.

Extended Data Figure 3 .Extended Data Figure 4 .Extended Data Figure 5 .

 345 . The green cells indicate a11.205, which are fated to become PPE-derived GnRH neurons. b, At the initial tailbud stage II, the lateral edges of the Six1/2 + cells begin to intercalate towards the midline. b9, The schematic shows a dotted circle where the future opening of the oral siphon forms. c, At the late tailbud I stage, theSix1/2 + cells have completed intercalation. The bright phalloidin signal in the centre of the pattern marks the apically localized actin of cells constricted towards the oral opening. The arrowheads indicate a PPE-derived cell fated to become a GnRH neuron. At this stage, the Six1/2 + cells are positioned on top of the anterior sensory vesicle over the ocellus. c9, The schematic shows a dotted circle where the oral opening forms. More anterior cells have undergone local cellular rearrangements. d-d0, Anterior lateral view of late tailbud I stage embryo electroporated with Six1/2. GFP and CNGA. mCherry. Asterisks indicate the ocellus. The arrowheads indicate a PPE-derived cell fated to become a GnRH neuron. d, Shows the Six1/2. GFP channel. d9, Shows the CNGA. mCherry channel (see Fig. 3c for larval expression). d0, Shows the merged image of Six1/2. GFP and CNGA. mCherry. Dorsal-ventral BMP patterning during PPE specification in Ciona intestinalis. a, Schematic of the anterior neural plate border during the gastrula-neurula transition. Cell lineage nomenclature is used. Chordin and Six1/2 are co-expressed in the lateral posterior derivatives of rows V and VI (dark blue; also see Fig. 1b). b, Ventral view of mid-gastrula stage embryo hybridized with a BMP2/4 mRNA probe and merged with the nuclear counterstain 49,6-diamidino-2-phenylindole (DAPI). c, Lateral view of an embryo during the gastrula-neurula transition hybridized with a BMP2/4 mRNA probe. Anterior is to the left. d, Dorsal view of an embryo during the gastrula-neurula transition hybridized with a Chordin mRNA probe and merged with a DAPI nuclear counterstain. Endogenous expression of newly described reporter genes in hatched Ciona larvae. a-d, Bright field anterior lateral views. a, Animal hybridized with a GnRH2 mRNA probe. White arrowhead indicates the position of the oral opening. b, Animal hybridized with an RXFP3 mRNA probe. c, Animal hybridized with a SOG/Chemokine-like mRNA probe. The heavily stained tunic was manually removed to reveal the expression signal. d, Animal hybridized with a CNGA mRNA probe. Red arrows mark areas of comparable expression throughout the panels in presumed PPE-derived neurons. An adjacent signal in the ocellus-associated photoreceptors makes it difficult to discriminate expression in PPE-derived neurons in panel d. GnRH expression requires BMP attenuation. a-d, Lateral view of a larva electroporated with GnRH. YFP and counterstained with phalloidin (violet). a, Larva co-electroporated with Dmrt. BMP2/4. Of the 200 larvae, 196 had no GnRH. YFP expression in aATENs and displayed mild to severe morphogenetic defects. Bracket shows mild anterior morphological defect. b, Larva co-electroporated with Dmrt. BMP5/7. Of the 200 larvae, 124 had GnRH. YFP expression in aATENs. c, Larva co-electroporated with Dmrt. BMPR1 CA . Of the 200 larvae, 197 larvae had no GnRH. YFP expression in aATENs. d, Larva co-electroporated with Dmrt. TGF-βR CA . Of the 200 larvae, 130 had GnRH. YFP expression in aATENs and most displayed severe anterior neural tube defects. Arrowheads in b and d indicate the position of GnRH expression in aATENs. Anterior is to the left in all images.

Figure 1 .

 1 Figure 1. Six1/2 expression requires BMP attenuation. a, Dorsal view of an embryo at early neurula stage electroporated with ZicL. GFP and hybridized with Six1/2 and Eya mRNA probes. b, Dorsal view of an embryo at mid-neurula stage co-electroporated with Six1/2. mCherry and Chordin. GFP after an additional division (16 Six1/2 + cells in total). c-f, Tailbud electroporated with Six1/2. mCherry and Dmrt. GFP; counterstained with phalloidin (grey). c, Tailbud co-electroporated with Dmrt. BMP2/4 (131 of 135 embryos showed no Six1/2. mCherry expression). d, Tailbud co-electroporated with Dmrt. BMP5/7 (120 of 135 embryos

Figure 2 .

 2 Figure 2. PPE-derived GnRH neurons express RXFP3. a, Tailbud electroporated with Six1/2. mCherry. Arrowhead indicates the posterior-most Six1/2 + cell that undergoes neurogenesis. b, Larva co-electroporated with GnRH. GFP and Six1/2. mCherry. Insets show Six1/2 expression (left) and GnRH expression (right) surrounding the oral opening (dotted line). The GnRH + neuron is marked with an asterisk. c, Dorsal view of larva electroporated with SOG/Chemokine-like. mCherry. The dotted line indicates the oral opening. d, Larva electroporated with RXFP3. mCherry. e, Phylogenetic placement among human rhodopsin-

Figure 3 .

 3 Figure 3. Ciliated GnRH neurons express functional CNGA. a-c, Larvae stained with phalloidin (grey). a, Dorsal view of a larva stained with an acetylated α-tubulin antibody; dotted line indicates the oral opening. b, Lateral view of a larva electroporated with GnRH. mCherry and stained with an acetylated α-tubulin antibody. c, Lateral view of a larva electroporated with CNGA. mCherry. Arrowheads indicate PPE-derived GnRH neurons in all panels. d, Ratiometric calcium imaging of HEK 293 cells transfected with Ciona CNGA, assayed by the intracellular Ca 2+ indicator Fura-2. The graph shows the efficiency of response of CNGA to the cyclic nucleotide analogues 8-Br-cGMP and 8-Br-cAMP in seven individual cells (coloured lines) over a 120-min period. Indicated are the time points (red arrows) and

  . The three nodal support values are (in order): ME bootstrap percentage, ML bootstrap percentage and BI posterior probability (only values. 50 are shown). Branch lengths are proportional to molecular change (amino acid substitutions per site) between nodes; see scale bar for measurement.Elevation of intracellular Ca 2+ concentration in cells expressing Ciona CNG channels in response to cyclic nucleotides. a, Intracellular Ca 2+ was visualized by Fura-2 ratiometric calcium imaging. Ca 2+ influx into HEK 293 cells transfected with CNGA and CNGC was induced by 10 mM 8-Br-cGMP and 1 mM 8-Br-cAMP. In contrast, no Ca 2+ influx was observed when HEK 293 cells transfected with CNGB were treated with either 10 mM 8-Br-cGMP or 1 mM 8-Br-cAMP. Coloured numbers in the 'Before' panels indicate cells that were subjected to a quantitative measurement of fluorescence. b, Dose-dependent response of HEK 293 cells transfected with CNGA to 8-Br-cAMP. Data were obtained from at least seven different cells in each of 16 different transfections. c, Dose-dependent response of HEK 293 cells transfected with CNGC to 8-Br-cAMP. Data were obtained from at least seven different cells in each of seven different transfections.

	Extended Data Figure 8.

Nature. Author manuscript; available in PMC 2016 September 01.

Acknowledgments

We thank Y. Miyamoto and M. Kotera for technical assistance and A. Stolfi for cloning Chordin. GFP. This work was supported by a grant from the National Institutes of Health (NS076542) and by Grants-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (25650118, 25290067) and from the Japan Space Forum (h160179). Portions of this study were facilitated by the National Bio-Resource Project of the Ministry of Education, Culture, Sports, Science and Technology in Japan. The work of C.H. in the laboratory of H. Yasuo was funded by the Agence Nationale de la Recherche (ANR-09-BLAN-0013-01). P.B.A. and A.N.K. were supported by predoctoral fellowships from the National Science Foundation and California Institute for Regenerative Medicine, respectively.

The coding sequence of SOG/Chemokine-like has been deposited in GenBank/EMBL/DDBJ under the accession KR902348.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.