Stephanie Hutin

Romain Blanc-Mathieu

Philippe Rieu

François Parcy

Xuelei Lai

Chloe Zubieta

Identification of Plant Transcription Factor DNA-Binding Sites Using seq-DAP-seq

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

Introduction

Proper gene regulation requires the finely orchestrated activity of thousands of transcription factors (TF) and TF complexes. Multicellular organisms dedicate 6-8% of their protein coding genes to the production of TFs and understanding the DNA-binding syntax of these factors is a fundamental challenge in biology (1). In vivo techniques such as chromatin immunoprecipitation followed by polymerase chain reaction (ChIP-PCR) or deep sequencing ChIP-seq have been commonly used to identify TFBSs in vivo, allowing the interrogation of specific potential TFBS or the identification of TFBS genome-wide, respectively (2, 3). New in situ techniques such as Cleavage Under Targets and Release Using Nuclease (CUT&RUN) (4-6) and Cleavage Under Targets and Tagmentation (CUT&Tag) (7) have further refined the identification of in vivo binding by investigating TF binding in the native chromatin environment without a cross-linking and DNA shearing step. CUT&RUN uses a protein A/G tethered microccocal nuclease to release bound antibody-TF-DNA complexes, avoiding disruption of the chromatin. CUT&Tag uses a similar strategy but with a protein A/G-Tn5 transposase preloaded with sequencing adaptors, allowing for release of antibody-TF-DNA complexes and adaptor ligation in a single step. These techniques have certain drawbacks including the need for antibodies specific to the TF of interest or the use of a transgenic line expressing a tagged version of the TF of interest, potential poor signal to noise ratios due to low endogenous TF expression levels and artefacts due to target TF overexpression if a transgene is used. In addition, as TFs often form multiprotein complexes including with other TFs, chromatin remodelers and scaffold proteins, identification of direct DNA binding and modeling of the TFBS of an individual TF or specific TF complex remains a challenge.

In vitro techniques such as protein binding microarrays (PBM) and systematic evolution of ligands by exponential enrichment (SELEX) and SELEX-seq use specific TFs or TFs complexes and either immobilized DNA or a random DNA library to identify TFBS (8, 9). In a PBM experiment, a TF or TF complex is expressed with an epitope tag such as FLAG-tag peptide (DYKDDDDK), Maltose Binding Protein (MBP) or Glutathione S-transferase (GST) for purification and identification of DNA binding events. The TF of interest is applied to a microarrays of immobilized dsDNA from a specific genome or random sequences, the microarray is washed to remove nonspecifically bound TFs and then a primary antibody is used to detect and quantify TF-DNA binding events (9). PBMs allow the modelling of TFBS of a specific TF or TF complex and the identification of putative binding sites in the genome but suffer from the use of short stretches of sequence of 10-12 base pairs that may miss longer TF recognition sites and require the use of a TFBS model to predict genomic binding sites. SELEX and SELEX-seq also use purified TF or TF complexes and query a random synthetic DNA library (8, [START_REF] Tuerk | Systematic evolution of ligands by exponential enrichment: RNA ligands to bacteriophage T4 DNA polymerase[END_REF][START_REF] Slattery | Cofactor binding evokes latent differences in DNA binding specificity between Hox proteins[END_REF][START_REF] Ellington | In vitro selection of RNA molecules that bind specific ligands[END_REF]. Successive rounds of DNA enrichment and amplification are performed allowing the generation of TFBS models or the identification of a single or limited number of high affinity binding motifs as a function of the number of rounds of selection performed. The TFBS models generated from SELEX experiments are not derived from genomic sequences and no direct mapping to the genome can be performed. DAP-seq and seq-DAP-seq use purified TF or TF complexes and fragmented genomic DNA (gDNA) to overcome some of the drawbacks of the in vivo and in vitro techniques described above [START_REF] O'malley | Cistrome and Epicistrome Features Shape the Regulatory DNA Landscape[END_REF](14)[START_REF] Li | Uncovering the "ZIP code" for bZIP dimers reveals novel motifs, regulatory rules and one billion years of cis-element evolution[END_REF][START_REF] Lai | The intervening domain is required for DNA-binding and functional identity of plant MADS transcription factors[END_REF][START_REF] Lai | Genome-wide binding of SEPALLATA3 and AGAMOUS complexes determined by sequential DNA-affinity purification sequencing[END_REF]. Recombinant expression of the epitope-tagged TF or TF complex is performed using a cell or cell-free expression system. One or more purification steps of the target TF or TF complex using multiple affinity tags may be required to ensure a single homogeneous species to be assayed. The purified TF complex of interest is immobilized and used to pull down specific binding sequences from a gDNA library pre-ligated with adaptors, the bound fragments are eluted, PCR amplified with barcode sequences and analyzed using next-generation sequencing, allowing for genome-wide binding-site mapping reflective of both local sequence context and DNA methylation status (13, 14). Alternatively, the gDNA library may be amplified prior to the DNA-binding step to remove any modifications such as methylation which may be present and may impact TF binding or site recognition (18). DAP-seq benefits from the use of purified specific TFs or TF complexes and the simplified in vitro system eliminates the complexity of different complexes containing the TF of interest that may form in vivo and alter DNA-binding specificity. In addition, by examining fragmented genomic DNA, the noise generated from non-specific binding is spread over the entire genome, thus eliminating the need of selection and amplification rounds required by SELEX-seq, for example. While DAP-seq and ampDAP-seq protocols have been well-described (13-15), we focus particular attention on the production and sequential pulldown steps required for working with heteromeric TF complexes (seq-DAP-seq) and provide the initial steps of an easy-to-follow data analysis pipeline.

Here we present a detailed schematic protocol designed for the seq-DAP-seq novice. The protocol includes (1) the preparation of a fragmented gDNA library from plant genomic DNA, ligation to Y-adaptors and optional PCR library amplification; (2) the production and purification of a heteromeric TF complex using multiple epitope tags and sequential pull downs;

(3) DNA binding to the TF complex and preparation for next generation sequencing and (4) a basic bioinformatics data analysis pipeline and on-line resources for processing DAP-seq sequencing data.

Materials

Equipment

Bioinformatics requirements

1. Basic knowledge of Bash language is needed to run the pipeline.

2. Dependencies are: FastQC [START_REF] Andrews | FastQC: A Quality Control Tool for High Throughput Sequence Data[END_REF], NGmerge (20), Bowtie 2 (21), Samtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] Bedtools [START_REF] Quinlan | BEDTools: The Swiss-Army Tool for Genome Feature Analysis[END_REF], MACS2 [START_REF] Gaspar | Improved peak-calling with MACS2[END_REF], MSPC, MEME Suite [START_REF] Bailey | MEME SUITE: tools for motif discovery and searching[END_REF], Python and R scripts available on the GitHub page provided (https://github.com/Bioinfo-LPCV-RDF/DAPseq_TFanalysis) and optionally SRA Toolkit if DAP-seq sequence reads are downloaded from the NCBI.

Methods

DAP-seq and ampDAP-seq input library preparation

The same genomic DNA library can be used for several individual DAP-seq and seq-DAP-seq experiments, increasing assay throughput and limiting library-to-library bias. Illumina adaptor barcoding sequences may be chosen for multiplexing and are added at the final PCR amplification step. For ampDAP, libraries with different barcodes can be prepared in parallel but should not be combined until final multiplexing. This DNA extraction and library preparation protocol is for gDNA from Arabidopsis thaliana Col-0 but genomic DNA from the target plant of choice may be used. The libraries may be prepared and stored at -20°C or -80°C for later use.

1.

Start the DNA library preparation with and without endogenous DNA modifications by extracting Arabidopsis Col-0 gDNA using a kit following the manufacture's manual or using standard CTAB protocols (26). Scale-up as required for the number of individually barcoded libraries that are desired. Each library preparation requires 5 µg of high quality gDNA. Perform library quantification using a Qubit kit. Note 1 2. Dilute 5 µg gDNA to a final volume of 125µl in EB into a microTUBE-130 AFA.

The protocol can be safely interrupted at this point

7.

Centrifuge the sample in a tabletop centrifuge at 4 °C for 30 minat 16,000 × g to pellet the DNA.

8.

Carefully remove the supernatant by pipetting without disturbing the DNA pellet.

Carefully inspect the tube for the DNA pellet which is usually visible as a clear/white pellet and take care not to touch the pellet with the pipette tip. 9. Add 1 mL of 70% ethanol to wash the DNA pellet.

10. Centrifuge the sample for 2 min at 16,000 × g at 4 °C. Carefully discard the supernatant by pipetting in order to remove as much of the remaining ethanol as possible.

11.

Air dry the DNA at room temperature by opening the tube on the benchtop for 5-10 min. Allowing the pellet to dry for too long may make it more difficult to resuspend.

12. Resuspend the DNA pellets in 34 μL of EB buffer by pipetting up and down 40x. Vortex to mix the sample.

18.

Repeat steps 6-11.

19.

Resuspend the DNA pellet in 32 μL of EB buffer by pipetting up and down 40x.

20.

If the pellet is still visible, place the sample at 37 °C in a water bath for 5 min to help resuspend the DNA and pipette up and down 10x.

[Figure 3 near here] 21. Transfer 32 µl sample to undergo an A-tailing reaction (Figure 3) into a clean PCR tube and add 5 µl 10x ThermoPol Buffer, 10 µl 1 mM ATP, 1 µl Taq DNA polymerase and mix gently.

22.

Centrifuge at room temperature for 5 s at 3,000g to collect any sample adhering to the walls of the tube.

23.

Incubate at 72 °C for 30 min. in a thermocycler.

24.

Transfer to a clean 1.5 ml reaction tube.

25.

Repeat steps 17-18.

26.

Resuspend DNA in 35 μl of EB. Place the sample at 37 °C in a water bath for 5 min to aid resuspension.

27.

Measure the DNA concentration using 2 µL DNA and a Qubit with dsDNA Qubit Assay kit.

[Figure 4 near here] 28. To prepare the Y-shape adaptor, mix equimolar quantities of Y-shaped adaptor A and adaptor strand B (Figure 4) at 30 µM concentration for each in a 1.5 ml reaction tube in EB.

29.

Heat in a thermocycler or heating block to 95°C for 15 min.

30.

Heat water to 95°C in a 500 ml beaker.

31.

Place the tube with the adaptors to be ligated in a foam floating tube rack and place in beaker. 35. Centrifuge at room temperature for 5 s at 3,000x g to collect any sample adhering to the sides of tube.

Let cool on

36.

Incubate at room temperature for 3 h.

Perform Ethanol precipitation (steps 17-18).

38.

Resuspend the DNA in 30 μl of EB. Place the sample at 37 °C for 5 min to aid resuspension.

39.

Measure the DNA concentration of 2 μl of samples of the Y-adaptor ligated sample using a Qubit kit and compare to the concentration before the ligation (step 27). This is the DAP-seq input library.

40.

Test 15 ng of the library before and after adaptor ligation with qPCRs. In a PCR tube, mix 15 ng of pre-ligated DNA (step 26) or DAP-seq input library with 1 µL TruSeq Illumina universal primer (25 μM), 1 µL of any index primer (25 μM) and top up to 10 µL with nuclease free water. Add 10 µL of 2x Syber Green Master Mix.

41.

Load the samples into a real time qPCR machine and run the following program: 98 °C for 2 min, 98 °C for 30 s followed by 98 °C for 15 s 60 °C for 30 s 72 °C for 2 min (30 cycles), 72 °C for 10 min. Compare non-ligated and ligated samples using Absolute Quantification option to verify successful ligation.

Adaptor ligated gDNA may be used for a DAP-seq experiment at this point (DAP-seq input library) if step 39 is successful. Note 6 and Note 7 3.1.1 ampDAP libraries without DNA modifications (Figure 5)-[Figure 5 near here] 1. Take 15 ng of gDNA sample after adaptor ligation (step 37).

Perform PCR amplification using TruSeq Illumina universal primer and the index

primer with distinct barcode sequence (example set of 12 unique index primers are given in materials list). Note 8 3. In a PCR tube, mix 15 ng of template DNA with 10 µL (5×) Phusion HF Buffer, 2.5 µL dNTPs, 1 µL TruSeq Illumina Universal primer (25 μM stock), 1 µL index primer with distinct index (25 μM stock), 1 µL and Phusion DNA Polymerase (2,000 U/ml) and top up to 50 µL with ddH2O. Perform for each indexed adaptor library desired.

4.

Load the samples into a thermocycler and run the following program: 95 °C for 2 min, (1 cycle); 98 °C for 30 s, 98 °C for 15 s, 60 °C for 30 s, and 72 °C for 2 min (12 cycles),

5.

Perform PCR clean-up using the Qiagen PCR clean-up kit according to the manufacturer's protocol.

6.

Quantify library using Qubit Fluorimeter with dsDNA HS Qubit Assay kit according to the manufacturer's protocol.

7.

Aliquot into 50 ng samples and freeze at -20 °C or -80 °C for long term storage. Note 9

The protocol can be safely interrupted at this point 3.2 TF expression, TF complex formation and sequential pull down -1. Clone each gene encoding the TF of interest into a pTnT vector with appropriate epitope tag (i.e. 3xFLAG or 5xMyc) to obtain the plasmids of interest.

2.

Plasmid DNA is amplified in E. coli DH5α and purified using a Qiagen miniprep kit or similar according to the manufacturer's protocol (Note 10). At least 2 µg of purified plasmid DNA is required per TnT reaction at a suitable concentration (~200 ng/µl).

Equimolar amounts of plasmid for each TF to be expressed should be used in a single reaction mix. DNA concentration may be measured on a microvolume spectrometer and purity determined by the A260/280 ratio (ratio should be > 1.8).

TnT expression systems such as the TnT SP6 High-yield Wheat Germ Protein

Expression System from Promega is used to express the targets proteins according to the manufacturer's protocol.

4.

Place 30 µl of TnT mix in a PCR tube. For single input reactions, 2-4 µg of total plasmid is used in a 50 µl final reaction volume. For double input reactions, 2 µg plasmid of each construct is used in a 50 µl reaction volume, top up to 50 µl with ddH2O. Note 11 5. Reactions are incubated for 2 h at 25 °C in a thermocycler to maintain a constant temperature.

6.

Reserve 2 µl of reaction mix for western blots to confirm protein expression with the appropriate anti-FLAG or anti-Myc antibody. Commercially available antibodies antic-Myc and anti-FLAG antibodies at a dilution of 5000 x are sufficient to determine expression of the full-length protein.

1. Prepare an SDS gel or use a commercially available one. The percentage of acrylamide should be chosen according to the size of the protein.

2. Place the gel in an adequate running equipment in 1x SDS running buffer.

3. Load the gel with your sample and a ladder.

4. Run a SDS-gel at 180 V until the proteins are well separated. An indication are the dyes in the loading dye.

5.

Cut 4 pieces of Whatman paper slightly bigger than your gel (e.g. 9x9 cm) and a membrane in the size of the gel (8x8 cm).

6. Wet all components with blotting buffer.

7.

Build the blot: In the cassette place on the -side a sponge tissue, then two Whatman papers followed by the membrane followed by the gel, then two Whatman paper. Pay attention, once the gel touches the membrane it should not me moved. Try to avoid air bubbles by wetting the components in between the stacks.

8. Roll over the stack with a 10 mL pipet to expel air bubbles.

9.

Place the second tissue (sponge) and close the cassette.

10. Place the cassette in the running chamber in the correct orientations -to -and + to +).

11.

Transfer over night to a membrane at 10 mA in a cold room.

12. Disassemble the blot. Handle the membrane with forceps.

13.

Transfer the membrane into a blot box and add 5 % Milk powder/1xPBS/0.1% Tween 20 for 1h.

14. Dilute the primary antibody (1:200, 1:500, 1:1000) in 1 mL 5% milk powder/1x PBS/0.1% Tween 20.

15.

Place the membrane onto a clean surface (e.g., clean bench top or glass surface).

16. Distribute the 1 mL of primary antibody solution onto the membrane taking care to fully cover the membrane and seal with a plastic wrap. 17. Incubate 1h at room temperature.

22.

Expose the membrane in an imager or using ECL-sensitive film.

23.

Heat the stripping buffer to 50˚C in a water bath.

24.

Incubate the membrane with the stripping buffer at 50˚C for 15 min two times.

25. Wash the membrane with 1x PBS/0.1% Tween 20 several times for 5 min until the stripping buffer is completely removed .

26.

Block the membrane and reuse to check for the second epitope tag.

27.

Repeat steps 19-22 with the second antibody.

For recombinant protein expression in bacteria see Note 12 and 13.

7.

Prepare the magnetic beads by resuspending them by inverting the storage tube 5x.

8.

Pipette 100 µl of magnetic beads (anti-FLAG or anti-Myc magnetic beads) in storage buffer into to a clean 1.5 ml reaction tube and add 200 µl of IP buffer. Note 14 9. Place the tube on a magnetic rack and pipette off the buffer once the solution has cleared and all the beads have been drawn to the magnet.

10.

Remove the tube from the magnetic rack and repeat the wash step 3x wash with 200 µl of IP buffer by gently pipetting up and down 5x.

11.

Add 100 µl of IP buffer and transfer the magnetic beads in a clean tube and store on ice.

[Figure 6 near here] 12. In a clean 1.5 ml reaction tube, 50 µl (Step 3.1) of the TnT reaction is combined with 130 µl cold IP buffer and 20 μl anti-FLAG magnetic beads in order to start the first pull down (Figure 6). Note 15 and 16 13. Gently mix by pipetting up and down 3x.

14.

Incubate 1 h at 4 °C on a rotating wheel or by gently tapping the tube every 5 min to resuspend the beads.

15.

Immobilize the anti-FLAG magnetic beads by placing the reaction tubes on a pre-chilled magnetic rack and carefully remove supernatant.

16. Wash three times with 100 μl IP buffer, gently pipetting up and own each time to resuspend the beads.

17.

Remove and discard supernatant.

18.

Elute the FLAG-tagged proteins or protein complexes by addition of 175 µl of IP buffer supplemented with 200 µg/ml 3xFLAG peptide.

19.

Immobilize the beads on the prechilled magnetic rack and carefully collect the supernatant.

20.

Place the supernatant containing the eluted TF complexes in a clean 1.5 ml reaction tube and add 20 µl of washed anti-c-Myc magnetic beads (step 11) (Figure 7). [Figure 7 near here]

21.

Incubate 1 h at 4°C on a rotating wheel or by gently tapping the tube every 5 min to resuspend the beads.

22.

Immobilize the beads on a magnetic stand and carefully pipette off the supernatant and discard it.

23.

Wash the beads 3x with 100 µl of IP buffer and resuspend in 50 µl volume of IP buffer.

The purified TF or TF-associated complex, while still bound on magnetic beads are used for seq-DAP-seq. 8 near here] 1. Resuspend the beads in 100 μl of IP buffer or optimized binding buffer.

Binding of DNA to immobilized TF complex-[Figure

2.

Add 50 ng DAP-seq input library pre-ligated with Illumina Y-adaptor (DAP input library step 38) or ampDAP input library (step 47) (Figure 8).

3.

Incubate at 4 °C for 90 min on a rotating wheel.

4.

Place the samples on a prechilled magnetic rack and remove the supernatant by pipetting.

5.

Wash the beads six times with 100 µL IP buffer or optimized binding buffer by gently pipetting up and down 2x while keeping the tube on ice.

6.

Transfer the buffer and beads to a new PCR tube and wash 2x with 100 µl IP buffer or optimized binding buffer by pipetting gently up and down 2x.

7.

Place the tube on the prechilled magnetic rack. Remove the supernatant and resuspend the beads in 30 μl of EB and transfer to a PCR tube.

8.

Heat the PCR tube in a thermocycler at 98 °C for 10 min.

9.

Resuspend beads and transfer to a 1.5 ml reaction tube.

10. Place the tube on a magnetic stand and collect the supernatant (Figure 8). Note 17 and 18

The protocol can be safely interrupted at this point samples may be frozen at -20 and stored

for later PCR amplification and sequencing.

DNA amplification, pooling and sequencing

1. In a PCR tube, mix 25 µL of the supernatant with 10 µL (5x) Phusion HF Buffer, 2.5 µL dNTPs, 1 µL TruSeq Illumina Universal primer (25 μM stock in EB), 1 µL index primer with distinct index (25 μM stock in EB), 1 µL and Phusion DNA Polymerase (2,000 U/ml), and top up to 50 µL with ddH20.

2.

Load the samples into a thermocycler and run the following program: 98 °C for 2 min, then 98 °C for 30 s, 98 °C for 15 s 60 °C for 30 s 72 °C for 2 min (20 cycles), 72 °C for 10 min.

3.

Purify the PCR products using AMPure XP beads. An optional qPCR verification step can be performed if the TF complex of interest has known binding sites. Note 19 and Figure 9 [Figure 9 near here] 4. Quantify DNA concentration using a Qubit fluorimeter and dsDNA HS assay kit. DNA concentration should be > 10 ng/ul.

5.

Measure the molarity of the libraries using the NEBNext Library Quant Kit for Illumina following the manufacturer's instructions. Note 20 and 21 6. Combine 5 μl of each library with different barcodes in equal molarity in a 1.5ml tube.

The rest of the reactions can be stored at -20 °C. Library concentration should be >10 ng/µl.

7.

Quantify DNA concentration using a Qubit fluorimeter and dsDNA HS assay kit.

8.

Run the sample on an Illumina sequencer such as MiSeq Instrument according to the manufacturer's specifications. Multiplexing of equimolar DNA libraries may be performed at this point.

9.

Use paired-end sequencing of 150 cycles for 10-20 million reads per lane. 10)-[Figure 10 near here] Processing of DAP-seq sequence reads to identify bound regions and determine TFs or TF complex binding sites can be decomposed in five main steps (Figure 10). A bioinformatics pipeline encapsulated in a Conda environment is provided at https://github.com/Bioinfo-LPCV-RDF/DAPseq_TFanalysis. The pipeline provides an example DAP-seq sequencing dataset for BIM2, an Arabidopsis thaliana TF, with sequencing data from an ampDAP-seq experiment with paired end reads with two replicates and an input ampDAP DNA library as control. An example of a successful and unsuccessful seq-DAP-seq experiment are shown in Figure 11.

Bioinformatic analysis (Figure

Explanations and main commands for each step are described below. Note 22 [Figure 11 near here] 1. Check the reads quality by using FastQC [START_REF] Andrews | FastQC: A Quality Control Tool for High Throughput Sequence Data[END_REF]. Note 23. Adapters are removed using NGmerge [START_REF] Gaspar | NGmerge: merging paired-end reads via novel empirically-derived models of sequencing errors[END_REF] with the -a option and reads are trimmed using Trimmomatic (27).

Commands are given below to run the bioinformatics pipeline with example or user provided data. TF_1 and TF_2 refer to forward and reverse reads, respectively. ). We do not recommend using a microvolume spectrometer for gDNA library quantification as the measured concentration may be higher due to contaminating molecules that absorb at A260. If using a spectrophotometer, the A260/A280 ratio should be > 1.8 for a high quality gDNA library.

2.

All ethanol precipitation steps can be replace by using Qiagen MinElute Reaction Cleanup kit. This minimizes the time for DNA purification/concentration to 10 mins. We have good results with kits for DNA purification/concentration steps.

3. DNA end repair and 3'dA tailing may be performed as one step using kits such as NEB End Repair/dA-Tailing Module. We have separated the steps to illustrate each step of the protocol.

4.

End-It™ commercial kits have all the required components for creating blunt ended DNA premixed and reduce pipetting errors and we recommend using these. The concentrations of the various component are given in parentheses.

5.

Steps 27-31 may be performed prior to or during End repair and A-tailing reactions to save time.

6. DNA quantification using a Qubit fluorimeter and dsDNA HS assay kit gives more accurate quantification than a microspectrophotometer.

7.

If PCR amplification of the adaptor ligated library requires more than 15 cycles to reach saturation, this may indicate a problem with the library and we recommend using a new gDNA extraction. DNA loss during wash steps and residual ethanol are likely causes of poor qPCR results.

8.

Care should be taken to avoid contamination of libraries or accidental pooling after bar code addition. A detailed example of single barcoding with one index primer and the universal primer is shown in Figure 12. [Figure 12 near here] 9. We have stored DAP and ampDAP libraries at -20°C for up to 1 year without an effect on downstream steps.

10.

Miniprep DNA extractions may be pooled to obtain sufficient plasmid DNA for cell free protein expression. Alternatively, phenol:chloroform:isoamyl alcohol (25:24:1) DNA extractions may be performed. The quantity and quality of plasmid DNA obtained from phenol chloroform extraction is generally better than that obtained from commercial miniprep kits.

11.

We have successfully produced up to four protein complexes in a single reaction mix.

Performing the protein expressions in a single TnT reaction may be important for the generation of stable heteromeric TF complexes. In our hands, expression of proteins from multiple plasmids in a single reaction volume results in better protein yields and active TF complexes when compared to separate reactions and subsequent mixing of the protein partners.

Stoichiometry of input plasmid should reflect final complex stoichiometry. However, in certain cases, plasmid concentrations may need to be optimized to produce active complexes, particularly if the proteins have different stabilities. These proteins can be produced in excess to ensure active complex formation. We recommend testing by EMSA with different plasmid concentrations.

12. Standard expression and purification of tagged proteins from E. coli or other expression systems may be used to obtain TF or TF complexes of interest. Affinity purification using tags including 6x-His, GST or MBP may be used for purification from cultures of 10 ml. For DAPseq applications 1 µg of purified protein is sufficient and is generally attainable from a 10 ml culture, although this may vary depending on the TF target. While protein tags (GST or MBP)

and His tags are often used for recombinant protein purification on inexpensive resins, we recommend including a 3x FLAG tags or 5x Myc tags in the construct for protein immobilization and DNA binding experiments due to the high affinity binding of 5x Myc and 3x FLAG tagged proteins to anti-FLAG and anti-Myc magnetic beads. We do not recommend the use of His tags and Ni-NTA magnetic beads due to nonspecific interactions with the resin and leaching/ protein loss during the DNA binding and washing procedure.

13.

Care should be taken with protein fusion tags. If the DNA-binding domain is N-terminal the tag should be located on the C-terminus and vice versa. Large tags may interfere with DNAbinding. Bulky protein tags such as MBP and GST may cause steric hindrance of the DNAbinding domain, resulting in poor DNA binding efficiency and care should be taken with tag location and/or inclusion of a spacer sequence between the tag and TF of interest.

14.

The buffer is critical for successful seq-DAP-seq. If an optimized buffer from EMSA or other experiments is known, we recommend using this buffer. For most TF complexes IP or EB buffer may be used without dithiothreitol or beta-mercaptoethanol as these reducing agents may negatively affect performance of anti-FLAG and anti-Myc magnetic beads. TCEP does not affect bead performance and may be used at concentrations up to 1 mM.

15.

All buffers and reagents should be at 4 °C and kept on ice.

16.

For sequential pull-down reactions, we prefer to use anti-FLAG magnetic beads and a FLAG peptide elution step followed by immobilization on anti-Myc magnetic beads due to the efficient elution of the proteins using 3x FLAG peptides and the high affinity and specificity of the anti-Myc magnetic beads.

17 Add 20 µl of 1x SDS-PAGE Protein Sample Buffer to the remaining beads for western blot verification of protein complex. Boil 20 µl of resuspended beads for 1 min, place on the magnetic stand and run 10 µl of supernatant on a 12% SDS PAGE for Western blot to control the immunoprecipitation.

18. These DNA-binding steps (1-9) should be performed in triplicate or quadruplicate for robust statistics.

19. AMPure beads are efficient for purification of DNA fragments and size selection. We have found that 1.2-1.5x volume of beads, 5 min. binding of DNA followed by 2 washes in 70% ethanol, 5-10 min. drying on the benchtop at room temperature and 30 µl elution (EB) give the best results. A slightly higher volume of beads is recommended by the manufacturer (1.8x).

20.

Optional size selection may be performed as follows. We have found that these steps decrease the amount of DNA available for sequencing and are often not necessary.

1. Place up to 60 µl of the pool into a new 1.5 ml reaction tube and add 12 µl of 6x gel loading dye.

2.

Load the sample across a large well of a 1% (wt/vol) agarose gel containing 0.005% (vol/vol) ethidium bromide or SybrSafe.

3. Load a 100-bp ladder in a separate lane.

4.

Run the sample at 100 V for approximately 20 min.

5.

Cut a ~200 to 400 bp DNA smear from the gel.

6.

Extract the DNA using the Qiagen Gel Extraction Kit following the manufacturers manual.

7.

Elution in 25 µl of EB.

Fragment size of 200-400 bp is selected due to optimal read length for Illumina MiSeq

Instrument and may vary depending on sequencing platform used. Longer reads up to 600 bp are possible but the genomic library should then be prepared with longer fragment size.

21.

If possible positive and negative qPCR controls should be performed prior to deep sequencing (if known binding sites are available from ChIP data, for example). If there is no difference in amplification of the positive vs. negative control, this is an indication that the seq-DAP-seq protocol has failed. We have not had a case where the NGS gave useable data if the positive control qPCR verification failed. We recommend optimizing the binding buffer for seq-DAP-seq and repeating the protocol. We have found that amplicon sizes of 150 bp and 5-10 positive and 5-10 control regions is sufficient for library verification. 1 µl of library diluted to 50-100 µl is sufficient for qPCR verification. Enrichment of 4-fold or higher (CT cycle differences of at least 2 for positive and control qPCR) is a good indication that the experiment worked (Figure 9). 36. This setting will remove leading/trailing low quality or N bases (below quality 3), scan the read with a 4-base wide sliding window, cutting when the average quality per base drops below 20 and drop reads which are less than 36 bases long after these steps. Quality score of 20 means that there is a probability of 1/100 that the base is incorrect. Trimmomatic has other options to trim reads with alternative strategies.

24.

This step identifies genomics regions bound by the factor. In the alignment files, an excess of read is expected at bound regions. The difficulty lies in differentiating regions that collected lot of reads due to TF binding from the ones that accumulate lot of reads due to random/biased sequencing. Results of the peak calling should be checked using a genome browser to sense the quality of the peak calling. For this coverage files (BAM or bedgraph) peaks coordinates files (BED format) for each replicate can be loaded on a genome browser (IGB or IGV) for visual inspection. We have found this to be a simple way to qualitatively assess the peak calling step.

25.

Since the tools used in this step are memory demanding it is preferable to run it on a subset of the data, for example on the 600 most significant peaks in the MSPC output as shown in the example here. This constitutes the training set for motif detection and model construction. The remaining sequences will be used as a positive test set in the model validation step.

26.

Meme-chip may be run several times varying the minimum and maximum motif width so as to tune the motif detection. One can start with wide range of width (4 to 20 bp) with options -meme-minw set to 4 and -meme-maxw set to 20. Then visualize the result to spot the most significant motif(s) and re-run with a size range that correspond to this motif to trim off edges with low information content values (below ~0. [START_REF] Bailey | MEME SUITE: tools for motif discovery and searching[END_REF]). An important option is the -meme-pal option that will look for palindromic motif only.

27.

A ROC curve close to the diagonal with an AUC near 0.5 means that the model randomly collects sequences, whereas a curves well above the diagonal (AUC >> 0.5) means that the models correctly identifies most sequences from the positive set. Protein binding to magnetic beads. For a two-protein transcriptional complex, sequential pull downs are performed to isolate a single species (green and yellow circles). In the first pulldown step, the TF complex is immobilized on magnetic beads via an epitope tag (i.e. 3xFLAG tag with anti-FLAG magnetic beads) and eluted.

15 .

 15 Microtube 130 AFA (Covaris S220) or similar 16. Sterile 1.5 ml reaction tubes (Eppendorf or similar) 17. Sterile thin-walled PCR tubes 2.2 Kits 1. DNeasy Pro Kit (Qiagen) 2. End-It™ DNA repair kit (Lucigen) or similar 3. AMPure XP beads (Beckman Coulter) or similar 4. Miniprep kit (Qiagen) or similar 5. PCR clean-up kit (Qiagen) or similar 6. Qubit dsDNA HS quantification kit (ThermoFisher) 7. NEBNext Library Quant Kit for Illumina or similar 8. Gel Extraction Kit (Qiagen) or similar 9. TnT® SP6 High-Yield Wheat Germ Protein Expression System (Promega) or similar 2.3 Reagents and materials 1. Taq polymerase (NEB Taq DNA polymerase or similar) 2. Taq polymerase buffer (NEB ThermoPol or similar) 3. 25 mM MgCl2 (in ddH2O) 4. 1 mM dATP (stock solution in ddH2O) 5. ddH2O 6. T4 DNA ligase 7. 20 µg/µL Glycogen 8. 7.5 M NH4CH3CO2 (ammonium acetate) 9. Phenol:chloroform:isoamyl alcohol (25:24:1)

3. 5 . 6 .

 56 Fragment the DNA (Figure 2). Use the target fraction size of 200-300 bp on a temperature controlled ultrasonicator with the program target 200 bp. (I.e. peak power 175.0, Duty Factor 10.0%, Cycle/Burst 200, duration 180 s, sample temp 4-8 °C for DNA fragmentation.) 4. Transfer the sample to a clean sterile 1.5 ml tube and performed ethanol precipitation of the DNA. Start an ethanol precipitation by adding 2.5 µL glycogen (20 µg/µL), 62.5 µl (0.5x the sample volume) of 7.5 M NH4OAc and top up the sample with 475 µl (2.5x sample + ammonium acetate volume) 100% ethanol. Vortex to mix the sample. Note 2 Place the tubes at -20 °C overnight to precipitate the DNA from the sample or at -80 °C for at least 1 hour.

13 . 15 . 16 . 17 .

 13151617 If the pellet is still visible, place the sample at 37 °C in a water bath for 5 min to help resuspend the DNA and pipette up and down 40x. 14. End-repair the fragmented DNA (Figure 2).Transfer 34 µl of the DNA to a clean PCR tube and add 16 µl of End-It Repair Kit Master Mix (5 µL 10× End-It buffer, 5 µL 2.5 mM dNTP (each) mix, 5 µL 10 mM ATP and 1 µL End-It enzyme mix (T4 polymerase., T4 kinase and Klenow fragment)) and mix by gently tapping the side of the tube. Note 3 and Note 4 Centrifuge at room temperature for 5 s at 3,000x g to collect any sample adhering to the sides of the tube. Incubate at room temperature for 45 min and stop the reaction by heating to 70 °C for 10 min in a thermocycler. Transfer the sample to a clean 1.5 ml tube. Perform an ethanol precipitation by adding 1 µL glycogen (20µg/µL), 25 µl (0.5 x the sample volume) of 7.5 M NH4OAc and top up the sample with 189 µl 100% ethanol.

  benchtop until RT. Note 5 33. Transfer 30 µl sample (step 26) to a clean PCR tube to be ligated to the Y-shape adaptor. 34. Add 10 µl Ligation Master Mix (5 µL T4 DNA Ligase buffer and 5 µL 10x T4 DNA ligase buffer) and 10 µL annealed Y-adaptors (step 31) in EB and mix gently by pipetting up and down 5 times.

18 . 19 . 20 . 21 .

 18192021 Place the membrane into a blot box and remove the excess antibody by washing the membrane 5 times, 5 min in 1x PBS/0.1% Tween 20. Dilute the secondary antibody (1:10000) in 10 mL 5% milk powder/1x PBS/0.1% Tween 20 and incubate 1h. This solution can be reused and stored at 4 °C or -20 °C for long term storage. Wash the membrane 5 times, 5 min in 1x PBS/0.1% Tween 20. Incubate the membrane with ECL Plus Western Blotting Detection reagent for 1 min.

24 5. 6 . 7 . 8 .

 24678 fastqc -o out_dir TF_1.fastq.gz TF_2.fastq.gz NGmerge -a -1 TF_1.fastq.gz -2 TF_2.fastq.gz -o out_dir trimmomatic PE TF_1_noAdpaters.fastq.gz -2 TF_2_noAdpaters.fastq.gz TF_1_paired.fq.gz TF_1_unpaired.fq.gz TF_2_paired.fq.gz TF_2_unpaired.fq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:36 2. Map the reads to the reference genome sequence by using command bowtie2 of the Bowtie 2 package, for all replicates and the control independently(21). Prior to mapping, the reference genome is first indexed with the command bowtie2-build. Resulting alignment files are in SAM format that are parsed to discard poorly mapped reads (more than three mismatches, mapping quality below 30 and suboptimal alignments) and to remove PCR duplicates before being converted to BAM file using Samtools commands. Example command line input is given below. bowtie2-build genome.fasta genome_index bowtie2 -x bowtie genome_index -1 TF_1_paired.fq -2 TF_2_paired.fq -S out.samdovetail samtools view -Sh out.sam | grep -e "^@" -e 'XM:i:[012][^0-9]' | awk '$1~/@/ || $5>30 {print $0}' | grep -v "XS:i:" > out.filtered.sam 3. To prepare for peak calling, the BAM files are converted to BED files using bedtools bamtoBed command (23, 28, 29). 4. Call the peaks for each replicate separately by comparing their alignment file with that of the control input DNA using MACS2 (24). Resulting peaks files are then compared for consistency using MSPC (30) to keep ones that are observed in all (-c option set to 100%) or a majority (-c set to 50%) of replicates. macs2 callpeak -t replicate.bed -c control.bed -B -f BEDPE -n out_dir -g genome_length --call-summits -q 0.00001 mspc_mk -i (peak_rep1.narrowPeak peak_rep2.narrowPeak) -r Tec -w 1e-4 -s 1e-8 -c 100% -o outdir -d 1 Note Analyze the bound regions to identify an enriched DNA sequence motif that likely corresponds to a transcription factor binding sites (TFBS). Several tools exist, the most popular being meme-chip from the MEME-suite package (25). Note 25 and 26. meme-chip -oc out_dir -nmeme 600 -meme-maxsize 600000 -meme-minw 4 -mememaxw 8 -meme-nmotifs 3 -dreme-m 0 -noecho training_peaks.fas Once a motif has been established, evaluate its ability to discriminate between bound regions (positive data) and unbound genomic regions or randomly shuffled positive sequences (negative test set). The letter-probability matrix (that gives the frequency of each of the four nucleotides at each position in the motif) generated by meme-chip in step 4 (located in meme/meme_out/meme_mini.txt file) can be used (providing some format conversion) by several tools to search DNA sequence for this motif. In our pipeline, script meme2pfm.sh operates this conversion. Run the score.py script to scan the positive and negative set fasta sequences. Run the ROC.R script to compute a ROC curve and AUC to determine the quality of the model. Note 27. For large scale genomic DNA extraction, cost effective methods such as CTAB-based extraction (31), phenol:chloroform:isoamyl alcohol (25:24:1) or non-organic solvent extraction may be preferred (26, 32

22 . 23 .

 2223 With the example data, the pipeline runs within an hour on a standard computer (6 CPU and 15 Gb memory). In the FastQC output, one should pay attention to the quality score at each position along the reads. The quality score usually decreases towards the 3' end of the reads. After adaptor sequences have been trimmed off (NGmerge), one can trim low quality reads with Trimmomatic with option LEADING: 3 TRAILING: 3 SLIDINGWINDOW: 4:20 MINLEN:

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Schematic of DNA end repair reaction. 5' phosphate is shown in red.

Figure 3 .

 3 Figure 3. Schematic of A-tailing reaction. Added 3' A is labelled.

Figure 4 .

 4 Figure 4. Y-adaptor annealing and ligation reaction. The Y-adaptor has an overhanging T

Figure 5 .

 5 Figure 5. AmpDAP library barcoding. Barcodes are added using a unique indexing primer.

Figure 6 .

 6 Figure 6. Protein binding to magnetic beads. For a two-protein transcriptional complex,

Figure 7 .

 7 Figure 7. Second pull down to isolate the hetero TF complex. A second pull down using a

Figure 8 .

 8 Figure 8. DNA binding and amplification. The DAP or ampDAP input library is mixed with

Figure 9 .

 9 Figure 9. Example of qPCR amplification quality control, green: positive sets; red: negative

Figure 10 .

 10 Figure 10. Overview of the bioinformatic analysis pipeline for DAP-seq sequencing reads. The

Figure 11 .

 11 Figure 11. Example of IGB snapshot from ampDAP-seq experiments. In black are successfully

Figure 12 .

 12 Figure 12. Detailed example of barcoding and DNA amplification. Primers are color-coded

5 Figure 9

 59 Figure 7


Sylvain Schmitt
email: sylvain.m.schmitt@gmail.com

Guillaume Salzet

Fabian Jörg Fischer

Isabelle Maréchaux

Jerome Chave

| Guillaume Salzet

| Fabian

Jörg Fischer

Aaron Ellison

rcontroll: An R interface for the individual-based forest dynamics simulator TROLL

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

| INTRODUC TI ON

Forests are often conceptualised as a dynamic collection of individual trees responding in different and complex ways to local interactions, environmental change and natural and anthropogenic disturbances [START_REF] Levin | Ecosystems and the biosphere as complex adaptive systems[END_REF][START_REF] Shugart | A theory of forest dynamics. The ecological implications of forest succession models. A Theory of Forest Dynamics, The Ecological Implications of Forest Succession Models[END_REF]. Forests, as major stores of carbon, home to considerable biodiversity and providers of key ecosystem services, also are at the nexus of the environmental challenges of this century. They are threatened, being exposed to a wide range of impacts, including commodity-driven deforestation, forestry, shifting agriculture and wildfires [START_REF] Curtis | Classifying drivers of global forest loss[END_REF], and to the consequences of climate change [START_REF] Hubau | Asynchronous carbon sink saturation in african and Amazonian tropical forests[END_REF][START_REF] Saatchi | Persistent effects of a severe drought on Amazonian forest canopy[END_REF].

To understand how the processes that govern forest dynamics interact and how forests will respond to anthropogenic pressure, it is vital to produce biologically informed scenarios and for this models are important tools [START_REF] Fisher | Vegetation demographics in earth system models: A review of progress and priorities[END_REF][START_REF] Shugart | Gap models and their individualbased relatives in the assessment of the consequences of global change[END_REF].

Modern forest IBMs are usually written in high-performance languages such as C++ or Fortran, with platform-specific compilation challenges. Until now, forest IBMs have therefore been limited to a small number of expert users, typically core developers. This contrasts with other modelling approaches, such as species distribution models (SDMs) which also require complex, low-level code, but which have been widely adopted, in part due to the wide availability of tools and R packages for researchers and students [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF][START_REF] Schmitt | ssdm: An r package to predict distribution of species richness and composition based on stacked species distribution models[END_REF]. Conversely, few forest IBMs are currently available as R packages-plant being a rare example [START_REF] Falster | plant: A package for modelling forest trait ecology and evolution[END_REF]-making it difficult to reproduce research and to use IBMs in applied settings, both in research and beyond.

Forest models create links between a range of data to simulate forest dynamics on spatial and temporal scales inaccessible to empirical studies [START_REF] Marechaux | Tackling unresolved questions in forest ecology: The past and future role of simulation models[END_REF]. In the diverse family of forest dynamic models, individual-based models (IBMs) are unique in simulating the growth and demography of each individual tree, and in the relative ease with which they integrate field data to simulate ecological and physiological processes. Remarkably, from a historical perspective, IBMs were also the first type of numerical forest models, developed in the early 1960s for the purpose of optimising the yield of plantations [START_REF] Ek | FOREST: A computer model for simulating the growth and reproduction of mixed species forest stands[END_REF] and soon after to extend modelling capacities to any forest type [START_REF] Shugart | Forest Succession Models[END_REF].

Individual-based forest simulators have since proved useful for studying species-rich tropical forests [START_REF] Kazmierczak | A neutral vs. nonneutral parametrizations of a physiological forest gap model[END_REF][START_REF] Köhler | The effects of tree species grouping in tropical rainforest modelling: Simulations with the individual-based model FORMIND[END_REF][START_REF] Purves | Predictive models of forest dynamics[END_REF][START_REF] Urban | Spatial applications of gap models[END_REF] and for addressing a broad range of fundamental and applied questions [START_REF] Marechaux | Tackling unresolved questions in forest ecology: The past and future role of simulation models[END_REF], such as the effect of biodiversity on forest functioning through virtual experiments [START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF][START_REF] Morin | Temporal stability in forest productivity increases with tree diversity due to asynchrony in species dynamics[END_REF][START_REF] Sakschewski | Resilience of Amazon forests emerges from plant trait diversity[END_REF][START_REF] Schmitt | Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment[END_REF], or the quantification of forest carbon stocks and fluxes [START_REF] Fischer | Simulating carbon stocks and fluxes of an African tropical montane Forest with an individualbased forest model[END_REF][START_REF] Fyllas | Solar radiation and functional traits explain the decline of forest primary productivity along a tropical elevation gradient[END_REF]. By leveraging remote-sensing products, IBMs have expanded both their spatial scale of application and range of parameterisation [START_REF] Fischer | Improving plant allometry by fusing forest models and remote sensing[END_REF][START_REF] Rödig | Estimating the carbon fluxes of forests with an individual-based forest model[END_REF][START_REF] Rödig | The importance of forest structure for carbon fluxes of the Amazon rainforest[END_REF][START_REF] Shugart | Computer and remote-sensing infrastructure to enhance large-scale testing of individual-based forest models[END_REF], and remote-sensing products have been benchmarked with the aid of virtual remotely sensed scenes [START_REF] Knapp | Linking lidar and forest modeling to assess biomass estimation across scales and disturbance states[END_REF].

To fill this gap, we present an R package named rcontroll, which integrates the TROLL individual-based model into R to simulate forest ecosystem processes and species dynamics. The TROLL model was originally developed for tropical forests, with the aim of simulating ecological succession from bare ground and the spatial patterns emerging from tree falls and gap dynamics [START_REF] Chave | Study of structural, successional and spatial patterns in tropical rain forests using TROLL, a spatially explicit forest model[END_REF].

Species-specific parameters are based on plant functional traits allowing the joint simulation of carbon and tree diversity [START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF]. TROLL has been applied to the study of tropical forest resilience [START_REF] Schmitt | Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment[END_REF], tree allometry and how it can be inferred from remote sensing [START_REF] Fischer | Improving plant allometry by fusing forest models and remote sensing[END_REF], and forest responses to wind-throw [START_REF] Rau | Transferability of an individual-and trait-based forest dynamics model: A test case across the tropics[END_REF].

| MATERIAL S AND ME THODS

| TROLL model

TROLL is a spatially explicit and individual-based forest model. We provide a short overview in the following section and a schematic diagram illustrating the functioning of TROLL (Figure 1a), but refer the reader to the original TROLL publications for further details [START_REF] Chave | Study of structural, successional and spatial patterns in tropical rain forests using TROLL, a spatially explicit forest model[END_REF][START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF], Appendix S5, Figure S1 therein for another schematic diagram).

Tree growth and competition in TROLL happen in a 3D space, which is discretised into voxels of one cubic metre. Within this space, trees establish, grow and interact with each other, both vertically and horizontally, with at most one tree per square metre. A simplified tree geometry is adopted: the tree trunk is represented by a cylinder, characterised by its diameter at breast height (DBH) and total height; its crown has a symmetrical geometric shape such as a cylinder, cone or elliptical surface, which is characterised by its radius and depth and adjustable via shape parameters. Each tree is attributed to a botanical species with corresponding species-specific

K E Y W O R D S
forest simulator, individual-based model, R package, spatially explicit, TROLL plant functional traits related to ecophysiological processes, including photosynthesis, carbon allocation and mortality, as well as tree architecture. Specifically, these traits include leaf mass per area LMA (g m -2 ), leaf nitrogen N (mg g -1 ), leaf phosphorus P (mg g -1 ) and wood density WD (g cm -3 ), as well as two parameters hmax and ah which describe allometric growth between trunk diameter dbh (m) and tree height h (m) according to a Michaelis-Menten equation: h = hmax × dbh (dbh + ah) . Horizontal growth of crown radius cr is modelled according to a power law cr = a × dbh b , with parameters a and b specified at the community level. Crown depth cd is modelled as a fixed fraction of tree height.

New trees are recruited into the community from a seed bank that is filled by an external seed rain. If the simulated stand contains mature trees, these will also contribute to the seed rain. Each model iteration, a random seed is selected from the seedbank and, provided that the environmental conditions are favourable for tree growth (the potential tree's light compensation point is exceeded), turned into a seedling of 1 cm DBH. Otherwise, the grid cell remains empty. For each tree thus established, TROLL simulates carbon assimilation, carbon allocation, trunk and crown growth, seed dispersal and tree death (Figure 1a). TROLL uses the C3 photosynthesis model of [START_REF] Farquhar | A biochemical model of photosynthetic CO 2 assimilation in leaves of C3 species[END_REF], a function of local temperature, light irradiance, vapour pressure deficit and atmospheric CO2 concentration. Model parameters, including the maximum rate of carboxylation, Vcmax and the maximum electron transport capacity, Jmax, are inferred from the species-specific plant functional traits N, P and LMA [START_REF] Domingues | Co-limitation of photosynthetic capacity by nitrogen and phosphorus in West Africa woodlands[END_REF].

Solar irradiance is calculated within each voxel as the fraction of irradiance transmitted by the voxels above following a Beer-Lambert law, while temperature and vapour pressure deficit are varying throughout the canopy according to empirical functions [START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF]. Carbon assimilation is calculated half-hourly for a representative day per month. Stomatal conductance responds to air vapour pressure deficit and is modelled according to [START_REF] Medlyn | Reconciling the optimal and empirical approaches to modelling stomatal conductance[END_REF]. Net assimilable carbon is calculated as gross assimilable carbon through photosynthesis minus respiration [START_REF] Atkin | Global variability in leaf respiration in relation to climate, plant functional types and leaf traits[END_REF] and is allocated to trunk and root growth, leaf production and reproduction based on empirical relationships and the biomass density of each component (LMA, WD).

If leaf area reaches its optimum, no further leaves are allocated, and the excess carbon is allocated to growth or a storage pool [START_REF] Fischer | Inferring the structure and dynamics of tropical rain forests with individual-based forest growth models (Doctoral dissertation)[END_REF]. Tree growth results in a change in tree height and crown size following DBH-driven allometries, which in turn influence the light environment at the next time step. In TROLL, tree mortality results from three distinct processes: (i) stochastic mortality, modelled as a function of a maximum background mortality rate and a linearly decreasing relationship with species-specific wood density; (ii) carbon starvation, if net assimilated carbon is negative over a consecutive period exceeding the leaf lifespan; and (iii) stochastic treefall events, assumed to depend on a tree height threshold.

The combination of these processes at the individual level leads to realistic community-level dynamics, including successional patterns in species composition, self-thinning and gradual saturation of stand biomass on decadal to centennial timescales even with minimal calibration [START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF][START_REF] Rau | Transferability of an individual-and trait-based forest dynamics model: A test case across the tropics[END_REF].

The minimal set of input files required for a TROLL run include (i) climate data for the focal location, (ii) functional traits for the list of species at the focal location and (iii) global parameters, that is parameters that do not depend on species identity. All these input files and simulation parameters can be easily configured and run using rcontroll. a typical personal computer). The rcontroll R package is a wrapper of TROLL that comprises the model and facilitates its access to users. rcontroll currently calls version 3.1.7 of TROLL using the Rcpp package [START_REF] Eddelbuettel | Rcpp: Seamless R and C++ integration[END_REF]. rcontroll includes functions that generate inputs for simulations and run simulations. Finally, it is possible to analyse the TROLL outputs through tables, figures and maps taking advantage of other R visualisation packages. rcontroll also offers the possibility to generate a virtual LIDAR point cloud that corresponds to a snapshot of the simulated forest.

| rcontroll workflow and usage

| Construction and manipulation of input files

As stated above, three types of input data are needed for a typi- and their uncertainties reported. Therefore, rcontroll users only need to input the species-specific trait data to run TROLL simulations, irrespective of the site. TROLL was originally developed for tropical and subtropical forests, so certain assumptions must be critically examined when applying it outside the tropics. The input files can be used to start a TROLL simulation run within the rcontroll environment (see below) or saved so that the TROLL simulation can be started as a command line tool.

| Simulations

The default option is to run a TROLL simulation using the troll function of the rcontroll package. The output is stored in a trollsim R class. For multiple runs, users can rely on the stack function, and the output is stored in the trollstack class. Both trollsim and trollstack values can be accessed using object attributes in the form of simple R objects (with @ in R). They consist of eight simulation attributes: (1) name, (2) path to saved files, (3) parameters, (4) inputs, (5) log, (6) initial and final state, (7) ecosystem output metrics and (8) species output metrics. The initial and final states are represented by a table with the spatial position, size and other relevant traits of all trees at the start and end of the simulation.

The ecosystem and species metrics are summaries of ecosystem processes and states, such as net primary production and aboveground biomass, and they are documented at species level and aggregated over the entire stand. Simulations can be saved using a user-defined path when run and later loaded as a simple simulation (load_output function) or a stack of simulations (load_stack function).

| Simulated airborne LIDAR scanning option

TROLL also has the capacity of generating point clouds from virtual aerial LIDAR scanning of simulated forest scenes. Within each cubic metre voxel of the simulated stand, points are generated probabilistically, with the probability depending both on the amount of light reaching the particular voxel and the amount of leaf matter intercepting light within the voxel. Extinction and interception of light are based on the Beer-Lambert law, but an effective extinction factor is used to account for differences between the near-infrared and visible light. The definition of the LIDAR parameters (generate_lidar function) is optional but allows the user to add a virtual aerial LIDAR scan for a time step of the TROLL simulation. When this option is enabled, the cloud of points from simulated aerial LIDAR scans is stored as LAS using the R package lidR [START_REF] Roussel | lidR: An r package for analysis of airborne laser scanning (ALS) data[END_REF] as a ninth attribute of the trollsim and trollstack objects.

| Manipulation of simulation outputs

rcontroll includes functions to manipulate simulation outputs. Simulation outputs can be retrieved directly from the trollsim or trollstack objects and summarised or plotted in the R environment with the print, summary and autoplot functions. The get_chm function allows users to retrieve canopy height models from aerial LIDAR point clouds (Figure 2). In addition, a rcontroll function is available to visualise TROLL simulations as an animated figure (autogif function, Figure 1b).

| APPLI C ATI ON S

| Sensitivity of TROLL parameters

As a first example, we performed a sensitivity analysis of 12 parameters (Table 1), the results of which could be used for effective model calibration. Calibration of TROLL consists of adjusting the parameter values so that the model behaves as closely as possible to the empirical data. We simulated an undisturbed forest plot located in French Guiana based on the inventory data from the Paracou research station (Gourlet-Fleury et al., 2004). We retained only the species provided as defaults in rcontroll. We sampled the parameter space using the Latin hypercube sampling (LHS) method, using the a priori parameter ranges (Table 1). LHS was modified to account for the expected correlation between parameters (Fischer, 2019) using the Huntington-Lyrintzis algorithm [START_REF] Huntington | Improvements to and limitations of Latin hypercube sampling[END_REF] implemented in the pse R package [START_REF] Chalom | pse: Parameter space exploration with Latin Hypercubes[END_REF]. Parameter combinations were replicated 10 times to assess stochastic variation 2041210x, 0, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14215 by CIRAD -DGDRS -DIST, Wiley Online Library on [02/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License (500 combinations × 10 replications leading to 5000 simulations).

Simulations were run from an initial bare ground for 600 years using the stack function in rcontroll, and they were subsequently summarised into four ecosystem metrics averaged over the last 100 years: the number of trees ≥10 cm DBH, the number of trees ≥30 cm DBH, above-ground biomass and gross primary production. To convert the discrete response into a continuous response function, covering the entire space of the 12 parameters, we derived a Gaussian process model for each ecosystem metric from the repeated samples using the hetGP R package [START_REF] Binois | HetGP: Heteroskedastic gaussian process modelling and sequential design in R[END_REF]. This surrogate method allows us to generate any number of additional simulations for any parameter combination and yields robust mean estimates of ecosystem metrics with TROLL parameters independent from the noise generated from replications.

Finally, a sensitivity analysis was carried out using the sensitivity R package to determine the relative influence of the 12 parameters on the mean estimates of the four ecosystem metrics.

We used Morris analysis, a qualitative approach to eliminate noninfluential parameters [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF]. Morris analysis is based on a one-step-at-a-time sampling scheme (i.e. at each run, only one input parameter is given a new value), which infers the effects of each parameter on the ecosystem metrics. The sensitivity factors can be compared globally, and the nonlinearity and interaction of the model is described qualitatively. The Morris analysis highlighted TROLL sensitivity to the parameters CR a of the allometric relationship between crown radius and trunk diameter, the fraction of net primary productivity allocated to the canopy f canopy and the apparent quantum yield for carbon fixation (Figure 3; [START_REF] Fischer | Improving plant allometry by fusing forest models and remote sensing[END_REF][START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF]. The Morris analysis also revealed the importance of nonlinear effects or interactions of the tested parameters on the four ecosystem metrics, which can be explored further quantitatively with variance-based analyses (e.g.

Sobol indices, see Supplementary Information SI2 and SI3). Thus, use case shows how rcontroll helps handle complex and simulationintensive experiments, fostering our understanding of the TROLL model behaviour, and hence of model predictions and uncertainties, as well as its transferability.

| Functional diversity and forest resilience

As a second example of rcontroll application, we revisited the question: does functional diversity improve the tropical forest resilience to disturbance, which has been previously addressed in a virtual experiment using TROLL [START_REF] Schmitt | Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment[END_REF]. We replicated a scaled-down version of this experiment to demonstrate how the rcontroll interface facilitates the implementation of a simulated experiment such as this one (Supplementary Information SI4 and SI5). For three levels of species richness (5, 10 and 20), we simulated 10 random assemblages of communities whose composition from the default TROLLv3_species dataset included in rcontroll (resulting in 30 simulations). This was done using the stack function. Simulations were initialised from bare ground and run for 600 years to reach a mature forest state. We then sampled half of the individual trees and removed the rest to simulate a random disturbance (simplified from [START_REF] Schmitt | Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment[END_REF].

After the disturbance stage, forest dynamics was simulated for another 600 years, again using the stack function. A control simulation was also conducted (with no disturbance). Taxonomic and functional diversity were assessed before and after disturbance using indices of species richness and functional diversity (FDiv, the volume of functional space occupied by the community; FEve, the regularity of the distribution of abundance within this volume; [START_REF] Villéger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF]. We computed changes in above-ground biomass, basal area, total number of stems, number of stems larger than 10 cm and 30 cm DBH, from which we calculated a resilience index (see [START_REF] Schmitt | Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment[END_REF]. The rcontroll package reproduces the analyses of Schmitt et al. ( 2020) with only about 80 lines of R code, increasing transparency and replicability (Supplementary Information SI4 and SI5; Figure 4).

F I G U R E 3 Sensitivity of TROLL to parameters using Morris analysis. Twelve parameters of TROLL were tested, with abbreviations found in Table 1.

The level of grey of the label indicates parameters related to the same processes (from dark grey to white, respectively: carbon allocation, crown allometry, leaf ecophysiology, mortality and reproduction). The x-axis represents the absolute mean of elementary effect μ*, and the y-axis represents the variance of elementary effect σ. A high μ* indicates a factor with a heavy overall influence on model outputs; a high σ indicates heterogeneity of the sensitivity across the parameter space, which could indicate nonlinearities or interactions with other parameters. A high μ* also tends to produce high σ, so σ should be interpreted relative to the ratio σ/μ*. The three lines dividing the space provide the type of relation of the focal parameter with other studied parameters.

| DISCUSS ION

The R package rcontroll is a wrapper for the TROLL individual-based forest simulator. The main goal is to increase the accessibility of the C++-coded model and to make it available to a wider community of researchers and students by facilitating the setting up and running of simulations. rcontroll also provides the community with tools to quickly analyse the outputs of TROLL in a standardised way. rcontroll is intended to facilitate research on complex issues in forest ecology and provide new use cases. This new wrapper of TROLL will also help increase transparency and reproducibility in simulation experiments.

We have provided two concrete use cases. In the first, rcontroll is used to perform a simulation-intensive sensitivity analysis of 12 parameters of the TROLL model, providing both a better understanding of the model's behaviour and uncertainties and opportunities for efficient model transferability with calibration across sites. In the second, we showed how rcontroll facilitates the replication of a previously published analysis [START_REF] Schmitt | Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment[END_REF] by simplifying the code while increasing transparency and reproducibility. The reproducibility of studies can be further increased by the use of other R packages, such as tidyverse [START_REF] Wickham | Welcome to the tidyverse[END_REF] essential for processing simulation results, analyses and graphs. rcontroll can be used to address a wide range of topical fundamental and applied ecological questions by means of modelling, such as the assessment of management scenarios [START_REF] Marechaux | Tackling unresolved questions in forest ecology: The past and future role of simulation models[END_REF].

The rcontroll package is tightly linked to TROLL, so a specific discussion of the latter is in order here. As with any model, TROLL is built upon assumptions and care should be taken to interpret its results carefully. We strongly recommend that users apply the model within its current range of application and discuss their results by taking into account the limitations of TROLL discussed in [START_REF] Maréchaux | An individual-based forest model to jointly simulate carbon and tree diversity in Amazonia: Description and applications[END_REF] and [START_REF] Rau | Transferability of an individual-and trait-based forest dynamics model: A test case across the tropics[END_REF]. TROLL has been developed and used for natural tropical forests. Although it is intended to be widely applicable and has been successfully extended to the subtropics [START_REF] Rau | Transferability of an individual-and trait-based forest dynamics model: A test case across the tropics[END_REF], certain assumptions need to be critically examined when applying it outside the trop- [START_REF] Schmitt | ssdm: An r package to predict distribution of species richness and composition based on stacked species distribution models[END_REF], to increase its user-friendliness, as well as the Canopy Constructor model [START_REF] Fischer | A simulation method to infer tree allometry and forest structure from airborne laser scanning and forest inventories[END_REF], which could be used to run simulations from empirical forest inventories. Finally, we aim to increase the interoperability of rcontroll with the rest of the R ecosystem by providing formatting functions that would make it easy to import forest inventories from existing repositories or models and to easily export simulation outputs to other tools designed for their analysis.

CO N FLI C T O F I NTER E S T S TATEM ENT

No conflict of interest.

F I G U R E 4 Trajectories of simulated above-ground biomass (AGB, Mg/ha) over time (year) for one community (with 20 species). The colours of the lines represent the types of trajectories, with the trajectory of the predisturbance simulation from bare soil in black, the trajectory after the disturbance event in red and the trajectory of the undisturbed control simulation from the mature state before disturbance in green.

  rcontroll provides user-friendly functions for setting up and analysing TROLL simulations. Users specify floristic community composition, ecological parameters and climatic conditions, and default values are available where possible. rcontroll provides default data for demonstration runs, a climate data generator from a global dataset, visualisation tools as well as a generator of virtual LIDAR point clouds. Here, we demonstrate the flexibility of rcontroll by testing the sensitivity of TROLL parameters and replicating analyses from a previously published study using reproducible code examples.

  TROLL is coded in C++, and it typically simulates hundreds of thousands of individuals over hundreds of years in several hectares in minutes (a simulation over 4 ha and 600 years lasts ca. 1 min on F I G U R E 1 TROLL forest simulator. (a) TROLL is a spatially explicit and individual-based forest model, which simulates carbon assimilation and allocation (light green), trunk and crown growth (light green), seed dispersal (brown) and tree deaths (red) and recruitments (dark green). (b) Output from a TROLL simulation using the autogif function in the rcontroll package. The image shows a vertical cut in the forest structure along the x-axis (in metres) with individual tree height (metres) on the y-axis. The tree colours indicate the identity of the species and can be changed using the ggplot2 grammar. The figure shows the forest structure after 625 years of a successional trajectory starting from bare ground, which can be seen dynamically over 200 years at https:// raw.githu buser conte nt.com/sylva insch mitt/rcont roll/main/inst/ figur es/troll.gif.

  cal TROLL simulation: (i) climate data, (ii) plant functional traits and (iii) global model parameters. Presimulation functions include global parameters definition (generate_parameters function) and climate data generation (generate_climate function). rcontroll also includes default data for species and climate inputs for a typical French Guiana rainforest site. The purpose of the generate_climate function with the help of the corresponding vignette is to create TROLL climate inputs from ERA5-Land (Muñoz-Sabater et al., 2021), a global climatic reanalysis dataset that is freely available (see Supplementary Information SI1). The ERA5-Land climate reanalysis is available at 9 km spatial resolution and hourly temporal resolution since 1950, and daily or monthly means are available

F

  Comparison of clouds of points obtained from a virtual airborne LIDAR scan of a forest scene simulated with TROLL (a) with a real airborne LIDAR scan from Nouragues station (b). The horizontal axes represent the x-axis and y-axis (in metres), and the vertical axis represents height (in metres). The thermal colour scale indicates the height of the points in the cloud, from 0 m in dark blue to 50 m in red. The LIDAR simulator is not a ray tracing algorithm, but an approximation of the real laser scans, and the comparison shows among other things simplified geometric shapes of the trees and no return loss due to oblique angles.TA B L E 1Parameters used for the sensitivity analysis. Parameters have been classified by associated processes. The a priori interval used for sensitivity analysis is given for each parameter with its source. Carbon allocation f wood : fraction of NPP allocated to wood growth (%) [0.01, 1.00] f wood + f canopy + f leaf = 1 f canopy : fraction of NPP allocated to canopy (%) [0.01, 1.00] Mortality m 0 : maximal basal mortality rate (events. year -1 ) [0.01, 0.05] Fischer (2019) wsg lim : wood specific gravity limiting mortality factor (dimensionless) [1.00, 1.20] Maréchaux and Chave (2017) V c : treefall stochastic threshold (dimensionless) [0.01, 0.15] Rau et al. (2022) Reproduction seedrain: Total number of reproduction opportunities coming from outside [100, 100,000] Maréchaux and Chave (2017) nbs 0 : local seed dispersed par mature tree [1, 1000] Crown allometry CR a : intercept of Log-Log regression to infer crown radius from DBH[1.5, 3] Fischer (2019) with a correlation factor: = 0.65 CR b : slope of Log-Log regression to infer crown radius from DBH [0.4, 0.8]

  ics. Some assumptions and limitations may change with TROLL versions (3.1.7 at the time of this manuscript), so users are encouraged to check future TROLL upgrades. Feedback on the TROLL model is welcome by the development team (JC, IM and FJF). As TROLL is being further developed, we plan to upgrade rcontroll concurrently with TROLL and to swiftly include any major new development. rcontroll includes a helper function that keeps track of the corresponding TROLL version (TROLL.version function). In future, we also plan to add a global user interface in shiny for rcontroll (e.g.

The

  rcontroll package is free and open source (version 0.1.0 with GPL3). It is available on the CRAN repository https://cran.r-proje ct.org/web/packa ges/rcont roll/index.html and can be installed in the R environment using the install.packages ('rcontroll') command. The project is hosted on GitHub (https://github.com/sylva insch mitt/rcont roll), which allows future users to openly contribute to the project. AUTH O R CO NTR I B UTI O N S Sylvain Schmitt and Guillaume Salzet conceived the ideas; Sylvain Schmitt led the package development; Sylvain Schmitt and Guillaume Salzet designed the examples and analysed model outputs; and Sylvain Schmitt led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication. ACK N O WLE D G E M ENTS This work was supported by the 'Amazonian Landscapes in Transition' project and by the 'Investissement d'Avenir' grants (CEBA, ref. ANR-10-LABX-25-01; TULIP, ref. ANR-10-LABX-0041; ANAEE-France: ANR-11-INBS-0001), all managed by the Agence Nationale de la Recherche, by CNES (BIOMASS-VALO TOSCA Project), and the LECOS Project funded by the INRAE-CNRS joint initiative.

2041210x, 0, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14215 by CIRAD -DGDRS -DIST, Wiley Online Library on [02/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Methods in Ecology and Evoluঞon SCHMITT et al.

Funding information

Grant/Award Number: ANAEE-France: ANR-11-INBS-0001, TULIP: ref. ANR-10-LABX-0041 and CEBA: ref. ANR-10-LABX-25-01; CNES (BIOMASS-VALO TOSCA Project); LECOS Project

DATA AVA I L A B I L I T Y S TAT E M E N T

All data are available on Zenodo [START_REF] Schmitt | sylvainschmitt/rcontroll: v0.1.0[END_REF].

O RCI D

PEER R E V I E W

The peer review history for this article is available at https:// w w w.w e b o f s c i e n c e . c o m /a p i /g a t e w a y/ w o s /p e e r-r e v i e w/10.1111/2041-210X.14215.

S U PP O RTI N G I N FO R M ATI O N

Additional supporting information can be found online in the Supporting Information section at the end of this article.

Supporting Information S1-S6: Supporting information SI1

describes the climate generation function. Supporting information SI2 further detail the sensitivity analysis of TROLL parameters with reproducible code in SI3 and supporting functions in SI6. Supporting information SI4 further detail the analysis of functional diversity on forest resilience with reproducible code in SI5.


Reynald Affeldt

Yves Bertot

Cyril Cohen

Pierre Roux

Kazuhiko Sakaguchi

Enrico Tassi

Porting Coq

Porting Coq Scripts to the Mathematical Components Library Version 2

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

1. MathComp users who have not been using the Canonical command should not see much difference compared to the past version upgrades of MathComp. Some identifiers that are now useless have been removed but this is documented in the changelog. For example bool_eqType might need to be replaced by bool : eqType or just bool. Also, the behavior of some rewritings might have changed, requiring explicit patterns. Typically, it might happen that rewriting with associativity lemmas requires the user to indicate whether it is supposed to happen on the left or the right-hand side of an equivalence relation, so that rewrite addrA might need to be rewritten rewrite [in LHS]addrA for example. See Sect. 4.4 for a concrete example. For such users, reading this tutorial through the end might not be necessary.

2.

The target readers are primarily MathComp users who have been instantiating structures using the Canonical command.

3. As for the few users who have been developing their own hierarchies of structures, this tutorial might be of little help and they rather need to refer to:

the original paper for an extensive introduction to HB commands [START_REF] Cohen | Hierarchy Builder: Algebraic hierarchies made easy in Coq with Elpi (system description)[END_REF],

the HB development for documentation and examples [START_REF]High level commands to declare a hierarchy based on packed classes[END_REF] (start with the README), various papers for more applications [START_REF] Affeldt | Porting the Mathematical Components library to Hierarchy Builder[END_REF] [2, Sect. 3] [3, Sect. 4],

already ported developments such as odd-order, multinomials, etc.

For the sake of concreteness, we illustrate the port of CompDecModal [START_REF] Doczkal | Completeness and decidability of modal logic calculi[END_REF] in Sect. [START_REF] Cohen | Hierarchy Builder: Algebraic hierarchies made easy in Coq with Elpi (system description)[END_REF]. Before that we review the basics of HB in Sect. 2 and review the documentation tools available for porting in Sect. 3.

Quick Reminder about the HB Vocabulary

The goal of this section is to briefly explain the three main commands introduced by HB: HB.mixin, HB.structure, and HB.instance. The knowledgeable reader can safely skip this section.

Let us consider the most basic scenario in generic terms. Here is the pattern to declare a structure Struct that sits at the bottom of a hierarchy. The interface of the structure goes into a mixin: HB.mixin Record isStruct params carrier := { ... properties about the carrier ... }

The structure itself is declared like a sigma-type:

#[short(type=structType)] HB.structure Definition Struct := {carrier of isStruct carrier} Note that HB is using Coq attributes to declare the type corresponding to a structure.

Here is the pattern to declare a new structure NewStruct that extends the existing structure Struct; note the of syntax.

HB.mixin Record NewStruct_from_Struct params carrier

of Struct params carrier := { ... more properties about the carrier ... }

In the case of the extended structure, the sigma-type makes appear the dependency to the parent structure; note the & syntax. This process results in the creation of the types structType and newStructType such that elements of the latter are also understood to be elements of the former.

Finally, the declaration of a mixin Struct is accompanied by the creation of a constructor Struct.Build which is used to instantiate a structure using the command: HB.instance Definition _ := Struct.Build params.

The command HB.instance should trigger the printing of several lines of information output such as module_type__canonical__struct_Struct is defined

The absence of this output often indicates failure of the HB.instance command.

Tools to Port MathComp Applications

Documentation

The following pieces of documentation are useful during the process of porting a MathComp application to MathComp 2:

The changelog is the primary source of information. See CHANGELOG.md [START_REF]Mathematical Components[END_REF].

Additionally, structures are documented in the headers of Coq scripts according to the following format: See for example the eqType structure defined in the file ssreflect/eqtype.v. See this wiki entry for more information about the documentation of scripts.

Optionally, the user can double-check the naming of identifiers and lemmas with the naming conventions explained in CONTRIBUTING.md [START_REF]Mathematical Components[END_REF].

HB Commands Useful to Explore an Existing Hierarchy

Besides the changelog and the headers of Coq scripts, the user can use HB commands to explore a hierarchy of mathematical structures.

Information about Structures with HB.about

Basic information about structures can be obtained via the command HB.about as in:

> HB.about eqType. HB: eqType is a structure (from "./ssreflect/eqtype.v", line 137) HB: eqType characterizing operations and axioms are:

-eqP -eq_op HB: eqtype.Equality is a factory for the following mixins:

-hasDecEq (* new, not from inheritance *) HB: eqtype.Equality inherits from: HB: eqtype.Equality is inherited by:

-SubEquality -choice.Choice ...

(The output message refers to a factory: this is a generalization of mixin.)

Graph of an HB Hierarchy It is also possible to explore a HB hierarchy using the command HB.graph. Inside a Coq file: HB.graph "hierarchy.dot".

From a terminal:

For example, Fig. 1 displays the immediate vicinity of eqType. -hasDecEq tells us that eqType instances can be built with hasDecEq.Build. (Note that by default HB.howto may not return all the available factories; it might be necessary to increase the depth search using a natural number as in HB.howto xyzType 5.)

To learn which parameters a xyz.Build constructor is expecting, one can use the HB.about command:

> HB.about hasDecEq.Build. HB: hasDecEq.Build is a factory constructor (from "./ssreflect/eqtype.v", line 135) HB: hasDecEq.Build requires its subject to be already equipped with: HB: hasDecEq.Build provides the following mixins:

-hasDecEq HB: arguments: hasDecEq.Build T [eq_op] eqP -T : Type -eq_op : rel T -eqP : Equality.axiom eq_op

The message indicates that hasDecEq.Build is expecting a type T, a predicate eq_op : rel T (implicit argument, as indicated by the square brackets) and a proof of Equality.axiom eq_op. One can thus instantiate an eqType on some type T with]

HB.instance Definition _ := hasDecEq.Build T proof_of_Equality_axiom.

or HB.instance Definition _ := @hasDecEq.Build T eq_op proof_of_Equality_axiom.

which should output a few lines among which (recall that the absence of this output often indicate an instantiation problem)1 :

module_T__canonical__eqtype_Equality is defined Discover Aliases and Feather Factories In addition to the structures and constructors listed by HB.about, the library defines some aliases (a.k.a. feather factories). These aliases are documented in the header comments. For instance, an eqType instance on some type T can be derived from some T' already equipped with an eqType structure, given a function f : T -> T' and a proof injf : injective f: HB.instance Definition _ := Equality.copy T (inj_type injf).

See eqType.v for inj_type.

Information about Instances with HB.about

Instances a type is already equipped with can be listed with HB.about, for instance:

> HB.about bool. HB: bool is canonically equipped with structures: -Order.BDistrLattice Order.BLattice Order.BPOrder (from "./ssreflect/order.v", line 6064) ... lists all the structures bool is already equipped with.

Porting a MathComp Development to MathComp 2

The basic strategy to port an existing MathComp development to MathComp 2 is to (1) install Math-Comp 2, (2) compile the existing Coq scripts, and (3) fix the errors one after the other. For the sake of concreteness, we explain the port of CompDecModal [START_REF] Doczkal | Completeness and decidability of modal logic calculi[END_REF]. This is a development with a moderate use of MathComp whose port involves fixing the instantiation of basic structures that most developments using MathComp are likely to use.

In the following, the offending commands appear in a gray areas command incompatible with MathComp 2

and their fixes are singly famed:

MathComp 2 fix for the command above

Import the HB Library

First thing first, any Coq file using HB must start with:

From HB Require Import structures.

Instantiation of Structures with MathComp 2

From the viewpoint of the MathComp user, the main change is the way mathematical structures are now instantiated. The fix can again be inferred from the changelog: > HB.instance Definition _ := [Choice of fset_type by <:]. HB_unnamed_factory_11 is defined choice_Choice__to__choice_hasChoice is defined HB_unnamed_mixin_14 is defined fset_fset_type__canonical__choice_Choice is defined fset_fset_type__canonical__choice_SubChoice is defined > HB.instance Definition _ (T : countType) := [Countable of fset_type T by <:]. T is declared HB_unnamed_factory_30 is defined choice_Countable__to__choice_hasChoice is defined choice_Countable__to__eqtype_hasDecEq is defined choice_Countable__to__choice_Choice_isCountable is defined HB_unnamed_mixin_34 is defined fset_fset_type__canonical__choice_Countable is defined fset_fset_type__canonical__choice_SubCountable is defined Note that, although HB does provide an #[hnf] attribute, it should not be necessary in general.

Finally, the last Canonical command causes a deprecation warning that needs to be addressed: Warning: Notation "[ subCountType of _ ]" is deprecated since mathcomp 2.0.0. Use SubCountable.clone instead.

[deprecated-notation,deprecated] fset_subCountType is defined

In fact, going back one step, it can be observed by the output the HB.instance commands that the instantiation of the Countable structure already triggered the instantiation of the SubCountable structure, rendering the last Canonical command harmful. It therefore needs to be deleted.

To sum up, here follows the complete fix:

Section FinSets.

Variable T : choiceType. In fact, we can go one step further. Instead of instantiating the Equality structure at line 5 and then the Choice structure at line 7, one can start by instantiating the Choice structure and get the Equality structure automatically, so that a shorter fix would be: This small example already illustrates the advantage of porting to MathComp 2.

Finitely Iterated Operators

The next series of compilation errors that occur when porting CompDecModal is about finitely iterated operators, which are likely to be used by most MathComp developments:

Canonical Structure fsetU_law (T : choiceType) := Monoid.Law (@fsetUA T) (@fset0U T) (@fsetU0 T). Canonical Structure fsetU_comlaw (T : choiceType) := Monoid.ComLaw (@fsetUC T).

The first compilation error indicates a change in the signature of a constructor:

> Canonical Structure fsetU_law (T : choiceType) := > Monoid.Law (@fsetUA T) (@fset0U T) (@fsetU0 T). Error:

In environment T : choiceType The term "fsetUA (T:=T)" has type "forall X Y Z : {fset T}, X `|`(Y `|`Z) = X `|`Y `|`Z" while it is expected to have type "Type".

We can use HB.about to discover the relevant structure (see Sect. [START_REF] Affeldt | Semantics of probabilistic programs using s-finite kernels in coq[END_REF] We have now enough information to fix the compilation error:

HB.instance Definition _ (T : choiceType) := Monoid.isLaw.Build {fset T} fset0 fsetU (@fsetUA T) (@fset0U T) (@fsetU0 T).

Note that we have made explicit the key (the operator fsetU) in the instantiation. Though not strictly necessary, this is good practice to document it on this occasion.

Next compilation error:

> Canonical Structure fsetU_comlaw (T : choiceType) := > Monoid.ComLaw (@fsetUC T). Error: In environment T : choiceType The term "fsetUC (T:=T)" has type "forall X Y : {fset T}, X `|`Y = Y `|`X" while it is expected to have type "Type". Note that we have added a type information to make the key (here, form) explicit as recommended when performing instantiations. This fix is sufficient because it generates the Choice instance along with the Countable instance. This kind of slow down is however rather exceptional.

More Instantiations of

There are a few more compilation errors but they are similar to the ones we explained so far.

Conclusion

This document illustrated the port of a typical MathComp development to MathComp 2. We reviewed the available tools (documentation and HB tactics) and went through a number of concrete sample errors and warnings that we explained and corrected.

The example of porting CompDecModal demonstrated that the use of HB improves the Coq scripts in terms of readability and even allows for improvements. The complete fix can be found online: branch (relevant commit), it required the edition of 10 files, 35 insertions, and 67 deletions, which arguably represents a moderate amount of work.

The migration process to MathComp 2 will surely generate more questions that the community is ready to answer via the math-comp streams on https://coq.zulipchat.com.

#

  [short(type=newStructType)] HB.structure NewStruct params := {carrier of NewStruct_from_Struct parames carrier & Struct params carrier}.

(

  ***************************************************************************

  ***************************************************************************)

Figure 1 :

 1 Figure 1: The vicinity of the structure eqType in MathComp

  The next two compilation errors are similarly due to the removal of choiceType and CountType, and to the change of the notations [choiceMixin of _ by <:] and [countMixin of _ by <:]: > Canonical Structure fset_choiceType := Eval hnf in ChoiceType _ [choiceMixin of fset_type by <:]. Error: The reference ChoiceType was not found in the current environment. > Canonical Structure fset_countType (T : countType) := > Eval hnf in CountType _ [countMixin of fset_type T by <:]. Error: The reference CountType was not found in the current environment.

>

  Canonical Structure fset_subCountType (T : countType) := > Eval hnf in [subCountType of fset_type T].

  Definition _ := [isSub for elements by fset_type_rect]. HB.instance Definition _ := [Equality of fset_type by <:]. Canonical Structure fset_predType := PredType (fun (X : fset_type) x => nosimpl x \in elements X). HB.instance Definition _ := [Choice of fset_type by <:]. End FinSets. HB.instance Definition _ (T : countType) := [Countable of fset_type T by <:].

  Definition _ := [isSub for elements by fset_type_rect]. HB.instance Definition _ := [Choice of fset_type by <:]. Canonical Structure fset_predType := PredType (fun (X : fset_type) x => nosimpl x \in elements X). End FinSets. HB.instance Definition _ (T : countType) := [Countable of fset_type T by <:].

  Most Canonical (or Canonical Structure) commands are replaced by HB.instance (see Sect. 2) and there are small changes to MathComp notations such [subType ...], etc. Regarding CompDecModal, the first offending set of commands is the following (file fset.v): Note that the instance need not be named and better not be since it is the job of HB to figure out instances automatically. It is important to check that HB displays more than one message as a response to HB.instance, otherwise this might indicate a failed instantiation. Next compilation error: This error is primarily due to the remove of the EqType constructor [6, Sect. 2.1]. In fact, most xyzType constructors from MathComp should not be necessary anymore. See the changelog. Similarly to the [subType for _ by _] notation above, the [eqMixin of _ by <:] has changed:

	> Canonical Structure fset_eqType := EqType _ [eqMixin of fset_type by <:].
	Error: The reference EqType was not found in the current environment.
	-in `eqtype.v.
	..
	+ notation `[eqMixin of T by <:]`, use `[Equality of T by <:].
	..
	The fix is therefore:
	> HB.instance Definition _ := [Equality of fset_type by <:].
	HB_unnamed_factory_8 is defined
	eqtype_Equality__to__eqtype_hasDecEq is defined
	HB_unnamed_mixin_10 is defined
	fset_fset_type__canonical__eqtype_Equality is defined
	Section FinSets. fset_fset_type__canonical__eqtype_SubEquality is defined
	Variable T : choiceType.
	...
	End FinSets.
	Canonical Structure fset_countType (T : countType) :=
	Eval hnf in CountType _ [countMixin of fset_type T by <:].
	Canonical Structure fset_subCountType (T : countType) :=
	Eval hnf in [subCountType of fset_type T].
	This error is due to a change of notation that is documented in the changelog. Search for the string, say,
	"[subType" in CHANGELOG.md:
	-in `eqtype.v.
	..
	+ notation `[subType for v by rec]`, use `[isSub for v by rec].
	..
	The fix is therefore the following:
	> HB.instance Definition _ := [isSub for elements by fset_type_rect].
	HB_unnamed_factory_3 is defined
	fset_fset_type__canonical__eqtype_SubType is defined

Canonical Structure fset_subType := [subType for elements by fset_type_rect]. Canonical Structure fset_eqType := EqType _ [eqMixin of fset_type by <:]. Canonical Structure fset_predType := PredType (fun (X : fset_type) x => nosimpl x \in elements X). Canonical Structure fset_choiceType := Eval hnf in ChoiceType _ [choiceMixin of fset_type by <:]. Let us consider compilation errors in order: > Canonical Structure fset_subType := [subType for elements by fset_type_rect]. Error: Syntax error: [reduce] expected after ':=' (in [def_body]).

  .2.1): We can now use HB.howto to discover how to build this structure (see Sect. 3.2.2): Finally, we can use HB.about to learn about the parameters of a constructot (see Sect. 3.2.2):

	> HB.about Monoid.Law.
	HB: Monoid.Law.type is a structure (from "./ssreflect/bigop.v", line 415)
	...
	> HB.howto Monoid.Law.type
	HB: solutions (use 'HB.about F.Build' to see the arguments of each factory F):
	-Monoid.isLaw
	-SemiGroup.isLaw; Monoid.isMonoidLaw
	> HB.about Monoid.isLaw.Build
	...
	HB: arguments: Monoid.isLaw.Build T idm op opA op1m opm1
	-T : Type
	-idm : T
	-op : T -> T -> T
	-opA : associative op
	-op1m : left_id idm op
	-opm1 : right_id idm op

  This is similar to above: use HB.about to learn about Monoid.ComLaw and use HB.howto to inquire about its construction.

	HB.instance Definition _ : isCountable form := PCanIsCountable formChoice.pickleP.
	> HB.howto Monoid.ComLaw.type.
	HB: solutions (use 'HB.about F.Build' to see the arguments of each factory F):
	-Monoid.isComLaw
	-SemiGroup.isComLaw; Monoid.isMonoidLaw
	-SemiGroup.isCommutativeLaw; Monoid.isLaw
	-SemiGroup.isLaw; SemiGroup.isCommutativeLaw; Monoid.isMonoidLaw

  Equality Structures The next compilation errors in the files Kstar_def.v, gen_def.v, and CTL_def.v are handled as already explained above. Failing Rewriting The next compilation "error" is in the file CTL/demo.v. It is actually a tactic whose execution has been substantially slowed down by the upgrade to MathComp 2: move => [p' y]. rewrite /MRel /Mstate (negbTE (root_internal _)) [_ && _]/= orbF.The offending rewrite is the one with orbF. It used to apply by default to the left-hand side of an <-> equivalence but now the user is required to indicate the rewrite location with a pattern to recover a reasonable execution time: move => [p' y]. rewrite /MRel /Mstate (negbTE (root_internal _)) [_ && _]/=. rewrite [in X in X <-> _]orbF.

We have also observed at the time of this writing that the output of the HB.instance command may not be visible by default with VSCoq.

PCanIsCountable : forall [T : countType] [sT : Type] [f : sT -> T] [f' : T -> option sT], pcancel f f' -> isCountable.axioms_ sTThe fix in this case is therefore: