High-affinity P2Y2 and low-affinity P2X7 receptor interaction modulates ATP-mediated calcium signaling in murine osteoblasts

The P2 purinergic receptor family implicated in many physiological processes, including neurotransmission, mechanical adaptation and inflammation, consists of ATP-gated non-specific cation channels P2XRs and G-protein coupled receptors P2YRs. Different cells, including bone forming osteoblasts, express multiple P2 receptors; however, how P2X and P2Y receptors interact in generating cellular responses to various doses of [ATP] remains poorly understood. Using primary bone marrow and compact bone derived osteoblasts and BMP2-expressing C2C12 osteoblastic cells, we demonstrated conserved features in the P2-mediated Ca2+ responses to ATP, including a transition of Ca2+ response signatures from transient at low [ATP] to oscillatory at moderate [ATP], and back to transient at high [ATP], and a non-monotonic changes in the response magnitudes which exhibited two troughs at 10−4 and 10−2 M [ATP]. We identified P2Y2 and P2X7 receptors as predominantly contributing to these responses and constructed a mathematical model of P2Y2R-induced inositol trisphosphate (IP3) mediated Ca2+ release coupled to a Markov model of P2X7R dynamics to study this system. Model predictions were validated using parental and CRISPR/Cas9-generated P2Y2 and P2Y7 knockouts in osteoblastic C2C12-BMP cells. Activation of P2Y2 by progressively increasing [ATP] induced a transition from transient to oscillatory to transient Ca2+ responses due to the biphasic nature of IP3Rs and the interaction of SERCA pumps with IP3Rs. At high [ATP], activation of P2X7R modulated the response magnitudes through an interplay between the biphasic nature of IP3Rs and the desensitization kinetics of P2X7Rs. Moreover, we found that P2Y2 activity may alter the kinetics of P2X7 towards favouring naïve state activation. Finally, we demonstrated the functional consequences of lacking P2Y2 or P2X7 in osteoblast mechanotransduction. This study thus provides important insights into the biophysical mechanisms underlying ATP-dependent Ca2+ response signatures, which are important in mediating bone mechanoadaptation.

1 Abstract P2 purinergic receptor family implicated in many physiological processes, including neurotransmission, mechanical adaptation and inflammation, consist of ATP-gated non-specific cation channels P2XRs and G-protein coupled receptors P2YRs. Different cells, including bone forming osteoblasts, express multiple P2 receptors; however, how P2X and P2Y receptors interact in generating cellular responses to various doses of [ATP] remains poorly understood. Using primary bone marrow and compact bone derived osteoblasts and BMP2-expressing C2C12 osteoblastic cells, we demonstrated conserved features in the P2-mediated Ca 2+ responses to ATP, including a transition of Ca 2+ response signatures from transient at low [ATP] to oscillatory at moderate [ATP], and back to transient at high [ATP], and a non-monotonic changes in the response magnitudes which exhibited two troughs at 10 -4 and 10 -2 M [ATP]. We identified P2Y2 and P2X7 receptors as predominantly contributing to these responses, and constructed a mathematical model of P2Y2R-induced inositol trisphosphate (IP3) mediated Ca 2+ release coupled to a Markov model of P2X7R dynamics to study this system. Model predictions were validated using parental and CRISPR/Cas9-generated P2Y2 and P2Y7 knockouts in osteoblastic C2C12-BMP cells. Activation of P2Y2 by progressively increasing [ATP] induced a transition from transient to oscillatory to transient Ca 2+ responses due to the biphasic nature of IP3Rs and the interaction of SERCA pumps with IP3Rs. At high [ATP], activation of P2X7R modulated the response magnitudes through an interplay between the biphasic nature of IP 3 Rs and the desensitization kinetics of P2X7Rs.
Moreover, we found that P2Y2 activity may alter the kinetics of P2X7 towards favouring naïve state activation. Finally, we demonstrated the functional consequences of lacking P2Y2 or P2X7 in osteoblast mechanitransduction. This study thus provides important insights into the biophysical mechanisms underlying ATP-dependent Ca 2+ response signatures, which are important in mediating bone mechanoadaptation.

ATP-mediated P2R Ca 2+ responses in murine osteoblasts
ATP-stimulated P2R Ca 2+ responses were assessed in three independent murine osteoblasts models: BMP2-transfected C2C12 osteoblastic cells (C2-OB), bone-marrow-derived osteoblasts (BM-OB), and compact-bone-derived osteoblasts (CB-OB). Osteoblasts were loaded with Ca 2+indicator dye Fura2, stimulated with varying doses of ATP, and changes in [Ca 2+ ] i were recorded using live cell fluorescent microscopy (Fig. 1). Qualitatively, the recorded Ca 2+ response timeseries signatures demonstrated a general trend of exhibiting transient single-peaked responses at low [ATP], multi-peaked oscillatory responses at mid-range [ATP], and relatively sustained single/multi-peaked response at high [ATP] (Fig. 1A). The Ca 2+ responses at each [ATP] were analysed by quantifying several parameters, including overall response magnitudes and activation times, as well as oscillatory fractions, magnitudes, periods and peaks (see Table S1 for definitions) (Mackay et al., 2016). Similar to previous study (S. Xing et al., 2016), we harmonized doseresponse profiles across osteoblast models by first aligning the responses along the dose-axis to match troughs/peaks, followed by rescaling the responses to a [0,1] interval (Fig. S1). Such alignment allowed us to account for (i) inconsistencies in ATP solution preparations between experiments and (ii) varied dose-sensitivities across cell lines. Calcium responses induced by low [ATP] (<10 -7 M) were consistently associated with low response magnitudes and slow activation kinetics (Fig. 1A, left two columns), with little to no oscillatory component (Fig. 1B). Increasing [ATP] further induced responses with faster activation kinetics and higher magnitudes (Fig. 1A, middle two columns). This also coincided with the emergence of high frequency oscillations (~10-20 s periods; Fig. 1B, Fig S1B, C) which peaked at ~10 -5 M ATP stimulation. Notably, the oscillatory peak did not coincide with the peak magnitude. Instead, as cells were stimulated with higher [ATP], the oscillatory component began to diminish, exhibiting lower frequency oscillations and fewer oscillatory peaks, while response magnitude continued to increase, peaking at ~2×10 -4 M [ATP] (Fig. 1A, right two columns, Fig. 1B). For [ATP] >2×10 -4 M, the response magnitude decreased with increasing [ATP] in all osteoblast lines. Thus, in all osteoblast models, the intracellular Ca 2+ response to ATP shifts with increase in [ATP] from a transient with a single narrow peak, to oscillatory and back to transient with a pronounced wide peak.  Fura2-loaded BMP2-transfected C2C12 osteoblastic cells (C2-OB), bone-marrow-derived osteoblasts (BM-OB), or compact-bone-derived osteoblasts (CB-OB) were stimulated with ATP (10 -8 to 10 -3 M), changes in [Ca 2+ ] i were recorded, and characteristic parameters of individual celllevel Ca 2+ responses were quantified. (A) Representative ATP-induced Ca 2+ signature responses for different osteoblastic lines. Recording duration: 120 s. (B) Activation time, magnitude and oscillatory characteristics of Ca 2+ responses in different osteoblastic cells were aligned to obtain consensus on dose-dependency behaviours (see Fig. S1 for intermediary alignment steps). Data are response means, normalized to peak dose-response. Solid curves: Loess curves fit to normalized response means. Vertical solid lines: peak magnitude; Vertical dashed lines: peak oscillatory activity. CI: confidence interval, M: molar concentration.

P2Y2 and P2X7 receptors orchestrate the ATP-mediated Ca 2+ responses
To examine which P2 receptors contribute to the ATP-induced Ca 2+ responses, we first examined their expression in osteoblastic cells of different origin. Among the P2Y family, P2ry2, P2ry4, P2ry12, and P2ry14 transcripts were detected in all osteoblastic cells by RT-qPCR ( Fig. 2A top).
Among the P2X family, P2rx4 and P2rx7 transcripts were the most abundantly expressed in all osteoblastic cells ( Fig. 2A bottom). These expression profiles suggest that P2Y2, 4, 12, 14 and P2X4, 7 are the predominant P2 receptor subtypes expressed in osteoblastic cells, among which P2Y2 and P2X7 were the most abundant transcripts. To confirm that P2Y2 and P2X7 receptors are functional, we stimulated Fura2-loaded C2-OB cells with ATP and receptor-specific agonists: the P2Y2-agonist UTP and P2X7-agonist bzATP (Fig. 2B). Consistent with previously characterized P2 receptor sensitivities (S. Xing et al., 2016), we found that the estimated EC50s were 1.0 µM for ATP, 2.8 µM for UTP and 26.4 µM for bzATP in C2-OB cells. Importantly, the oscillatory Ca 2+ responses evoked by 10 -6 M ATP were recapitulated following stimulation with 10 -6 M UTP, and similarly, the sustained responses evoked by 1 mM ATP were observed following 1 mM bzATP stimulation, suggesting that P2Y2 and P2X7 receptors dominate the responses to lower and higher [ATP], respectively (Fig. 2C). Using CRISPR-Cas9 double-nickase constructs, we generated clonal C2-OB cells harboring mutations in P2ry2 (P2ry2Δ) or P2rx7 (P2rx7Δ) (Fig.   2D) to further investigate the independent contribution of P2Y2 and P2X7 to the P2-mediated Ca 2+ responses (presented in subsequent sections).   Li & Rinzel, 1994) with a Markov model of P2X7R kinetics adapted from (Khadra et al., 2013) (Fig. 3). The cell was divided into two compartments ( Fig. 3A) through the receptor channels generated by a 12-state Markov P2X7R sub-model ( Fig. 3B) (Khadra et al., 2013). The P2X7R submodel assumes that each receptor has three ATP binding sites, two of which must be occupied for the receptor to be open, and that each state represents the fraction of receptors with a given number of occupied ATP-binding sites (Fig. 3B, solid circles).
The closed, , and desensitized, , states are non-conducting, whereas the open states ( = 1 -4) possess the same conductance 7 . The states were divided into three rows corresponding to desensitized (Fig. 3B, top row), naïve (Fig. 3B, middle row) and sensitized or primed (Fig. 3B, bottom row) states, respectively. The naïve row is comprised of the states 1 , 2 , 1 , 2 that have not been exposed to ATP for a prolonged period of time, whereas the sensitized and desensitized rows are comprised of the states 3 , 4 , 3 , 4 ( 1 , 2 , 3 , 4 ) that have been previously exposed to ATP. The forward and backward transitions along each row describes ATP binding and unbinding, respectively, whereas downward and upward transitions between the rows represent receptor sensitization (middle to bottom row), desensitization (middle to top row) or recovery (bottom/top to middle row). The rate of desensitization increases as more ATP molecules bind to P2X7R and the open probability along the sensitized row is larger than that for the naïve row.
Combining the two submodels together produced the following model: where and represent the fraction of free Ca 2+ in the cytosol and ER, respectively, as a result of buffering with 0  < << 1, is the ratio between cytosolic and ER volume, is the maximum rate of IP 3 production by P2Y2 in response to ATP, is the half-maximum production of IP 3 through P2Y2R and is the rate of IP 3 degradation.  Using the model parameters provided in Table 1, we simulated Ca 2+ responses to different [ATP] in three specific cases: (i) in naïve cells expressing both P2Y2R and P2X7R using the full model ( Fig. 4A, blue curves), (ii) in cells that do not express P2Y2R, using the submodel for P2X7 component only (Fig. 4A, grey curves), and (iii) in cells that do not express P2X7R using the submodel for P2Y2 component only (Fig. 4A, yellow curves). The simulated Ca 2+ responses were compared with those obtained experimentally in WT (Fig. 4B, blue curves), P2ry2Δ (Fig. 4B, grey curves) and P2rx7Δ (Fig. 4B, yellow curves) C2-OB cells. As shown, the responses to low [ATP] were predominantly P2Y2-mediated, while the response to high [ATP] were jointly mediated by P2Y2 and P2X7. Notably, the characteristic two-peaked response to 10 -3 M observed in experimental recordings (Fig. 1A, Fig. 4B) was predicted by the full model (Fig 4A) of WT cells, but was abolished in P2ry2Δ and P2rx7Δ recordings and in simulations of P2X7 and P2Y2 submodels. These data strongly support the interaction between P2Y2 and P2X7 receptors in generating this unique signaling feature.
One interesting aspect of the recordings and simulations displayed in    with IP3Rs (Y. Li & Rinzel, 1994); this results in a high magnitude semi-persistent response ( Fig.   5D, blue curve). The aforementioned mechanism suggests that the heterogeneity in Ca 2+ response profiles at a given [ATP] observed experimentally ( Fig. 1A)   P2rx7Δ C2-OB cells stimulated with ATP ( Fig. 6A, red circles) and in the model lacking P2X7Rs ( Fig. 6A, red curve), strongly implicating P2X7 in this phenomenon. Therefore, we next plotted the total Ca 2+ flux through P2X7Rs, estimated as the area under the P2X7R flux curve, as well as the maximal Ca 2+ fluxes through P2X7Rs and IP 3 Rs predicted by the model (Fig. 6B-E). When P2X7R-mediated Ca 2+ entry became evident at 10 -5 M ATP (Fig. 6B, light grey region), the maximum flux through IP 3 Rs in the full model (Fig. 6C, black curve) dropped below that of the P2X7R-lacking submodel (Fig. 6C, red curve) due to the biphasic dependence of IP 3 Rs on [Ca 2+ ] i , resulting in the first trough in the dose-response (Fig. 6A, light grey region). At high 10 -2 M [ATP], on the other hand, the time required for the P2X7R flux to decay to half of its maximum (t 1/2 ) decreased ( Fig. 6D). As a result, despite the maximum P2X7R flux continually increasing (Fig.   6E), the Ca 2+ entering through P2X7R began to decrease at elevated [ATP] (Fig. 6B, dark grey), resulting in the second trough in the dose-response (Fig. 6A, dark grey). Taken together, these simulations indicate that the non-monotonic Ca 2+ dose-response to ATP is driven by an interplay between the biphasic nature of IP 3 Rs and the desensitization kinetics of P2X7Rs.  Next, we examined why the Ca 2+ response to high [ATP] is dramatically affected in P2ry2Δ cells ( Fig. 4B and Fig. 7A). We examined the magnitudes of Ca 2+ responses to 10 -2 M [ATP] in WT P2ry2Δ cells, and found that in the absence of P2Y2R, response magnitudes exhibited a distinct bimodal distribution with one cluster of responses similar to those in WT, and another one with much higher response magnitudes (Fig. 7B). We hypothesized that the bimodality in the P2X7Rmediated responses is due to P2X7Rs being in the naïve or sensitized initial state (Yan et al., 2010).
To verify this, we used the P2X7 submodel to simulate the Ca 2+ response in two different scenarios: initiating the P2X7-simulations from the naïve closed state 1 (i.e., 1 (0) = 1, closely matching the mean of the right mode of the distribution (Fig. 7B). These data suggest that P2Y2 activation may alter the kinetics of P2X7 towards favouring naïve state activation.

Functional contributions of P2Y2 and P2X7 to mechanotransductive signaling
To investigate the potential functional consequences of the complex interactions between P2Y2 and P2X7, we examined how the absence of each of these receptors affects ATP-mediated mechanotransduction. We have previously shown that mechanical stimulation of a single "primary" osteoblast with a glass micropipette leads to the release of 10 -5 to 10 -4 M ATP into the pericellular space, which then diffuses to stimulate neighbouring non-mechanically perturbed "secondary" cells (Mikolajewicz et al., 2019;Mikolajewicz, Zimmermann, et al., 2018). Here, we mechanically stimulated a single fura2-loaded osteoblast from parental C2-OB, or clones deficient in P2Y2, P2ry2Δ, or P2X7, P2rx7Δ and recorded [Ca 2+ ] i responses in the primary and secondary cells (Fig. 8A). We found that while in P2rx7Δ cells the response was qualitatively similar to WT, in P2ry2Δ cells secondary responses were abolished (Fig. 8B, C). Quantitatively, the primary response was unaffected in P2ry2Δ cells, but exhibited higher magnitude and faster decay in P2rx7Δ cells (Fig. 8C, D). The suppression of the secondary response in P2ry2Δ cells was evident by the reduced signaling radius ( = 9 × 10 -3 ), fractions of responding cells ( = 5 × 10 -11 ), response magnitudes ( = 2 × 10 -4 ) and areas under the curves (AUC; = 0.02) (Fig. 8B-D, yellow). In contrast, in P2rx7Δ cells, the signaling radius and fraction of responding secondary cells was unaffected; however, the response magnitudes and areas under the curves of secondary cells were significantly higher in P2rx7Δ cells compared to parental C2-OB cells (Fig. 8C, D, grey). These data demonstrate that P2Y2 receptor is critical for the secondary responses, consistent with its high sensitivity to ATP. In addition, the contribution of P2X7 to Ca 2+ responses is evident even though extracellular ATP in these experiments remained below [ATP] required to stimulate P2X7. Taken together, these data strongly support the importance of an interplay between P2Y2R and P2X7R.   (Khadra et al., 2012;Yan et al., 2011), suggesting that Ca 2+ release from the ER through the P2Y2 pathways may underly the altered kinetics of P2X7. Taken together, our study demonstrates multiple points of interactions between P2Y2 and P2X7 receptors, which are not only activated at very different ranges of ATP concentration, but also belong to different classes of receptors.
Finally, our study demonstrates that the absence of either P2Y2 or P2X7 has significant implications on ATP-mediated mechanotransduction. We used an experimental setup in which the mechanical stimulation of a single osteoblasts generates a micro-injury in its cell membrane, leading to a release of ATP that signals to neighboring (secondary) cells (Mikolajewicz, Zimmermann, et al., 2018). First, we showed that in the absence of P2Y2, the transmission of ATP signal to neighboring cells is effectively interrupted. These findings are consistent with previous reports that osteoblasts from P2Y2 -/mice exhibited dramatic reduction in fluid flow-induced Ca 2+ responses even though the ATP release was similar to WT osteoblasts (Y. Xing et al., 2014).
Second, we have found that in the absence of P2X7, both primary and secondary responses are significantly altered. While local ATP concentrations at the site of micro-injury may support the involvement of P2X7 in generating the Ca 2+ response of the primary cell (Mikolajewicz et al., 2019;Mikolajewicz, Zimmermann, et al., 2018), the observed changes in the secondary responses are surprising, since we have previously shown that the amount of ATP released in these experiments is below the concentrations required for P2X7 activation (Mikolajewicz, Zimmermann, et al., 2018). Nevertheless, this observation is consistent with previously suggested alterations in mechanotransductive signaling in P2X7 deficient mice (Ke et al., 2003;D. Zeng et al., 2019). Thus our study supports the important role of P2 receptor network in generating a mechanotransductive signal that conveys complex information to neighbouring cells In conclusion, this study provided a complex mechanism of interdependency between the action of high affinity G-protein coupled receptor P2Y2 and a low affinity ligand gated ion channel P2X7.
Using a combination of experimental studies in osteoblastic cells with the full compliment of P2 receptors, as well as osteoblasts deficient in P2Y2 or P2X7, and mathematical modeling of P2Y2Rmediated Ca 2+ release coupled to a Markov model of P2X7R dynamics, allowed us to explore the intricate details of the subcellular signaling induced by ATP in bone forming osteoblasts. The conclusions drawn demonstrated causative links between the exposure to mechanical force, early ATP-mediated signaling, and mechanoadaptive response of bone tissue.

Materials and Methods
Software. transferred into puromycin-free media, allowed 3 days for recovery, and re-plated in 96 well plates at a ~1 cell per well density. After 3 weeks of expansion, half of each single-cell colony was replated in 96-well plates and the other half was collected for genomic DNA extraction using DNeasy kit. Genomic DNA for each single-cell colony was amplified by touchdown PCR using primer sets designed to flank the genomic region targeted by Cas9 (Table S2), and amplicons were separated on a gel to screen for clones with evident band shifts. Selected clones were subsequently validated by immunoblot analysis, and termed P2ry2Δ and P2rx7Δ cells, for P2ry2 and P2rx7, respectively.

Quantitative real-time polymerase chain reaction (qRT-PCR). Total RNA was isolated using
RNeasy kit and QIAshredder columns and reverse transcribed using cDNA reverse transcription kit. Real time qPCR was performed using QuantStudio 7 Flex PCR System, with SYBR Green or TaqMan Master Mix. Primer sequences are provided in Table S2 and cycling conditions in Table   S3.

Intracellular Ca 2+ recordings and analysis.
Cells were plated on glass-bottom 35 mm dishes or 48-well plates (MatTek Corporation), for single-cell mechanical stimulation and agonist application experiments, respectively. Cell were loaded with Fura2-AM for 30 min, acclimatized in physiological solution (PS) for 10 min on the stage of an inverted fluorescence microscope (Nikon T2000), and imaged as described previously (Mikolajewicz, Zimmermann, et al., 2018).
Immunoblotting. Cell lysates were extracted in RIPA lysis buffer and samples were prepared and subject to SDS-PAGE on a 10% (w/v) acrylamide gel as previously described (Mikolajewicz, Zimmermann, et al., 2018). Blotted nitrocellulose membranes were incubated with primary antibodies overnight (1:1000 dilution, 5% BSA in TBST, 4°C) and secondary antibodies were applied for 1 h (1:1000 dilution, 5% BSA in TBST, rt) prior to visualization with chemiluminescence system.

Mechanical-stimulation.
Single osteoblastic cells were stimulated by local membrane indentation with a glass micropipette using a FemtoJet microinjector NI2 (Eppendorf Inc.), as previously described (Mikolajewicz, Zimmermann, et al., 2018). fluxes considered in this model, as described below.

Statistical
(1) Plasma Membrane Ca 2+ Leak ( ). A constant inward leak across cell membrane to ensure that total [Ca 2+ ] within the cell remained positive. It was assumed to be constant (see Table 1).
(2) IP 3 R Ca 2+ Flux ( 3 ): The Ca 2+ flux through IP 3 Rs, given by where 3 is the maximum rate of Ca 2+ release by the IP 3 R and 3 is the IP 3 R open probability, defined by which follows a Hodgkin-Huxley gating formalism adopted by De Young and Keizer (De Young & Keizer, 1992) and later simplified by Li and Rinzel (Y. Li & Rinzel, 1994). In this simplification, the activation by IP 3 (defined by ∞ ) and [Ca 2+ ] i (defined by ∞ ) through binding to the receptor were assumed to be instantaneous, given by whereas the inactivation by [Ca 2+ ] i (defined by the gating variable ℎ), also through binding, was assumed to occur at a much slower time scale governed dynamically by where ℎ = In the study by De Young and Keizer (De Young & Keizer, 1992) the values of ( = 1, 2, 3, 5) were fit to experimental data (Watras & Ehrlich, 1991 where is the maximal rate of Ca 2+ leak from the ER. (4) Ca 2+ ATPase Activity ( and ). Ca 2+ removal by PMCA and Ca 2+ re-uptake into the ER by SERCA described by Hill functions (Chen et al., 2014;Croisier et al., 2013), given by where = , , is the maximal pumping rate and is the affinity of the pump to bind to Ca 2+ .
(5) P2X7R Ca 2+ Flux ( 2 7 ). Ca 2+ flux through P2X7Rs. A 12-state Markov model (Khadra et al., 2013) was initially used to compute the current, given by where 7 is the conductance of the P2X7R open states , = 1-4 ( Khadra et al., 2013). With emerging evidence suggesting that P2XRs do not dilate (M. Li, Toombes, Silberberg, & Swartz, 2015;Mackay et al., 2017), the maximum conductance of open ( 1 and 2 ) and sensitized/primed ( 3 and 4 ) states in this P2X7R sub-model were assumed to be equal. We also assumed that the rate of desensitization increased with ATP binding ( 2( 2) < 2( 1) < 2( 2) ) and that the open probability is higher in the desensitized row (i.e., 7 > , = 2, 4, 6). These modifications kept the time series simulations of the P2X7R model generally unchanged. To obtain the overall Ca 2+ flux through these channels, we then used the formalism from Zeng et al. (S. Zeng, Li, Zeng, & Chen, 2009) to convert ionic current to flux, scaled by a fraction that represents the average Ca 2+ flux (Egan & Khakh, 2004) relative to that for Na + and K + . The latter was necessary as P2X7Rs are non-selective cation channels. Using the above description, the following expression was used to describe this flux  (Mikolajewicz, Smith, Komarova, & Khadra, 2021). These simulations can be run by solving the function "fullmodel.m" using the ordinary differential equation solver ode15s. Figure    Transition rates between naïve (Fig. 3B, middle row) and sensitized (Fig. 3B,  0.8 Transition rates between naïve (Fig. 3B, middle row) and desensitized (Fig. 3B,