Long-term Mechanical Loading is Required for the Formation of 3D Bioprinted Functional Osteocyte Bone Organoids using Human Mesenchymal Stem Cells

Jianhua Zhang Swiss Federal Institute of Technology in Zurich Julia Griesbach Swiss Federal Institute of Technology in Zurich Marsel Ganeyev Swiss Federal Institute of Technology in Zurich Anna-Katharina Zehnder Swiss Federal Institute of Technology in Zurich Peng Zeng ScopeM Gian Nutal Schädli Swiss Federal Institute of Technology in Zurich de Leeuw Anke Swiss Federal Institute of Technology in Zurich Marina Rubert Swiss Federal Institute of Technology in Zurich Ralph Mueller (  ram@ethz.ch ) Swiss Federal Institute of Technology in Zurich https://orcid.org/0000-0002-5811-7725


Introduction
Organoids are self-organized three-dimensional (3D) in vitro tissue cultures that are derived from stem cells. Organoids offer unique possibilities for modeling and studying organogenesis and human pathologies. They open new avenues for drug discovery, diagnostics and regenerative medicine 1 . Up to now, bioengineered organoids are small and cannot be grown to tissue scale. Therefore, they lack the architectural features of native organs and cellular diversity that would allow the emergence of higherlevel tissue and organ characteristics 2,3 . Functional bone organoid formation presents challenges including development of a multicellular (osteoblast-osteocyte) tissue through self-organization and differentiation of stem cell and mineral formation with a 3D hybrid hierarchical architecture at the macroscopic scale. 3D bioprinting is a process based on additive manufacturing that shows promise for creating complex composite tissue constructs through the precise layer-by-layer placement of living cells and biomaterials 4,5,6 . Recently, Jonathan et al. 7 have shown how 3D bioprinting and organoid technology can be merged to control organoid formation from millimeter to centimeter scales. Organoidforming stem cells, such as human mesenchymal stem cells (hMSCs), were bioprinted directly into extracellular matrices (ECM) conducive to spontaneous self-organization. hMSCs are multipotent cells capable of differentiation along various tissue-speci c pathways, which have emerged as an attractive cell type for 3D bioprinting and organoid formation 8 . The readily available nature of hMSCs paired with an ability to undergo osteogenesis and adopt an osteoblast-or osteocyte-like phenotype makes these cells especially useful for bone tissue engineering. In our previous study, hMSCs-laden graphene oxide (GO)/alginate/gelatin composite bioink has been shown to keep good scaffold delity along the entire cell culture time, promoting osteogenic differentiation and ECM mineralization 9 . However, acquisition of an osteoblast-osteocyte co-culture characterized by lacuna-canaliculi network formation and mineral formation closely resembling those of native bone tissue has not been yet achieved.
In vitro mechanical loading, which is closely mimicking the native mechanical environment of the bone tissue, may be one strategy for optimizing osteogenesis and overcoming the above functional limitations.
The rationale for this strategy is given by the vital role of mechanical loading in bone development, bone remodeling and fracture healing 10 . In the development of healthy bone, each of these loading effects (compression or shear stress) has been shown to in uence various osteogenic responses of bone-derived cells (e.g., osteoblasts and osteocytes) and their progenitors (hMSCs). Osteocytes are considered to be mechanosensitive cells that sense shear stresses caused by the load-induced movement of the interstitial uid within the lacuna-canalicular system 11 . The stimulation of osteocytes leads to the osteogenic differentiation of hMSCs, the recruitment of osteoblasts to the bone surface, and the subsequent formation of mineralized bone matrix. Although the biomechanical modulation of in vivo bone tissue has been well studied, the precise role of mechanical stimulation in in vitro bone organoid formation is not clearly understood. To our knowledge, no literature has demonstrated that matured osteocyte phenotype cells were differentiated from hMSCs and the formation of 3D lacuna-canaliculi network in a 3D bioprinted bone organoid in vitro. Bioreactors are considered to be one of the most promising approaches for the application of mechanical loading on 3D cell-laden constructs for engineering bone tissue with sterile conditions, gas exchange and nutritional supply in vitro 12 . Externally applied mechanical stimuli, such as compression and shear stress, have been reported to promote osteogenic differentiation in vitro 13,14 . However, most of the studies were limited to short-term loading duration and focused on initial changes in gene expression 15 . These studies have not yet established a link between these transient gene expressions and functional outcomes related to the long-term ECM mineralization and cell selforganization of the forming osteocyte organoid.
In this work, we updated previous compression bioreactor designs 16 to fabricate culture systems compatible with 3D bioprinting, in situ microcomputed tomography (micro-CT) monitoring and cyclic mechanical loading. In addition to being X-ray translucent, the compression bioreactor could provide supplemental nutrition and gas exchange under sterile conditions and enables easy transfer. These properties are bene cial to cell culture and bone tissue engineering 16 . The 3D bioprinted hMSCs-laden scaffolds were cultured in the compression bioreactor system, which can apply individual cyclic mechanical loading through the compression of the 3D constructs. This study aims to form a 3Dbioprinted functional osteocyte bone organoid by applying long-term mechanical loading and changing the loading initiation time on 3D-bioprinted hMSCs-laden GO composite constructs over 8 weeks of culture. We hypothesize that the cells respond to the mechanical stimulus, directly sensing the scaffold deformation due to compression loading, which initiates a cascade of events that lead to 3D functional osteocyte bone organoid formation. Cell viability, DNA content, osteoblast/osteocyte-related gene expression, organoid stiffness, local Young's moduli, histology and immunohistochemistry staining were systematically investigated. The in uence of mechanical loading on bone-like tissue and lacunacanaliculi network formation in 3D bioprinted cell-laden scaffolds was visualized using in situ micro-CT and focused ion beam scanning electron microscopy (FIB-SEM).
2. Results And Discussion 2.1 3D bioprinting, compression bioreactor assembly and mechanical loading The components of the compression bioreactor and the assembled compression bioreactor are shown in Figure S1 AB. Compared to previous bioreactor systems 17 , our bioreactor could provide individual mechanical loading for construct and in situ visualization of mineral formation over time in micro-CT scanning. This study aimed to investigate the in uence of mechanical loading and loading initiation time on engineering 3D functional osteocyte bone organoids using 3D bioprinted hMSCs-laden GO composite constructs. For this purpose, we conducted cell culture studies with hMSCs. The 3D hMSCs-laden scaffold was bioprinted on the platform, and alginate chains were crosslinked with Ca 2+ after bioprinting ( Figure 1 AB), which was similar to our previous works 18, 19 . The scaffold morphology, including pore size, lament diameter and scaffold overall area, were similar to our previous study 9 . The platform with the 3D bioprinted hMSCs-laden scaffold was successfully assembled in the compression bioreactor, and the scaffold was cultured with 5 ml of osteogenic medium ( Figure 1C). Six bioreactors were kept in a rack constructed in-house and cultured in the incubator (5% CO 2 , 37°C) ( Figure 1D). The compression bioreactor was mounted into the developed mechanical stimulation unit (MSU) 16 to perform vertical cyclic mechanical loading on the 3D bioprinted hMSCs-laden scaffold. Mechanical compression is an important component of the load seen in bone. Our developed MSU provided highly accurate and precise displacement and load detection 16 . The load detection resolution of the system was 0.01 N. The forcetravel curve of the 3D bioprinted cell-laden scaffold is shown in the toe region and the linear region, exhibiting a typical stress-strain curve for the mechanical testing of soft tissue ( Figure 1F). The preload is an important parameter, which should be higher than the background force and avoid destroying the scaffold. The preload was set at 0.07 N in the linear region to reveal the scaffold surface and record the absolute length of the scaffold in our setup. Previous loading protocols varied with different bioreactor setups, materials and target engineered tissues 17 . Markus et al. 20 developed a pin-on-ball bioreactor system to apply mechanical loading to 3D cell-scaffold constructs. Cyclic axial compression was performed at a frequency of 0.1 Hz, with 5% compressive strain, and the ball oscillated at 0.6 Hz with an amplitude of ± 60° to mimic the native cartilage. In our pilot studies, the displacement plot lower than 10 Hz could illustrate a clear sine wave as the effective displacement. Therefore, the loading protocol was 0.07 N preload, 1% strain, 5 Hz, 5 min per loading time, 5 times per week, which were optimized to keep the scaffold intact for our study. The compression bioreactor was combined with nondestructive micro-CT monitoring to visualize the effect of mechanical loading on bone-like tissue formation.

Cell behavior and biochemical analysis
Numerous cellular biochemical responses to mechanical loading are transient, indicating the cell's ability to adapt its behavior to a new mechanical environment 21 . We sought to test the hypothesis that a cell can respond at the cellular (cell viability and morphological) and molecular (cell proliferation and osteogenic differentiation) levels to adapt to external mechanical loading in the short (day 9) and medium timepoints (day 31). The cell morphological responses to mechanical loading was also investigated at the end timepoints (day 56). Cell viability remained over 85% for all groups on day 9 and day 31. Predominant green uorescence and low red uorescence demonstrated the dominant population of live cells in Figure 2A-F. Figure 2G shows that cell viability was not signi cantly different between groups, and a small decrease in cell viability was observed in ML01 and ML21 on day 31. These results indicated that mechanical loading had a minimal effect on cell survival, as reported previously 15 . There were no differences in cell morphology (dendrite length, percentage of dendritic cells or percentage of spherical cells) at day 9 among the groups ( Figure 2H-J). Nevertheless, differences in the cell morphology among groups were observed from day 9 to day 31. Cell spreading behavior increased from day 9 to day 31 in all groups. This phenomenon was most pronounced in the mechanical loading groups, as the dendrite lengths of ML01 and ML21 were signi cantly higher on day 31 than on day 9. The percentage of dendritic cells was signi cantly increased and the percentage of spherical cells was signi cantly decreased in all groups from day 9 to day 31. More speci cally, the highest dendrite length and percentage of dendritic cells were 49.6 ± 3.2 µm and 86.7 ± 5.5% in ML01 at day 31, which was signi cantly higher than 23.53 ± 4.2 µm and 42.8 ± 4.0% in NL. The percentage of spherical cells in NL was signi cantly higher than that in ML01 on day 31. The dendrite length and percentage of the dendritic cells in ML21 were not signi cantly different than those in ML01 on day 31 but were higher than those in NL. The same results were also shown in the actin staining from the sections in the top region of NL, ML01 and ML21 at day 56 culturing in the bioreactors. ML01 had more spreading osteocyte-like cells with longer dendrite length compared to NL and ML21 groups ( Figure S2 A-B). Two single cells were connected to form the cellular network in loading groups, while single cells did not connect in NL groups at day 56 ( Figure S2 D-F). Similar 3D cell spreading behavior enhanced by dynamic cyclic mechanical stimulation was observed for NIH-3T3 broblasts and hMSCs in previous reports 15 , suggesting that mechanical loading had a positive in uence on cell morphology.
The in uence of mechanical loading on cell proliferation in 3D bioprinted cell-laden scaffolds was investigated by analyzing DNA content on day 9 and day 31 ( Figure 2K). On day 9, the DNA content in ML01 was signi cantly higher than that in NL and ML21, which indicated that mechanical loading increased cell proliferation in the early culture period. Yoshikawa et al. 22 have also shown that the DNA content signi cantly increases in the stretch-stimulated group on day 8 when rat bone marrow cells were cultured on Petriperm TM dishes. An approximately 70% decrease in the average DNA content from day 9 to day 31 within the same groups were found after culture in osteogenic media for 31 days. The overall DNA content decrease was similar to that in our previous results 19 , which may have resulted from two factors: cell washout and the osteogenic differentiation of hMSCs. Consistent with previous reports, we observed that some encapsulated cells had emigrated from the scaffolds and migrated on the platform due to the movement of medium within the bioreactor 15 . Previous studies also found increased cell death and decreased cell number with increased osteogenic differentiation of hMSCs, which could have occurred in this system 23 . Soe et al. 15 reported a 36-42% decrease in the average DNA content from day 7 to day 21 during the osteogenic differentiation of hMSCs. Similar to previous mechanical loading studies in 3D cell-laden constructs 24 , our results showed that there were no signi cant differences in DNA content among the different groups at day 31. These results indicate that mechanical loading and loading initiation time did not impact cell proliferation after 31 days of culture.
Alkaline phosphatase (ALP) is known as an early biochemical marker for the osteoblast phenotype. ALP is considered an essential factor in assessing bone differentiation and mineralization. Here, ALP activity was measured on day 9 and day 31 after culture in the compression bioreactor. Figure 2L shows that ALP activity was not signi cantly different among groups on day 9. However, ML01 had the highest ALP activity (11.8 ± 1.3 nmol/ng), which was signi cantly higher than that of NL (5.3 ± 0.7 nmol/ng) and ML21 (5.6 ± 0.5 nmol/ng) on day 31. Based on the results obtained on cell morphology, mechanical loading enhanced the cell spreading behavior and the longest dendrite length was observed in ML01 ( Figure 2H). Here, we suggest that higher cell spreading observed in ML01 might enhance differentiation of hMSCs towards the osteoblastic lineage. The increases in ALP activity in ML01 on day 31 may indicate both more osteogenic protein production of mature osteoblastic cells and the higher induction of differentiation in osteoprogenitors. Meanwhile, ALP activity was signi cantly increased from day 9 to day 31 in the same group. An overall increase in ALP activity can be due to an increase in the production level of the osteogenic protein and/or a rise in the numbers of alkaline phosphatase-positive cells during the osteogenic differentiation of hMSCs. Interestingly, there was no difference between ML21 and NL on day 31, which indicates that the in uence of mechanical loading on osteogenic differentiation may be related to loading time. In vitro mechanical stimulation has been shown previously to upregulate ALP protein expression and activity within hMSCs culture 25 . These results support the contention that mechanical stimulation promotes differentiation to osteoblasts and enhances bone formation.

In situ micro-CT monitoring of mineral formation and maturation
Bone adapts to mechanical loading, leading to an overall increase in bone mass and bone strength through a positive shift in the balance between bone formation and bone resorption in vivo. In vitro mineral formation and maturation in response to external mechanical loading were investigated through in situ micro-CT monitoring in our study. Figure 3A shows the weekly mineral formation in 3D reconstructed micro-CT images of the NL, ML01 and ML21 groups in horizontal view. Compared to NL, mechanical loading enhanced the even distribution of mineralized ECM deposition throughout the scaffold in ML01 and ML21. More homogeneous mineral distribution was observed in ML01 and ML21 in the top view from day 28 to day 56. From the vertical view in Figure S3, mineralized tissue showed more uniform growth along the scaffold surface than inside the scaffold. Previous studies have shown that ow perfusion induced tissue modeling with the formation of pore-like structures in the scaffolds and enhanced the distribution of cells and matrix throughout the scaffolds in vitro 26 . Our results indicated that compression mechanical loading also improved the distribution of mineral formation, resulting in homogeneously mineralized scaffolds ( Figure 3A).  Figure S4 C) groups is presented. The low-density volume decreased and the high-density volume increased for all groups over time, which indicated that the mineral matured and became stiffer. The quanti cation of total mineral volume over time in the NL, ML01 and ML21 groups is shown in Figure 3B. On day 35, the total mineral volume was signi cantly higher in ML01 than in NL. Afterward, there were no signi cant differences between the two groups on days 42, 49 and 56. The mineral volume in ML21 and NL was similar on day 21. No signi cant differences were found on days 21, 28, 35, 42, and 49, while there was a signi cant difference between ML21 and NL on day 56. Interestingly, ML21 had a higher total mineral volume than ML01 on days 42, 49, and 56, but no signi cant differences were found. The peak mineral volumes for all groups are shown in Figure S4 E, when the threshold value was over 135.27 mg HA/cm 3 . In the rst four weeks, there was no mineral formation in any group. Afterward, the mineral volumes increased, especially in the ML01 group. ML01 had the highest peak mineral volume on days 42, 49, and 56, which were signi cantly higher than those of NL and ML21. The mineral formation rates and the total mineral volume in NL, ML01 and ML21 between the two timepoints are shown in Figure 3C. The mineral formation rate was signi cantly higher in ML01 than in ML21 and NL on days 21-28, after which it decreased. ML21 had a signi cantly higher mineral formation rate on days 35-42 and 42-49 than ML01. No signi cant differences in the mineral formation rate were found between ML21 and NL. Organoid mineral density (OMD) based on the total mineral volume was measured at different timepoints. Figure 3D shows that the OMD increased, which indicated that the maturation of mineral increased over time. No signi cant differences were observed after 42 days of culture in the compression bioreactor in any groups. ML01 had the highest OMD values on day 49 and day 56, which were signi cantly higher than those of NL and ML21. The OMD in ML21 was similar to that in NL, and no signi cant differences were found between NL and ML21. Our results indicated that mechanical loading enhanced mineral maturation with higher OMD and peak mineral volume but not total mineral volume in ML01 on day 56. Mauney et al. 27 showed that mineralized matrix production increased under mechanical stimulation with 10 nM dexamethasone for 8 and 16 days when hMSCs were cultured on 3D partially demineralized bone scaffolds. In vivo mechanical stimulation has been shown previously to play an integral role in bone formation by stimulating osteoprogenitor cells in the marrow stroma to differentiate into osteoblasts at cortical bone surfaces. Our results are the rst to show the long-term in uence of mechanical loading on the mineral formation and maturation in vitro for 3D bioprinted cell-laden scaffolds. These results lead us to the hypothesis, that organoid mineral density is more important than mineral volume regarding the mineral properties in micro-CT evaluation.

Scaffold stiffness
Furthermore, the mechanical adaptation of the 3D bioprinted hMSCs-laden scaffolds in response to external mechanical loading was investigated. The median curves of force-displacement changing in different groups after culture in osteogenic medium for 56 days are shown in Figure 4A. The forces were dramatically increased in all groups in the rst 100 µm of displacement, and the slope and force in ML01 were much higher than those in ML21 and NL. As shown in Figure 4B, the stiffness at day 1 was 75.7 ± 10 N/m, and stiffness was signi cantly increased in all groups after 56 days of culture in the compression bioreactors. The stiffness in ML01 was 2190 ± 727 N/m, which was signi cantly higher than that in NL (415 ± 156 N/m) and ML21 (569 ± 200 N/m) on day 56. Additionally, the average Young's modulus of the mineral from the micro-CT images were calculated by Finite Element (FE) stimulation using the conversion function: The average Young's modulus of the ML01 group was 445 ± 181 kPa and signi cantly higher than that in both the NL (154 ± 70 kPa) and the ML21 (111 ± 76 kPa) on day 56 ( Figure 4C). Several studies have shown that dynamic mechanical loading enhances the maturation of MSCs-laden constructs, resulting in higher mechanical and functional properties 28 . There was no signi cant difference between NL and ML21 in both stiffness and average Young's modulus. These results suggest that the mechanotransduction pathways initiated by dynamic compression may differ fundamentally between osteoblasts and undifferentiated or osteogenically differentiated MSCs. To our knowledge, this is the rst study to demonstrate improved mechanical properties of 3D bioprinted MSCs-based engineered bone organoid through the long-term application of dynamic compressive loading. Although we achieved a 5fold improvement in stiffness and a 3-fold improvement in Young's modulus over 8 weeks in the ML01 group compared to the NL group, these values remained lower than those of native bone. Future studies will optimize these parameters over longer culture durations to further explain load-induced increases in mechanical properties.
Structural integrity of the bone is affected by the composition and organization of the constituents within this structure. Hydroxyapatite, which is the mineral component, accounts for more than 50% of the bone volume and the rest is composed of collagen and water 29 . To evaluate the bone health and particularly for diagnosing osteoporosis, assessment of the amount of the mineral (bone mineral content) or the density of the mineral (bone mineral density) with imaging is considered as the gold standard. Many studies have investigated the correlation between the elastic modulus and bone mineral density, to form the reference relations which can estimate the mechanical properties from the image-based measurements of the bone mineral density. In our study, Figure 4D shows that there was no correlation between organoid stiffness and total mineral volume. Instead, organoid stiffness had a strong power law correlation (R 2 =0.98) with organoid mineral density (OMD) ( Figure 4E). The power value is 10.97. These results indicated that organoid stiffness in the 3D bioprinted osteocyte bone organoid was positively related to OMD. This relationship is similar to trabecular and compact bone, where modulus is related to density through a power law relationship at the bulk length-scale 30, 31 . Carter et al. 30 demonstrated that the compressive stiffness of bone is proportional to the square of the apparent density and to the strain rate raised to the 0.06 power. Two factors may contribute to the different power values relationship in stiffness and density between our 3D osteocyte bone organoid and bone specimen. The rst factor is given to the different units of OMD and apparent density used between studies; while in our study OMD was presented in mg HA/cm 3 , the apparent density is given in mg/cm 3 in most of the literature. The second factor is due to the difference in strain rate. Mathieu et al. 31 demonstrated that the modulus and collapse stress increased linearly with strain rate in 3D poly( L -lactic acid)/hydroxyapatite scaffolds.
Those results demonstrated that our mineralized bone organoid exhibited similar power correlation behavior to that of natural bone at the bulk length-scale. To verify the hypothesis that OMD and mechanical properties also correlate to tissue level. FE Young's modulus calculated by Function (1) was converted to micro-FE stiffness and the results showed a strong correlation with experimental stiffness results (R 2 = 0.94) ( Figure 4F). It indicates that small increases in OMD can lead to substantial increases in FE Young's modulus and micro-FE stiffness. Experimental mechanical property could be estimated by FE stimulation using the Function (1). Scha er et al. 32 showed a power law correlation with similar exponent in cortical bone, although with a lower preceding factor. In overall, these results con rmed the power law correlation of OMD with Young's modulus, and the bigger in uence of OMD on scaffold mechanics compared to scaffold morphometry in mineral volume. With the ability to convert local OMD to Young's modulus using FE stimulation, these ndings will provide in the future a powerful new tool for non-destructive assessment of local mechanical feedback. This will help estimate the mechanical properties from the micro-CT image measurements of the mineral density in our organoids and bone tissues.

Osteogenesis-related gene expression
The expression of several osteogenic gene markers was characterized to provide a better understanding of the events involved during the in vitro osteogenic differentiation of hMSCs under mechanical stimulation. We evaluated the osteogenic differentiation stage of hMSCs in terms of early (COL1A2, ALPL, RUNX2) and late (BGLAP) stage osteoblast-related genes and osteocyte-related genes (PDPN, PHEX, DMP1, SOST) by RT-qPCR analysis after 56 days of culture in compression bioreactors. COL1A2, ALPL and RUNX2 are important genes in the early osteogenic differentiation pathway of progenitor cells and can be upregulated by mechanical stimulation in 2D and 3D environments 15,33 . The NL group had the highest COL1A2 and ALPL gene expression, which were signi cantly higher than ML01 and ML21. No differences were found in COL1A2 and ALPL gene expression between ML01 and ML21 (Figure 5A B). A similar trend could be found in RUNX2 gene expression: expression in the NL group was signi cantly higher than in ML21, but no difference was found between NL and ML01 ( Figure 5C). The osteoblastrelated gene expression results indicated that the largest number of early-stage osteoblasts was present in the NL group. Zhang et al. 34 demonstrated that mechanical strain regulated RUNX2 activation and favored osteoblast differentiation through activation of the phosphorylation level of extracellular regulated protein kinase (ERK)1/2 signaling pathway. However, most of those mechanical loading experiments were performed in short terms, such as on day 7 and day 14.
The opposite trend was found in late-stage osteoblast-related gene expression. ML01 had the highest BGLAP gene expression ( Figure 5D), which was signi cantly higher than that in NL and ML21. No signi cant difference was found between NL and ML21. The gene expression of BGLAP indicates that more late-stage osteoblasts were present in the mechanical loading groups, especially in ML01. Our results are in line with other related studies showing that osteocalcin expression is sensitive to mechanical loading within hMSCs in vitro 27 . During the bone development process, matured osteoblasts will be embedded within the bone matrix and differentiate into osteocytes. Osteocytes are dendritic cells that express multiple genes involved in mineral metabolism, including PDPN, PHEX, DMP1 and SOST. For osteocyte-related gene expression, PDPN and PHEX have been implicated in the initial stages of osteocyte differentiation. Figure 5E F shows that PDPN and PHEX gene expression was similar in each group and exhibited no signi cant difference. However, mechanical loading upregulated osteocyte-related (DMP1 and SOST) gene expression. Only two in four samples in the NL group exhibited DMP1 and SOST gene expression. All samples in the ML01 and ML21 groups upregulated DMP1 and SOST gene expression. DMP1 gene expression was 20-fold and 7.2-fold higher in ML01 and ML21 than in NL, respectively. The DMP1 gene expression in ML01 was signi cantly higher than that in ML21. A similar trend was observed for SOST gene expression. SOST gene expression was found to be 7.8-fold higher in the ML01 group and 4.7-fold higher in the ML21 group than in the NL group. DMP1 protein is localized in the lacunar and canalicular walls of osteocytes, and the gene responds rapidly to mechanical stimulation 35 . In vivo studies have shown that DMP1 expression increases in response to a mechanical stimulus 36, 37 . However, mechanical stimulation has been shown to induce downregulation of sclerostin gene and protein expression and promote the Wnt signaling pathway and osteoblast formation 38 . Our results showed that early osteoblast-related gene expression in the osteogenic differentiation of hMSCs was decreased and the expression of late-stage osteoblast-and osteocyte-related genes were increased after long-term mechanical loading. However, the detailed mechanism of mechanical loading effects on the osteogenic differentiation of hMSCs is still unclear.

Histology staining and collagen I analysis
The adaptive response of bone to mechanical loading requires increased synthesis and turnover of matrix proteins and ECM mineralization, with special emphasis on collagen I. In this study, we have evaluated the synthesis and ECM mineralization under different mechanical loading conditions. Alizarin Red S staining showed the mineral was stained red, and no differences were observed between ML01 and ML21. More yellowish minerals were found in NL, which indicated less mineralization in one section ( Figure S5 A). H&E staining showed a uniform cell distribution in sections from the NL, ML01 and ML21 groups (Figure S5 B). Figure 6D shows representative 2D original micro-CT slides without thresholding from the top region of the NL, ML01 and ML21 scaffolds. The mineral was visualized in white and the unmineralized tissue or culture medium in black, showing that ML01 had more mineralization and more homogeneous mineralization than NL in one slide. The 2D micro-CT image was consistent with the Alizarin Red S staining ( Figure S5A) and 3D micro-CT results ( Figure 3A). As shown in Figure S5 C, Sirius Red-stained sections appeared pink-red, and all areas were stained; no differences were observed in the NL, ML01 and ML21 groups under bright eld microscopy. Observing the same section with polarized light allowed not only clear identi cation of collagen networks but also semiquantitative appreciation of the collagen I network during maturation. Because of the birefringent properties of collagen I 39 , we observed different alignments of the collagen bers using different polarized light orientations, as shown with arrows in Figure 6A. Figure 6A show that thicker and longer collagen bers were observed in the ML01 and ML21 groups than in the NL group. Collagen bers at different maturation levels appear in various colors, including red, orange, yellow and green, which indicate differences in collagen bers during maturation, as analyzed by ImageJ. More mature (red) collagen bers were found in the ML01 and ML21 groups, while more immature collagen bers (green) were found in the NL group ( Figure 6B). The quanti cation results ( Figure 6C) showed that the red mature collagen was signi cantly higher in the ML01 (63.3 ± 3.4%) and ML21 (67.4 ± 4.2%) groups than in the NL (29.0 ± 12.6%) group. No differences were observed in the percentage of mature collagen between ML01 and ML21. After reconstruction using CT-re, the collagen I network is shown in Figure 6E. Collagen I bers are shown in different colors to distinguish stages of maturation. Three examples of the quanti cation of collagen I ber width, length, straightness and angle by CT-re reconstruction in NL, ML01 and ML21 groups are shown in Figure S6. The mechanical loading groups had signi cantly higher ber width and ber length than the non-loading group. ML21 had the longest ber width and length, while ML01 had the highest ber number, which was signi cantly higher than that in NL and ML21. Our results indicated that the level of collagen I maturation and ber parameters ( ber width and length) were increased in the mechanically loaded groups. Similar results showed that collagen content as assayed by Sirius Red was signi cantly (2-fold) higher on days 15 and 20 in loaded samples than in static controls in which osteoblasts were seeded on polyurethane scaffolds 40 . ECM turnover is in uenced by physical activity in tendon and skeletal muscle, and both collagen I synthesis and maturation increase with mechanical loading 41 . These changes may modify the mechanical properties and the viscoelastic characteristics of the scaffold and likely increase its stiffness in our study. Some researchers have demonstrated that physiological signaling pathways (mitogenactivated protein kinase) and growth factors (transforming growth factor-β, insulin-like growth factor) play an important regulatory role in collagen I synthesis and turnover concerning mechanical loading 41 .

Osteoblast and osteocyte immunohistochemistry staining
To identify the presence of osteoblasts and osteocytes in our 3D bioprinted hMSC-based osteogenic organoid after culture in the compression bioreactor for 56 days, osteocalcin (OC), podoplanin and sclerostin immunochemistry stainings were performed. The images ( Figure S7 A-C) for NL, ML01 and ML21 show the nuclei in cyan and OC in red. More OC staining was present in the mechanical loading groups than in the non-loading group. The osteocalcin integrated intensity by area in the NL, ML01 and ML21 groups after culture in osteogenic medium for 56 days was 587 ± 224, 3140 ± 1174.77 and 1913 ± 518, respectively (Figure S7 D). ML21 and ML01 had signi cantly higher OC integrated intensity than NL. These results are consistent with the BGLAP gene expression, which indicates that mechanical loading enhances the development of mature osteoblasts with increasing osteocalcin expression at both gene and protein levels.
In vitro, functional osteocytes formation through the full differentiation from MSCs has not yet been demonstrated for human cells. Here, we showed that mechanical loading induced 3D osteogenic organoid formation from hMSCs to functional osteocyte differentiation. First, the immunohistochemistry control staining data are shown in Figure S8, including cell nuclei, actin, primary antibody and secondary antibody control staining. Nuclei and actin staining showed the presence of cells, but total darkness was shown in the primary and secondary antibody controls under the different wavelengths. Figure S8 indicates that immunohistochemistry staining was performed only by combining the primary and secondary antibodies, and the staining signal was from the target protein staining. Figure 7 (low magni cation) and Figure S9 (high magni cation) show the cell nuclei, actin, podoplanin, sclerostin and merged staining images in the NL, ML01 and ML21 groups on day 56 of culture in the compression bioreactor. Podoplanin is considered a critical transmembrane glycoprotein in the early stage of osteocytogenesis. All groups exhibited podoplanin protein expression, and the mechanical loading groups showed higher podoplanin expression than the non-loading group in Figure 7. Sclerostin is a product of the SOST gene, which is expressed by mature osteocytes embedded in mineralized bone. In line with the gene expression results, both mechanical loading groups exhibited sclerostin protein expression. ML01 showed much higher sclerostin expression than ML21, and sclerostin was surrounded by osteocytes. The NL group showed no sclerostin protein expression at day 56. Our results indicated that mechanical loading enhanced osteocyte differentiation and maturation from osteoblasts in vitro. In vivo studies showed that osteocytes regulate bone mass in response to mechanical loading by decreasing levels of sclerostin, which acts by inhibiting the Wnt/β-catenin pathway in osteoblasts and increasing osteogenesis 42 . The sclerostin protein expression in ML01 may explain why the total mineral volume was not signi cantly different between ML01 and NL after day 42, as the sclerostin expression in ML01 may activate the Wnt/β-catenin pathway in osteoblasts and inhibit osteogenesis while improving mineralized scaffold density.

Osteocyte morphology and lacuna-canaliculi formation
Osteocytes are the most abundant cells in bone and are embedded in sites in the bone mineral matrix called lacunae. They form an extended 3D network whose processes interconnecting the cell bodies reside in thin canals, the canaliculi. Together with the osteocyte lacunae, the canaliculi form the lacunarcanalicular network (LCN). Denser mineral and higher osteocyte gene and protein expression were con rmed in ML01 scaffolds compared to the other groups. Here, we investigated osteocyte morphology and LCN formation in NL, ML01 and ML21 after culture in a compression bioreactor. Video S1 showed that a 3D cellular network was formed in the ML01 scaffold, and low-magni cation of the video showed the 3D cellular network morphology on day 31. High-magni cation images showed how two single cells connect on day 56 in Video S2. The results indicated that mechanical loading induced 3D cellular connection. Single-osteocyte morphology is shown in Figure S10, which shows numerous slender processes. These processes were often branched, showing similarity to osteocyte morphology in human bone tissue. Osteocytes occupied small pores in the mineralized tissue in ML01, as shown by the combination of cell nuclei/actin staining in confocal microscopy and backscattered electron (BSE) images in SEM in Figure 8, S11, S12 and S13. While most of cells were on the surface of the mineralized section in NL and ML21 groups, which did not locate at the pore structure ( Figure S11). The 2D osteocyte organization in the ML01 section is shown in Figure S12 Figure 8 and S13 show that a single osteocyte was located in a lacuna in the high magni cation view. Those results indicated lacunae were formed in ML01 after 56 days cultured in the compression bioreactor, while no lacunae formation was found in NL and ML21 groups.
Calcium, phosphate and oxygen were enriched around lacunae, which con rmed that osteocytes were embedded in the mineralized tissue in Figure 8E-G, Figure S12 G-I and Figure S13 G-I. The Ca/P ratio was 1.1 ( Figure 8H), which was close to the 1.67 Ca/P ratio of hydroxyapatite in bone tissue. These results con rmed that mechanical loading induced lacuna formation in our 3D bioprinted centimeter-scale osteocyte organoid. In vivo, osteocytes are connected with neighboring cells through dendritic cell processes. The processes within the mineral matrix reside in thin canals, which are called canaliculi. The presence of canaliculi in our 3D bioprinted osteocyte bone organoid was con rmed by FIB-SEM in Figure   8D1-D4 and Figure S14. Figure 8D1-D4 are the SEM images of the FIB cross-sections in the same positions, which are shown as dashed lines in Figure 8D. Osteocytes were present in the lacunae and surrounded by calcium phosphate nanoparticles ( Figure 8D1). FIB cross-sections were cut along the process direction, and canaliculi are shown in Figure 8D2-D4. The canaliculi were surrounded by calcium phosphate nanoparticles. The results indicated that mechanical loading induced canalicular formation along the process direction. The SEM images in Figure S14 show mineral deposition in ML01 scaffolds, including nanoparticles and the formation of small and large crystals of calcium phosphate. Overall, osteocytes, LCNs and minerals were formed in our 3D bioprinted osteocyte organoid (ML01) after culture in the compression bioreactor for 56 days. Mechanical loading plays a vital role in bone osteocyte organoid formation from 3D bioprinted hMSC-laden scaffolds.

YAP staining
Yes-associated protein (YAP)/transcriptional coactivator with PDZ-binding motif (TAZ) activity was found to be profoundly regulated by both ECM elasticity and cell geometry and was instrumental for cells to carry out the corresponding biological responses to mechanical inputs 43,44 . Cells can be exposed to different types of mechanical cues, which may regulate YAP/TAZ expression. To determine the loading effect on YAP protein expression in cells, YAP immunochemistry staining was performed on the sections in the top regions of the scaffolds with different mechanical loading conditions. Figure 9 and Figure S15 show the YAP staining in different groups at high and low magni cation. Mechanical loading upregulated YAP protein expression, which was more pronouncedly expressed in ML01 as compared to ML21.
Although also expressed, YAP was lower expressed was in NL and the majority of YAP protein was expressed intercellularly. Recent studies showed that shear stress facilitates MSC differentiation into osteogenic cells by sustaining TAZ 45 . Cyclic stretching is su cient to promote the formation of a tensile cytoskeleton and thus nuclear YAP in broblasts when seeded on a soft substratum 46 , or in cells undergoing contact inhibition 47 . Our compression bioreactor system may provide a new approach to investigate the in uence of compression mechanical loading on YAP/TAZ activity in the 3D ECM environment. Mechanical loading could regulate and increase YAP expression, although the detailed interaction mechanism by which YAP expression affects 3D osteogenic organoid formation is still unclear.
In summary, we have presented a compression bioreactor that could be combined with a 3D bioprinting and micro-CT monitoring system to provide mechanical loading for 3D functional osteocyte bone organoid formation in vitro. Our results showed that mechanical loading and loading initiation time play critical roles in the formation of self-organizing 3D bioprinted osteocyte bone organoid. The 3D osteogenic organoid exhibited mature collagen I bers, osteoblast/osteocyte phenotypes, and LCN embedded in dense mineral with good mechanical properties, which opens new avenues for therapeutic applications. Our compression bioreactor system could provide mechanical loading not only for osteogenic organoid formation but also for the in vitro formation of different organoids, such as kidney, cartilage, and heart.

Bioinks preparation and 3D bioprinting of cell-laden GO composite scaffolds
Details of the materials are provided in the supporting information. GO/alginate/gelatin (0.1%/0.8%/4.1% w/v) ink solution was prepared by using sodium alginate, gelatin and GO (5 mg/ml), according to our previous method and presented in the supporting information. hMSCs (Lonza, Walkersville, MD, USA) were isolated from human bone marrow aspirate and characterized as described previously 19 . P3 hMSCs were expanded in expansion medium (DMEM, 10% FBS, 1% P/S/F, 1% NEAA and 1 ng/ml bFGF) under standard cell culture conditions (37°C, 5% CO 2 ) for 7 days until approximately 80% con uence. To prepare the bioink solution, cells were resuspended in control medium (DMEM, 10% FBS, 1% P/S/F) at a concentration of 5 million cells per 100 µl. The cell suspension was mixed with 1 ml of GO/alginate/gelatin ink solution, and the bioink was thoroughly mixed with sterile magnetic sh at 200 rpm/min for 1 min at room temperature.
A microextrusion-based two-syringe cell bioprinter (INKREDIBLE + , CELLINK, USA) was used to print the 3D cell-laden GO composite scaffolds. The scaffold model (10 mm × 10 mm × 2.5 mm), printing parameters and printing processes were as described in our previous work 19 . Scaffolds were bioprinted on doublesided tape (3M TM , USA) on the platform of the compression bioreactor, which was autoclaved before the experiment. After bioprinting, the scaffold on the platform was crosslinked with 2% (w/v) CaCl 2 (VWR, Dietikon, Switzerland) in DMEM for 10 min, and then the scaffold was washed three times with DMEM.

Compression bioreactor assembly and mechanical loading
In this work, we updated a new platform and piston design to culture our 3D bioprinted cell-laden scaffolds and apply mechanical loading. The compression bioreactor could be assembled with all the components, including platform, piston, ring, inner wall, outer wall, base, and spring, shown in Figure S1 A and assembled as shown in Figure S1 B. Before assembly, all the components were autoclaved at 121°C for 45 min. On day 1 after bioprinting, the scaffold on the platform was assembled with the remaining components and cultured with 5 ml of osteogenic medium.
A mechanical stimulation unit (MSU) was custom-made to induce load onto the scaffold. The compression bioreactor could be mounted into MSU-like cartridges (see Figure 1E). The apparatus was controlled via an in-house program on LabView (National Instruments, Austin, TX). The MSU was placed inside a sterile laminar ow hood next to the incubator. The culture compression bioreactors were individually removed from the incubator and mounted to the MSU for mechanical compression. To de ne the preload between the scaffold and piston for scaffold detection and avoid scaffold destruction during loading, force-displacement curves were measured with three different scaffolds in MSU. The preload was 0.07 N, which was in the linear region of the force-displacement curve. The loading region was in the 1% strain in the force-displacement curve, according to the force research on the preload between the piston and scaffold. The loading frequency was 5 Hz with 1% sinusoidal strain (see Figure S1 C). The loading time was 5 min per day and 5 times per week. The strain, frequency and loading time were input into LabView to control the mechanical loading. When mechanical loading was applied ( Figure 1G), the scaffold was cultured on the platform between the base and piston in the culture position, and then the base was lifted slowly. The scaffold touched the piston surface, and LabView recorded the scaffold absolute length when forces reached the preload in the contact position. The base moved up and down in 5 Hz and applied mechanical loading to the scaffold in the loading position based on the loading protocol. We carried out a total of three independent experiments ( Figure 1H). Group 1 was the control group, and non-loading was applied to the 3D cell-laden scaffolds cultured in compression bioreactors. Group 2 received mechanical loading on the scaffolds from day 1. Group 3 was a static culture in the compression bioreactors for the rst 21 days and was subjected to mechanical loading on the scaffolds from day 21. The three different groups were labeled NL, ML01, and ML21. During the experiment, the osteogenic medium was changed three times per week.

Cell viability and cell morphology
Cell viability in the 3D bioprinted cell-laden GO composite scaffolds with different loading conditions was assessed using a LIVE/DEAD® Viability/Cytotoxicity assay after 9 days and 31 days of culture in osteogenic medium. Confocal microscopy (Visitron Spinning Disk, Nikon Eclipse T1) was used to visualize the living (green) and dead (red) cells. Cell viability was calculated as the percentage of living cells among total cells. The details of cell viability and cell morphology analysis are presented in the supporting information.

DNA content and ALP activity
Cell proliferation and osteogenic activity were determined by measuring DNA content and alkaline phosphatase (ALP) activity. After culturing in osteogenic medium for 9 and 31 days, scaffolds in the NL, ML01 and ML21 groups were carefully removed from the platforms and washed twice in PBS. For DNA content and ALP activity measurements, the scaffolds were disintegrated in 0.7 ml of 0.2% (v/v) Triton X-100 and 5 mM MgCl 2 solution using two steel beads and a minibead beater (Biospec Products, Bartlesville, OK). ALP activity measurements were carried out directly after the scaffolds had been disintegrated. The measurements were performed using a commercially available assay based on the conversion of p-nitrophenyl phosphate to p-nitrophenol, according to the protocol of the manufacturer. DNA contents were measured after 48 h of incubation at room temperature when scaffolds had been disintegrated. The Quant-iT TM PicoGreen assay (Life Technologies, Switzerland) was carried out from the supernatant fraction and followed the manufacturer's instructions. The DNA content was calculated with three samples in every group and averaged. The calculated ALP activity was normalized by the DNA content measured for each scaffold.

In situ micro-CT monitoring of mineral formation and maturation
In situ micro-CT images of all groups were taken on days 7, 14, 21, 28, 35, 42, 49 and 56 of culture in a compression bioreactor with 5 ml of osteogenic medium to monitor mineral formation as described previously 19 . The bioreactors were removed from the rack and transferred to a micro-CT 40 (SCANCO Medical AG, Brüttisellen, Switzerland) for imaging in a xed position. After each scan, the bioreactors were returned to the rack, and placed in the incubator (37°C, 5% CO 2 ). The methods of micro-CT image reconstruction, mineral histogram distribution, total and mineral volume, and mineral scaffold density (MSD) are provided in the supporting information.

Scaffold mechanics
The mechanical properties were characterized by measurement of the stiffness. The stiffness was assessed on a custom-made mechanical stimulation unit (MSU) at room temperature. An uncon ned uniaxial compression test was performed under displacement control to measure the force-displacement curve, with a preload of 0.07 N and a speed of 4 µm s −1 (corresponding strain rate of 0.002 s −1 ) until 500 µm deformation of the scaffold was reached. 3D cell-laden scaffolds (n = 4) in different groups were tested after culturing in osteogenic medium and applying mechanical loading for 56 days. The stiffness was calculated from the linear region of the force-displacement curve according to the following formula: where ΔF is the force change in the linear region, and ΔL is the displacement change in the linear region. The stiffness was calculated in several small linear regions in the force-displacement curve, and the highest value was selected as the nal value for the sample. The correlation between stiffness and mineral volume and MSD was analyzed by origin 2018 (OriginLab, USA) with a power law 30 in the following allometric formula: Y = a + bX c .
The correlation between local MSD and Young's modulus was investigated using the linear elastic micro-Finite Element (micro-FE) solver ParOSol 48 . To overcome inadequate boundary conditions due to uneven surfaces of the scaffolds as well as to reduce computing time, suitable subsets of the micro-CT images were used for the investigation. Correlation of subset density with scaffold density and scaffold stiffness showed subsets in form of slices with 1/6 height from the bottom third are representative and were thus chosen for further investigation. Subsets' local density was converted to Young's modulus using multiple conversion functions 49 for the subsequent micro-FE analysis simulating a 1% uniaxial compression. The subset stiffness was calculated from the resulting force and subsequently correlated with experimental stiffness. Consequently, the conversion function with the strongest linear correlation was adjusted to t the experimental results, and extrapolated to full height by factor 6. Micro-CT images were then converted to Young's modulus and the mean value of the mineralized tissue was calculated.

Osteogenesis-related gene expression
After culturing in the osteogenic medium for 56 days, compression bioreactors were carefully disassembled, and scaffolds were removed from the platform. The total RNA amount in the 3D cell-laden  Table S1. The C t values obtained for each sample were normalized to those for the housekeeping gene and presented as the fold-changes in gene expression of the loading groups versus the non-loading group. Data were analyzed using the comparative Ct method (2 −ΔΔCt ).

Crysection, histology staining and collagen I analysis
After 56 days of culture in the compression bioreactors with osteogenic medium, scaffolds were washed twice with PBS and xed in 4% formaldehyde in 10 mM CaCl 2 and 0.15 M NaCl solution for 2 h. Then, samples were placed in 10% sucrose with 10 mM CaCl 2 for 2 h at room temperature. Next, the samples were placed in 30% sucrose with 10 mM CaCl 2 for 16 h at room temperature. Finally, the samples were embedded in optimum cutting temperature compound (OCT, VWR) and snap-frozen on dry ice. Samples were sectioned using Kawamoto's cryo lm type 2C (SECTION-LAB Co. Ltd., Japan) to 10 µm thickness using a cryotome (CryoStar NX70, Thermo Scienti c) and according to Kawamoto's protocol 50 . Before staining, sections were xed on microscope slides using 1% chitosan adhesive. For this purpose, two drops of chitosan solution were deposited on the slide, and the sections were placed on the slide and kept in a running fume hood to allow the chitosan to dry 50 . Alizarin Red S, H&E and Sirius Red staining are described in the supporting information. Sections were imaged on an automated slide scanner (Panoramic 250 Flash II, 3Dhistech, Hungary) in the bright eld channel with a magni cation of 20 ×. The polarized images of Sirius Red staining sections were taken on a Wide eld Zeiss Polarization Microscope (Zeiss, Germany) with a magni cation of 10 x under 0°-90°a nd 35°-125°, relying on the birefringence properties of collagen I 51 . Collagen I bers were classi ed from mature to immature with different thresholding values and shown in red, orange, yellow and green in ImageJ. The percentage of the red color of mature collagen I ber was calculated by ImageJ 52 . Collagen I ber parameters, including ber width, ber length, ber straightness, and angle, were analyzed by CT-re reconstruction 53 .

Immunohistochemistry staining
Osteocalcin staining procedure and analysis are described in the supporting information. FIB/SEM: The same actin staining sections were studied using FIB (Helios 5 UX, ThermoScienti c, the Netherland) at the same position as osteocytes (actin staining) and lacunae (SEM) under high vacuum conditions. To protect and lable the surface, a 1 µm thick layer of W was deposited over the ROI, at 30 kV with a beam current of 2.4 nA. The cross-sections of the ROI were rst roughly milled at 30 kV with a FIB beam current of 9.1 nA, followed by a ne polishing at 30 kV with a FIB beam current of 0.44 nA. Serial images were then acquired using the through lens detector (TLD) in backscattered electron (BSE) mode with the following settings: Acceleration voltage 2 kV, Beam current 0.1 nA, Pixel dwell time 10 µs, voxel size: 30 × 30 × 30 nm.

Statistical analysis
GraphPad Prism 8 was used to do the statistical analysis for the obtained data. The comparison of data for NL, ML01 and ML21 groups at one timepoint were done using a two-way ANOVA test together with pair-wise comparison, followed by Tukey corrections. The comparison of data for NL, ML01 and ML21 groups at different timepoints were done using a two-way ANOVA test together with Tukey's multiple comparisons test. The statistical analysis for FE average Young's modulus was using one-way ANOVA together with Newman-Keuls Multiple Comparison Test. * P < 0.05 were considered statistically signi cant, # P < 0.05 when compared between different time points in the same group. Supplementary Materials Table S1. Primer sequences used for the analysis of relative gene expression by RT-qPCR.      scaffolds. *P < 0.05, *** P <0.001, data are shown as mean ± standard deviation (N=4).        Video S1. Low-magni cation images showed that a 3D cellular network was formed in the ML01 scaffold on day 31.
Video S2. High-magni cation images showed how two single cells connect on day 56 in ML01 group.    groups. (F) Correlation of FE stiffness to experimental stiffness of all samples in the NL, ML01 and ML21 groups at day 56. * P < 0.05, ** P < 0.01, data are shown as the mean ± standard deviation (n = 4). # represents a signi cant difference (P < 0.05) compared to the scaffold before loading on day 1.