S.-F Kammerer-Jacquet

Angelique Brunot

Mathilde Lefort

Sahar Bayat

Benoit Peyronnet

Gregory Verhoest

Romain Mathieu

Alexandra Lespagnol

Jean Mosser

Brigitte Laguerre

Solène-Florence Kammerer-Jacquet

Alain Ravaud

Jean-Christophe Bernhard

Frantz Dupuis

Mokrane Yacoub

Marc-Antoine Belaud-Rotureau

Karim Bensalah

Nathalie Rioux- Leclercq

Solene-Florence Kammerer-Jacquet

Metastatic clear cell renal cell carcinoma with a long-term response to sunitinib: a distinct phenotype independently associated with low PD-L1 expression

Keywords: renal cell carcinoma, sunitinib, PD-L1, long-term responders

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype of renal cancer (70%). [START_REF] Srigley | The International Society of Urological Pathology (ISUP) Vancouver Classification of Renal Neoplasia[END_REF] Along with hereditary ccRCC, sporadic ccRCC is also frequently characterized by an alteration of the VHL gene, a tumour suppressor gene, leading to the transcription of genes regulated by HIF such as VEGF, which triggers the angiogenic process. [START_REF] Gossage | VHL, the story of a tumour suppressor gene[END_REF] In 20% of patients, ccRCC is diagnosed at metastatic stage and 30% of the remaining patients will further develop metastases detected during follow-up. With an approximately 50% risk of metastasis, the prognosis of ccRCC is poor, and the mortality rate is 40% at 5 years. [START_REF] Ljungberg | The epidemiology of renal cell carcinoma[END_REF] Anti-angiogenic therapies have significantly improved the prognosis of patients with metastatic clear cell renal cell carcinoma (ccRCC). [START_REF] Motzer | Targeted therapy for metastatic renal cell carcinoma[END_REF] Sunitinib, a tyrosine kinase inhibitor (TKI) targeting VEGFR, is currently an approved first-line treatment option for patients and is the most commonly administered treatment worldwide. [START_REF] Motzer | Sunitinib versus interferon alfa in metastatic renal-cell carcinoma[END_REF][START_REF] Miller | Patterns of Care in Patients with Metastatic Renal Cell Carcinoma Among a U.S. Payer Population with Commercial or Medicare Advantage Membership[END_REF] However, up to 30% patients experience progressive disease (PD) upon response evaluation criteria in solid tumours (RECIST) at their first evaluation and most patients ultimately develop PD. [START_REF] Bergers | Modes of resistance to anti-angiogenic therapy[END_REF] Although several resistance mechanisms have been described, few molecular markers of sensitivity or primary resistance with a prognostic impact have been identified. [START_REF] Garcia-Donas | Single nucleotide polymorphism associations with response and toxic effects in patients with advanced renal-cell carcinoma treated with first-line sunitinib: a multicentre, observational, prospective study[END_REF] Another treatment approach is based on targeted immunotherapy using checkpoint inhibitors because ccRCC is considered an immunogenic tumour with high numbers of mononuclear immune cells such as tumour-infiltrating lymphocytes (TIL). [START_REF] Leite | PD-L1 expression in renal cell carcinoma clear cell type is related to unfavorable prognosis[END_REF][START_REF] Thompson | Implications of B7-H1 expression in clear cell carcinoma of the kidney for prognostication and therapy[END_REF][START_REF] Thompson | PD-1 is expressed by tumor-infiltrating immune cells and is associated with poor outcome for patients with renal cell carcinoma[END_REF] PD-L1 is a transmembrane protein which binds to its co-stimulatory receptor, PD-1 (B7-1), expressed by activated TILs, as a means to down-regulate antitumour immune responses by promoting TIL apoptosis and thus favour tumour progression. Recent clinical trial (Checkmate 025) demonstrated the superiority of anti-PD-1 (nivolumab) over everolimus in second-line treatment. [START_REF] Motzer | Nivolumab versus Everolimus in Advanced Renal-Cell Carcinoma[END_REF] Following the results of Checkmate 214, the combination of nivolumab plus ipilimumab is now recommended as first-line therapy for patients with intermediate or poor prognosis according to International Metastatic Renal Cell Carcinoma Database Consortium (IMDC) model. [START_REF] Powles | Updated European Association of Urology Guidelines Recommendations for the Treatment of First-line Metastatic Clear Cell Renal Cancer[END_REF][START_REF] Heng | External validation and comparison with other models of the International Metastatic Renal-Cell Carcinoma Database Consortium prognostic model: a population-based study[END_REF] Although numerous studies have focused on primary refractory patients defined by progression within the first 3 months, few have studied LTR. [START_REF] Heng | Primary anti-vascular endothelial growth factor (VEGF)-refractory metastatic renal cell carcinoma: clinical characteristics, risk factors, and subsequent therapy[END_REF] [START_REF] Molina | Long-term response to sunitinib therapy for metastatic renal cell carcinoma[END_REF] With the advent of new treatment in ccRCC, it is crucial to better identify them as these patients are the most likely to benefit from sunitinib. [START_REF] Powles | Updated European Association of Urology Guidelines Recommendations for the Treatment of First-line Metastatic Clear Cell Renal Cancer[END_REF] In this study, we aimed to describe the pathological and immunohistochemical phenotype and VHL status of ccRCC and the clinical outcome of patients according to their long-term responder status.

MATERIAL AND METHODS

Patient selection and classification

Primary ccRCC-specimens were collected from patients undergoing nephrectomy in two French University Hospitals (Rennes and Bordeaux) from 1997 to 2013 with metastases either at diagnosis or during the follow up. For inclusion in the study, patients with metastatic ccRCC received sunitinib (50mg/day, four weeks-on/two weeks-off) as first-line treatment (prior cytokine therapy was allowed) and completed at least one 28-day cycle of sunitinib, and undergone their first CTscan assessment with Response Evaluation Criteria In Solid Tumours (RECIST 1.1). [START_REF] Eisenhauer | New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1)[END_REF] Drug schedule and dose-reduction policy complied with local practice guidelines. Follow-up chest/abdomen CT-scans were performed every 2 cycles of treatment (3 months). LTR were defined by at least 18 months of treatment without any progression. [START_REF] Molina | Long-term response to sunitinib therapy for metastatic renal cell carcinoma[END_REF] For each patient, the following clinical and pathologic information was gathered (using data base UroCCR): age, sex, the six factors in the IMDC model (anaemia, thrombocytosis, neutrophilia, Karnofsky performance status <80, <1 year from diagnosis to first-line targeted therapy and hypercalcemia) before sunitinib introduction, pTNM stage at nephrectomy, tumour size, and nucleolar ISUP grade. [START_REF] Ko | The International Metastatic Renal Cell Carcinoma Database Consortium model as a prognostic tool in patients with metastatic renal cell carcinoma previously treated with first-line targeted therapy: a population-based study[END_REF][START_REF] Delahunt | The International Society of Urological Pathology (ISUP) grading system for renal cell carcinoma and other prognostic parameters[END_REF] Histopathologic assessment was performed by three experienced pathologists (SFKJ, MY and NRL). For each patient, frozen ccRCC were available. Informed consent was signed from each patient and institutional review board approval was obtained for this study (CNIL declaration receipt 1812601v0).

Immunohistochemical study

For each primary tumour, a representative slide of the tumour with the highest nucleolar grade and the corresponding paraffin block was selected. Four µm-thick whole tissue sections were cut and mounted on glass slides (Superfrost+, Menzel Glazer). The preparations were dried for 1 hour at 58°C, and then overnight at 37°C. The sections were deparaffinized with toluene

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 6
and rehydrated with ethanol. The preparations were pretreated and immunostained using Ventana Benchmark XT. VEGFA (Anti-VEGF antibody, sc-152, dilution 1/100 ; Santa Cruz Biotechnology, Santa Cruz, CA, USA), CAIX (Anti-CAIX antibody, ab15086, dilution 1/1500, Abcam, Cambridge, UK), c-MET (Anti-Total c-MET, SP44, Rabbit Monoclonal Primary Antibody, ready-diluted, Ventana, Roche, Switzerland), PD-L1 (Anti-PD-L1 antibody, clone 130021, dilution 1/200, RD System, Minneapolis, USA) and PD1 (anti-PD-1 antibody, clone NAT105, dilution 1:50; Abcam, Cambridge, UK) expressions were assessed by immunohistochemistry as previously described. [START_REF] Rioux-Leclercq | Plasma level and tissue expression of vascular endothelial growth factor in renal cell carcinoma: a prospective study of 50 cases[END_REF][START_REF] Patard | Low CAIX expression and absence of VHL gene mutation are associated with tumor aggressiveness and poor survival of clear cell renal cell carcinoma. International journal of cancer[END_REF][START_REF] Shin | Clinicopathologic Analysis of PD-L1 and PD-L2 Expression in Renal Cell Carcinoma: Association with Oncogenic Proteins Status[END_REF] The reactivity of antibodies was revealed with HRP-labeled polymer conjugated secondary antibodies using diaminobenzidine (DAB) as chromogen (Sigma-Aldrich, France). Negative controls were performed by omitting the primary antibody. The tumour expression for each antibody was independently evaluated (SFKJ and NRL), without knowledge of the case. The cut-off for positive cases was 30% of tumour cells for VEGF and 85% for CAIX as previously described. [START_REF] Rioux-Leclercq | Plasma level and tissue expression of vascular endothelial growth factor in renal cell carcinoma: a prospective study of 50 cases[END_REF][START_REF] Patard | Low CAIX expression and absence of VHL gene mutation are associated with tumor aggressiveness and poor survival of clear cell renal cell carcinoma. International journal of cancer[END_REF] For PD-L1 and MET, absent (0), weak (1), moderate (2) and strong expression (3) were reported and cases were then subdivided into negative (0-1) or positive (2-3) subgroups. [START_REF] Shin | Clinicopathologic Analysis of PD-L1 and PD-L2 Expression in Renal Cell Carcinoma: Association with Oncogenic Proteins Status[END_REF][START_REF] Kammerer-Jacquet | Correlation of c-MET Expression with PD-L1 Expression in Metastatic Clear Cell Renal Cell Carcinoma Treated by Sunitinib First-Line Therapy[END_REF] For PD-1, immunostaining density was evaluated in tumour infiltrating lymphocytes and was semiquantified as absent, rare, moderate or dense as previously reported. [START_REF] Schmidt | PD-1 and PD-L1 Expression in NSCLC Indicate a Favorable Prognosis in Defined Subgroups[END_REF]

VHL status

Next generation sequencing

For VHL gene, the entire coding sequence and exon-intron junctions of exons 1, 2 and 3 were analyzed. Genomic DNA was extracted using Magtration System 12GC (Bionobis) according to the manufacturer's instruction. Genomic DNA from all samples was quantitated with the Quan-iT PicoGreen dsDNA assay kit (Thermo Fisher Scientific). DNA target preparation and enrichment were performed by amplification using the Access Array® system (Fluidigm, San Francisco, USA). A 10-nucleotide "barcode" tag, specific to each sample and Illuminaspecific sequencing adaptors were attached using secondary PCR. Purified products were then pooled and sequenced on the Illumina MiSeq NGS instrument (Illumina Inc., San Diego, California).

Multiplex Ligation-dependent Probe Amplification

Methylation-Specific-MLPA (MS-MLPA) was used to detect CpG methylation islands in VHL gene promotor. The SALSA MS-MLPA kit ME001B Tumour suppressor-1 allows detecting aberrant methylation of CpG-Islands located in the promoter region of the VHL gene. [START_REF] Jeuken | MS-MLPA: an attractive alternative laboratory assay for robust, reliable, and semiquantitative detection of MGMT promoter hypermethylation in gliomas[END_REF] The unmethylated DNA will not generate a signal, and a normal probe signal will be

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 7 detected if the site is methylated.

Statistical analysis

The phenotype of ccRCC (histologic and immunohistochemical features, VHL status in long term responders) was compared with other patients. Associations were analyzed with χ², Fisher and Mann-Whitney tests. For logistic regression, we used a backward stepwise selection with p<0.05 inclusion criteria. The rank of elimination was obtained when a variable was removed from the equation, and the odds ratio, 95% CI, and p-value for the removed variables were obtained on the removal step. For clinical outcome, we represented progression-free survival (PFS) from the sunitinib introduction to progression and overall survival (OS) from sunitinib introduction to death using Kaplan-Meier curves. We also represented PFS starting from 18 months for LTR and OS from progression to death. All pvalues were 2-sided, and p-values less than 0.05 were considered statistically significant. All statistical analyses were performed using Stata 14.1 (College Station, TX) software.

RESULTS

Patient and tumour characteristics

Patient and tumour characteristics are presented in Table 1. The study retrospectively included 90 consecutive metastatic patients with primary ccRCC. The mean follow-up period was 25 months (range 1-76 months) from sunitinib introduction. The median PFS and OS were 10 and 22 months, respectively. Eighty-five patients (94.4%) experienced progression and 71 (78.9%) died from their cancer. Men were most represented (62.2% versus 37.8%). According to the IMDC model, most patients were in the intermediate group (47.7%). Locally advanced tumours were mostly represented (stage pT3-T4 in 75.5% of patients), and showed a high nucleolar ISUP grade 3 or 4 (n=88, 92.3%). Metastases were present at the initial diagnosis for 55.6% of patients. More than two-third of patients presented multiple metastatic sites (n=69; 76.7%). The most common metastatic sites were lung (75.6%), bone (53.3%) and liver 2). They more frequently experienced a complete or partial response at the first radiologic evaluation (p=0.035). Corresponding ccRCC were associated with less nucleolar ISUP grade 4 (p=0.037) and less hilar fat infiltration (p=0.006). They were also associated with low expression of PD-L1 by immunohistochemistry (p=0.02), Figure 4.

No association with VHL status was identified. All the significant variables were included for logistic regression except liver metastases that could interfere with Heng score criteria. The 3.

DISCUSSION

To our knowledge, this is the first study to assess the phenotype of ccRCC in LTR along with a long term clinical follow-up. Previous studies identified clinic-biological criteria associated with LTR.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

9 Molina et al. described for the first time the patients who benefit the most from sunitinib at clinical level [START_REF] Molina | Long-term response to sunitinib therapy for metastatic renal cell carcinoma[END_REF] . This retrospective study (n=186) mainly included ccRCC but also other subtypes representing 12% of their cohort even though their carcinogenesis was different.

Moreover, the majority of them (52.1%) received sunitinib in combination with gefitinib, bevacizumab or everolimus whose targets were not the same. In their study, favourable Memorial Sloan Kettering Cancer Center (MSKCC) score was associated with LTR.

A recent study by Escudier et al. proposed a pooled analysis of clinical trials including retrospectively 5714 patients. [START_REF] Tannir | Long-Term Response to Sunitinib Treatment in Metastatic Renal Cell Carcinoma: A Pooled Analysis of Clinical Trials[END_REF] Among them, 897 (15.7%) patients were LTR. They were associated with higher percentage of early tumour shrinkage at the first scan, white race, favourable criteria of IMDC mode, clear cell histology, no liver metastasis, body mass index

≥ 25 kg/m².
Similarly, in our study, LTR exhibited good or intermediate prognosis (IMDC model related to MSKCC), which remained unsurprisingly significant after multivariate analysis. Moreover, they had significantly fewer liver metastases than other patients. Indeed, liver metastases were previously associated with poor prognosis in patients with metastatic ccRCC and the absence of liver metastases was already correlated to LTR. [START_REF] Tannir | Long-Term Response to Sunitinib Treatment in Metastatic Renal Cell Carcinoma: A Pooled Analysis of Clinical Trials[END_REF][START_REF] Mckay | Impact of bone and liver metastases on patients with renal cell carcinoma treated with targeted therapy[END_REF] McKay et al previously described that their impact was attributed to alteration of sunitinib metabolism and the liver microenvironment that could favour an aggressive phenotype. [START_REF] Mckay | Impact of bone and liver metastases on patients with renal cell carcinoma treated with targeted therapy[END_REF] One limitation of our study is the selection of 18 months as a cut-off to define LTR but was reproducible with previous studies. [START_REF] Molina | Long-term response to sunitinib therapy for metastatic renal cell carcinoma[END_REF][START_REF] Tannir | Long-Term Response to Sunitinib Treatment in Metastatic Renal Cell Carcinoma: A Pooled Analysis of Clinical Trials[END_REF] Another limitation is our sample size impairing multivariate analysis. However, contrary to previous studies, we included only clear cell histology and performed a pathological and immunohistochemical study with VHL status on primary tumours along with clinical data.

At molecular level, VHL was inactivated by mutually exclusive mutation or promoter methylation in the majority of cases without any association with the sunitinib response. As previously described, no association between VHL status and outcome was observed. [START_REF] Rini | Clinical response to therapy targeted at vascular endothelial growth factor in metastatic renal cell carcinoma: impact of patient characteristics and Von Hippel-Lindau gene status[END_REF][START_REF] Choueiri | von Hippel-Lindau gene status and response to vascular endothelial growth factor targeted therapy for metastatic clear cell renal cell carcinoma[END_REF] VHL inactivation, considered as an archetypical tumour-initiating event in ccRCC carcinogenesis, failed to identify sunitinib responders, as other mechanisms probably interfere.

We report a distinct pathological and immunohistochemical phenotype of the primary tumours of LTR. As a matter of fact, ccRCC were particularly associated with a lower nucleolar ISUP grade and a less frequent infiltration of hilar fat. They were independently expression with poor outcome in patients with metastatic ccRCC. [START_REF] Choueiri | Correlation of PD-L1 Tumor Expression and[END_REF] Recent update on metastatic ccRCC treatment based on recent Checkmate 214 results recommends a dichotomized approach for first line therapy according to IMDC model. [START_REF] Powles | Updated European Association of Urology Guidelines Recommendations for the Treatment of First-line Metastatic Clear Cell Renal Cancer[END_REF] In intermediate-and poor risk patients, targeted immunotherapies are now indicated whereas in favourable-risk patients, antiangiogenic therapy remains the standard of treatment. In our study, some LTR belonged to intermediate-and poor risk and could have been treated by immunotherapies according to the new standards; however, they clearly benefited from sunitinib. The high expression of PD-L1 associated with better response to immunotherapy in Checkmate 214 trial, could suggest, in LTR associated with low expression of PD-L1, less interest in targeted immunotherapies and reinforce its potential use as predictive biomarker.

Interestingly, the OS was also found to be prolonged in these patients, even from progression.

LTR demonstrated an increased overall survival that was not only explained by their response to sunitinib first-line treatment, but also to second-line targeted therapy. The hypothesis could be that these patients who are more than good responders to sunitinib may be likely to do well under other targeted therapy. Further studies could help clarify whether such responses are specific to sunitinib or reflect underlying favourable biology.

In conclusion, LTR showed a prolonged OS even from progression. Their primary tumours demonstrated a different phenotype with PD-L1 low expression suggesting a potentially lower impact of targeted immunotherapy in these patients.

  characterized them without assessing the phenotype of the tumours.

( 23 .

 23 Figures 1, 2 and 3. LTR had a median PFS of 28 months versus 4 months for other patients (p<0.001). Similarly, LTR had a median OS of 49 months versus 14 months for other patients (p<0.001). From progression, LTR still had a difference of survival, median OS of 21 versus 7 months (p=0.29).

only 2

 2 factors that remained significantly associated with LTR were good or intermediate risk (IMDC model) (p=0.014, OR (95% CI)=5.29 (1.39, 20.01)) and low PD-L1 expression (p=0.029, OR (95% CI)=3.145 (1.22, 8.31)), Table

  expression of PD-L1 by immunohistochemistry. First demonstrated in our study, this is consistent with the results of Choueiri et al. who correlated PD-L1

  from COMPARZ, a Randomized Controlled Trial. Clin Cancer Res. 2015;21:1071-1077. M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Figure 1 :

 1 Figure 1: Kaplan Meier curve representing progression-free survival (PFS) and overall survival (OS)PFS: Median survival: 4 vs 28 months (p<0.001), OS: median survival: 14 vs 49 months (p<0.001)

Figure 2 :

 2 Figure 2: Kaplan Meier curve representing progression free survival of long-term responders (n=28) with time zero starting at 18 months.

Figure 3 :

 3 Figure 3: Kaplan Meier curves representing overall survival (OS) from progression. OS: median survival: 7 vs 21 months (p=0.029)

Figure 4 :

 4 Figure 4: PD-L1 expression A-Absence of PDL1 expression, immunohistochemistry (IHC) x100 B-Low expression of PDL1, IHC x100 C-Moderate expression of PDL1, IHC x100 D-High expression of PDL1, IHC x100

Table 1 : Patient and tumour characteristics of 90 metastatic ccRCC .

 1

	Tables
	Characteristics	No	%
	Age (years)
	Mean		61,3
	Range		37-85
	Sexe
	Male	56	62,2%
	Female	34	37,8%
	IMDC
	Favourable	20	22,2%
	Intermediate	43	47,8%
	Poor	27	30,0%
	T stage
	T1	14	15,6%
	T2	8	8,9%
	T3	63	70,0%
	T4	5	5,6%
	N stage
	N0	75	83,3%
	N1-N2	15	16,7%
	M stage
	M0	40	44,4%
	M1	50	55,6%
	Tumor size (cm)
	Mean		9.6
	Range		2-9.5
	ISUP nucleolar grade
	Grade 2	7	7,8%
	Grade 3	33	36,7%
	Grade 4	50	55,6%
	RECIST 1
	Complete response	1	1,1%
	Partial response	27	30,0%
	Stable disease	34	37,8%
	Progression disease	28	31,1%

Table 2 : Comparison between long-term responders and other patients.

 2

	Variables	Long-term	%	Other patients	%	p-value
		responders		(n=62)
		(n=28)
	Clinical and radiological data
	Male	20	71,4%	36	58,1% 0,226 †
	Age >65	8	28,6%	29	46,8% 0,104 †
	Good/intermediate risk (IMDC)	25	89,3%	38	61,3% 0,007 †
	Distant lymph node metastasis	11	39,3%	28	45,2% 0,603 †
	Pulmonary metastasis	20	71,4%	48	77,4% 0,54 †
	Bone metastasis	16	57,1%	32	51,6% 0,626 †
	Liver metastasis	2	7,1%	19	30,6% 0,015 †
	Cerebral metastasis	4	14,3%	15	24,2% 0,286 †
	Multiple metastasis	21	75,0%	48	77,4% 0,802 †
	RECIST 1 (CR-PR vs SD-PD)	13	46,4%	15	24,2% 0,035 †
	Pathological analysis
	Size >7cm	16	57,1%	42	67,7% 0,331 †
	Nucleolar ISUP grade 4	11	39,3%	39	62,9% 0,037 †
	Sarcomatoid component	3	10,7%	16	25,8% 0,104 †
	Tumor necrosis	19	67,9%	52	83,9% 0,085 †
	Microvascular invasion	11	39,3%	30	48,4% 0,422 †
	T3-T4 stage	19	67,9%	49	79,0% 0,253 †
	Hilar fat infiltration	7	25,0%	35	56,5% 0,006 †
	Peri-renal fat infiltration	14	50,0%	31	50,0%	1 †
	Venal invasion	14	50,0%	27	43,5% 0,569 †
	N1-N2 stage	3	10,7%	12	19,4% 0,375 ‡
	M1 stage	15	53,6%	35	56,5% 0,799 †
	Immunohistochemistry
	CAIX >85%	10	35,7%	18	29,0% 0,526 †
	VEGFA >30%	12	42,9%	38	61,3% 0,103 †
	MET 2-3 intensity	18	64,3%	44	71,0% 0,526 †
	PD-L1 2-3 intensity	16	57,1%	50	80,6% 0,02 †
	PD-1 moderate or dense	18	64,3%	48	77,4% 0,192 †
	VHL status
	VHL mutation	20	71,4%	44	71,0% 0,964 †
	Promoter methylation	2	7,1%	8	12,9% 0,718 ‡
	VHL inactivation	22	78,6%	52	83,9% 0,543 †
	†, χ² test; ‡, Fisher exact test
	RECIST 1: First RECIST evaluation, CR, complete response; PR, partial response; SD, stable
	disease; PD, progression disease

Table 3 : Logistic regression with rank of elimination and intermediate results

 3

	Variables	Rank	OR	95% CI	p-value
	Good/intermediate prognosis (Heng)	-	5.292	1.399, 20.01	0.014
	RECIST 1 (CR-PR vs SD-PD)	3	0.389	0.139, 1.083	0.071
	Nucleolar ISUP grade 4	1	0.509	0.182, 1.419	0.197
	Hilar fat infiltration		2	0.401	0.136, 1.184	0.098
	PD-L1 low intensity		-	3.145	1.122, 8.812	0.029
	Supplementary table
	Variables	Long-term responders	%	Short-term responders	%	p-value
		(n=28)			(n=62)
	Heng score
	Favourable	10	35,7%		10	16,1%	0.039 †
	Intermediate	15	53,6%		28	45,2%	0.46 †
	Poor	3	10,7%		24	38,7%

ACKNOWLEDGMENTS

The authors would like to thank the Ligue Contre le Cancer, CORECT, Rennes Hospital and the French Institute of Cancer (INCa) for their financial support.

The authors would also like to thank the Rennes Hospital Biological Resources Center (BB-0033-00056, http://www.crbsante-rennes.com/) for patient sample management as well as Pascale Bellaud and Roselyne Viel from the Histopathology platform H2P2-BIOSIT at the Faculty of Medicine in Rennes for their technical support.

The authors want to thank Nathalie Costet from Team 9, IRSET and Chloé Rousseau from CIC Inserm 1414 for their contribution to statistical analysis.

AUTHORS' CONTRIBUTION

DISCLOSURE OF POTENTIAL CONFLICTS OF INTEREST

The authors declare that they have no conflict of interest.

RESEARCH INVOLVING HUMAN PARTICIPANTS AND/OR ANIMALS

All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. For this type of study formal consent is not required.

INFORMED CONSENT

Informed consent was obtained from all individual participants included in the study.


K Oh

L Joanny

F Gouttefangeas

B Fabre

V Dorcet

B Lassalle-Kaiser

A Vacher

C Mériadec

S Ababou-Girard

G Loget

Black Silicon Photoanodes Entirely Prepared with Abundant Materials by Low-Cost Wet Methods

Keywords: solar energy conversion, oxygen evolution reaction, nickel, silicon, iron

black silicon (BSi) substrate modified with catalytic metal nanoparticles (NPs). The overall fabrication strategy is cost-efficient as it only requires the use of abundant materials and simple wet procedures such as electrochemical etching and electrodeposition and does not involve the use of buried homojunction and protection layer. Importantly, these results demonstrate that electrodeposited transition metal NPs can stabilize structured Si photoelectrodes without the need for a protection layer.

The depletion of fossil fuels and the increasing environmental pollution imply an urgent development of alternative technologies to foster the energy transition. [START_REF] Lewis | Powering the Planet: Chemical Challenges in Solar Energy Utilization[END_REF] In this frame, the largescale use of solar energy is particularly promising because this is the most abundant renewable energy on Earth. [START_REF] Abbott | Keeping the Energy Debate Clean: How Do We Supply the World's Energy Needs?[END_REF] Solar energy research went through a frantic period after the oil crisis in the early seventies and, in this period, water-splitting photoelectrochemical cells (PECs) started to gain attention. [START_REF] Fujishima | Electrochemical Photolysis of Water at a Semiconductor Electrode[END_REF] These devices are based on light-absorbing semiconductor (SC) materials interfaced with a liquid phase. Under illumination, the charge carriers generated in the SC phase are driven to the solid-liquid interface where they react to generate products in solution. [START_REF] Gerischer | Solar Photoelectrolysis with Semiconductor Electrodes[END_REF] In the case of water-splitting PECs, photogenerated holes (h + ) and electrons (e -) respectively oxidize and reduce water to generate O 2 and H 2 , the latter being a highly promising clean energy carrier. [START_REF] Lewis | Powering the Planet: Chemical Challenges in Solar Energy Utilization[END_REF] An effective utilization of sunlight is required for PEC systems, which can be ideally done by employing small band gap and/or structured SCs. [START_REF] Zhu | Nanostructured Photon Management for High Performance Solar Cells[END_REF] Si is a highly abundant low band gap SC (E g = 1.1 eV) that would be an excellent candidate for manufacturing PECs, especially because this material can be readily processed at the industrial level. [START_REF] Sun | Enabling Silicon for Solar-Fuel Production[END_REF] However, Si is difficult to be employed as a photoelectrode due to its notorious instability. This issue is even more pronounced when Si A c c e p t e d M a n u s c r i p t is used as a photoanode in alkaline media because of: i) photocorrosion which induces fast electrochemical deactivation by generating an insulating SiO x layer [START_REF] Bae | Strategies for Stable Water Splitting via Protected Photoelectrodes[END_REF] and ii) alkaline etching (AE) which causes the irreversible dissolution of the material. [START_REF] Zhang | Electrochemistry of Silicon and Its Oxide[END_REF] The main strategy reported so far for manufacturing stable Si-based oxygen evolution reaction (OER) photoanodes consisted in first modifying the Si surface with a transparent and conformal protection layer such as TiO 2 , [START_REF] Scheuermann | Design Principles for Maximizing Photovoltage in Metal-Oxide-Protected Water-Splitting Photoanodes[END_REF][START_REF] Hu | Amorphous TiO 2 Coatings Stabilize Si, GaAs, and GaP Photoanodes for Efficient Water Oxidation[END_REF] Al 2 O 3 , 11 SiO 2 [START_REF] Ji | Localized Dielectric Breakdown and Antireflection Coating in Metal-Oxide-Semiconductor Photoelectrodes[END_REF] and then applying a co-catalytic (cocat) coating in order to improve oxygen evolution reaction OER kinetics. Recently, unexpected reports from our group [START_REF] Loget | Dispersed Ni Nanoparticles Stabilize Silicon Photoanodes for Efficient and Inexpensive Sunlight-Assisted Water Oxidation[END_REF][START_REF] Oh | Elucidating the Performance and Unexpected Stability of Partially Coated Water-Splitting Silicon Photoanodes[END_REF] and others [START_REF] Hill | An Electrodeposited Inhomogeneous Metalinsulator-semiconductor Junction for Efficient Photoelectrochemical Water Oxidation[END_REF][START_REF] Xu | Silicon Photoanodes Partially Covered by Ni@Ni(OH) 2 Core-shell Particles for Photoelectrochemical Water Oxidation[END_REF][START_REF] Lee | Tailored NiO x /Ni Cocatalysts on Silicon for Highly Efficient Water Splitting Photoanodes via Pulsed Electrodeposition[END_REF] have revealed that heterogeneous nanoscale pinched-off n-Si/metal (metal = Ni or Co) junctions, prepared only by a simple aqueous-phase electrodeposition process, can be employed to manufacture Si photoanodes exhibiting high photocurrent densities and stability without a protection layer. In addition, we have shown lately that these photoanodes can be implemented in complete monolithic water-splitting PEC to produce H 2 under simulated sunlight. [START_REF] Kasemthaveechok | A General Concept for Solar Water-Splitting Monolithic Photoelectrochemical Cells Based on Earth-Abundant Materials and a Low-Cost Photovoltaic Panel[END_REF] The transposition of this concept to structured n-Si surfaces has been reported only once, [START_REF] Tung | Tunable Electrodeposition of Ni Electrocatalysts onto Si Microwires Array for Photoelectrochemical Water Oxidation[END_REF] however, with a Si structuration process based on lithography and vacuum-based physical methods. [START_REF] Tung | Tunable Electrodeposition of Ni Electrocatalysts onto Si Microwires Array for Photoelectrochemical Water Oxidation[END_REF] In addition, there has been recently a strong interest in producing photoelectrode based on black silicon (BSi) (that is, highly absorbing structured Si). [START_REF] Liu | Black Silicon: Fabrication Methods, Properties and Solar Energy Applications[END_REF][START_REF] Oh | Nanoporous Black Silicon Photocathode for H 2 Production by Photoelectrochemical Water Splitting[END_REF] These systems required, so far, the use of a TiO 2 protection layer deposited by atomic layer deposition (ALD). [START_REF] Yu | Enhanced Photoelectrochemical Efficiency and Stability Using a Conformal TiO 2 Film on a Black Silicon Photoanode[END_REF] In this communication, we report the first example of protection layer-free BSi photoanode, prepared by a simple and cost-efficient electrodeposition strategy. The fabrication strategy is described in Scheme 1. First, n-type Si(100) wafers were structured by employing a rapid (<10 min) two-step method previously introduced by our group to fabricate a ~10 m-thick porous Si layer (Figure 1a) composed of vertically-aligned macropores bearing A c c e p t e d M a n u s c r i p t spikes at their extremities (Figure 1b). As previously reported, [START_REF] Loget | Enhancing Light Trapping of Macroporous Silicon by Alkaline Etching: Application for the Fabrication of Black Si Nanospike Arrays[END_REF] these structures efficiently trap the incident light, as shown by their pitch-black visual appearance and their low reflectance (vide infra), we refer herein to them as BSi. Next, we modified the BSi surfaces with a catalytic and stabilizing coating. Toward this goal, we first electroplated Ni 0 on the freshly-hydrogenated BSi using the same method as the one previously employed on planar Si. [START_REF] Loget | Dispersed Ni Nanoparticles Stabilize Silicon Photoanodes for Efficient and Inexpensive Sunlight-Assisted Water Oxidation[END_REF][START_REF] Oh | Elucidating the Performance and Unexpected Stability of Partially Coated Water-Splitting Silicon Photoanodes[END_REF] This electrodeposition leads to the homogeneous decoration of the porous BSi layer with randomly-dispersed Ni nanoparticles (NPs). The role of these NPs role is to generate an array of robust nanoscale Schottky junctions, promoting an efficient hole transfer across the Si/Ni and Ni/electrolyte interface as well as a high photovoltage for sunlight-assisted OER, caused by the pinch-off effect. [START_REF] Loget | Water Oxidation with Inhomogeneous Metal-Silicon Interfaces[END_REF][START_REF] Laskowski | Junction Behavior of N-Si Photoanodes Protected by Thin Ni Elucidated from Dual Working Electrode Photoelectrochemistry[END_REF][START_REF] Tung | The Physics and Chemistry of the Schottky Barrier Height[END_REF] In addition, considering that NiFe alloys are currently among the best OER electrocatalysts, [START_REF] Mccrory | Benchmarking Hydrogen Evolving Reaction and Oxygen Evolving Reaction Electrocatalysts for Solar Water Splitting Devices[END_REF] and that Fe incorporation improves OER kinetics in Si/Ni planar photoanodes, [START_REF] Oh | Elucidating the Performance and Unexpected Stability of Partially Coated Water-Splitting Silicon Photoanodes[END_REF] we have decided to coat the Ni NPs with a NiFe layer by a second aqueous electrodeposition step. To do so, the BSi surfaces modified by the Ni NPs were first thermally oxidized after the first electrodeposition step by placing them into an oven at 90 °C for 1 h, in order to generate a SiO x layer onto the unprotected Si surface (red color in Scheme 1). The presence of this electrically-passivating layer ensures the electrodeposition of the NiFe coating only on the Ni NPs by preventing the deposition on the BSi, as shown in Scheme 1 (note that the electrodeposition was performed employing an aqueous electrolyte containing Ni and Fe salts at a 10/1 molar ratio, see SI and Figure S4 for more details). Scanning and transmission electron microscopy (SEM and TEM) analyses were performed before and after electrodeposition in order to characterize the coatings. Cross-section SEM (Figure 1c), confirmed the presence of randomly-dispersed Ni/NiFe NPs on the walls of the BSi pores and spikes. In addition, TEM allowed a clear morphological characterization of the Si/Ni nanoscale Schottky junctions (Figure 1d). As shown in these figures, the metal particles were hemispherical with diameters below 100 nm.

The chemical composition of the modified surfaces was determined by employing several methods. First, scanning transmission electron microscopy coupled with energy-dispersive X-ray spectroscopy (STEM-EDS) mapping (Figure 2a-d), performed on a BSi area containing 2 isolated NPs (Figure 2a), confirmed the presence of Si in the substrate (red signal in Figure 2b) and Ni in the NPs (blue signal in Figure 2c). However, Fe mapping shown in Figure 2d was not fully conclusive: whereas its presence could be discerned for the largest NP, it was not possible to observe it for the smallest NP. X-ray photoelectron spectroscopy (XPS, survey spectra are shown in Figure S1) confirmed the presence of the two metals after surface modification, as revealed by the appearance of two Ni 2p 3/2 peaks at 851.9 and 855.8 eV (Figure 2e), corresponding respectively to Ni 0 and oxidized Ni (generated during thermal annealing before NiFe deposition), as well as the Fe 2p 3/2 peak (Figure 2f) at 710.8 eV. The analysis of these two regions allowed to estimate a Ni/Fe ratio of 13 ±5, a stoichiometry that fits well with those reported for highly-efficient OER-active Ni 1-x Fe x OOH materials. This oxyhydroxide species is the known OER-active phase for NiFe mixed oxides, that is electrochemically generated at the OER potential. [START_REF] Stevens | Measurement Techniques for the Study of Thin Film Heterogeneous Water Oxidation Electrocatalysts[END_REF] The successful electrodeposition of NiFe on the surface was further confirmed by X-ray absorption near edge structure (XANES, Figure 2g), which showed no signal at the Fe K-edge after the first electrodeposition step but clearly exhibited an edge at E0 = 7133 eV after the second electrodeposition step. These results confirm that our preparation method is effective to modify BSi surface by Ni/NiFe NPs. Next, we study the optical and photoelectrochemical properties of these surfaces.

A c c e p t e d M a n u s c r i p t

The total reflectance of several samples was measured using an integrating sphere in order to quantify their light-trapping ability (Figure 3a). The comparison between the reflectance spectrum obtained for planar Si (purple curve) and that measured for BSi (green curve) shows a considerable reflectance decrease caused by the multiple internal reflections of the incident light within the porous structure. In addition, it is interesting to note that the modification with the Ni/NiFe NPs (red curve) did not affect significantly the reflectance of the BSi, that remained below 10% in the visible range. The surface was tested for OER in 1 M NaOH (measured pH = 13.6), as shown by the cyclic voltammograms (CVs) of Figure 3b. In the dark, no photoresponse could be obtained (black curve). In contrast, high photocurrents were recorded under illumination with simulated sunlight (AM 1.5G, 100 mW cm -2 , red curve). TEM experiments, performed on a photoanode that was employed for OER during 50 CV cycles, demonstrated the good adhesion of the NPs on the BSi surface (Figure 1d). The photoelectrochemical measurements obtained for 10 independently-prepared electrodes (Table S1 contrast, non-modified BSi surfaces failed to promote OER, as shown by the CV obtained under illumination with such a surface (green curve in Figure 3b). Photocurrent spectroscopy, recorded with a modified BSi surface at +2 V, revealed that the incident photon-to-current efficiency (IPCE, red curve in Figure 3b) is comprised between 50 and 68% in the visible range. The integrated value of the photocurrent spectrum (purple curve in Figure 3b) obtained by using the IPCE and the AM 1.5G solar spectrum is in good agreement (22.8 mA cm -2 ) with the value obtained by cyclic voltammetry at the same potential, confirming the validity of our measurements. An optical detection system was employed to measure the quantity of produced O 2 during a 30 min-long preparative electrolysis (red curve in Figure 3d). This measurement, which was performed under illumination, confirmed that O 2 was produced when the modified photoanode was biased at +2 V (t = 15 min), which ceased immediately when the polarization was stopped (t = 45 min). The calculated production rate is 220.6 mol h -1 cm -2 . Based on the electrical charge consumed during this experiment (the chronoamperogram (CA) is shown in Figure S2), we calculated the theoretical amount of O 2 , that is plotted as purple disks in Figure 3d. Taking into account the experimental and the theoretical numbers of mole of O 2 , we calculated a Faradaic efficiency () of 96%. The stability of a photoanode was tested by a prolonged electrolysis test at +2 V in 1 M NaOH (Figure 3e), which revealed that OER could be performed over 16 h with a rather slow degradation of the electrode, as shown by the fact that it maintained a photocurrent higher than 75% of its initial value at the end of the test.

The performance of our photoanodes is compared with that of other reported systems based on structured n-type Si in Table S2. Note that in this table, we have deliberately omitted the electrodes based on a Si homojunction (e.g. np + -Si), that require additional fabrication steps and fall outside the scope of the present article. From that table, it can be observed that our method is the only one that does not imply the use of a conformal protection layer (such as ALD TiO 2 ) and that is entirely based on wet methods. Interestingly, we can also notice that our operation time is the highest that has been reported so far for protection-free structured photoanodes and that the performance we obtained in terms of and operation time is higher than that reported for 𝑗 1.23𝑉

BSi using an ALD-deposited TiO 2 protection layer. [START_REF] Yu | Enhanced Photoelectrochemical Efficiency and Stability Using a Conformal TiO 2 Film on a Black Silicon Photoanode[END_REF] In conclusion, we have reported a new strategy based on simple, low-cost wet methods to manufacture BSi-based photoanodes. The light-trapping n-type BSi substrate is prepared by a ACS Paragon Plus Environment

ACS Applied Energy Materials

A c c e p t e d M a n u s c r i p t fast (<10 min) double-step technique, and then covered by Ni/NiFe NPs using aqueous electrodeposition. In this system, the Ni cores generate an effective pinched-off heterogeneous junction with n-Si [START_REF] Loget | Water Oxidation with Inhomogeneous Metal-Silicon Interfaces[END_REF] and the NiFe layer promote efficient OER catalysis. These photoanodes can operate for 16 h, producing O 2 at high rates with a quasi-quantitative Faradaic efficiency. To date, this is the first example reported of buried homojunction-free and protection layer-free BSi only based on simple wet methods. In addition, these results demonstrate that the stabilization of n-Si by electroplated Ni NPs, which was reported for planar Si(100), [START_REF] Loget | Dispersed Ni Nanoparticles Stabilize Silicon Photoanodes for Efficient and Inexpensive Sunlight-Assisted Water Oxidation[END_REF][START_REF] Oh | Elucidating the Performance and Unexpected Stability of Partially Coated Water-Splitting Silicon Photoanodes[END_REF] is also effective on structured n-Si. We are confident that these results will inspire other groups to develop original designs for photoanodes that will be employed for cost-effective water splitting PECs.

  A c c e p t e d M a n u s c r i p t 1d) as well as the metal NPs (inset of Figure

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. Chemical characterization of the BSi photoanodes. a) Bright field STEM picture

Figure 3 .Scheme 1 .

 31 Figure 3. Optical and photoelectrochemical characterization of the BSi photoanodes. a) Total

  s c r i p t

Funding Sources

This work is funded by Agence Nationale de la Recherche, ANR (project EASi-NANO, ANR-16-CE09-0001-01).

ASSOCIATED CONTENT Supporting Information. Experimental section, XPS survey spectra, supplementary figures, and tables are provided in supporting information. The following files are available free of charge.

AUTHOR INFORMATION

Corresponding Author *gabriel.loget@univ-rennes1.fr

Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

ACS Paragon Plus Environment ACS Applied Energy Materials


Victor Cazalis
email: victor.cazalis@laposte.net

Soumaya Belghali

Ana S L Rodrigues

Paul Caplat

Willson Gaul

Using a large-scale biodiversity monitoring dataset to test the effectiveness of protected areas at conserving North-American breeding birds

Keywords: biodiversity conservation, biodiversity monitoring, protected areas effectiveness, birds, North-American Breeding Bird Survey, conservation effectiveness

Protected areas currently cover about 15% of the global land area, and constitute one of the main tools in biodiversity conservation. Quantifying their effectiveness at protecting species from local decline or extinction involves comparing protected with counterfactual unprotected sites representing "what would have happened to protected sites had they not been protected". Most studies are based on pairwise comparisons, using neighbour sites to protected areas as counterfactuals, but this choice is often subjective and may be prone to biases. An alternative is to use large-scale biodiversity monitoring datasets, whereby the effect of protected areas is analysed statistically by controlling for landscape differences between protected and unprotected sites, allowing a more targeted and clearly defined measure of the protected areas effect. Here we use the North American Breeding Bird Survey dataset as a case study to investigate the effectiveness of protected areas at conserving bird assemblages. We analysed the effect of protected areas on species richness, on assemblage-level abundance, and on the abundance of individual species by modelling how these metrics relate to the proportion of each site that is protected, while controlling for local habitat, altitude, productivity and for spatial autocorrelation. At the assemblage level, we found almost no relationship between protection and species richness or overall abundance. At the species level, we found that forest species are present in significantly higher abundances within protected forest sites, compared with unprotected forests, with the opposite effect for species that favour open habitats. Hence, even though protected forest assemblages are not richer than those of unprotected forests, they are more typical of this habitat. We also found some evidence that species that avoid human activities tend to be favoured by protection, but found no such effect for regionally declining species. Our results highlight the complexity of assessing protected areas effecti veness, and the necessity of clearly defining the metrics of effectiveness and the controls used in such assessments.

Introduction

The increasing human footprint on natural ecosystems is leading to major declines in species ' populations [START_REF] Mcrae | The Living Planet index[END_REF] and has already resulted in thousands of extinctions (IUCN, 2018), to such an extent that [START_REF] Ceballos | Biological annihilation via the ongoing sixth mass extinc tion signaled by vertebrate population losses and declines[END_REF] characterised current times as a period of "biodiversity annihilation". Habitat loss and degradation are the most important pressures on biodiversity [START_REF] Vié | Wildlife in a changing world: an analysis of the 2008 IUCN red list of threatened species[END_REF][START_REF] Balmford | Trends in the state of nature and their implications for human well -being: Trends in the state of nature[END_REF], as a result of anthropogenic activities such as agriculture, urbanisation, industry, transport and recreation [START_REF] Foley | Global consequences of land use[END_REF]. The most evident response to these threats is to establish areas with restricted, or even no human activities, i.e., to create protected areas (PAs). Modern PAs have their origins in the 19 th century and currently represent the most important conservation tool, with about 15% of the global land area already protected to some extent, and coverage planned to reach 17% by 2020 (UNEP-WCMC IUCN, 2016).

Understanding the extent to which PAs are effective as biodiversity conservation tools is thus fundamental for guiding future conservation efforts. Accordingly, there is a substantial literature on PA effectiveness: as of the 1 st October 2018, 260 publications in the Web of Science included in their title "protected AND area* AND effective*". However, within this literature there are disparate approaches to the concept of "effectiveness".

A first set of studies questions whether PAs are effective at representing species or ecosystems, using gap analyses for measuring the overlap between PAs and the distributions of species or of ecosystem types (e.g., [START_REF] Rodrigues | Effectiveness of the global protected area network in representing species diversity[END_REF][START_REF] Brooks | Coverage provided by the global protected-area system: is it enough?[END_REF]. These studies do not directly quantify the effectiveness of PAs at conserving biodiversity, but the extent to which species or ecosystems are buffered from human impacts under the assumption that PAs are highly effective in doing so. A second set of studies focuses on the means employed locally by PA managers in order to protect biodiversity, for example in terms of staff or money (e.g., [START_REF] Leverington | A global analysis of protected area management effectiveness[END_REF]. These analyses do not directly measure PA effectiveness in reducing human impacts, but rather the resources allocated to this purpose. A third type of studies quantifies the effectiveness of PAs at preventing the conversion of natural ecosystems, typically by comparing land use change (e.g., deforestation rates) in protected versus unprotected areas (e.g., [START_REF] Nelson | Protected area effectiveness in reducing tropical defo restation[END_REF][START_REF] Andam | Measuring the effectiveness of protected area networks in reducing deforestation[END_REF]. These studies quantify PA effects at the habitat or ecosystem level, rather than at the species level. Finally, a set of analyses focuses on measuring the effect of PAs on species themselves, either on the diversity of assemblages or on the abundance of individual species, typically by contrasting protected versus unprotected sites. This fourth approach to PA effectiveness is the focus of the present study.

The effectiveness of PAs at conserving species can be assessed by comparing population trends inside and outside PAs (e.g. [START_REF] Gamero | Tracking progress toward EU biodiversity strategy targets: EU policy effects in preserving its common farmland birds[END_REF][START_REF] Devictor | Can common species benefit from protected areas?[END_REF][START_REF] Pellissier | Assessing the Natura 2000 network with a common breeding birds survey: Assessing Natura 2000 with a common breeding birds survey[END_REF]. Indeed, if PAs are effective, populations inside these areas are expected to be better buffered from threats and thus to decline less, or even to increase more, than those outside. Trends however can be misleading, because they are calculated in relation to a reference date (that seldom precedes all anthropogenic impacts) and because they are measured as percentages (which emphasise changes in small numbers). In this study, we focus instead on measures of PA effectiveness that assess current state, namely by contrasting population abundances and species diversity inside versus outside PAs (e.g. [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF][START_REF] Kerbiriou | Common bats are more abundant within Natura 2000 areas[END_REF][START_REF] Devictor | Can common species benefit from protected areas?[END_REF]. These measures combine two types of effects: the effectiveness at selecting as PAs sites of higher-than-ave ra ge conservation interest (i.e. differences that existed at the time of PA designation); and effective ne ss at maintaining species richness and abundance within existing PAs (i.e. differences established subsequently to PA designation).

Three recent meta-analyses investigated the effects of PAs on the state of species abundance and/or diversity by synthesising studies that made pairwise comparisons between protected and unprotected sites [START_REF] Geldmann | Effectiveness of terrestrial protected areas in reducing habitat loss and population declines[END_REF][START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF][START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF]. The underlying studies used in these meta-analyses did not necessarily aim to measure PA effectiveness; more often they investigated the effects of anthropogenic pressure, using PAs as benchmarks (e.g. [START_REF] Sinclair | Protected areas as biodiversity benchmarks for human impact: agriculture and the Serengeti avifauna[END_REF][START_REF] Bihn | Do secondary forests act as refuges for old growth forest animals? Recovery of ant diversity in the Atlantic forest of Brazil[END_REF]Wunderle et al., 2006, all integrated in Coetzee metaanalysis). In the meta-analyses, unprotected sites were treated as counterfactuals to the protected sites (i.e., by assuming that the latter would be in a similar condition to the former if it had not been protected), measuring the effect of protection as the observed difference between the two. These pairwise comparisons often contrast neighbouring sites, which presents the advantage of ensuring that both have broadly similar environmental characteristics (e.g. same climate). However, they do not necessarily take into account the fact that PAs tend to be biased in their location towards higher altitudes and lower productivity areas [START_REF] Joppa | High and far: biases in the location of protected areas[END_REF]. To reduce these biases, [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF] controlled for the differences in altitude, slope and agricultura l suitability. Controlling for these factors means that their results are less influenced by PAs' location biases and, therefore, that they reflect more strongly the effects of protection itself. Another potential confounding effect in pairwise comparisons of neighbouring sites arises from the leakage effect, whereby the human activities that would have taken place inside a PA are displaced to areas around it, artificially inflating the perceived effectiveness of PAs [START_REF] Ewers | Estimates of reserve effectiveness are confounded by leakage[END_REF]. This effect is difficult to control for, but should be reduced if the counterfactua l sites are not immediately adjacent to the PAs.

An important decision when choosing a suitable spatial counterfactual to a PA, one that strongly affects the definition and thus the measure of PA effectiveness, is whether to control for habitat type or not. On the one hand, not controlling for habitat can lead to comparing sites that are not expected to have similar biodiversity regardless of their protection status (e.g. protected grasslands vs unprotected forests). On the other hand, controlling for habitat type can result in an overlooking of the effects that PAs have on biodiversity by preventing habitat changes (e.g. deforestation or urbanisation). For instance, given a hypothetical PA covering a natural grassland, possible counterfactuals include an unprotected site of similar habitat (i.e., natural grassland), an unprotected site with a different type of natural habitat (e.g., forest, wetland), or an unprotected site with human-modified habitat (e.g., extensive pasture, herbaceous cropland, urban area). The choice of counterfactual is certain to have a major impact on the differences observed, and thus on the measure of PA effectiveness, but it is not necessarily obvious which option is the best counterfactual. In theory, it is the site that best represents "what would have happened to the PA in the absence of protection"; in practice, this is not necessarily easily determined. All three metaanalyses [START_REF] Geldmann | Effectiveness of terrestrial protected areas in reducing habitat loss and population declines[END_REF][START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF][START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF] include comparisons where habitat has not been controlled for, meaning that the counterfactual's habitat may be different or similar to the protected site's habitat. Additionally, a subset of [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF]'s analyses focuses on comparisons between protected and unprotected sites with matched habitats. In the latter, the measure of PA effectiveness concerns protection from habitat degradation rather than protection from habitat conversion.

Another key consideration in analysing PA effectiveness is the biodiversity metrics applied to comparing protected and unprotected sites. The three meta-analyses employed a diversity of metrics, some at the level of species' assemblages, some focused on individual species. [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF] used only assemblage-level metrics and found higher species richness and overall abundance inside PAs than outside, but no difference in rarefaction-based richness (i.e. number of species for a given number of individuals) nor in the proportion of endemic species. When matching sites with similar habitats, species richness was only higher in young and small PAs than in unprotected sites (no difference between other protected and unprotected sites), suggesting that the effect of PAs in preventing habitat degradation was light. Conversely, [START_REF] Geldmann | Effectiveness of terrestrial protected areas in reducing habitat loss and population declines[END_REF] considered only species-level metrics (presence, abundance, nest survival) and found contrasted but mainly positive effects of PAs. Finally, [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF] considered both levels; at the assemblage level, they found higher species richness and overall abundance in protected than in unprotected sites; at the species level, they found that individual species abundances were typically higher inside PAs.

An alternative to measuring PA effectiveness through pairwise comparisons is to use statistica l models in which covariates control for differences between protected and unprotected sites. This approach requires access to a large dataset on the spatial distribution of biodiversity, but reduces the subjectivity in the choice of counterfactuals, by making explicit which variables are controlled for, and the measure of effectiveness being investigated. For example, [START_REF] Devictor | Can common species benefit from protected areas?[END_REF] applied this approach to survey data on common birds across France to find a positive effect of PAs on bird abundances for half of the species investigated, especially declining species . [START_REF] Duckworth | Effectiveness of protected areas for bird conservation depends on guild[END_REF], working on bird abundance data collected across South Africa, found that PA coverage was positively correlated with occupancy of frugivorous, insectivoro us, vegivores and predator species, and negatively correlated with occupancy of granivorous species.

In the present study, we quantify the effectiveness of Protected Areas at protecting birds by taking advantage of a large dataset of bird counts across a near-continental areathe North American Breeding Bird Survey [START_REF] Pardieck | North American Breeding Bird Survey Dataset 1966 -2016[END_REF]. Controlling statistically for altitude and productivity in order to reduce the effect of PA location biases, we estimated PA effectiveness on two levels of biodiversity: on species' assemblages, through indices of richness and summed abundance; and on individual common species, by estimating the effect of PAs on species' abundance. At the assemblage level, we expected to find higher species diversity and abundance inside PAs. Indeed, as human activities are causing species population declines and local extinctions [START_REF] Ceballos | Biological annihilation via the ongoing sixth mass extinc tion signaled by vertebrate population losses and declines[END_REF], and as PAs are expected to buffer against these activities, this should predictably lead to overall higher species richness and higher total abundance inside PAs, as found by [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF] and [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF]. At the individual species level, we expected higher abundances within PAs. However, given differences in species' habitat requirements, this result cannot be expected to hold universally (i.e., species are not all expected to be more abundant in all PAs). For example, we expected protected forests to have a positive effect on forest species, but not on grassland species. Hence, we controlled in our analyses for broad vegetation structure (forest, shrub, herbaceous), by investigating separately the effects of PAs dominated by a particular vegetation structure on species with different habitat requireme nts. Additionally, we expected species with overall declining populations (thus more affected by anthropogenic activities), and species that avoid human presence (more sensitive to human disturbance) to present higher abundances inside PAs.

Methods

In this study, we use the term "PA effectiveness" to refer to the difference in species diversity or abundance between protected and unprotected sites. This difference combines the effects of PA site selection (i.e., differences existing prior to the implementation of PA, for example if they are implemented in sites with higher-than-average richness or abundances) and the effects of protection itself (i.e., difference that arise after PA designation, if PAs effectively reduce populatio n depletions and species local extinctions).

Bird data

We used data from the North-American Breeding Bird Survey (BBS), a long-term volunteer-based monitoring scheme in Canada, the USA, and Mexico [START_REF] Pardieck | North American Breeding Bird Survey Dataset 1966 -2016[END_REF], version 2016.0. Our study region encompasses solely Canada and the USA, as few Mexican routes are monitored. This program is based on the annual monitoring of thousands of routes, each approximately 40 km long, during the bird breeding season. Each route is split into 50 stops; at each stop, the observer counts every bird heard or seen during three minutes, before moving to the next stop.

Given the length of BBS routes, they often intersect multiple land use types (e.g. forested, urban, agriculture; each with different bird assemblages), and they are rarely wholly contained within protected sites (most of the routes that cross PAs do so only in small fractions of their length). As a result, whole BBS routes are not particularly suited sampling units for investigating how PAs affect bird species. We chose instead to focus on small sections of routes -sequences of five stops, covering about 4 km -in order to obtain field-sampling units that are more homogeneous in terms of land use types and whose bird assemblages can be more directly linked to local landscape characteristics, especially protection. For each route, we only used the first sequence of five stops, because the only precisely georeferenced point we had access to was the first stop of each route. Indeed, even if in principle additional stops are spaced about 0.8 km from each other, in practice this distance can vary, making the location of additional stops in each route progressively more imprecise. Henceforth, and for simplicity, we use the term "routes" to refer to these initial sections of five stops rather than to entire BBS routes.

We excluded aquatic and nocturnal taxa, which are not well detected by this diurnal road-based monitoring scheme, as well as hybrid individuals. We also excluded seven non-native species, as they are not the focus of conservation efforts. The main dataset we analysed included 400 species in total. To test if removing non-native species can bias analyses (e.g, because they replace native species) we also ran analyses including these species (results are presented in Appendix S5).

Following [START_REF] Kendall | First-time observer effects in the North American breeding bird survey[END_REF]'s recommendations, we removed from the dataset the first year of participation of each BBS observer, to reduce bias due to differences in observer experience. To do so, we extracted the observers' identifying number from the "Weather" file of the dataset [START_REF] Pardieck | North American Breeding Bird Survey Dataset 1966 -2016[END_REF] and calculated for each observer the first year of data collection reported in the dataset. We then removed every observation made by this observer this given year. We also removed double counts, which can be either due to the presence of two observers or an observer sampling several times in one year, by excluding observations for which the 'RPID' code (i.e. Run Protocol type) was 102, 103, 104, 203 or 502. We then focused on routes sampled at least 5 years between 2007 and 2016, obtaining a set of 3,046 routes. For routes sampled more than five years, we analysed only five (randomly selected) years of data, thus ensuring a consistent sampling effort across all routes. We fixed this arbitrary threshold of five years as a trade-off between obtaining high data quantity (number of routes analysed) and data quality (number of species detected per route, which increases with the number of years sampled). For each species, counts were summed across the five points and the five years, to obtain a single value per species per route, which we used as a measure of abundance. We acknowledge that these values correspond only to detected birds rather than true abundances. Detection is known to vary between habitats, depending on vegetation structure [START_REF] Pacifici | Effects of vegetation and background noise on the detection process in auditory avian point-count surveys[END_REF]. This could lead to a difference of detection probability (and thus of perceived abundance) between protected and unprotected sites if vegetation structure differs; controlling for vegetation structure in our analyses should reduce this bias.

Landscape data

For each route, we analysed the properties of the landscape within a 500 m buffer around the route's 4 km track (total area ca. 5 km 2 ). Given that 500 m corresponds broadly to the bird detection radius of the BBS [START_REF] Sauer | The North American breeding bird survey, results and analysis 1966 -2015[END_REF], we considered this a suitable description of the environment affecting the composition of birds detected by the BBS. [START_REF] Small | Local landscape composition and configuration around North American Breeding Bird Survey routes: Ecological Archives E093-215[END_REF] found that the immediate landscape composition (buffer of 0.4 km) of BBS routes was similar to the largescale landscape composition (buffer of 10 km), so we do not expect this choice to strongly affect the results.

A Protected Area is defined by the IUCN as "a clearly defined geographical space, recognised, dedicated and managed, through legal or other effective means, to achieve the long term conservation of nature with associated ecosystem services and cultural values" (UNEP-WCMC IUCN, 2016). PAs are categorised by the IUCN within seven categories based on their protection level. At one extreme, Ia PAs are "strictly protected areas set aside to protect biodiversity […], where human visitation, use and impacts are strictly controlled and limited to ensure protection of the conservation values". At the other extreme, VI PAs "conserve ecosystems and habitats together with associated cultural values and traditional natural resource management systems [and] are generally large, with most of the area in a natural condition, where a proportion is under sustainab le natural resource management and where low-level non-industrial use of natural resources compatible with nature conservation is seen as one of the main aims of the area" (UNEP-WCMC and IUCN, 2018). We used data on locations (polygon shapefile) and IUCN categories of PAs collated in the World Database on Protected Areas (UNEP-WCMC and IUCN, 2018). According to the Word Database of Protected Areas methodology to calculate area covered by PAs (UNEP-WCMC and IUCN, 2019), we excluded "Man and Biosphere" reserves and PAs for which implementation was not finalised, keeping only PAs with a status "designated", "inscribed" or "established". In addition, we also removed PAs that were not spatialized (no polygon associated). Using QGis (QGIS Development Team, 2017), we calculated the proportion of area inside each route's buffer that falls within a PA (all IUCN-categories combined, and dissolved to avoid doublecounting of areas under multiple PA designations). We have also run analyses considering stricter PAs only (categories I-IV), as the effectiveness can vary with protection level [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF][START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF].

We characterised each route according to four additional landscape variables, using QGis (QGIS Development Team, 2017): net primary productivity, altitude, human footprint, and type of vegetation structure. The first three are continuous variables, available as raster files, and we obtained a value per route by calculating the mean value across all pixels that overlap the respective buffer. We calculated net primary productivity as the mean during spring months (March to June) between 2004 and 2015 according to the monthly Net Primary Productivity Terra/Modis (NASA, 2017; resolution 0.1 degree, about 62 km 2 at 45°N). Altitude was obtained from the GLOBE Digital Elevation Model (National Geophysical Data Center, 1999; resolution 0.008 degree, about 0.40 km 2 at 45°N). Human footprint was derived from the 2009 Global terrestrial Human Footprint map [START_REF] Venter | Global terrestrial Human Footprint maps for 1993 and 2009[END_REF]; resolution 0.01× 0.008 degrees, about 0.50 km 2 at 45°N). We defined the vegetation structure as a categorical variable with three types: forest, shrub and herbaceous. We used the Climate Change Initiative -Land Cover layer, using 2011 values as this is the central year of our sampling period (ESA, 2015, resolution 300m) and reclassified the land cover classes into the three vegetation structure types: forest from land cover classes 50-90 and 160-170 (N=1,282 routes; 97 protected by 50% at least); shrub, 120-122 (N=298, 30 protected); herbaceous, 130-153 (which includes croplands; N=1,214; 19 protected). We then obtained the main vegetation structure type for each route as the dominant one in the buffer. Routes dominated by other land use classes (mosaic, [30-40, 100-110 and 180]; bare areas [200][201][202]; water bodies [210]; urban [190], other [220]) were not analysed because they were too scarce. The 2,794 routes used in analyses are mapped in Appendix S1.

Statistical analyses

We estimated the effect of PAs on each of two assemblage-level indices (species richness and summed abundance) and on the abundance of individual species using General Additive Models (GAMs). Models all had identical structures, with the response variable modelled as a function of a one-way interaction between the proportion of PAs inside the buffer and the type of vegetatio n structure. We added smoothed terms controlling for productivity and altitude, as well as longitude and latitude in order to correct for spatial autocorrelation. The general structure of the GAMs was: Response ~ PA * vegetation + s(productivity, altitude, longitude, latitude)

Assemblage-level analyses

For each route, and across all 400 bird species analysed, we calculated two assemblage indices, in each case using the cumulative number of species or individuals seen across the 5 stops, over 5 years: species richness (μ = 28.5 ± 9.4 species); summed bird abundance across all species (μ= 248 ± 90 individuals). We then used a GAM to model each of these two assemblage variables against the above-mentioned covariates, assuming a Gaussian distribution for richness and a negative binomial distribution for abundance.

Species-level analyses

We excluded the rarest species from this analysis in order to have enough statistical power, keeping only the 133 species observed in at least 150 routes. Under this threshold, numerous species were too rarely detected within protected routes, leading to aberrant estimates of PA effectiveness (either highly positive or highly negative). For each species, we only analysed routes that fall within the species' distribution within our study area. We defined this distribution as the 90 % spatial kernel of the routes where the species was observed, obtained using the 'adehabitat' R package [START_REF] Calenge | The package "adehabitat" for the R software: A tool for the analysis of space and habitat use by animals[END_REF]. We treated all routes inside the kernel where the species was not observed as having zero abundance.

We modelled each species' abundance using a GAM as mentioned above, with a negative binomia l distribution. We then calculated for each species a "PA effect" (PAE), measured as the differe nce in predicted abundance between a fully protected (i.e. 100 % protected) and an unprotected route (i.e. 0 % protected) with all control variables fixed to their median values. We calculated PAE separately for each of the three types of vegetation structure, to obtain for each species a value of PAEFor for routes dominated by forest, PAEShrub for shrub routes, and PAEHerb for herbaceous routes.

For each type of vegetation structure, we studied PAE values in order to understand the factors explaining which species are favoured or not by PAs. To do so, we used a linear model and a phylogenetic linear model with similar structures using species-level covariates. We considered three covariates: species' habitat preference, population trend, and human-affinity. We used species' main habitat compiled in 11 classes by [START_REF] Barnagaud | Temporal changes in bird functional diversity across the United States[END_REF] see Fig. 2). We used species' population trends in North America between 1966 and 2015, calculated for each species by [START_REF] Sauer | The North American breeding bird survey, results and analysis 1966 -2015[END_REF] from the BBS data (negative numbers for declining species, positive for increasing species). We winsorized these values, folding down the 2.5% extreme values on each side, bringing estimates to a Gaussian distribution. Finally, we estimated for each species a humanaffinity index, as the median human footprint of the routes where the species was observed, weighted by species' abundance on the route. This index was calculated across all 3,046 routes prior to the exclusion of routes based on habitat types (i.e. also including routes dominated by mosaic, bare areas, water bodies and urban land cover) to be more representative of the diversity of habitats where species are present.

To account for phylogenetic autocorrelation, we ran a Brownian motion model implemented in the 'phylolm' R package [START_REF] Ho | A linear-time algorithm for Gaussian and non-Gaussian trait evolution models[END_REF]. To obtain the bird phylogeny, we selected randomly 100 phylogenetic Hackett backbone trees over 10,000 from [START_REF] Jetz | The global diversity of birds in space a nd time[END_REF] and calculated a maximum clade credibility tree using Tree Annotator from Mr Bayes [START_REF] Drummond | Bayesian Phylogenetics with BEAUti and the BEAST 1.7[END_REF] with no burnin, and node heights calculated with the median. Confidence intervals of the phylogenetic model were estimated with the 'phylolm' function, using a bootstrap with 100 bootstrap replicates [START_REF] Ives | Phylogenetic logistic regression for binary dependent variables[END_REF].

Results

Assemblage-level analyses

At the assemblage level, species richness and summed abundance differed significantly between vegetation structure types (respectively P=0.013, P=4.10 -10 ), underlying the importance of accounting for habitat differences when studying PA effects.

Neither species richness nor summed abundance were significantly affected by the proportion of PAs in models that did not control for vegetation structure (respectively -0.15 ± 0.72, P=0.835 and -0.046 ± 0.031, P=0.143). In models controlling for vegetation structure, species richness did not vary significantly with protection within forest or within shrub routes (respectively -1.39 ± 0.90, P=0.121 and -0.305 ± 1.581, P=0.847) but increased with protection for herbaceous routes (4.35 ± 1.90, P=0.022). Summed abundance lightly decreased with protection within forest routes (-0.084 ± 0.039, P=0.030) but did not vary with protection within shrub or within herbaceous routes (respectively 0.082 ± 0.069, P=0.232 and 0.084 ± 0.082, P=0.307).

Species-level analyses

According to the linear model, values of PAEFor -the predicted difference in a given species' abundance between fully protected versus unprotected forest routesdiffered significa ntly depending on the species' main habitat. Hence, species that have forest as their main habitat showed higher abundances in protected than unprotected forests (Table 1, Fig. 2). This differe nce was significant for the three main forest habitat preferences (i.e. mixed, deciduous, conifer) but not for the general forest category, which only represented five species (Table 1). Species favouring open habitats were significantly less abundant in protected forests than in unprotected forests (Table 1, Fig. 2). We found no significant PA effect within forest routes for species favouring other habitat types. Species' population trends between 1966 and 2015 did not significantly explain PAEFor (Table 1). In contrast, species' human-affinity tended to be negative ly correlated with PAEFor (i.e., we found higher effects of PAs for species with lower affinity to humans in forested routes; Table 1, Fig. 3). This trend was still present when we considered only forest species, but was not significant either (green dots in Fig. 3; see Supporting Information in Appendix S2 for additional test).

None of these effect was significant in the phylogenetic model (Table 1 and Fig. 2), which presented large confidence intervals.

The effect of PAs within shrub routes (PAEShrub) was not affected by species' main habitat and was only significantly negatively correlated with human-affinity (Supporting information, Table S1). PAEHerb was only significantly affected by habitat preferences, being negative for coniferforest species (Supporting information, Table S2). These results, however, need to be interpreted taking into account that shrub or herbaceous protected routes were rare in our dataset: on average, each species' kernel included only 13 shrub and 9 herbaceous routes protected by 50% or more, in contrast with 50 protected forest routes (Fig. 1; see Appendix S1 and S7 in Supporting Information).

Results, both at the assemblage and at the species levels, were similar but less significant when we considered only PAs of stricter management, as defined by IUCN categories I-IV [START_REF] Dudley | Guidelines for applying protected area management categories[END_REF]; see Supporting Information, Appendix S6). For shrub and herbaceous routes, the number of protected routes was even smaller than when all PAs were considered, leading to outlier results. Including non-native species in the analyses had little effect on the results at the assemblage level, the only difference being that the effect of protection within herbaceous routes became no longer significant (Appendix S5). At the species level, adding the three non-native species detected on more than 150 routes did not change the results (Appendix S5).

Table 1: Model summaries regarding the estimated effect of Protected Areas on species within forest routes (PAEFor): linear model and phylogenetic linear model with Brownian motion model. The top part gives estimates and P-values for all covariates, the bottom part gives estimates and P-values for all species' habitat preferences, with trend and human-affinity fixed to zero. N corresponds to the number of species in each case.

* P-values for the habitat variable as a whole could not be obtained, as Anova tables are not implemented in the 'phylolm' package.

Linear model Figure 3: Effect of human-affinity (higher for species found preferably in areas of higher human footprint) on species' responses to Protected Areas within forest routes (PAEFor, above zero for species whose abundance in forest routes is higher in protected rather that in unprotected areas). Forest species (green) are species whose main habitat is "forest", "conifer forest", "mixed forest" or "deciduous forest"; non-forest species (brown) are all other species. Lines correspond to the effect of human-affinity on PAEFor for deciduous forest species (green) and semi-open species (brown), predicted by the linear model and their 95% confidence intervals in dashed lines.

Discussion

We compared the effect of PA coverage on bird species diversity, using assemblage indices (species richness, summed abundance) and individual species' abundances.

At the assemblage level, we found very little effect of protection, only a small increase in richness in herbaceous routes and a small decrease in overall abundance in protected routes. In one sense, this is not surprising, as several large-scale studies found that assemblage metricsparticular ly species richnessare relatively resilient to disturbance through species substitution [START_REF] Dornelas | Assemblage time series reveal biodiversity change but not systematic loss[END_REF][START_REF] Supp | Species -level and community-level responses to disturbance: a crosscommunity analysis[END_REF]. Moreover, areas with low human-induced disturbance can have higher species richness than a pristine area, as predicted by the intermediate disturbance hypothes is [START_REF] Roxburgh | The intermediate disturbance hypothesis: patch dynamics and mechanisms of species coexistence[END_REF]. Accordingly, [START_REF] Hiley | Impacts of habitat change and protected areas on alpha and beta diversity of Mexican birds[END_REF] found lower alpha avian diversity in Mexican PAs than in unprotected areas. However, our results contrast with previous studies investigating this question such as [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF] or [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF], who found a positive effect of PAs on species richness and on summed abundance, including in North America [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF]. These two studies being meta-analyses, it is possible that a publication bias against studies showing negative or null effects of PAs (discussed by [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF] artificia lly increased the difference they measured. This is even more so given that the underlying studies of the meta-analyses were often designed to measure the effect of anthropogenic pressures, using PAs as benchmarks, rather than to measure the effectiveness of PAs (e.g. [START_REF] Sinclair | Protected areas as biodiversity benchmarks for human impact: agriculture and the Serengeti avifauna[END_REF][START_REF] Bihn | Do secondary forests act as refuges for old growth forest animals? Recovery of ant diversity in the Atlantic forest of Brazil[END_REF]Wunderle et al., 2006, all used in Coetzee meta-analysis), and may thus have focused on particularly intact protected sites and/or in highly degraded non-protected sites.

Conversely, our results are not necessarily generalizable to other regions or other taxa, for example if North American birds are less sensitive to human activities than other taxa in North America or than birds in other regions, or if there is less contrast in human impacts in protected versus unprotected areas in North America than elsewhere. In our study, the observed lack of differe nce between protected and unprotected sites in terms of richness and abundance may also be an artefact of differences in species' detectability [START_REF] Boulinier | Estimating species richness: the importance of heterogeneity in species detectability[END_REF], if PAs protect mainly species that are difficult to detect. This detection problem should not affect our result at the species level.

Even when overall species richness and abundance are similar, PAs may nonetheless have an effect on avian assemblages if different species respond differently to protection. We were only able to investigate this in depth for routes whose vegetation was dominated by forest, for which there was adequate sampling quality in PAs (Fig. 1). We found that among forest routes, PAs have an overall positive effect on species' abundance, but only for those species with forest as their main habitat.

In contrast, abundances of species favouring open habitats are negatively correlated with protection in forests. Forest PAs thus maintain a more forest-typical bird assemblage than comparable unprotected forests. These effects were significant with the linear model, but not the phylogenetic model. This suggests that much of the effect attributed to habitat preferences under the linear model relies on phylogenetic relatedness among species, which is not surprising as bird habitat preferences and phylogeny are correlated.

Contrary to our prediction, we did not find that species with low human-affinity (i.e., species that avoid human-impacted areas) are significantly favoured by forest PAs, even if there was as nonsignificant positive effect. Also contrary to our expectation, and to previous results for common French birds [START_REF] Devictor | Can common species benefit from protected areas?[END_REF], we found no correlation between species' population trends over the past 50 years and PAEFor. This may reflect the fact that our model included only relative ly common species (i.e., observed on at least 100 routes in the studied years). It is thus possible that the most human-averse and endangered species are favoured by PAs, but that we could not measure it.

Our results suggest that PAs in herbaceous areas have a negative impact on conifer-forest species and on those with low human-affinity, whereas the effect of PAs in shrub routes was negative ly correlated with human-affinity. Given the scarcity of protected routes within both of these vegetation structure types, we do not consider these results particularly robust or informative of the effectiveness of PAs, but they nonetheless emphasise the biases of BBS routes against shrub areas and herbaceous PAs (Appendix S7).

Given that PAs located in forests are not expected to favour the same species as PAs located in grasslands or shrub lands, we controlled for vegetation structure in our analyses of PA effects . However, this control masked the effect PAs may have had in preventing changes in vegetatio n structure (and associated changes in bird assemblages). For instance, given the vegetation structure categorisation we applied, the counterfactual for a protected forest was an unprotected forest, which does not take into account the possibility that the PA may have prevented the forest from being cleared. In other words, our approach does not measure the effect PAs can have on species diversity by preventing habitat destruction (that modifies vegetation structure type). Instead, it only measures the effects PAs can have in preventing habitat degradation modifying the vegetation structure type), for example from natural forest to exploited forest, or from natural grassland to croplands.

Pairwise comparisons of protected versus unprotected sites, and thus the meta analyses by [START_REF] Geldmann | Effectiveness of terrestrial protected areas in reducing habitat loss and population declines[END_REF], [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF] and [START_REF] Gray | Local biodiversity is higher inside than outside terrestrial protected areas worldwide[END_REF], can integrate the combined effects of habitat destruction and habitat degradation on species diversity, given that the counterfactual chosen may have a different habitat structure from the protected site (e.g., a protected forest compared with an unprotected cropland). Nonetheless, because these metaanalyses build from underlying studies with a diversity of criteria in the choice of the counterfactuals, it is not straightforward to interpret the effectiveness values obtained. For instance, as discussed before, numerous underlying studies compared a highly degraded site with a protected site used as benchmark, with the purpose of estimating the impact of anthropogenic degradation, which can lead to an overestimate of PA effectiveness. Other studies aimed at estimating PA effectiveness directly (e.g. [START_REF] Wasiolka | Comparing biodiversity between protected savanna and adjacent nonprotected farmland in the southern Kalahari[END_REF][START_REF] Lee | The importance of protected areas for the forest and endemic avifauna of Sulawesi (Indonesia)[END_REF], but their choice of counterfactual was subjectively based on what authors considered likely to have happened to the protected site had it not been protected [START_REF] Coetzee | Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta -analysis[END_REF]. Finally, some other studies used in meta-analyses were not particularly interested in differences between protected and unprotected sites, with protection used only as a covariate to potentially explain variation around the signal the authors were interested in (e.g. [START_REF] Naidoo | Species richness and community composition of songbi rds in a tropical forest-agricultural landscape[END_REF][START_REF] Mccarthy | Assessing variation in wildlife biodiversity in the Tien Shan mountains of Kyrgyzstan using ancillary camera -trap photos[END_REF]. So even though our approach does not allow us to measure the full effects of PAs, the differences we measured between protected and unprotected sites are defined statistically depending on the covariates included, which allowed to define more clearly how we measured PA effectiveness. A main advantage of using large biodiversity monitoring datasets (such as bird monitoring schemes) rather than pairwise comparisons is thus the possibility of applying a well-defined and repeatable control.

More broadly, our results emphasise that it is impossible to clearly measure the effectiveness of PAs in conserving species diversity without defining precisely what is expected from them. In this study, we measured PAs effectiveness as the difference in bird diversity between protected and unprotected sites, controlling for main landscape differences. This definition combines conservation ability to protect richest areas and to reduce effectively human pressures in these areas. If PAs are expected to present higher diversity in terms of assemblage metrics (species richness or summed abundance), then we found no evidence in our analyses that PAs are effective.

If PAs are expected to protect all species' populations, then we did not find they were effective either, as for about half of the 133 species studied here we found a negative effect of forest PAs. However, our results show that North-American forest PAs present higher abundances in forest species when compared with unprotected forest sites, and lower abundances of species favouring open habitats. That this result holds even though we found no significant difference in total abundance suggests that bird assemblages in protected forests are more forest-typical than those in unprotected forests. Our results thus indicate that forest PAs in North America are contributing to prevent forest habitat degradation, and associated losses in the abundance of forest specialis t species. BBS routes do not currently cover sufficiently well other habitats besides forest to allow us to investigate whether the same result applies to PAs with a different vegetation structure, but datasets with a bigger proportion of sampling points inside PAs, across all habitats, would bring further light into this question.

Figure 1 :

 1 Figure1: Estimated Protected Areas effect per species (PAE) (represented on a log scale in both negative and positive values), against PAs sampling quality, per vegetation structure type of the routes. PAs sampling quality is the number of routes within the species' kernel with at least 50% of the buffer area covered by Protected Areas. Each point in the plot corresponds to a species, and each species can be represented by up to three points, one for each vegetation structure type.

Figure 2 :

 2 Figure 2: Estimated effect of Protected Areas on species within forest routes (PAEFor) per species' main habitat, estimated with a linear model (blue) and a phylogenetic linear model with Brownian motion (orange). Estimates were all calculated with species' population trend and human-affinity fixed to zero. Error bars represent 95% CI; dot sizes are proportional to the number of species in each habitat group. Stars indicate significant effects for the particular model, for the particular species' main habitat (P: 0.05 < * < 0.01 < ** < 0.001 < ***). Habitat types are ordered fro m the highest to the lowest PAEFor values under the phylogenetic model.

Acknowledgements

We thank Jean-Yves Barnagaud for his insightful comments and suggestions concerning the analyses. We are grateful to the thousands of U.S. and Canadian participants who annually perform and coordinate the Breeding Bird Survey. This preprint has been peer-reviewed and recommended by Peer Community In Ecology (https://dx.doi.org/10.24072/pci.ecology.100018)

Data accessibility

All data used in the study (bird abundances, bird traits and landscape covariates) are public ly available, and can be obtained from the sources listed in the references.

Supplementary material

Script and codes are available online: https://www.biorxiv.org/content/10.1101/433037v4.supplementary-material

Conflict of interest disclosure

The authors of this preprint declare that they have no financial conflict of interest with the content of this article.