Length of day (LOD) varies throughout Earth history, dictated in the long term by the retreat of the Moon and growth of Earth's core. Reconstruction of LOD in deep time has greatly benefited from analyses of the orbital cycle from geological records. The Precambrian database, however, is very limited. Here we present a cyclostratigraphic study on c. 1.1 Ga laminated argillaceous limestones of the Nanfen Formation from the North China craton. Milankovitch cycles are successfully recognized from the consistent cyclical variations of the high-resolution magnetic susceptibility data measured on drillcores from the Nanfen Formation. Sedimentary cycles with wavelengths of 12.45, 10.72 and 9.96 m in three drillcores are interpreted as 405 kyr eccentricity cycles. An obliquity period of 20.9–21.4 kyr and a precession period of 12.7–16.0 kyr are also identified. The Earth–Moon distance and LOD are thus estimated at 3.43 ± 0.04 × 108 m and 18.94 ± 0.39 h, respectively. Our new results together with the data previously published reveal an important discontinuity between 1.9 and 1.4 Ga in the evolutionary trend of the LOD. We suggest that Earth's core grew rapidly during the interval, which significantly changed the moment of inertia of the Earth–Moon system, and tidal dissipation again ruled the LOD variations after that stage.

Supplementary material: A description of the calculation method, magnetic susceptibility data, anisotropy of magnetic susceptibility data and supplementary figures are available at https://doi.org/10.6084/m9.figshare.c.6086544

Length of day (LOD), a function of Earth's rotation rate, varies throughout Earth history. Theoretically, whereas tidal dissipation increases the LOD, the superimposed effect of Earth's core growth could decelerate or reverse that trend by reducing Earth's moment of inertia (Runcorn 1964, 1965; Denis et al. 2011; Hinnov 2013, 2018). The reconstruction of LOD changes is thus of great significance for understanding the evolution of the Earth–Moon system and Earth's core growth. LOD values in deep time can be estimated in a quantitative manner through analyses of palaeontological and sedimentary indicators and the records of orbital forced palaeoclimate changes (e.g. Berger and Loutre 1994; Hinnov 2013; Meyers and Malinverno 2018; Laskar 2020). Numerical models have also tentatively established solutions for variations of the LOD (Berger and Loutre 1994; Laskar et al. 2004; Waltham 2015), which are, to first order, consistent with the geological records in the Phanerozoic (Laskar 2020). However, the data for the LOD in Precambrian time are rather rare, hindering an understanding of its evolution in earlier Earth history (Williams 2000; Zhang et al. 2015; Hinnov 2018; Meyers and Malinverno 2018; Laskar 2020).

Here we present a cyclostratigraphic study using high-resolution magnetic susceptibility (MS) measurements to investigate sedimentary cyclic patterns in the laminated argillaceous limestones of the c. 1.1 Ga Nanfen Fm in three drillcores from the North China craton (NCC). The data provide new constraints on Earth's obliquity and precession periods. This new LOD constraint, together with the data previously published, reveals a discontinuity between c. 1.9 and c. 1.4 Ga in the evolutionary trend of the LOD, which may indicate a period of rapid growth of Earth's core.

The NCC is bounded to the north and west by the late Paleozoic Central Asian Orogenic Belt and the early Paleozoic Qilian Orogenic Belt, respectively, and to the south and east by the Mesozoic Qinling–Dabie Orogenic Belt and the Su–Lu metamorphic belt, respectively (Fig. 1). The basement of the NCC consists of high-grade Archean to Paleoproterozoic metamorphic rocks. After the final stabilization of the NCC in the Paleoproterozoic (e.g. Zhao et al. 2005, 2012), the basement rocks were unconformably overlain by post-1.85 Ga little-deformed volcano-sedimentary successions of Mesoproterozoic to Neoproterozoic age and Phanerozoic cover.

Regional stratigraphy

The studied Nanfen Fm is exposed in the Benxi region in eastern Liaoning Province, northeastern NCC. Late Mesoproterozoic to early Neoproterozoic strata in the region include the Xihe Gp in the lower part of the succession and the Kangjia Fm in the upper part (Fig. 1b). The Xihe Group is subdivided into three formations: in ascending order, the Diaoyutai Fm, Nanfen Fm and Qiaotou Fm (LBGMR 1989). The Nanfen Fm conformably overlies the quartz sandstone of the Diaoyutai Fm and in turn disconformably underlies the Qiaotou Fm (Fig. 2c). The Nanfen Fm, deposited in shelf to littoral environments, is subdivided into three lithological members (LBGMR 1989). The Lower Member is mainly composed of pale blue laminated argillaceous limestones (Figs 3a, b and 4a, b) with greenish-grey to red siltstones and sandy shales near the bottom and purplish red limestones near the top. The Middle Member is dominated by purplish red thin-bedded calcareous mudstones intercalated with pale blue limestones, whereas the Upper Member is characterized by yellowish-green shales and sandstones (Fig. 2c).

Age and palaeomagnetism of the Nanfen Formation

The Diaoyutai, Nanfen and Qiaotou formations have no direct radiometric age constraints thus far. Together, a youngest detrital zircon age peak (1136 ± 11 Ma, n = 33 analyses) from the Diaoyutai Fm (Zhao et al. 2020) and the c. 945–920 Ma mafic sills intruding the Qiaotou Fm (Zhang et al. 2016; Zhao et al. 2020) constrain the depositional age of the Nanfen Fm between c. 1136 and c. 945 Ma (Zhao et al. 2020; Zhang et al. 2021; Fig. 2c). New geochronological results using the authigenic monazite Pb–Pb method from the Xinxing Fm in the Xuzhou region and the Liulaobei Fm in the Huainan region in the southeastern NCC (Zhang et al. 2022), the correlative strata of the Nanfen Fm (Su 2016; Zhao et al. 2020), suggest that the minimum depositional age of the Nanfen Fm is c. 1086 Ma (Zhang et al. 2022), also supporting the depositional age of the Nanfen Fm at c. 1.1 Ga.

The NCC has been controversially placed close to Siberia, Baltica, São Francisco or India at the margin of the supercontinent Rodinia during the late Mesoproterozoic to early Neoproterozoic (e.g. Li et al. 2008; Peng et al. 2011; Zhao et al. 2018; Liu et al. 2020; Merdith et al. 2021). Most of these models are mainly based on geological correlation, but have yet to be tested directly by palaeomagnetism. The palaeomagnetic directions of the pale blue limestone in the Lower Member of the Nanfen Fm are characterized by notably steep inclinations (c. 85°), indicating that the Benxi region of the NCC was then located at high palaeolatitudes (c. 80°; Fig. 2a; Zhao et al. 2020). In the Middle Member of the Nanfen Fm, however, the palaeomagnetic directions are characterized by moderate inclinations (c. 60°, palaeolatitude at c. 41°). A primary origin interpretation for these remanent magnetizations is supported by the presence of magnetic reversals and their distinctiveness from the younger poles of the NCC (Zhao et al. 2020). Palaeomagnetically constrained palaeogeographical reconstructions proposed an NCC–northwestern Laurentia connection in the supercontinent Rodinia during c. 1110–775 Ma on the basis of apparent polar wander path (APWP) comparison and geological evidence (Fig. 2a; Fu et al. 2015; Zhao et al. 2020; Ding et al. 2021; Zhang et al. 2021). In that model, the poles of the Lower and Middle members of the Nanfen Fm overlap with the well-sampled and well-dated c. 1110–1085 Ma portion of the Laurentian APWP, which also indicates that the pale blue limestone unit in the Lower Member is probably c. 1100 Ma in age, which is consistent with the age constraints based on the geochronological analyses mentioned above.

Drillcore records

The present study focuses on the pale blue and purplish red limestones in the Lower Member of the Nanfen Fm in three drillcores from the southern Benxi region, namely ZK301, ZK001 and ZK201 (Fig. 2b and c). The ZK001 and ZK201 wells were drilled in 2016, and the ZK301 well was drilled in 2017. The drillcore hole configuration was designed to be perpendicular to the geoid and the azimuth was not oriented. The drillcores have been stored at the Liaoning Eighth Geological Brigade Limited Liability Company. Because of the gentle regional dip of strata (c. 5°) (Fig. 3) and the high recovery (Table 1), the recovered stratum thickness is approximately equal to the drillcore depth. ZK301 has 100% recovery, covering the upper part of the Lower Member, the pale blue laminated argillaceous limestones (12.80–80.10 m in depth) and purplish red limestones (0.00–12.80 m) (Table 2; Fig. 5). ZK001 has 99% recovery, and penetrated 110.06 m of surface soil (0–4.60 m) and pale blue laminated argillaceous limestones (4.60–110.06 m) (Table 2; Figs 4a, b and 5). Drillcore ZK201 has 95% recovery, covering the lower part of the Lower Member (2.50–85.03 m) and the top of the Diaoyutai Fm (85.03–90.64 m). The pale blue laminated argillaceous limestone unit is from 2.50 to 31.69 m in ZK201 (Table 2; Fig. 5).

The three drillcores have covered the entire pale blue laminated argillaceous limestone portion. There is no obvious unconformity or erosional surface observed in the outcrops or the drillcores of the Nanfen Fm. Considering the short distance between the well locations, the stable bedding attitude and the lack of fault structures in the drilling area, the correlation between drillcores has proved straightforward (Fig. 2; Table 1). Based on the correlation of lithological marker beds and the drillcore parameters including well location, elevation, drillcore depth, drillcore interval of the pale blue limestones and regional bedding attitude (Table 1), the stratigraphic overlap ranges in drillcore depth of the pale blue limestone units between ZK301 and ZK001 and between ZK001 and ZK201 are estimated at c. 3 and c. 25 m, respectively. These expected overlaps are consistent with the results of the MS series analyses (see the next section; Table 1; Fig. 5).

Standard 2.5 cm diameter cylindrical core samples were drilled perpendicular to the axis of each drillcore (Fig. 3c and d), at an average sampling spacing of c. 7.5 cm, yielding a total of 2640 samples covering the entire pale blue limestone unit (Fig. 5). Each sample was cut into two 2.2 cm long specimens using a nonmagnetic dual-blade rock saw. One specimen from each sample was selected for rock magnetic measurements. All the specimens were subjected to MS measurements and 963 of them were selected at every c. 20 cm for anisotropy of magnetic susceptibility (AMS) measurements (Supplementary material Tables S1 and S2). Both the MS and AMS were measured using an AGICO MFK1-FA Kappa-bridge (with a sensitivity of 2 × 10−8 SI, <10−11 m3 kg−1) at the Palaeomagnetism and Environmental Magnetism Laboratory of the China University of Geosciences (Beijing). The volume susceptibility data were normalized by mass for further processing and analysis.

MS is one of the most successful proxies used in cyclostratigraphic analyses (Kodama and Hinnov 2015). MS variations are widely used to determine the astronomical cycles preserved in Proterozoic sedimentary rocks (e.g. Zhang et al. 2015; Bao et al. 2018; Mitchell et al. 2021). Several lines of evidence support the MS data for the Nanfen Fm in the drillcores being reliable indicators of depositional conditions and suitable for cyclostratigraphic analyses. First, the argillaceous limestones of the Nanfen Fm preserve thin laminations without obvious post-depositional disruption (Figs 3 and 4a, b). Second, no mineral indicative of significant metamorphism is identified through petrographic micrograph observations. Third, the AMS measurements, both for the palaeomagnetic oriented samples from outcrops (Fig. 4c and d; Zhao et al. 2020) and those from the drillcores (Fig. 4e and f), strongly suggest that the rocks of the Nanfen Fm still preserve the sedimentary magnetic fabric in the region. Fourth, the palaeomagnetic signals isolated from the Nanfen Fm have proved primary in origin (Zhao et al. 2020).

The MS series were linearly interpolated, and then the data for ZK301 and ZK001 were detrended by subtracting smooth curves using Robust Locally Weighted Regression in the software KaleidaGraph™ (Cleveland 1979). The maximum and minimum of MS in ZK301 are 9.914 × 10−8 m3 kg−1 and 7.679 × 10−9 m3 kg−1, respectively, with an average value of 3.453 × 10−8 m3 kg−1. The average sampling resolution is 8 cm and the minimum sampling resolution is 3 m. The MS series is linearly interpolated every 5 cm and detrended with a smoothing factor of f = 0.35 (Fig. 6). The maximum MS of ZK001 is 1.0616 × 10−7 m3 kg−1 and the minimum MS value is 1.118 × 10−8 m3 kg−1, with an average value of 3.9975 × 10−8 m3 kg−1. The average sampling resolution is 7.4 cm and the minimum sampling resolution is 3 cm. The MS series are linearly interpolated and then detrended with an f = 0.20 (Fig. 6). The maximum and minimum of MS in ZK201 are 7.8188 × 10−8 m3 kg−1 and −3.3822 × 10−9 m3 kg−1, respectively. The average MS value and the average sampling resolution are 4.0069 × 10−8 m3 kg−1 and 7.5 cm, respectively. The MS series of ZK201 are linearly interpolated every 5 cm. The data in the three drillcores are finally prewhitened and normalized.

Periodic cycles in the stratigraphic domain

Multitaper method (MTM) power spectra are conducted using the Matlab software with red noise modelling at 85, 90, 95 and 99% confidence levels (Mann and Lees 1996; Fig. 7). Comparing the significant frequency ratios in a spectrum with Milankovitch frequency ratios is a common method to test whether the observed cycles in strata were formed under orbital forcing (Weedon 2003). Wavelet analyses are used in cyclostratigraphy to check the stability of the depositional rate and orbital forcing through the studied strata (Grinsted et al. 2004).

MTM power spectral analysis (Fig. 7) of the pre-processed MS series from ZK301 reveals cycles with 12.45 m, 60 cm, 45 cm and 31 cm wavelengths above 99% confidence level. MS series from ZK001 reveals 10.72 m, 47 cm and 21 cm wavelengths above 99% confidence level. For the MS series of ZK201, MTM power spectral analysis reveals cycles with wavelengths of 9.96, 0.44, 0.32 and 0.21 m above 95% confidence level. Wavelet scalograms of ZK301, ZK001 and ZK201 reveal relatively continuous cycles with periods of comparable spectral bands with power spectra (Fig. 7).

The following periodic features are notable in the spectra. (1) The wavelengths of 12.45 m in ZK301, 10.72 m in ZK001 and 9.96 m in ZK201 (signal A), and the wavelengths of 0.60 m in ZK301, 0.47 m in ZK001 and 0.44 m in ZK201 (signal B), are very significant in the three drillcores. Also, the wavelength ratios of signal A and signal B are fairly constant, at 20.75, 22.80 and 22.60, in the three drillcores ZK301, ZK001 and ZK201, respectively. (2) In the spectra of the ZK301 and ZK201 MS series, the wavelength of 0.45 m in ZK301 and 0.32 m in ZK201 (signal C) has been identified, and the wavelength ratios of signal B and signal C in the two cores are 1.33 and 1.37, respectively. The spectra of three drillcores have a similar hierarchy of wavelengths.

Identification of the orbital periodic signals

Milankovitch cycles of 405 kyr eccentricity are thought to be the most stable during Earth history (Laskar et al. 2004), whereas both obliquity and precession periods varied (Berger and Loutre 1994; Laskar et al. 2004; Waltham 2015). Obliquity and precession periods at c. 1.1 Ga are estimated to be 22.27 kyr and 13.65–15.94 kyr, respectively, from the numerical model (Table 3; Waltham 2015).

If signal A represents the c. 405 kyr eccentricity period, then signal B and signal C are within the range of the obliquity and precession periods, respectively. The power spectra of the tuned MS series show periodic signals of 20.9–21.4 kyr and 12.7–16.0 kyr, which are consistent with the obliquity and precession periods at c. 1.1 Ga, respectively (Fig. 8; Table 3; Waltham 2015). Wavelet analysis results further support the stability of the interpreted astronomical signals (Fig. 8). Bandpass filters of the presumed short eccentricity signal show amplitude modulation of the bandpass filters of the presumed precession signal (Supplementary material Fig. S1), further proving that the periods are Milankovitch cycles.

Complete depositional records of the laminated pale blue limestones in the Lower Member of the Nanfen Fm are recovered from the three drillcores. The duration of the pale blue limestone unit in cores ZK301, ZK001 and ZK201 is 3.167 myr, 4.865 myr and 1.511 myr, respectively (Fig. 9). Overlapping strata between ZK301 and ZK001, and between ZK001 and ZK201, are identified through the correlation of lithological marker beds, the drillcore locations, elevation, drillcore depth and bedding attitude. Field comparison can be further calibrated by smooth curves, filtering curves and box plots (Table 1; Fig. 5). Consistency of smooth curves in the base of the drillcore ZK301 and in the top of ZK001 suggests that 77.35–80.10 m in ZK301 corresponds to 4.62–7.42 m in ZK001, also meaning that 79.40 m of ZK301 corresponds to 6.77 m of ZK001 and the stratigraphic overlap range in core depth is c. 2.8 m (Fig. 5). Moreover, on the basis of the consistency of the box plot and filtering curve of the eccentricity cycle of 405 kyr, 85.82–110.06 m in the lower part of ZK001 corresponds to 3.83–27.63 m in the upper part of ZK201. Datum lines are set on the peak of E11 of ZK001 (108.10 m) and on the peak of E2 of ZK201 (26.18 m) (Fig. 5). The stratigraphic overlap range in the drillcore depth is c. 25.1 m.

Taking the overlapped intervals into account, the duration of deposition of the pale blue limestone is calculated at 8.116 myr in the region. The average sedimentation rate of the pale blue limestone unit is estimated at c. 2.1 cm ka−1, which is within the range of carbonate sedimentation rate (0.5–20 cm ka−1, Fairchild et al. 2016). Correlation coefficient (COCO) analyses (Li et al. 2018, 2019) on the three drillcores show significant peaks at c. 2 cm ka−1, consistent with the average sediment rate revealed in MTM spectral analyses and the wavelet spectra (Supplementary material Fig. S2).

Based on the obliquity and precession periods of 20.9–21.4 and 12.7–16.0 kyr derived from the Nanfen Fm, the Earth–Moon distance and LOD at c. 1.1 Ga are constrained at 3.43 ± 0.04 × 108 m and 18.94 ± 0.39 h, respectively (Fig. 10), through the use of Poisson equations (see the Supplementary material; Berger and Loutre 1994; Laskar et al. 2004; Laskar et al. 2011). The new data are consistent with the estimated value of W15 (Waltham 2015), supporting that the W15 model is probably effective for the period since at least 1.1 Ga (Fig. 11). Furthermore, the cyclostratigraphic results from the rhythmites of the Xiamaling Fm in the NCC yield an LOD of 18.68 ± 0.25 h at c. 1.4 Ga (Zhang et al. 2015; Meyers and Malinverno 2018) and the coeval strata of the Roper Group in Australia yield similar results (Mitchell et al. 2021). The confidence interval of LOD revealed from the Xiamaling Fm partially overlaps with the W15 model, suggesting that the W15 model is probably valid at c. 1.4 Ga. New cyclostratigraphic results from the Gaoyuzhuang Fm in the NCC reveal an LOD of 19.16 ± 0.7 h at c. 1.56 Ga (Liu et al. 2022). The range of LOD is significantly higher than the predicted value of the W15 model (Fig. 11), implying that the W15 model may be invalid at 1.56 Ga. The new LOD data seem not to conform to the hypothesis that Earth entered a resonant state keeping a stable LOD of c. 21 h throughout the so-called Boring Billion, between c. 2.0 and 0.6 Ga (Bartlett and Stevenson 2016; Klatt et al. 2021).

The data for LOD derived from the Paleoproterozoic (c. 2.6–1.9 Ga) geological records, however, are inconsistent with the W15 model prediction, with a 2.2 h difference (Fig. 11; Table 4). The cyclostratigraphic datum at c. 2.65 Ga is from carbonate cycles of the Cheshire Fm in Zimbabwe, for which an average precession period of c. 11.6 kyr is presented (Hoffmann et al. 2004). Tidalites datapoints at c. 2.5 Ga are from laminae sequences of the banded iron formation (BIF) in the Weeli Wolli Fm of Australia (Walker and Zahnle 1986; Williams 1989) and the cyclostratigraphic datapoint at c. 2.5 Ga is from BIF of the Dales Gorge Member of the Brockman Iron Fm in Australia (de Oliveira Carvalho Rodrigues et al. 2019). The stromatolite datapoints of c. 2.0 Ga and c. 1.88 Ga are from the Great Slave Supergroup, Gunflint Fm of Canada and the Biwabik Fm of America (Pannella 1972a, 41b; Fralick et al. 2002). We speculate that the offset between geological data and the W15 model for c. 2.6–1.9 Ga is probably caused by Earth interior geodynamic processes, probably either inner core nucleation or the further segregation of Earth's (outer) core.

Most estimates for the age of inner core nucleation (ICN) are during the Proterozoic based on analyses of the palaeointensity variation of the Earth's magnetic field (e.g. Biggin et al. 2015; Bono et al. 2019; Veikkolainen and Pesonen 2021). The density difference between the inner core (12.394 × 103 kg m−3) and the outer core (10.746 × 103 kg m−3) is small (Denis et al. 1997), and the influence of ICN on LOD is estimated at only c. 0.02 h using the simple Earth model (Denis et al. 2011), demonstrating that ICN has very little influence on the evolutionary trend of LOD.

We thus next consider further segregation of Earth's outer core. The total change of LOD caused by Earth's core growth can be estimated through use of a simple Earth model (equations 13; Denis et al. 2006, 2011), and is up to 3.6–3.8 h. This calculation is consistent with the previous estimates (Runcorn 1964, 1965; Birch 1965):
LODLOD=II
(1)
I=8π15[ρcr5+ρ(R5r5)]
(2)
ρ=M(4/3)πρcr3(4/3)π(R3r3)
(3)
where I is the moment of inertia of Earth, representing the rest of Earth's density, LOD′ and I′ are changes of LOD and I, R and r are the radii of Earth and Earth's core, respectively, with density of Earth's core ρ = 10.817 × 103 kg m−3 and mass of Earth M = 5.9737 × 1024 kg. If the size of Earth's core at 1.4 Ga is as large as it is at present, the average c. 2.2 h offset at 2.6–1.9 Ga means 1.0 × 106 m growth of the radius of Earth's core.

As a consequence, the c. 2.2 h offset of LOD may indicate a rapid growth of Earth's core between c. 1.9 and c. 1.4 Ga, and growth in the size of the radius from c. 2.5 × 106 to c. 3.5 × 106 m. Because the W15 model based on tidal dissipation fits well with geological data from c. 1.4 Ga to the present, the core growth may have not played a significant role in the LOD variations. The variation of LOD at 2.6–1.9 Ga significantly has a similar trend to the LOD curve in W15 (Fig. 11), implying that tidal dissipation also dominated LOD variations in the Paleoproterozoic. We thus suggest that the Earth's core growth significantly changed the moment of inertia of the Earth–Moon system between c. 1.9 and c. 1.4 Ga, whereas tidal dissipation dominated the LOD variations at 2.6–1.9 Ga and from c. 1.4 Ga to the present.

  1. Milankovitch cycles are identified from the c. 1.1 Ga high-latitudinal Nanfen carbonate successions exposed in the NCC, indicating Earth's obliquity and precession periods of 20.9–21.4 kyr and 12.7–16.0 kyr, respectively, at this late Mesoproterozoic age.

  2. The Earth–Moon distance and the LOD at c. 1.1 Ga are constrained at 3.43 ± 0.04 × 108 m and 18.94 ± 0.39 h, respectively.

  3. Earth's core has probably undergone a rapid growth between c. 1.9 and c. 1.4 Ga, which can explain the offset between proxy data and the W15 model. The radius increment is estimated to be up to 1.0 × 106 m.

We appreciate discussions with G. E. Williams and K. P. Kodama. Helpful reviews from R. Mitchell and an anonymous reviewer, and editorial comments from W. Xiao improved the paper. The authors thank K. Yang, J. Ding and W. Ren for assistance in the field. We thank Z. Wang, W. Wang, G. Lyu and C. Gong for their support in the sampling process.

XB: conceptualization (equal), data curation (equal), formal analysis (lead), visualization (lead), writing – original draft (lead), writing – review & editing (equal); HZ: data curation (equal), investigation (lead), resources (equal), writing – original draft (equal), writing – review & editing (equal); SZ: conceptualization (lead), project administration (lead), resources (equal), supervision (lead), writing – original draft (equal), writing – review & editing (lead); XL: investigation (equal), resources (equal); WT: investigation (equal), resources (equal); CL: investigation (equal); HW: writing – review & editing (supporting); HL: writing – review & editing (supporting); TY: writing – review & editing (supporting)

This work was funded by the National Natural Science Foundation of China (grants 41888101 and 41830215) and Chinese ‘111’ project B20011.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All data generated or analysed during this study are included in this published article (and its supplementary information files).

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)