David Pallier

Vincent Le Cam

Sébastien Pillement

Energy-efficient GPS synchronization for wireless nodes

Keywords: Clock synchronization, WSN, Energy efficiency, FPGA, GPS, Kalman filter, Structural Health Monitoring

published or not. The documents may come

Energy-efficient GPS synchronization for wireless nodes

I. INTRODUCTION

M ONITORING the health of structures often implies the deployment of many sensor nodes to sample physical parameters such as ambient temperature, constraints, vibrations, acoustic waves, electrical fields, etc. Posted to a central supervisor, in delayed or in real-time, those samples are fed to Structural Health Monitoring (SHM) algorithms to assess the infrastructure health. As most of these parameters are timedependent and sampled from different nodes, the system needs to have a common time reference. The following examples underline the need for synchronization and the accuracy needed for SHM applications:

1) The natural frequency range of vibrations for structural diagnostics (modal analysis) is typically comprised between 0 and 100 Hz. Thus, accelerations are usually sampled at 1000 Hz. To ensure good modal reconstruction, the synchronization error between nodes must be lower than 1 ms. 2) Acoustic wave propagation is often used to localize the crack origin in a steel bar or a wire break in bridge cables. The method based on Time Difference Of Arrival (TDOA) is very dependent on the sensor's ability to timestamp acoustic wavefronts. With a typical velocity of 5000 m/s and a localization accuracy expected to +/-15 cm, synchronization error between two nodes must be lower than 30 microseconds.

3) The TDOA method is also used to monitor the localization of lightning on high electrical voltage lines. With

This work was supported by the WISE program of Pays-de-la-Loire region D. Pallier and V. Le Cam, Univ. Gustave Eiffel, Inria, COSYS-SII, I4S, F-44344 Bouguenais, France S. Pillement, Univ Nantes, CNRS, IETR UMR 6164, F-44000 Nantes, France a speed of 200 000 km/s and an expected localization of +/-50 m, synchronization error between two nodes must be under 250 nanoseconds. Time-counting on an electronic device, such as a wireless node, is mainly based on the use of an oscillator. Unfortunately, oscillators are highly sensitive to environmental conditions such as temperature, acceleration, or voltage supply stability [START_REF] Vig | Introduction to quartz frequency standards[END_REF]. Consequently, the frequency of the oscillator drifts from its theoretical value over time. Without a synchronization protocol, each node will have its own local time and the offset between the clocks of the nodes increases over time. In most of the cases, this synchronization protocol relies on timestamps exchange between the nodes. The synchronization accuracy that can be obtained in the network depends on the accuracy of the timestamps exchange, the synchronization update rate, and the stability of local oscillators.

A wireless sensor network (WSN) is a cost-effective way to deploy sensing nodes onto already existing large structures (bridges, wind power plants, railway lines, etc). This is why WSNs are more and more used for SHM applications. However, synchronization between nodes is harder in wireless networks due to non-deterministic propagation delays. Classical wireless synchronization protocols such as [START_REF] Elson | Fine-grained network time synchronization using reference broadcasts[END_REF], [START_REF] Ganeriwal | Timingsync protocol for sensor networks[END_REF] and [START_REF] Maróti | The flooding time synchronization protocol[END_REF] cannot achieve sub-microsecond accuracy, especially in multihop scenarios. In this context the Precision Time Protocol (PTP) has been extended over 802.11 communications with hardware timestamping [START_REF] Exel | Clock synchronization in ieee 802.11 wireless lans using physical layer timestamps[END_REF] or synchronization over UWB communications [START_REF] Beluch | A subnanosecond synchronized mac-phy cross-layer design for wireless sensor networks[END_REF] and can achieve sub nanosecond accuracy. While these protocols can offer good synchronization accuracy for the monitoring of small structures, they cannot be used for large structures like bridges, railways or power-lines, as they need close and iterative synchronization process. For large outdoor systems (which are most of SHM applications) a satel-lite navigation system, like GPS, can be used to synchronize each node to UTC as in [START_REF] Le Cam | Embedded subspace-based modal analysis and uncertainty quantification on wireless sensor platform pegase[END_REF]. This method is suitable for large infrastructures and allows timestamping accuracy up to tens of nanoseconds. The main limitation of this technique is the energy cost of a GPS receiver. Since energy consumption needs to be minimized in wireless systems, GPS solution needs to be optimized according to the required timestamping accuracy.

In this context, we developed a new synchronization scheme based on GPS holdover. Namely, GPS is used for precision purpose but it will be turned off to save power. During the off-state, the local time of the node is estimated using a clock model. The synchronization accuracy obtained with this scheme depends on the chosen holdover parameters: from tens of ns to 1 ms. To further improve the estimation of the clock model parameters, the GPS signal can be filtered when in on-state. The contributions of this work are: (1) the design of a node with dedicated hardware for timestamping, (2) the accuracy and energy efficiency evaluations of our synchronization scheme for different clock models, and (3) the study of the impact of frequency filtering using a Kalman filter.

The paper is organized as follows. The related works are summarized in Section II. Node architecture, timestamps computation, the use of GPS for periodic re-synchronization and the experimental setup are described in Section III. Section IV introduces the clock models, and the Kalman filter. The results are discussed in Section V, while Section VI resumes our contributions and outlines future perspectives.

II. RELATED WORKS

Since synchronization in distributed systems was first formalized in [START_REF] Lamport | Time, clocks, and the ordering of events in a distributed system[END_REF] and later extended in [START_REF] Srikanth | Optimal clock synchronization[END_REF], numerous solutions have been developed to tackle the problem of synchronization in networks. The Network Time Protocol [START_REF] Mills | Internet time synchronization: the network time protocol[END_REF] and the Precision Time Protocol [START_REF] Lee | Standard for a precision clock synchronization protocol for networked measurement and control systems[END_REF] are the most used network synchronization protocols over the wired Ethernet. It is possible to obtain microsecond synchronization accuracy with these protocols, and an extension of the PTP [START_REF] Serrano | The white rabbit project[END_REF] aims to achieve sub-nanosecond accuracy.

Classical synchronization protocols for WSNs are based on the exchange of dedicated timestamped RF beacons. These protocols can be divided into two categories: the sender to receiver and the receiver to receiver protocols. In the sender to receiver scheme, a node sends its time to a receiver that timestamp the arrival of the packet using its local time. Then, the receiver synchronizes itself by estimating the offset between its local time and the time of the sender. To estimate the propagation delay a two-way exchange has been developed in [START_REF] Ganeriwal | Timingsync protocol for sensor networks[END_REF]. The same principle is used in [START_REF] Maróti | The flooding time synchronization protocol[END_REF] but on a mesh topology instead of a tree and with the estimation of the packet processing overhead. Other works [START_REF] Van Greunen | Lightweight time synchronization for sensor networks[END_REF] [START_REF] Römer | Time synchronization in ad hoc networks[END_REF] [15] also rely on the same principle but are optimized for specific nodes with constrained resources like PicoRadios [START_REF] Rabaey | Picoradios for wireless sensor networks: the next challenge in ultra-low power design[END_REF]. In the receiver to receiver scheme [START_REF] Elson | Fine-grained network time synchronization using reference broadcasts[END_REF] [17] a third party node broadcast a beacon that doesn't contain a timestamp. All the receivers timestamp the arrival of the beacon with their local clock and then exchange their observations pairwise to calculate their offsets. The performances of all these protocols are limited by the propagation delays and the packet processing overhead. The interested reader can refer to [START_REF] Djenouri | Synchronization protocols and implementation issues in wireless sensor networks: A review[END_REF] that outlines the existing implementations of these techniques and their limitations. While some implementations of these protocols can achieve microsecond accuracy most of them are in the tens of microsecond range. Moreover, the accuracy decreases with the number of hops which can be problematic on large structures such as bridges, railways, or power-lines. Signal processing techniques have been described in [START_REF] Wu | Clock synchronization of wireless sensor networks[END_REF] to deal with non-deterministic propagation delays but, to the best of our knowledge, these techniques have not been implemented on real wireless nodes.

Other works focus on the implementation of clock synchronization in wireless communication norms. ZigBee is the most used protocol in the Low Rate Wireless Area Personal Network (LR WPAN), which adds network and application layers on top of the 802.15.4 protocol. In [START_REF] Cox | Time synchronization for zigbee networks[END_REF] the Flooding Time Synchronization Protocol (FTSP) is used to synchronize nodes over ZigBee to obtain a timestamping accuracy of 61 µs. The use of UWB for the physical layer was added to IEEE 802.15.4 as an amendment in 2007 and merged in IEEE 802.15.4-2011. In [START_REF] Beluch | A subnanosecond synchronized mac-phy cross-layer design for wireless sensor networks[END_REF] UWB is used to synchronize the clocks of a pair of chips on the same PCB with an accuracy of 374 ps (standard error of 677 ps). Several implementations of synchronization protocols over IEEE 802.11 are compared in [START_REF] Mahmood | Clock synchronization over ieee 802.11-a survey of methodologies and protocols[END_REF]. IEEE 802.11 compliant synchronization solutions can achieve µs clock synchronization accuracy [START_REF] Mahmood | Performance of ieee 802.11's timing advertisement against synctsf for wireless clock synchronization[END_REF] and non compliant solutions using PTP (IEEE 1588) can achieve nanosecond accuracy [START_REF] Exel | Clock synchronization in ieee 802.11 wireless lans using physical layer timestamps[END_REF] or sub-nanosecond accuracy with a hardware timestamping mechanism. While these synchronization schemes can offer a good synchronization accuracy on small structures, they are not suited for large structures, where the nodes can be distant of a few hundred meters to several kilometers. More recently, Low-Power Wide-Area Network (LPWAN) communications have gained interest for industrial applications. Based on LoRa, Sigfox, or NB-IoT, these protocols allow for communications over a few kilometers in urban areas. In [START_REF] Gao | A multichannel low-power wide-area network with high-accuracy synchronization ability for machine vibration monitoring[END_REF] LoRa and NB-Iot are used for node communications in the context of machine vibrations monitoring. Nodes are synchronized every 10 s using LoRa with an accuracy of 5 µs.

Another approach to wireless node synchronization is to use a satellite navigation system, like the GPS, to make each node synchronized to UTC. This solution has been implemented in [START_REF] Le Cam | Embedded subspace-based modal analysis and uncertainty quantification on wireless sensor platform pegase[END_REF] and [START_REF] Le Cam | Wireless sensors synchronization: an accurate and deterministic gps-based algorithm[END_REF]. The clock of the node is steered by the PPS signal coming out of its GPS receiver. To achieve an accuracy of tens of nanosecond, dedicated hardware is used with high-speed counters. The use of hardware induces deterministic timestamping which removes software overhead and thus improves the timing accuracy. While this approach is suitable for accurate synchronization, a GPS receiver requires a lot of energy (up to 100 mW). This cost can be affordable for applications requiring a high synchronization accuracy, but for less critical applications that require less accuracy, the energy consumption needs to be optimized. In this work, we propose an adaptive synchronization protocol based on automatic on/off switching of the GPS receiver to minimize the energy consumption under a defined accuracy constraint.

Adaptable synchronization accuracy has already been described in [START_REF] Palchaudhuri | Adaptive clock synchronization in sensor networks[END_REF] as an extension of RBS [START_REF] Elson | Fine-grained network time synchronization using reference broadcasts[END_REF]. In [START_REF] Hamilton | Aces: adaptive clock estimation and synchronization using kalman filtering[END_REF] a Kalman filter is used to track the offset of a local clock from timestamped message exchange and an algorithm is developed to find the optimal message exchange rate for a parametric accuracy. Our work differs from these solutions as it is designed on top of GPS synchronization (although any GNSS receiver can be used) instead of beacons exchange and relies on dedicated timestamping hardware to achieve up to tens of nanoseconds synchronization accuracy.

III. WIRELESS NODE

A. Node Architecture

The node design presented in this paper is similar to the one presented in [START_REF] Le Cam | Embedded subspace-based modal analysis and uncertainty quantification on wireless sensor platform pegase[END_REF]. It includes a dedicated hardware timestamping unit. The timestamping unit is connected to a GPS receiver and has inputs connected to digital signals that have to be timestamped (see Figure 1). The node is synchronized with the PPS signal and the epoch transmitted on a serial link from the GPS receiver. A rising edge on the PPS signal corresponds to the beginning of a UTC second. The epoch corresponds to the number of seconds elapsed since January 6th, 1980. It is used to compute the UTC date at the next rising edge of the PPS signal. The detection of a rising or a falling edge on a digital input of the node is called an event. The events are timestamped with high-speed counters in the timestamping unit and are stored until they are read by the main control unit of the node. This node architecture has several advantages:

• No processing overhead: The use of dedicated hardware on FPGA suppress processing overhead at event detection. Since the delay between the rising edge on the input and the timestamping is deterministic and not affected by software overhead due to process scheduling and interrupt latency, the timestamping accuracy is improved. • Parallel high granularity timestamping: On most microcontroller units, the number of parallel counters able to timestamp independent events is limited. Another limiting factor is the time granularity fixed by the quartz oscillator frequency. With a custom design, the number of parallel counters is only limited by the FPGA size. The frequency is fixed by the PLL included in the FPGA and can be higher than the quartz oscillator frequency, allowing a higher granularity. • Asynchronous data fetching and main control unit independence: Since timestamping and data processing are separated, the main control unit can fetch data (on digital inputs) asynchronously on the timestamping unit. Thus any operating system and any microcontroller can be used as long as the microcontroller is equipped with a serial link.

In this work, we focus on the timestamp unit design. For the main control of the node, without loss of generality and for quick prototyping purposes, we used a Raspberry Pi. minus the measured ∆[n] for this particular event. This time is positive when the event is detected after the true PPS edge and negative when detected before. Hence:

t[n] = t pps [n] + θ[n] -∆[n] (1)
C. Timestamping unit Architecture

The above timlestamp computation is implemented on a Spartan 6 FPGA and is shown in Fig. 1. Since this design aims at turning off the GPS, the clock synchronization is different from the one used in [START_REF] Le Cam | Embedded subspace-based modal analysis and uncertainty quantification on wireless sensor platform pegase[END_REF]. In this design, the clock counter is not reset by the PPS but instead is bounded by N 0 . This integer corresponds to the number of clock beats counted over one second at the nominal frequency f 0 of the quartz oscillator. In the meantime, the local time offset θ is recorded to convert the local timestamp to UTC. The fake PPS signal is generated inside the FPGA, by dividing the local clock frequency to obtain a 1Hz square signal. The offset corresponds to the differential between this signal and the true PPS signal from the GPS receiver. Tracking the local time offset has been chosen over the control of the clock as it allows for both offline and online synchronization. This enables us to perform the timestamp conversion, from local time to UTC, on the fly or post facto if the offset is recorded alongside the events. This post facto synchronization allows us to test different synchronization schemes from the same events dataset for evaluation purposes.

The values θ and ∆ are not directly measured inside the timestamp unit as they are integers corresponding to counter values and not timing values. Let N θ and N ∆ be the counter values corresponding to the offset counter and the event counter, respectively. These integers are converted into time values using the true frequency of the oscillator. Equation 1 becomes:

t[n] = t pps [n] + 1 f [n] (N θ [n] -N ∆ [n]) (2)
With N f , N θ and N ∆ the counter values from the timestamp unit, while t pps is obtained from the GPS receiver.

The number of local clock beats counted during a full period of the true PPS signal is N f . Let T GP S be the time elapsed between two rising edges of the true PPS signal. If the true PPS signal is perfect then T GP S is one second. The true frequency of the local oscillator can be computed as follows:

f [n] = N f [n] T GP S [n] (3)
giving:

t[n] = t pps [n] + T GP S N f [n] (N θ [n] -N ∆ [n]) (4)
If the true PPS signal is perfect, the only residual error is due to the granularity resolution as N f , N θ and N ∆ are integer values. The frequency output of the PLL inside the FPGA was set to 240MHz giving a counter resolution of 4.17 ns.

D. Periodic re-synchronization from GPS

As described in [START_REF][END_REF], GPS synchronization allows for a timing accuracy under 60ns. However, it entails a substantial energy cost as the GPS receiver is constantly on. When a lower accuracy is needed, GPS can be used periodically to re-synchronize the node. Between two synchronizations, the GPS receiver is off to save energy. The accuracy that can be obtained depends on the re-synchronization period. To be able to deliver the PPS and the epoch, the GPS receiver needs to decode parts of the navigation message broadcasted by the GPS satellites. In normal conditions the GPS receiver has to:

• Generate PPS: In order to generate the PPS signal, the GPS receiver needs to track the "in-view" satellites, compute its time and position (called fix) and then lock its clock. According to [START_REF][END_REF] the time to first fix is under 1 second when the ephemeris are still valid. From experimentations, we observed less than 1 second to get a fix and less than 4 seconds to lock the clock. • Refresh ephemeris: In order to track the in-view satellites the receiver needs to compute its position, therefore the ephemeris needs to be updated. In this work, we update them at their typical update rate of two hours. • Acquire a complete navigation message: Since the ephemeris are already acquired every two hours the complete navigation message is acquired once a day. In between the above phases, we propose to turn on and off the GPS to save energy. The off phase is called holdover and correspond to the time between two consecutive GPS time synchronizations. Figure 3 shows the on and off phases of the GPS receiver for a 24 hours period. The T_NAV is 25 minutes long every 24 hours. The T_EPH is 1 minute long every 2 hours. The GPS is on during k on seconds and in holdover during k -k on seconds. In this mode, the GPS receiver is switched off (except its memory to keep the ephemeris). Equation 5gives r, which represents the on/off ratio of the GPS receiver over 24 hours as a function of the parameters k, k on , T N AV and T EHP . The minimum k on is 5 seconds to get at least one locked pulse on the PPS signal, after this first delay the receiver delivers one pulse per second. Two pulses are required to get N θ , N ∆ and N f . The main difficulty is then to track the offset of the fake PPS, while the GPS is in holdover mode.

Fig. 3: Synchronization process from GPS

IV. OFFSET TRACKING A. Clock model

The offset needs to be predicted during the holdover to have accurate timestamping. To do so a clock model is used. Let C(t) be the node local time as a function of the real time.

C(t) = t 0 ω(τ )dτ + θ 0 (6)
With ω the clock rate and θ 0 the offset at the origin. A perfect clock would yield C(t) = t. The C(t) = ω.t + θ model is often used [2] [4]. In this model, the clock rate is constant but not equal to one. This is due to the fact that most of the oscillators are not functioning at their nominal frequency, but at a frequency supposed to be constant. While the latter is known to be false for crystal quartz oscillator, the variations of the frequency over a short time can be small enough to be considered constant since ambient temperature moves slowly. Fig. 4 represents the local time as a function of the real time for a perfect clock, a clock with a constant rate and a clock with a varying rate. The oscillator used in this design is an SPXO with an overall frequency stability of +/-20 ppm. Typically, the frequency stability of these oscillators is in the ppb range for sampling intervals of one second [START_REF] Mahmood | Clock synchronization over ieee 802.11-a survey of methodologies and protocols[END_REF]. The offset of a clock at time t is the difference between its local time C(t) and the real time t. Thus the local time offset θ(t) can be expressed as:

θ(t) = C(t) -t (7) = t 0 ω(τ )dτ + θ 0 -t (8)
Since the offset is sampled every second in the node, a discrete time model is needed. Equation 8 becomes:

θ[n] = n i=1 ω[i].τ [i] - n i=1 τ [i] + θ 0 (9)
With τ representing the duration of a supposed second for the node. This duration depends on the frequency of the node's oscillator and can be expressed as:

τ [n] = N 0 f [n] = N 0 N f [n] .T GP S [n] (10)
The clock rate ω[n] can be computed from the number of local clock beats counted during T GP S [n]:

ω[n] = N f [n] N 0 (11)
Therefore:

θ[n] = n i=1 T GP S [i] - n i=1 τ [i] + θ 0 (12)
If the PPS is considered perfectly synchronized to UTC, T GP S = 1s and thus:

θ[n] = n - n i=1 N 0 N f [i] + θ 0 ( 13
)
In holdover mode, τ is not known since the GPS receiver is off. The simplest method to estimate it is to use its last value that was observed when the GPS was on. Since the error on the offset is the sum of the error on τ , this error will grow as τ drift away from its last observed value. In that case, the synchronization accuracy will depend on the holdover time and the drift rate of τ . This model will now be referred to as the Constant Skew Clock Model (CSCM). Another model, that will be referred to as the Linear Skew Clock Model (LSCM), is to consider that τ is evolving linearly during the holdover mainly due to slow variations of the ambient temperature. In this model, the clock skew is considered as dynamic (see Fig. 4) and approximated by a piecewise linear function. The clock skew can be computed as:

τ [n + 1] = τ [n] + u ( 14
)
with u the common difference that can be evaluated for every holdover as:

u = τ [n] -τ [n -k] k ( 15
)
Under the hypothesis that u changes slowly, the last seen u can be used to predict τ during the holdover.

B. Kalman filter

In the previous section, the PPS signal was considered as perfect. In reality, this signal from the GPS receiver is subject to jitter, which creates noise in the frequency and offset measurement. According to [START_REF][END_REF], the RMS of the timing accuracy is 30ns. While this noise can be neglected when the GPS receiver is always on, it is not the case in the holdover mode. Indeed, the absolute offset error is the sum of the offset errors over time. The offset error growth, during holdover, depends on the accuracy of the measured frequency of the local oscillator at the beginning of the holdover.

To filter the noise on the PPS signal when the GPS is on, a Kalman filter was implemented. This filter uses the good short term accuracy of the local clock to filter the noisy PPS observations. The behavior of the fake PPS period τ during the short GPS on-time period is modeled by a random walk process:

τ [n + 1] = τ [n] + v[n] with v[n] ∼ N (0, σ 2 v ) (16)
Let z be the observations of τ from the noisy PPS signal:

z[n] = τ [n] + w[n] with w[n] ∼ N (0, σ 2 w ) (17)
The Kalman filtering process is based on the classical following equations [START_REF] Reid | Estimation ii[END_REF]:

• Prediction: τn+1|n = τn|n (18) p n+1|n = p n|n + q (19)
• Update:

ỹn = z[n] -τn|n-1 (20) s n = p n|n-1 + r (21) k n = p n|n-1 .s -1 n (22) τn|n = τn|n-1 + k n ỹn (23) p n|n = (1 -k n ).p n|n-1 (24)
With q = σ 2 v and r = σ 2 w . This filter is stable [START_REF] Souza | Riccati equations in optimal filtering of nonstabilizable systems having singular state transition matrices[END_REF]. The state variable is initialized with the observation. The variables r, q and p 0 are chosen empirically.

V. RESULTS AND DISCUSSION

To evaluate the synchronization algorithms detailed in the previous section, an experimental setup has been developed. Two wireless nodes based on the architecture detailed in section III were implemented on raspberry pis and spartan 6 development boards with their digital inputs connected to the same output of a low-frequency signal generator. This generator generates rising edges (i.e. events) every 2.34 sec that are timestamped by the nodes.

The nodes were placed indoor and the GPS antennas on the roof of the building. The fixed mode was used on the GPS receivers [START_REF][END_REF]. This GPS configuration is available on certain timing receivers from Ublox. It allows the user to define the GPS antenna location in the receiver to enhance the local time estimation and thus deliver a better PPS signal. The nodes recorded events during more than 2 days. Over 2 × 10 5 frequency and offset samples were recorded for each node as well as up to 8.5 × 10 4 event samples.

A. Clock models comparison

First, we compared the clock models. In Fig. 5 the estimated τ with the CSCM and the LSCM are plotted against the observed τ on a part of the dataset and for k = 1000s. The LSCM seems to be closer to the observations than the CSCM. However, Figure 6, which shows the same variables on another part of the dataset, the CSCM seems closer to the observations. The RMSE of the LSCM and the CSCM for the estimation of τ and θ on the full dataset and for k ranging from 10 to 3000, are summarized in Table I. The RMSE for τ and θ were computed on the difference between the predicted τ and θ and the GPS observations. We can see from this Table that the LSCM outperforms the CSCM for both τ and θ, when k > 1000. The PPS noise explained in Subsection IV-B can be seen in the observations in both figures. The fact that the LSCM yields poorer results than the CSCM for k 1000 is most probably due to this PPS noise. The impact of the Kalman filter on the offset tracking performance was evaluated. As seen in Section IV-B the Kalman filter is only used when the GPS is "on" to filter the PPS noise. Since the performance of the filter depends on the "on" mode duration, the effect of this duration was studied. To do so we ran the Kalman filter for different k on and for different k. Then the on/off ratio r was calculated according to Eq. 5. Figure 7 shows the r against the RMSE of the offset for different values of k on . The RMSE increases as the ratio decrease since a longer holdover yields a poorer synchronization accuracy but a better ratio. From these results we chose k on = 5 seconds since it seems to be the optimal duration for a RMSE ranging from 10ns to 50 µs. A lower k on slightly improves the ratio for longer holdover but worsen it for shorter ones.

Fig. 7: Optimal k on depending on the ratio r Figure 8 shows τ for k = 1000s and for the same part of the dataset than in Fig. 6. As can be seen the filtered LSCM seems to fit the observations better than the CSCM. Fig. 9 shows a comparison of the filtered and non-filtered models. The RMSE of the offset is plotted against k. It can be seen that both models can be used for the filtered and non-filtered cases depending on the holdover duration. It can also be seen that the Kalman Filter models yield better results than the nonfiltered models. For a holdover duration of less than 134s, the filtered CSCM produces better results, for a holdover duration over 134s the hypothesis of a slowing changing u seems valid and the filtered LSCM yields better results.

C. Timestamps

Then we compared the timestamping errors between the different models. The timestamps are computed according to Eq. 2. The timestamp errors are the differences between the timestamps computed with the GPS always "on" and the timestamps computed with the periodic GPS extinctions. Fig. 10 shows the RMSE of the timestamps errors as a function of k for both nodes. As expected the filtered CSCM yields Finally, we computed the GPS "on"/"off" ratio for all the models. These ratios are plotted against the RMSE on a semilog scale, in Fig. 11. It can be seen that for all models the ratio decreases exponentially as the RMSE increases and that they all converge to the limit of 2.5%. This limit corresponds to the minimal GPS ratio according to Eq. 5. The Kalman filter improves the ratio by a factor 3 under the microsecond accuracy despite a longer "on" state. The filtered CSCM yields a maximum improvement of 7% on the ratio over the filtered LSCM under 200s and the gain of the filtered LSCM above is marginal (0.5% max) as it happens for a RMSE greater than 300 ns. The results are similar for both nodes.

D. Discussion

Our synchronization solution allows for a RMSE under 10 ns with the GPS always "on". The RMSE stays under 20 ns on both nodes with a GPS "off" 60% of the time. For a GPS "off" 80% of the time and 95% of the time, the RMSEs are under 50 ns and 420 ns, respectively. As the computations were conducted offline (i.e a posteriori) to test different models, the GPS receiver was only virtually off but our study showed that the GPS receiver energy consumption could be divided However, the data showed that the distribution of the timestamps errors is not exactly normal as it has heavy tails. It was observed that the absolute maximum error is approximately ten times bigger than the RMSE for all our experiments. But even with an accuracy limit ten times bigger than the RMSE, the GPS can be switched off 60% of the time for a maximum error under 200 ns or 95% of the time for a maximum error under 4.2 µs. Table II summarizes the node to node theoretical range and the accuracy of existing synchronization protocols. TPSN, RBS, or FTSP with Berkley motes can offer µs synchronization accuracy on small networks with low complexity hardware. Standard IEEE 802.11 synchronization protocols offer the same synchronization accuracy, but subnanosecond accuracy can be obtained with the use of PTP. However, the range is still limited to local networks. Subnanosecond accuracy can also be obtained with UWB for an even smaller range. For long-range and wide networks, LoRa has been used to obtain 5 µs accuracy with synchronization updates every 10s. In comparison, the presented GPS based solution has a maximum accuracy of 10 ns (RMSE) and an accuracy of 4.2 µs for synchronization updates every 100 s.

While dedicated hardware is important for multiple channel high accuracy timestamping, the spartan 6 development board used in this work is not the best target to implement the design as it draws almost as much current as the GPS receiver. But it is only a proof of concept and the timestamp unit can be implemented on a dedicated ASIC or a smaller FPGA or a CPLD. We already implemented this unit on a smaller FPGA (lattice ICE40 [START_REF]Lattice ice40 datasheet[END_REF]) with a slightly lower resolution for a power consumption of 35 mW. For comparison, the ublox receiver needs 150mW when fully ON and the active GPS antenna draws another 10 to 60 mW (depending on the antenna used).

VI. CONCLUDING REMARKS

Synchronization is a challenging problem for wireless nodes, especially for applications requiring good timestamping accuracy across wide areas. Taking into consideration these applications we designed a node with a simple timestamping added in [START_REF] Palchaudhuri | Adaptive clock synchronization in sensor networks[END_REF] FTSP [START_REF] Maróti | The flooding time synchronization protocol[END_REF] ≈ 100 to 300 m (Berkley motes [START_REF] Ruiz-Sandoval | Sensor development using berkeley mote platform[END_REF]) TABLE II: Range and accuracy of synchronization protocols dedicated hardware that uses GPS for synchronization. In order to adapt the energy cost of the GPS receiver to the timestamping accuracy needed by the applications, we developed a synchronization scheme based on periodic GPS extinction. We found that both proposed models can be used depending on the timestamping accuracy needed. For a timestamp RMSE shorter than 300 ns, the constant skew model shows longer holdovers and for a RMSE above the 300ns, the linear model shows longer holdovers. However, it was shown that since most of the GPS ratio gains occur before 300 ns, the linear model yields minor improvements in terms of energy consumption. As expected the use of a Kalman filter to reduce the GPS noise improves the estimation of the local oscillator frequency and thus allows for longer holdover. But we also showed that this holdover length gain compensates for the longer required "on" state duration and overall the filtered models yield better GPS "on"/"off" ratios. We demonstrated that our solution can significantly reduce the energy cost of the GPS receiver depending on the required synchronization accuracy.

In future works, "online" implementations of our solution will be studied to evaluate the gains in energy consumption. Another interesting direction will be the automatic interruption of the holdover according to ambient temperature measurements. This could improve our synchronization scheme, especially for long holdovers, since it was observed that quick temperatures changes can lower the prediction accuracy. Finally, the use of time to digital converters will be investigated as it could allow us to lower the clock frequency inside the FPGA while maintaining, or even improving, the granularity.

Fig. 1 :

 1 Fig. 1: Timestamp unit

Fig. 2 :

 2 Fig. 2: Timing diagram The date t[n] of a recorded event occurring during the n th local second is computed based on the date of the true PPS rising edge t pps [n], that occurred during the fake PPS period and the elapsed time between the event and this rising edge. The elapsed time corresponds to the measured offset θ[n] minus the measured ∆[n] for this particular event. This time is positive when the event is detected after the true PPS edge and negative when detected before. Hence:

Fig. 4 :

 4 Fig. 4: Constant and dynamic skew of the local time over UTC

Fig. 5 :Fig. 6 :

 56 Fig. 5: Observed and predicted τ on the first part of dataset

Fig. 8 :

 8 Fig. 8: Observed and predicted τ using a Kalman filter

Fig. 9 :

 9 Fig. 9: Offset RMSE for the different models

Fig. 10 :

 10 Fig. 10: Timestamps RMSE on both nodes

Fig. 11 :

 11 Fig. 11: GPS ratio as a function of timestamps RMSE

TABLE I :

 I RMSE of the CSCM and LSCM

	k	CSCM	LSCM
		RMSE τ	RMSE θ	RMSE τ	RMSE θ
	10	1.3394E-08	5.8209E-08 1.8527E-08	7.8107E-08
	25	1.3626E-08	1.4230E-07 1.8483E-08	1.9750E-07
	50	1.3656E-08	2.8072E-07 1.8467E-08	3.9577E-07
	100	1.3855E-08	5.7021E-07	1.8680E-08	8.0479E-0
	250	1.4234E-08	1.4609E-06	1.8604E-08	2.0049E-06
	500	1.5341E-08	3.1602E-06	1.8715E-08	4.0681E-06
	1000	1.9235E-08	7.9718E-06	2.0007E-08	8.6390E-06
	2000	2.9425E-08	2.3955E-05	2.5034E-08	2.0925E-05
	3000	3.9110E-08	4.6289E-05	3.2004E-08	3.7230E-05


Mogens Aoxing Liu

Sandø Lund

Didier Boichard

Xiaowei Mao

Emre Karaman

Sebastien Fritz

Gert Pedersen Aamand

Yachun Wang
email: wangyachun@cau.edu.cn

Guosheng Su
email: guosheng.su@mbg.au.dk

Mogens Sandø Lund

imputation for sequencing variants preselected to a customized low-density chip

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.

of additional SNPs preselected through association analyses and bioinformatics analyses in major dairy breeds of Denmark-Finland-Sweden (DFS) and France (FRA). A previous study in genotyped animals of Danish Jersey showed that the prediction reliability of milk and protein improved by including DFS and/or FRA SNPs [START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF] . In Denmark, this customized low-density chip has started to be used to genotype young bulls and heifers, which were used to be genotyped by the standard low-density chip. However, there are still many animals that have already been genotyped by the standard SNP chip, e.g., the 54 K chip or the standard low-density chip, and thus lack of information on preselected sequencing SNPs (DFS and FRA SNPs). Because these animals are older than the animals genotyped with the customized SNP chip and may have phenotypes, they are especially important to assess the effects of the preselected sequencing SNPs.

Missing genotypes of an individual can be imputed by using the information from the family or/and the population [START_REF] Weigel | Prediction of unobserved single nucleotide polymorphism genotypes of Jersey cattle using reference panels and population-based imputation algorithms[END_REF] . Many studies have been performed to impute SNPs from a lower density SNP panel to a higher density SNP panel, with the majority focusing on SNPs from the standard SNP chips in which SNPs were evenly distributed on the genome and had relatively high minor allele frequency (MAF) [START_REF] Boichard | Design of a bovine low-density SNP array optimized for imputation[END_REF][START_REF] Hozé | High-density marker imputation accuracy in sixteen French cattle breeds[END_REF][START_REF] Su | Genomic relationships based on X chromosome markers and accuracy of genomic predictions with and without X chromosome markers[END_REF][START_REF] Su | Comparison of genomic predictions using medium-density (∼ 54,000) and high-density (∼ 777,000) single nucleotide polymorphism marker panels in Nordic Holstein and Red Dairy Cattle populations[END_REF] . Unlike the SNPs in the standard SNP chip, the preselected sequencing SNPs are unevenly distributed on the genome and some of the SNPs could have low MAF [START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF][START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF] . Brøndum et al. [START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF] showed that the correlation between imputed and observed genotypes from the high-density chip to WGS was 0.90, 0.89, and 0.87 for Nordic Holsteins, Danish Jersey, and Nordic Red, with a combined reference population of 116 Nordic Holsteins, 15 Danish Jersey, 23 Nordic Red, and 16 Brown Swiss. For WGS data, the relatively low imputation accuracy was also one of the reasons caused limited improvements in genomic prediction compared with using a high-density SNP chip [START_REF] Pérez-Enciso | Sequence-vs. chip-assisted genomic selection: accurate biological information is advised[END_REF] . Instead of using the WGS data as a reference panel, using a customized low-density chip could achieve a much larger reference population but with lower linkage disequilibrium between preselected sequencing SNPs and the other SNPs in the chip. By using a customized low-density SNP chip as a reference panel, Sanchez et al. [START_REF] Sanchez | Confirmation of candidate causative variants on milk composition and cheesemaking properties in Montbéliarde cows[END_REF] imputed 194 candidate variants for milk protein and fatty acid composition in Montbéliarde cows and showed that the mean squared correlation between imputed and true genotypes was 91.6% for SNPs with MAF higher than 1%.

The imputation accuracy could be highly influenced by the MAF of SNPs [START_REF] Ma | Comparison of different methods for imputing genome-wide marker genotypes in Swedish and Finnish Red Cattle[END_REF][START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF] . Ma et al. [START_REF] Ma | Comparison of different methods for imputing genome-wide marker genotypes in Swedish and Finnish Red Cattle[END_REF] showed that the lower the MAF, the higher the concordance rate but the lower the correlation between imputed and observed genotypes, when imputing from 3 K to 54 K or from 54 K to the high-density chip in a combined population of Swedish and Finnish Red Cattle, no matter which imputation algorithm was applied. For the imputation of preselected sequencing SNPs using a customized low-density chip as a reference panel, the reference population usually comprises young animals. This situation could be rare when imputing SNPs in the standard SNP chip or to WGS, where the animals with the target SNP panel are usually key animals such as ancestors [START_REF] Calus | Evaluation of measures of correctness of genotype imputation in the context of genomic prediction: a review of livestock applications[END_REF] . The relationship between reference and target populations can largely impact the imputation accuracy [START_REF] Mao | Imputation of markers on the bovine X chromosome[END_REF] . Therefore, the imputation accuracy of preselected sequencing SNPs and the influence of MAF on the imputation accuracy could differ for different groups of animals.

The integration of preselected DFS and/or FRA SNPs covered by the EuroGenomics customized SNP chip could benefit genomic prediction [START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF] . In the present study, the primary objective was to investigate the imputation accuracy for DFS + FRA SNPs in both contemporary animals and old bulls, using the customized low-density chip as a reference panel. Besides, we investigated the pattern of imputation accuracy in relation to MAF for the preselected DFS + FRA SNPs.

Methods

ethics statement. All data were recorded as part of the routine genomic evaluation, no animal handling or experiment was performed specially for the present study.

Data.

A total of 3,745 Danish Jersey bulls, 1,168 US Jersey bulls, and 28,678 Danish Jersey cows were available in the present study. The Danish Jersey bulls were mainly genotyped with the Illumina Bovine SNP50 chip (54 K, Illumina, Inc). The US Jersey bulls, which had close relationships with Danish Jersey [START_REF] Su | Sharing reference data and including cows in the reference population improve genomic predictions in Danish Jersey[END_REF] , and currently used in routine genomic evaluation of Danish Jersey, were genotyped with either the 54 K chip or the GeneSeek Genomic Profiler HD chip (HD, GeneSeek, Neogen Corp.). For US Jersey bulls, only SNPs in the 54 K chip were kept for analyses. Danish Jersey cows were genotyped with the 54 K chip (3%), the standard Illumina Bovine LD chip (standard LD, Illumina, Inc.) (49%), or the customized Illumina Bovine LD chip [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF] (customized LD, Illumina, Inc.) (48%).

The customized LD chip included SNPs in the standard LD chip, and more importantly, additional SNPs preselected from the WGS data in major dairy breeds of Denmark-Finland-Sweden (DFS) and France (FRA) [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF] . For DFS SNPs, a total of 1,715 sequencing SNPs were preselected through association analyses based on imputed WGS data in Nordic Holsteins, Nordic Red, and Danish Jersey [START_REF] Brøndum | Quantitative trait loci markers derived from whole genome sequence data increases the reliability of genomic prediction[END_REF][START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF] . For FRA SNPs, a total of 4,325 SNPs were preselected through both association analyses and bioinformatics analyses [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF] . In bioinformatics analyses, SNPs from regulatory regions of genes, breakpoints of structural variations, and SNPs with strong predicted effects using Variant Effect Predictor 23 (e.g., frameshift, stop gain, splicing site, and non-synonymous substitution) were preselected [START_REF] Boichard | Experience from large scale use of the EuroGenomics custom SNP chip in cattle[END_REF] .

The SNPs on sex chromosome or with an unknown position, and with monomorphism or multiple alleles were removed from analyses using Plink [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF] . Ultimately, 53,921 autosomal SNPs (including the SNPs of the 54 K chip and the preselected sequencing SNPs) from 33,591 individuals were available for further analyses.

Genotype imputation.

Genotype imputation was carried out using a family and population-based approach implemented in a commercial version of FImpute [START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF] , where no limitation was set for the reference population size. The FImpute has been proven to be able to provide accurate imputation with low computational requirements in dairy cattle [START_REF] Ma | Comparison of different methods for imputing genome-wide marker genotypes in Swedish and Finnish Red Cattle[END_REF][START_REF] Mao | Imputation of markers on the bovine X chromosome[END_REF] . In this study, the pedigree information was provided for all analyses. The pedigree was constructed by tracing back genotyped animals as many generations as possible, consisting of 6,102 bulls and 66,466 cows. The SNP-wise imputation accuracy was measured as correlation coefficient and concordance rate. The correlation was calculated as the Pearson correlation between observed and imputed genotypes (coded as 0, 1, or 2). The concordance rate was calculated as the proportion of correctly imputed genotypes to all imputed genotypes. In the present study, the imputation accuracy of preselected sequencing SNPs was investigated in both contemporary animals and old bulls.

Imputation of preselected sequencing SNPs for contemporary animals.

In the present study, the customized LD chip was used as a reference panel for the imputation of DFS and FRA SNPs. The animals genotyped with the customerized LD chip were relatively young compared with animals genotyped with other versions of SNP chips. For the imputation of DFS and FRA SNPs, a random subset of animals genotyped with the customized LD chip were used as the target population, while all remaining animals genotyped with the customized LD chip were used as the reference population. Therefore, the target population for the imputation of DFS and FRA SNPs was the contemporary animals of the reference population. For contemporary animals, three imputation scenarios were compared regarding the imputation strategies of preselected sequencing SNPs. In scenario one, all genotyped individuals were imputed to the standard LD chip together with preselected DFS and FRA SNPs (LD + DFS + FRA). Here, for animals genotyped with the 54 K chip, only SNPs which appeared in the standard LD chip were kept for analyses. In scenario two, all genotyped individuals were directly imputed to 54 K together with the preselected DFS and FRA SNPs (one-step 54 K + DFS + FRA). In scenario three, individuals genotyped with different LD chips were firstly imputed to 54 K, and then individuals genotyped with 54 K or with imputed 54 K were further imputed to 54 K together with preselected DFS and FRA SNPs (two-step 54 K + DFS + FRA).

In all three scenarios, the imputation of DFS and FRA SNPs was assessed by masking genotypes of DFS and FRA SNPs for a random subset of 2,837 individuals (around 20%) from 14,372 individuals genotyped with the customized LD chip. Besides, for scenarios one-step and two-step 54 K + DFS + FRA, the imputation from the standard LD chip to the 54 K chip was assessed by masking genotypes for a random subset of 986 individuals (around 20%) from 4,913 individuals genotyped with the 54 K chip to become the standard LD chip. The accuracy of imputation was calculated by comparing the imputed genotypes with the observed genotypes for the masked SNPs. The proportion of animals being masked to assess the imputation accuracy was subjective, but a ratio about 80:20 is often used for dividing a whole dataset into training and test datasets. In the present study, considering the number of individuals being available for each SNP group, 20% animals were masked in order to have a good balance between validation and reference population sizes.

Imputation of preselected sequencing SNPs for old bulls.

In the present study, most old bulls were genotyped with the 54 K chip and therefore lack of the information of DFS + FRA SNPs. To mimic the imputation of preselected sequencing SNPs (DFS + FRA SNPs) in old bulls, 500 SNPs selected from the standard LD chip were taken as preselected sequencing SNPs. Regarding the number of SNPs being selected to validate the imputation accuracy of preselected sequencing SNPs (500 in the present study), it was based on the considerations that the number of SNPs could be sufficient to control the sampling error in calculating imputation accuracy and moving this number of SNPs could not break the linkage disequilibrium structures among the remaining SNPs in the standard LD chip. Here, the reference population for the imputation of these 500 SNPs in old bulls comprised all genotyped cows. In total, 864 bulls were sires, grandsires, or great-grandsires of the reference population, while the remaining 4,157 bulls had more distant relationship with the reference population. To further investigate the influence of the number of genotyped progenies on the imputation of preselected sequencing SNPs in old bulls, the imputation accuracy average over these 500 SNPs was calculated for individuals with different number of genotyped daughters.

The distributions of MAF for SNPs in the standard LD chip were different from those of the preselected sequencing SNPs, which could potentially influence the assessment. Thus, we applied two strategies to select these 500 SNPs from the standard LD chip. One strategy was to randomly select 500 SNPs from the standard LD chip without considering MAF (Random500). In this case, the distribution of MAF for the 500 SNPs was similar to that of the standard LD chip. The other strategy was to select SNPs based on the MAF distribution of preselected sequencing SNPs (MAF500), where the number of SNPs being selected for each MAF class was equal to 500 multiplied the proportion of preselected sequencing SNPs in that MAF class (0.01-0.05, 0.05-0.10, 0.10-0.20, 0.20-0.30, 0.30-0.40, and 0.40-0.50). In this case, the distribution of MAF of these 500 selected SNPs was similar to that of preselected sequencing SNPs (DFS + FRA SNPs). Then, a two-step imputation was used to impute the selected 500 SNPs in old bulls. Animals genotyped with different LD chips were firstly imputed to the reduced 54 K (54 K minus the 500 SNPs), and then to the full 54 K (including the 500 SNPs). We repeated these procedures five times.

Results

imputation of preselected sequencing Snps for contemporary animals. For contemporary animals, imputation accuracy for DFS SNPs, FRA SNPs, and for the imputation from the standard LD chip to the 54 K chip are presented in Table 1. For different imputation scenarios, ranks of imputation accuracy were consistent regarding correlations and concordance rates. Among different SNP sets, imputation accuracy was highest for the imputation from the standard LD chip to the 54 K chip and lowest for FRA SNPs. For example, for SNPs with MAF higher than 0.01 in the two-step imputation to 54 K + DFS + FRA, average correlations ranged from 88.0% for FRA SNPs to 96.2% for the imputation from the standard LD chip to the 54 K chip, and average concordance rates ranged from 94.2% for FRA SNPs to 98.0% for the imputation from the standard LD chip to the 54 K chip. For different imputation scenarios, the highest accuracy was obtained using two-step 54 K + DFS + FRA and the lowest using one-step 54 K + DFS + FRA. For example, for DFS and FRA SNPs with MAF higher than 0.01, average correlations ranged from 88.2% when using one-step 54 K + DFS + FRA to 89.2% when using two-step 54 K + DFS + FRA, while average concordance rates ranged from 94.7% when using one-step 54 K + DFS + FRA to 95.1% when using two-step 54 K + DFS + FRA.

The pattern of imputation accuracy in relation to MAF is presented in Table 1 and Fig. 1. Generally, SNPs with lower MAF tended to have lower correlations but higher concordance rates, which were observed within all three SNP sets. For example, in the scenario two-step 54 K + DFS + FRA, compared with no restriction on MAF, removing SNPs with MAF lower than 0.01 increased correlations by 0.8, 2.3, and 2.3 percentage points while decreased concordance rates by 0.3, 0.6, and 1.8 percentage points, for the imputation from the standard LD chip to the 54 K chip, the imputation of DFS SNPs, and the imputation of FRA SNPs, respectively. imputation of preselected sequencing Snps for old bulls. The proportions of SNPs in different MAF groups for SNPs in the standard LD chip, and DFS + FRA SNPs are presented in Table 2. Generally, more DFS and FRA SNPs had low MAF than SNPs in the standard LD chip. For example, the number of SNPs with MAF within 0.01~0.05 in DFS and FRA SNPs was more than twice as that in the standard LD chip. The distributions of MAF for the 500 selected SNPs in Random500 and MAF500 were similar to those for SNPs in the standard LD chip and DFS + FRA SNPs, respectively.

In old bulls, imputation accuracy for MAF500 and Random500 are presented in Table 3. All results presented here were specific for SNPs with MAF higher than 0.01. The imputation accuracy for MAF500 and Random500 was similar and relatively high. For example, the correlation was 97.3% for MAF500 and 97.7% for Random500. In all cases, imputation accuracy of the 500 SNPs, which was used to mimic the imputation accuracy of preselected sequencing SNPs in old bulls, was higher than that of the preselected sequencing SNPs in contemporary animals.

All genotyped cows were from 588 paternal half-sib families (sire families), with family size up to 2,182 in the largest half-sib family. Among these 588 sires of genotyped cows, 257 were genotyped. Imputation accuracy for bulls with a different number of genotyped daughters in MAF500 and Random500 are presented in Table 3. Compared with the imputation for contemporary animals, imputation accuracy was high even for bulls without genotyped daughters, e.g., correlations higher than 97.2% and concordance rates higher than 98.4%. Generally, imputation accuracy in terms of both correlation and concordance rate increased when bulls had more genotyped daughters, but the increase slowed down after bulls had more than 100 genotyped daughters. For example, for the scenario MAF500, correlation increased 2.5 percentage point on average when the number of daughters increased from 0 to 100, but only increased 0.1 percentage point on average when the number of daughters increased from 100 to 200.

Discussion

The benefits of integrating sequencing SNPs preselected from association analyses and bioinformatics analyses have been well documented in genomic prediction [START_REF] Brøndum | Quantitative trait loci markers derived from whole genome sequence data increases the reliability of genomic prediction[END_REF][START_REF] Van Den Berg | Using sequence variants in linkage disequilibrium with causative mutations to improve across breed prediction in dairy cattle: a simulation study[END_REF][START_REF] Liu | Improvement of genomic prediction by integrating additional single nucleotide polymorphisms selected from imputed whole genome sequencing data[END_REF] . However, the magnitude of improvement in prediction reliability could be highly related to the imputation accuracy of these preselected sequencing SNPs [START_REF] Zhang | Impact of rare and low-frequency sequence variants on reliability of genomic prediction in dairy cattle[END_REF][START_REF] Pausch | Evaluation of the accuracy of imputed sequence variant genotypes and their utility for causal variant detection in cattle[END_REF] . By using a customized LD chip as a reference panel, the present study investigated the imputation accuracy of sequencing SNPs preselected from association analyses and bioinformatics analyses in both contemporary animals and old bulls.

In contemporary animals, imputation accuracy for preselected sequencing SNPs was much lower than that for the imputation from the standard LD chip to the 54 K chip, although the number of animals used as the reference population for imputing preselected sequencing SNPs was almost three times as that for the imputation from the standard LD chip to the 54 K chip. The accurate imputation from the standard LD chip to the 54 K chip has been well documented in various dairy cattle populations [START_REF] Boichard | Design of a bovine low-density SNP array optimized for imputation[END_REF][START_REF] Su | Genomic relationships based on X chromosome markers and accuracy of genomic predictions with and without X chromosome markers[END_REF] . One of the reasons for the high accuracy of imputation from the standard LD chip to the 54 K chip is that the standard LD chip was explicitly designed for the imputation to the 54 K chip, where both the SNPs in the standard LD chip and the SNPs to be imputed had relatively high MAF [START_REF] Boichard | Design of a bovine low-density SNP array optimized for imputation[END_REF] . The preselected sequencing SNPs, however, had lower MAF than SNPs in the standard SNP chip. In the present study, 14.7% of SNPs from the 54 K chip had MAF lower than 0.1, whereas 20.9% of DFS SNPs and 30.8% of FRA SNPs had MAF lower than 0.1. The poor imputation accuracy of SNPs with MAF lower than 0.1 has been widely observed in previous studies [START_REF] Daetwyler | Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle[END_REF][START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF][START_REF] Van Den Berg | Imputation to whole-genome sequence using multiple pig populations and its use in genome-wide association studies[END_REF] . Bolormaa et al. [START_REF] Bolormaa | Accuracy of imputation to whole-genome sequence in sheep[END_REF] showed that missense SNPs had lower MAF and lower imputation accuracy than intergenic and intronic SNPs in sheep. In addition, the difference in MAF distributions between SNPs in the standard LD chip and the preselected sequencing SNPs also indicated that there would be weak linkage disequilibrium between these two groups of SNPs. Previous studies have shown that the linkage disequilibrium between SNPs to be imputed and their closest SNPs could influence the imputation accuracy [START_REF] Sargolzaei | A new approach for efficient genotype imputation using information from relatives[END_REF][START_REF] Calus | Evaluation of measures of correctness of genotype imputation in the context of genomic prediction: a review of livestock applications[END_REF] . Besides, within preselected sequencing SNPs, imputation accuracy was lower for FRA SNPs compared with DFS SNPs. One explanation for this could be that FRA SNPs had lower MAF than DFS SNPs. In addition, in order to be significant from association analyses, the genotyping and imputation qualities of these SNPs are usually good [START_REF] Gibson | Rare and common variants: twenty arguments[END_REF] . Among different imputation strategies investigated in contemporary animals of the reference population, differences in imputation accuracy were observed for preselected sequencing SNPs. The two-step 54 K + DFS + FRA yielded higher imputation accuracy for DFS + FRA SNPs than the one-step LD + DFS + FRA. In two-step 54 K + DFS + FRA, the first step was to impute from the standard LD chip to the 54 K chip, and the second step was to impute from imputed 54 K to 54 K + DFS + FRA. Given the high accuracy for the first step, the two-step 54 K + DFS + FRA scenario could be considered as imputation from 54 K to 54 K + DFS + FRA. Regarding the imputation of DFS + FRA SNPs, the two-step 54 K + DFS + FRA scenario had a denser reference panel than the LD + DFS + FRA scenario, and therefore, a higher imputation accuracy. When imputing to 54 K + DFS + FRA, the two-step strategy yielded higher imputation accuracy than the one-step strategy which imputed missing genotypes of animals genotyped with the standard LD, the customized LD, and the 54 K chips simultaneously. The performance of a two-step imputation largely impacted by the performance of its first step [START_REF] Korkuc | Finding the optimal imputation strategy for small cattle populations[END_REF] . In the two-step imputation used in the present study, the first step (imputation from the standard LD chip to the 54 K chip) achieved high imputation accuracy in terms of both correlation and concordance rates (higher than 95%), which laid a good foundation for the second step (from imputed 54 K to 54 K + DFS + FRA). The result was consistent with the previous studies. In Montbéliarde cows [START_REF] Sanchez | Confirmation of candidate causative variants on milk composition and cheesemaking properties in Montbéliarde cows[END_REF] , a similar two-step imputation strategy was used to impute preselected sequencing SNPs when using a customized low-density SNP chip as a reference panel. In other dairy cattle populations, a multiple-step imputation strategy was superior to a one-step imputation strategy when imputing from 3 K and 6 K chips to the HD chip where animals were firstly imputed to the 54 K chip and then to the HD chip [START_REF] Vanraden | Genomic imputation and evaluation using high-density Holstein genotypes[END_REF] , or from the 54 K chip to WGS where animals were firstly imputed to the HD chip and then to WGS [START_REF] Van Binsbergen | Accuracy of imputation to whole-genome sequence data in Holstein Friesian cattle[END_REF] . Besides, when using the one-step imputation strategy, the one-step LD + DFS + FRA slightly outperformed the one-step 54 K + DFS + FRA for imputing genotypes of DFS + FRA SNPs. A previous study in Holstein Friesian cattle showed that the imputation from the 54 K chip to WGS, where one half of the reference animals were genotyped with WGS and the remaining half were genotyped with the HD chip, outperformed that using all reference animals with WGS [START_REF] Van Binsbergen | Accuracy of imputation to whole-genome sequence data in Holstein Friesian cattle[END_REF] . The authors of that study explained that the imputation algorithm cannot select the correct haplotype when there were multiple possible matches between sequencing haplotypes and the 54 K haplotypes, whereas less possible matches when the HD chip was added in between [START_REF] Van Binsbergen | Accuracy of imputation to whole-genome sequence data in Holstein Friesian cattle[END_REF] . In the present study, more haplotypes could be phased in the data set of the one-step imputation to 54 K + DFS + FRA than that to LD + DFS + FRA. Due to a larger number of missing SNPs to be imputed, the difficulty to match to the correct haplotype when there were multiple matches available could be one of the reasons why one-step 54 K + DFS + FRA was slightly inferior to one-step LD + DFS + FRA regarding the imputation of DFS + FRA SNPs.

In the present study, both correlation and concordance rates were applied to measure the imputation accuracy. For contemporary animals, SNPs with lower MAF tended to have lower correlations but higher concordance rates. This was consistent with that reported in Ma et al. [START_REF] Ma | Comparison of different methods for imputing genome-wide marker genotypes in Swedish and Finnish Red Cattle[END_REF] . The reason was that the SNP with a lower MAF tended to be imputed to the major allele, and thus, a relatively high concordance rate could be achieved since most major alleles can be correctly imputed [START_REF] Ma | Comparison of different methods for imputing genome-wide marker genotypes in Swedish and Finnish Red Cattle[END_REF] . Another study performed by Mulder et al. [START_REF] Mulder | Imputation of genotypes with low-density chips and its effect on reliability of direct genomic values in Dutch Holstein cattle[END_REF] showed that the concordance rate was highly influenced by the MAF, while the correlation was expected to be linearly related to the accuracy of genomic prediction. Calus et al. [START_REF] Calus | Evaluation of measures of correctness of genotype imputation in the context of genomic prediction: a review of livestock applications[END_REF] strongly suggested to use correlation as a measure of imputation accuracy, especially for sequencing SNPs with low MAF.

A previous study in Nordic Holsteins, however, reported a very limited improvement of reliability from integrating imputed sequencing SNPs with low MAF into genomic prediction [START_REF] Zhang | Impact of rare and low-frequency sequence variants on reliability of genomic prediction in dairy cattle[END_REF] , which could be caused by the poor imputation of low MAF SNPs. The preselected sequencing SNPs generally had lower MAF than SNPs in the standard SNP chip (Table 2), more efforts on the imputation of sequencing SNPs with low MAF could be necessary in order to benefit from preselected sequencing SNPs. Or at least, filtering SNPs by imputation accuracy should be performed before further processes in genomic prediction or association studies [START_REF] Van Den Berg | Imputation to whole-genome sequence using multiple pig populations and its use in genome-wide association studies[END_REF] . In a previous study of Danish Jersey 7 , reliabilities of genomic prediction increased for milk (e.g. from 0.320 for 54 K to 0.424 for 54 K + DFS + FRA in Danish Jersey bulls) and protein (e.g. from 0.268 for 54 K to 0.306 for 54 K + DFS + FRA in Danish Jersey bulls) after integrating imputed DFS and FRA SNPs with imputation accuracy, measured by both correlation and concordance rate, higher than 0.8. Alternatively, dosage scores instead of integer genotype codes could be used as explanatory variables in genomic prediction and association studies in order to account for the uncertainty of imputation, especially when SNPs has low imputation accuracy [START_REF] Pausch | Evaluation of the accuracy of imputed sequence variant genotypes and their utility for causal variant detection in cattle[END_REF][START_REF] Van Den Berg | Imputation to whole-genome sequence using multiple pig populations and its use in genome-wide association studies[END_REF] .

Previous studies showed that the closer the relationship between reference and target populations, the higher the imputation accuracy [START_REF] Weigel | Prediction of unobserved single nucleotide polymorphism genotypes of Jersey cattle using reference panels and population-based imputation algorithms[END_REF][START_REF] Zhang | Marker imputation with low-density marker panels in Dutch Holstein cattle[END_REF] . In the present study, the imputation for preselected sequencing SNPs in old bulls was investigated by imputing 500 SNPs of the standard LD chip. The imputation accuracies of these 500 selected SNPs, even for bulls without daughters, were higher than those for DFS and FRA SNPs in contemporary animals of the reference population. For old bulls, the missing genotypes of these 500 SNPs could be accurately imputed since most haplotypes of old bulls were available in genotyped cows, especially when old bulls had daughters, granddaughters, or great-granddaughters in the reference population. Furthermore, for SNPs with low MAF, both correlations and concordance rates in old bulls were quite high. These results indicated that the genotypes of old bulls could be accurately imputed by using information from multiple genotyped progenies. Thus, we suggested that it would be unnecessary to re-genotype old bulls for preselected sequencing SNPs with the customized LD chip.

In the present study, a customized LD chip was used as a reference panel to impute preselected sequencing SNPs. An alternative strategy is to impute all WGS SNPs using the WGS data as a reference panel. The reference population could be a small number of key animals selected from pedigree relationships or haplotype diversities or WGS data from relevant consortium such as the 1000 genomes bull project [START_REF] Druet | Toward genomic prediction from whole-genome sequence data: impact of sequencing design on genotype imputation and accuracy of predictions[END_REF] . Brøndum et al. [START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF] showed that the correlation between imputed and observed genotypes from the HD chip to WGS was 0.89 for Danish Jersey with a combined reference population. Compared with imputing preselected sequencing SNPs using a customized SNP chip as a reference panel, using the WGS data as a reference panel would be challenged by the limited reference population size, the intensive computational burden, and the difference in diversity and frequency of phased haplotypes between reference and target populations if only a small number of animals being sequenced [START_REF] Butty | Optimizing Selection of the Reference Population for Genotype Imputation From Array to Sequence Variants[END_REF] . The present study showed that the imputation accuracy using a customized LD chip as a reference panel was much higher than that based on WGS imputation reported in previous studies [START_REF] Daetwyler | Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle[END_REF][START_REF] Brøndum | Strategies for imputation to whole genome sequence using a single or multi-breed reference population in cattle[END_REF][START_REF] Van Binsbergen | Accuracy of imputation to whole-genome sequence data in Holstein Friesian cattle[END_REF] . However, the imputation accuracy for preselected sequencing SNPs in contemporary animals was still not very high (correlation lower than 92% for SNPs with MAF higher than 0.01). It would be a good strategy to directly genotype these preselected sequencing SNPs, such as integrating the preselected sequencing SNPs into a standard LD chip or a standard 54 K chip.

conclusions

For contemporary animals, the imputation first to the 54 K chip and then to 54 K + DFS + FRA SNPs achieved highest imputation accuracy and SNPs with lower MAF tended to have lower correlations but higher concordance rates. Given the limited imputation accuracy especially for SNPs with low MAF, it would be a good strategy to directly genotype the preselected sequencing SNPs using a customized SNP chip. For old bulls, high imputation accuracy was achieved for SNPs with all MAF ranges, and therefore it would be unnecessary to re-genotype preselected sequencing SNPs.

a54K=

  imputation from the standard LD chip to 54 K; DFS SNPs = SNPs selected from major dairy breeds in Denmark-Finland-Sweden; FRA SNPs = SNPs selected from major dairy breeds in France. b MAF = minor allele frequency. c LD + DFS + FRA = imputation to standard LD chip together with preselected DFS and FRA SNPs. d One-step 54 K + DFS + FRA = directly imputed to 54 K together with preselected DFS and FRA SNPs. e Two-step 54 K + DFS + FRA = first imputed to 54 K and then to 54 K together with preselected DFS and FRA SNPs.

Figure 1 .

 1 Figure 1. The accuracy of imputation using the two-step imputation in relation to MAF, for imputation from the standard LD chip to 54 K, and for imputation of DFS SNPs 2 and FRA SNPs 3 in contemporary animals.

aMAF500=

  randomly removed 500 SNPs from the standard LD chip according to minor allele frequency. b Random500 = randomly removed 500 SNPs from the standard LD chip without considering minor allele frequency.

Table 1 .

 1 Number of SNPs to be imputed and imputation accuracy (%) for contemporary animals.

				Correlation			Concordance rate
	SNP a	N	MAF b	LD + DFS + FRA c	One-step 54 K + DFS + FRA d	Two-step 54 K + DFS + FRA e	LD + DFS + FRA	One-step 54 K + DFS + FRA	Two-step 54 K + DFS + FRA
		40,260 >0.00	-	94.4	95.4	-	98.0	98.3
	54 K	34,342 > 0.01	-	95.4	96.2	-	97.7	98.0
		30,494 > 0.05	-	95.9	96.6	-	97.5	97.9
		1,681	>0.00	89.3	89.1	89.5	97.6	97.4	97.6
	DFS	1,329	> 0.01	91.4	90.8	91.8	97.0	96.8	97.0
		1,215	> 0.05	94.0	93.4	94.3	96.9	96.7	97.0
		4,144	>0.00	86.1	85.0	86.5	95.9	95.7	96.0
	FRA	2,837	> 0.01	87.3	87.0	88.0	94.0	93.7	94.2
		2,469	> 0.05	89.0	88.7	89.5	93.5	93.2	93.8

Table 2 .

 2 Percentages of SNPs (%) in different MAF a groups for SNPs in a standard LD chip or DFS SNPs plus FRA SNPs b . a MAF = minor allele frequency. b DFS SNPs plus FRA SNPs = SNPs selected from major dairy breeds in Denmark-Finland-Sweden and France.

	MAF	Standard LD chip	DFS + FRA
	0.01-0.05	4			9
	0.05-0.10	6			10
	0.10-0.20	15			22
	0.20-0.30	21			22
	0.30-0.40	26			18
	0.40-0.50	28			19
		No. of	No. of			Concordance
	Scenario	bulls	daughters Correlation	rate
	MAF500 a	5,021	-	97.3 (0.3)	98.6 (0.1)
		4,764	0	97.2 (0.2)	98.5 (0.1)
		144	>10	98.5 (0.2)	99.6 (0.0)
		45	>50	99.2 (0.2)	99.8 (0.0)
		26	>100	99.7 (0.1)	99.9 (0.0)
		16	>200	99.8 (0.1)	99.9 (0.0)
	Random500 b 5,021	-	97.7 (0.1)	98.4 (0.0)
		4,764	0	97.6 (0.1)	98.4 (0.1)
		144	>10	99.1 (0.2)	99.7 (0.0)
		45	>50	99.5 (0.1)	99.8 (0.0)
		26	>100	99.8 (0.0)	99.9 (0.0)
		16	>200	99.9 (0.0)	99.9 (0.0)

Table 3 .

 3 Imputation accuracy (%) for 500 SNPs randomly removed from the standard LD chip in bulls.

Scientific RepoRtS |(2020) 10:9524 | https://doi.org/10.1038/s41598-020-66523-7

Acknowledgements

We thank Viking Genetics (Assentoft, Denmark), the Danish Cattle Federation, and Nordic Cattle Genetic Evaluation (Aarhus, Denmark) for providing genotype and pedigree. We thank Bernt Guldbrandtsen (Aarhus University, Aarhus, Denmark), Goutam Sahana (Aarhus University, Aarhus, Denmark), and Ulrik Sander Nielsen (SEGES, Aarhus, Denmark) for generous helps on data preparation. This work is supported by the "Genomics in herds" research project financed by VikingGenetics (Randers, Denmark) and Nordic Cattle Genetic Evaluation (Aarhus, Denmark), and the "Genomic Selection in Animals and Plants (GenSAP)" research project financed by the Danish Council of Strategic Research (Aarhus, Denmark). Aoxing Liu acknowledges the scholarship provided by the China Scholarship Council (CSC). Yachun Wang is financially supported by the Program for Changjiang Scholar and Innovation Research Team in University (IRT_15R62).

Data availability

The datasets used in the present study are available from the corresponding author on reasonable request.

Author contributions

Aoxing Liu conceived the study, analyzed the data, and wrote the manuscript. Guosheng Su conceived and supervised the study, made substantial contributions to the interpretation of the results, and revised the manuscript. Mogens Sandø Lund supervised the study and made substantial contributions to the interpretation of the results. Didier Boichard provided the information of FRA SNPs and revised the manuscript. Xiaowei Mao www.nature.com/scientificreports www.nature.com/scientificreports/ provided pipelines of imputation, participated in the discussion during data analyses, and revised the manuscript. Emre Karaman participated in the discussion during data analyses and revised the manuscript. Sebastien Fritz provided the information of FRA SNPs. Yachun Wwang supervised the study and revised the manuscript. Gert Pedersen Aamand gave constructive comments and suggestions during data analyses. All the authors read and approved the manuscript.

competing interests

The authors declare no competing interests.