Biomanufacturing of organ-specific tissues with high cellular density and embedded vascular channels

 advances.sciencemag.org  09/06/2019 17:58:59 

The ability to construct whole organs for therapeutic use remains a daunting challenge, requiring billions of cells to be rapidly organized into functional microarchitected units that are supplied with nutrients via pervasive vascular channels (1). Without a readily perfusable circulatory network, engineered human tissues are limited to several hundred micrometers in thickness (210). This constraint arises due to the delay between implantation and anastomosis with host vasculature, which necessitates a reliance on the diffusive transport of oxygen and nutrients to maintain cell viability (2, 8, 9). Although avascular tissue grafts may provide a measurable improvement in organ function upon implantation (3, 8, 9), the de novo biomanufacturing of three-dimensional (3D) grafts and, ultimately, full-scale organs will inevitably require a perfusable vascular network. While 3D vascularized tissues (~1 cm thick) have recently been fabricated via multimaterial 3D bioprinting (11, 12) and stereolithography (13), they lack the requisite cellular density and microstructural complexity needed to achieve physiologically relevant levels of function. Engineered tissues composed of individual cells suspended in extracellular matrices (ECMs), i.e., so-called cells in gels, typically contain at least one to two orders of magnitude lower cell density than those observed in vivo (1, 14).

Recent advances in the self-assembly of human embryonic stem cells and induced pluripotent stem cells (iPSCs) have led to the development of organoids that have several characteristics akin in their in vivo organ counterparts (1517). Most organoid protocols [e.g., cerebral (16), kidney (1719), and cardiac (20) organoids] begin by generating embryoid bodies (EBs) composed of many iPSCs placed into microwells and cultured under static conditions, whose differentiation can be directed into mini organs of interest. We posit that these organoids may serve as ideal organ building blocks (OBBs) for biomanufacturing patient- and organ-specific tissues with the desired cellular density, composition, microarchitecture, and function provided that a perfusable network of vascular channels can be introduced within these living matrices.

Embedded 3D printing offers one viable strategy for achieving this goal. Lewis et al. first demonstrated this method by writing a viscoelastic, sacrificial (or fugitive) ink within acellular hydrogel (21) and silicone (22) matrices. After printing, these matrices were cured, and the sacrificial ink was removed, leaving behind a 3D network of interconnected channels. Building on this advance, other researchers developed synthetic (23) and biopolymer (24, 25) matrices that exhibit a self-healing, viscoplastic response that simplifies the patterning of complex 3D architectures. However, to date, these methods have only been used to construct acellular (23, 24) or sparsely cellular (25) matrices.

Here, we report a biomanufacturing method that relies on sacrificial writing into functional tissue (SWIFT) composed of a living OBB matrix to generate organ-specific tissues with high cell density, maturation, and desired functionality (Fig. 1). First, we create hundreds of thousands of iPSC-derived OBBs in the form of EBs, organoids, or multicellular spheroids. Next, these OBBs are placed into a mold and compacted via centrifugation to form a living OBB matrix. We then rapidly pattern a sacrificial ink within this matrix via embedded 3D printing, which upon removal yields perfusable channels in the form of single or branching conduits. Last, we demonstrate that these bulk vascularized tissues function and mature when perfused over long durations.

Fig. 1 Sacrificial writing into functional tissue (SWIFT).

(A) Step-by-step illustration of the SWIFT process. (B) (i) Large-scale microwell culture of approximately (ii) 2.5 ml of EB-based OBBs, compacted to form an (iii) OBB tissue matrix composed of approximately half a billion cells. Scale bar, 300 m (i). Scale bar, is 200 m (iii). (C) Time-lapse of sacrificial ink (red) writing via embedded 3D printing within an EB matrix observed from beneath the reservoir. (D) Front view of a vertical line of sacrificial ink printed within an EB matrix. Scale bars, 1 mm (C and D). (E) Examples of the SWIFT process for different OBB-based matrices composed of the following: (i) EBs, (ii) cerebral organoids, and (iii) cardiac spheroids. Row 1: Individual OBBs with characteristic markers. Rows 2 and 3: Cross sections [as indicated in (D) by the dashed line] of immunostained slices and bright-field images, respectively, of the OBB types. Scale bars, 50 m (top row) and 500 m (middle and bottom rows). (F) Generation of a helical (vascular) feature in an EB matrix via SWIFT: (i) CAD representation of the system and (ii) corresponding image of sacrificial ink writing within an EB matrix, and (iii) image sequence acquired during embedded 3D printing of a sacrificial ink (left), sacrificial ink evacuation upon incubation (middle), and tissue perfusion using media (dyed blue) through the printed helical vascular channels.

We grow iPSCs in adherent culture and transfer them into large-scale microwell arrays to form a large volume (> 1 ml) of EBs for SWIFT.After harvesting, these EBs can either be used directly as OBBs (Fig. 1A) or further differentiated into specific organoids of interest. For each SWIFT construct, we cultured and harvested up to 400,000 OBBs, which are subsequently mixed with a tailored ECM solution composed of collagen I and Matrigel held at 0 to 4C. This ECM solution is ideally suited for creating perfusable vascular channels because its fluid-like behavior at low temperatures facilitates casting and embedded printing. However, once the SWIFT construct is placed in an incubator at 37C, the ECM solution undergoes gelation and stiffens markedly to facilitate sacrificial ink removal and subsequent tissue perfusion (3, 4). The cold OBB-ECM slurry is compacted via centrifugation to produce a living tissue matrix (~2.5 ml in volume) that contains nearly half a billion cells at a cell density of ~200 million cells/ml (Fig. 1B). The cold, compacted OBB-ECM slurry provides sufficient support for embedded 3D printing of a sacrificial gelatin ink, which serves as a template for the vascular channel(s) of interest. During SWIFT, we directly observed yielding of tissue matrix roughly 1 mm ahead of the translating cylindrical nozzle and self-healing in its wake (Fig. 1C, fig. S1, and movie S1). Because of its self-healing behavior, the translating nozzle does not generate defects, such as crevasses, within the living tissue matrix. After SWIFT, the tissue construct is warmed to 37C. The sacrificial gelatin ink melts and is removed, leaving behind a network of tubular channels embedded within the tissue construct. The resulting tissue is immediately connected to an external pump and perfused in oxygenated media to maintain cell viability.

SWIFT can be used to embed vascular channels into living matrices composed of a broad range of OBBs. For example, we printed patent vascular channels into three distinct OBB matrices: compacted EBs, cerebral organoids (after 21 days of differentiation), and beating cardiac spheroids that contain iPSC-derived cardiomyocytes mixed with primary cardiac fibroblasts (Fig. 1, D and E). SWIFT does not adversely affect their complex microarchitecture, as evidenced by the presence of intact ventricle-like structures within individual cerebral organoids in these matrices (fig. S2, A and B). Once warmed, the printed tissue is supported by the surrounding ECM (fig. S2C). Given its free-form nature, one can use SWIFT to pattern nearly arbitrary vascular networks, as highlighted by the 3D helical structure printed and perfused within an EB-based matrix (Fig. 1F and movie S2).

The SWIFT biomanufacturing method requires that the OBB-based matrices initially display the requisite self-healing, viscoplastic response to enable sacrificial ink writing followed by a subsequent stiffening of the matrix to facilitate ink removal (Fig. 2) (26). A crucial component of our OBB matrix design is the engineered ECM, whose stiffness is highly temperature dependent. Following the protocol described above, a cold OBB-ECM slurry is created composed of a nearly monodisperse population of EBs (mean diameter of 212 m and a coefficient of variation of 4.7%) (Fig. 2A). Upon compaction, the resulting tissue matrix exhibited a strong shear thinning behavior with a shear yield stress of y ~10 Pa (Fig. 2B and fig. S3). When maintained at 2C, the tissue matrix exhibits a solid-like response, i.e., its plateau storage modulus (G2, ~200 Pa) exceeds its loss modulus (G, ~50 Pa). The sacrificial gelatin ink is designed, such that its shear yield stress (y, ~400 Pa) is at least an order of magnitude higher than that of the matrix (Fig. 2B) (22, 26) to enable the desired cylindrical filaments to be patterned in an omnidirectional manner. During printing, the fluid-like ECM within the tissue matrix contributes minimally to its rheological behavior (fig. S3, B and C). In the absence of densely packed OBBs within the matrix, the cold ECM solution (held at 2C) exhibited a purely viscous response (, ~1 Pa) that is not suitable for embedded 3D printing (21). However, upon printing and warming the tissue matrix to 37C, its plateau shear elastic modulus (G2) increases by roughly three orders of magnitude to a value of 280 kPa (Fig. 2C and fig. S3D), as the ECM undergoes thermal gelation to produce a fibrillar gel that occupies the interstitial space between densely packed OBBs. As stated previously, the pronounced increase in tissue matrix stiffness allows it to retain its structural integrity during removal of the sacrificial ink and support pressure-driven flow during subsequent perfusion. As controls, we measured the rheological properties of compacted EBs alone (no ECM) and found that its stiffness substantially increased upon heating from 2 to 37C. By contrast, the pure ECM rapidly stiffened by roughly two orders of magnitude to a value of 110 Pa under those same conditions (fig. S3E). Hence, by tailoring the tissue matrix composition and sacrificial ink rheology, we developed a system that allows vascular channels to be printed in any arbitrary direction via SWIFT (Fig. 2D). To demonstrate this, we patterned embedded channels with smoothly varying diameters between 400 m and 1 mm using a 250-m nozzle, through which the sacrificial ink is deposited at a constant volumetric flow rate and varying print speed (Fig. 2E). While there is no upper limit on channel diameter, those below 400 m could not be printed with high fidelity since their size approaches the characteristic diameter (d) of the OBBs, in this case ~200 m. While smaller EBs would enable higher-feature resolution during SWIFT, the volume of the EB matrix generated decreases as d3 for a fixed number of microwell arrays. For the print speeds evaluated (0.5 to 4 mm/s), the shear stress imparted by the nozzle does not adversely affect either the cell viability or matrix integrity, as shown in Fig. 2Eii and fig. S4. The embedded channel features can also be seamlessly connected to form branching vascular networks (fig. S5). We note that the viscoplastic nature of the tissue matrix also allows dense OBB-based inks to be directly printed within it akin to the sacrificial ink. As a simple example, we printed an ink composed of compacted fluorescent proteinexpressing EBs into a tissue matrix composed of wild-type (nonfluorescent) EBs in the form of a helical pattern (fig. S6).

Fig. 2 Living matrix and ink rheology for SWIFT.

(A) Size distribution (n = 413 EBs) of EBs used to form EB matrices. (B) (i) Apparent viscosity as a function of shear rate and (ii) shear storage (closed markers) and loss moduli (open markers) as a function of shear stress of the EB matrix and sacrificial gelatin ink. (C) Temperature effects on the plateau storage moduli (or loss modulus indicated by an asterisk) of the EB matrix and the sacrificial gelatin ink. (D) SWIFT printing of (i) horizontal and (ii) vertical features (vascular templates) embedded at print speeds of 0.5, 1, 2, and 4 mm/s. (E) Effect of print speed on the lumen diameter shown as (i) bright-field and (ii) viability staining images in the context of vertically printed channels and (iii) lumen (channel) diameter as a function of print speed for vascular templates embedded via horizontal and vertical SWIFT printing. Error bars indicate SD (n = 4). Scale bars, 2 mm (D and E). EthD-1, ethidium homodimer-1.

To determine cell viability, we manufactured a perfusable EB-based tissue in a perfusion chamber that enables circulation of oxygenated media through the embedded vascular channels and around the periphery of the tissue (Fig. 3A and fig. S7). As a control, we cast a nonperfusable EB-based tissue (referred to as no channel) and only flowed media around its periphery in the same chamber (Fig. 3A, top). We find that the control tissue (cross section, 4.2 mm by 4.2 mm) developed a necrotic core within 12 hours with viable, metabolically active cells (per calcein-acetoxymethyl (AM) ester staining) restricted to a narrow band (~0.8 mm thick) at the tissue periphery (Fig. 3B and fig. S8). By contrast, luminal perfusion of the perfusable EB-based tissue enhanced cell viability throughout the bulk tissue (Fig. 3B). Their high cell density made it difficult to image the core of individual EBs within these matrices (see fig. S4). Hence, we reported a tissue viability based on the ratio of a thresholded calcein signal (viable cells) to a thresholded nuclear signal (total cells) normalized by the value derived for the full cross section of an avascular tissue composed of the same EBs (Fig. 3Biv). Notably, similar cell viability values are observed initially (t = 0 hours) between the perfused and nonperfused (control) tissues (fig. S8) even for EBs located adjacent to the printed lumen, indicating that the SWIFT process has no adverse effects (fig. S4). We also find that tissues subjected to perfusion of normoxic (21% O2, 5% CO2) or hyperoxic (95% O2, 5% CO2) media (Fig. 3Biv) exhibited similar cell viability values.

Fig. 3 Perfusable EB tissue fabricated by SWIFT.

(A) Perfusion system used to assess tissue viability following SWIFT printing. (B) Viability staining and analysis following 12 hours of culture of a tissue featuring (i) no channel or perfused with either (ii) normoxic (21% O2) or (iii) with hyperoxygenated (95% O2) media along with the corresponding quantification of the (iv) normalized viability. Scale bars, 500 m. Error bars indicate SD (n = 4). The dashed lines highlight viability regions that arise from external perfusion. The core-only region corresponds to the area located within the innermost line. (C) An image sequence showing the embedded 3D printing of a branched, hierarchical vascular network within a compacted EB-based tissue matrix connected to inlet and outlet tubes, seen entering the tissue from the left and right. Scale bar, 10 mm. (D) Image of the perfusable tissue construct after 12 hours of perfusion (top) and fluorescent image of LIVE/DEAD (green/red) cell viability stains at various sections through the tissue (bottom). The dashed line represents the equivalent viability depth for an avascular control perfused only from the outside, see (Bi). Scale bars, 1 mm (E and F). (E) SWIFT printing of a bifurcating channel for lumen endothelialization. (F) Evacuated channel (highlighted by the white dashed lines) undergoing the perfusion of HUVEC cells. Scale bar, 1 mm. (G and H) Formation of a VECad-positive monolayer of an HUVEC endothelium. Scale bars, 500 m (G) and 50 m (H).

To scale up the SWIFT biomanufacturing process, we created a custom-shaped mold (4 mm thick) with a single inlet and outlet for perfusion (fig. S9) and filled it with 2.5 ml of our compacted EB-ECM slurry composed of ~400,000 EBs and a total of 0.5 billion cells. We then printed a branched, hierarchical channel network that is designed to both distribute flow and maintain a constant wall shear stress throughout the tissue matrix (Fig. 3C and movie S3). After printing, the tissue construct is warmed, the sacrificial ink is removed, and the resulting channels are perfused with hyperoxygenated medium (95% O2, 5% CO2) at a flow rate of 250 l/min. After 12 hours of constant perfusion, the embedded vascular channels remained patent, and the tissue contracts slightly (Fig. 3D, top) as adjacent OBBs fuse together and undergo remodeling. Similar to the smaller tissues (0.2 ml in volume) shown in Fig. 3Bii and iii, cells within these bulk perfusable tissues (2.5 ml in volume) remain viable throughout most of their 4-mm thickness, as indicated by the fluorescence of calcein-AM (Fig. 3D, bottom).

We introduced human umbilical vein endothelial cells (HUVECs) to form endothelial-lined channels within these SWIFT constructs. Specifically, the endothelial cells are seeded by first perfusing a suspension of HUVECs (107 cells/ml) through an embedded, bifurcating vascular network (Fig. 3, E to H, and movie S4) within a SWIFT construct. The cells were allowed to adhere to the luminal surface. Following 20 hours of perfusion, we sliced the opaque tissue longitudinally and directly imaged the endothelialized lumen. While endothelial cells do form a characteristic cobblestone pattern on some portions of the luminal surface, as indicated by vascular endothelial cadherin (VECad) staining (Fig. 3, G and H), their coverage was incomplete. Our ongoing efforts are focused on producing confluent endothelium within these printed vascular channels.

To further demonstrate SWIFT, we created a functional, perfusable tissue containing human iPSC-derived cardiac OBBs (Fig. 4). Because of its scalability, we adapted a growth factorfree suspension culture protocol that starts from EBs and modulates the Wnt pathway via addition of small molecules to the cell culture medium (27), a strategy that yields cardiomyocytes with a high efficiency (Fig. 4, A to D) (28). Using this protocol, we generated up to 6 ml (upon compaction) of beating cardiac OBBs (movie S5), containing 79 6% cardiomyocytes (cardiac troponin T-positive, cTnT+) and 19 6% stromal cells (cTnT, Vimentin+) (Fig. 4D and fig. S10). These OBBs were mixed with our engineered ECM and human neonatal dermal fibroblasts (HNDFs) to form a slurry that was subsequently compacted via centrifugation yielding a cardiac tissue matrix (fig. S11). Once compacted, the cardiac matrix had an estimated cell density of 240 million cells/ml and a cardiomyocyte density of 180 million cells/ml (Fig. 4E). When formed into a small disc-shaped mold, this cardiac OBB construct initially beats asynchronously, with calcium waves propagating sporadically and chaotically after a single day (day 1) of perfusion (fig. S12 and movie S6). However, after 7 days in culture (day 7), the tissue construct beats spontaneously and synchronously with calcium waves propagating rhythmically and rapidly throughout the tissue (fig. S12 and movie S6).

Fig. 4 Perfusable cardiac tissue fabricated by SWIFT.

(A) Cardiac organoid differentiation protocol. (B) Cardiac troponin T and 42,6-diamidino-2-phenylindole (DAPI) staining in a single cardiac OBB at day 9. Scale bar, 50 m. (C) Cardiac troponin T, -actinin, and DAPI staining in a single cardiac OBB at day 9. Scale bar, 10 m. (D) Cardiac spheroid composition in iPSC-derived cardiac OBB. Cardiomyocytes (CM) are identified as cardiac troponin T-positive (cTnT+) and stromal-like cells (strom.) as cTnT/Vimentin+. (E) Cellular density in compacted cardiac OBB tissue. (F) An image sequence showing the embedding, evacuation, and perfusion of branched vascular channels within a cardiac tissue matrix (tissue dimensions: top width, 6 mm; bottom width, 4.2 mm; depth, 4.2 mm; and height, 12 mm). Scale bars, 2 mm. (G) Viability staining of a SWIFT cardiac tissue (cross section) after 24 hours of perfusion. Scale bar, 500 m. (H) cTnT, -actinin, and DAPI staining in a SWIFT cardiac tissue after 8 days of perfusion that shows evidence of sarcomeric remodeling (arrowheads). Scale bar, 10 m. (I) Vertical displacement of the anchoring flexible prongs due to spontaneous cardiomyocyte contraction showing increasing amplitude over time. On day 8, 2 mM calcium is added to the medium to increase cardiomyocyte contractility (d8 + Ca). (J) Comparison of anchor displacement pattern between spontaneous contraction and electrical pacing (1 and 2 Hz) of SWIFT cardiac tissues. (K) Spontaneous contraction pattern before and after administration of 10 M isoproterenol. (L) Average contraction frequency under isoproterenol treatment. (M) Spontaneous contraction pattern before and after administration of 1 mM 1-heptanol. (N) Maximum peak-to-peak contraction amplitude under 1-heptanol treatment. (O) 3D CAD model of a normal human heart, including a segment of the left anterior descending (LAD) artery and a diagonal branch, downloaded from the National Institutes of Health 3D Print Exchange (additional septal and diagonal branches were added manually, pink). (P) A 1:2 scale polydimethylsiloxane mold is formed using the 3D computed tomography data, and the LAD artery together with diagonal and septal (arrowheads) branches are embedded into a septal-anterior wall wedge [yellow section in (O)] of the cardiac tissue matrix via SWIFT. Scale bar, 5 mm.

Next, we embedded, evacuated, and perfused a branching channel architecture within this cardiac OBB matrix (Fig. 4F and movie S7). To generate multiple samples for analysis, we produced simple tissue constructs that each contain a single printed channel, which is perfused by a luminal flow of 500 l/min, resulting in a moderate shear stress of ~0.8 dyn/cm2. Each T-shaped tissue is connected to an inlet tube at the bottom and was mounted at the top via a pair of anchoring prongs that serve to both stabilize the tissue and allow one to monitor its contraction and beating over time (fig. S13, A and B). After 1 day of perfusion, the bulk tissue remained viable (Fig. 4G). Over an 8-day period, the cardiac tissue developed a pervasive sarcomeric architecture (Fig. 4H and fig. S13, C and D). Over this same period, the overall tissue contractility increased by more than 20 times, as visualized by the amplitude of anchoring prong motion (Fig. 4I, fig. S13B, and movie S8). A ~40% enhancement in beating capability is obtained upon adding 2 mM calcium to the medium. Beating synchronization also increased, as suggested by coordinated contraction patterns on opposite sides of the tissue section at day 8 (fig. S14). In addition to spontaneous contraction, the tissue responded to paced stimulation via a pair of platinum electrodes situated on either side of the tissue (Fig. 4J and movie S8).

To further characterize the functional performance of cardiac tissue constructs produced by SWIFT biomanufacturing, we perfused isoproterenol, a nonselective agonist for the -adrenoreceptor, through an embedded channel at a concentration of 105 M, which resulted in a doubling of the spontaneous beating frequency (Fig. 4, K and L). We also perfused 1-heptanol, a gap junction blocker, at a concentration of 1 mM and observed a reduction in their contractile amplitude (Fig. 4, M and N). As a final demonstration, we printed an arterial vascular network geometry within a cardiac OBB matrix using patient-specific, cardiac structural data available from the National Institutes of Health (NIH) 3D Print Exchange (Fig. 4, O and P). We used a 3D printed mold to replicate a wedge of the myocardium at a 1:2 scale and, upon filling the mold with a compacted cardiac OBB matrix, used SWIFT to replicate the geometry of the left anterior descending (LAD) coronary artery, which passes anteriorly and then septally and lies adjacent to the outer edge of the cardiac tissue.

The SWIFT biomanufacturing method uses densely cellular matrices composed of iPSC-derived OBBs to produce autologous tissues composed of patient-specific cells. This method enables one to create perfusable organ-specific tissues of arbitrary volume and shape in a scalable manner. We estimate that large volumes (100 ml) of perfusable OBB-based tissues could be patterned with an embedded vascular channel network. However, successful translation of these organ-specific tissues for therapeutic applications will require additional advances. For example, current protocols for iPSC-derived organoids yield OBBs that lack sufficient cell maturation and microvascular network formation. Hence, only a modest contractility (~1% strain) is observed for SWIFT cardiac tissues compared to the contractility (~20% strain) observed for adult cardiac tissue (29). We are currently focused on generating anisotropic cardiac OBBs with the aim of further enhancing tissue function, including contractility. In addition, we envision adopting new protocols, such as those recently reported for renal organoids (19), that offer a pathway to creating more mature, microvascularized OBBs. Alternatively, one could enhance tissue function by using living matrices composed of primary cell spheroids harvested from adult tissues for SWIFT.The field of organ engineering is rapidly evolving. For example, other researchers recently reported the embedded 3D printing of a large heart-shaped construct using the traditional cells in gel approach. Although this work received considerable attention, those tissues were not perfused long term in vitro, and their cell viability and contractile function were not reported (30).

In summary, we have demonstrated a new biomanufacturing method, known as SWIFT, that uses iPSC-derived OBB tissue matrices that exhibit the requisite cell density, microarchitecture, and function approaching that of native tissues. SWIFT can be implemented with a wide range of OBBs, including EBs, differentiated organoids, and multicellular spheroids. Our method opens new avenues for creating personalized organ-specific tissues with embedded vascular channels for therapeutic applications.

Microwell array plates were used to generate large numbers (>105) of uniform EBs via forced aggregation. The microwells, in the form of inverted pyramids with 400-m base width, were arranged in a six-well plate format and manufactured via polydimethylsiloxane (PDMS) molding. The following process was used to facilitate demolding and ensure casting compatibility between polymers. Briefly, SYLGARD 184 (Dow Corning) was mixed 10:1 with its cross-linker and poured into a commercially available microwell plate, degassed in a vacuum chamber, and cured at 80C for 2 hours to generate a negative cast of the microwells. Next, Ecoflex 00-30 (Smooth-On Inc.) was used to create a positive mold by following a similar method of pouring, degassing, and curing as for the SYLGARD 184. After removing the cured Ecoflex, a two-part polyurethane (Smooth-On Inc.) was poured onto the Ecoflex positive to create a polyurethane-negative mold. Next, ~80 g of SYLGARD 184 was mixed with 8 g of its cross-linker in a planetary mixer, poured into a polyurethane negative mold, degassed for 10 min in a vacuum, and allowed to cure at room temperature for 24 to 48 hours. The cured SYLGARD 184 was then placed in the oven at 80C for >2 hours to ensure complete cross-linking before removal from the polyurethane mold.

To prepare the microwell arrays for cell culture, the arrays were immersed in isopropyl alcohol for 1 hour and then immersed in water and autoclaved for 1 hour. Immediately before preparing the cells, the autoclaved water was drained from the wells, and then, they were placed into a sterile one-well tray (Omnitray, Falcon). Next, 1.5 ml of AggreWell Rinsing Solution (ARS; STEMCELL Technologies) was added per well before centrifuging the wells at 2000g for 5 min to remove trapped air bubbles from the wells. The ARS was aspirated, and the wells were rinsed twice with 2 ml of Dulbeccos modified Eagles medium (DMEM). Last, the DMEM was replaced with 3 ml of EB culture medium (EBCM) supplemented with 10 M ROCK inhibitor (ROCKi; Y27635; EMD Millipore), and the microwell arrays were maintained at room temperature. To prepare EBCM, mTeSR1 was mixed 1:50 with a polyvinyl alcohol (PVA) stock solution to generate a solution of 4 mg/ml. Penicillin-streptomycin was added to 100 U/ml, and the solution was sterile-filtered and stored at 4C for up to 2 weeks. The PVA stock solution was prepared by adding 20 g of PVA (30 to 70 kDa; Sigma-Aldrich) to 80 ml of deionized water while stirring at room temperature. Under continuous stirring, the solution temperature was increased to 80C to fully dissolve PVA in water, after which the PVA solution was stored at 4C until ready for use.

To form EBs, PGP1 iPSCs were grown to 60 to 80% confluency in Matrigel-coated 225-cm2 T225 flasks (Falcon). One T225 of PGP1 cells was typically used to seed 24 wells (four six-well plates), each containing ~4200 microwells. To lift off the iPSCs from the substrate, the cells were first rinsed in PBS without calcium or magnesium and incubated in Accutase (Innovative Cell Technologies) for 15 min at 37C to generate a mostly single-cell suspension. The cells were added to prewarmed DMEM/F12 and Hepes and centrifuged at 220g for 5 min. The supernatant was removed, and the cells were resuspended in EBCM and10 M ROCKi. Cells were counted using a LIVE/DEAD imaging system, and 1 ml of the cell suspension is seeded into each well at a cell density of approximately 500 cells per microwell. The well plates were then centrifuged at 100g for 3 min to compact the cells. After 24 hours, the cell media was replaced with 4 ml of fresh EBCM to remove the ROCKi. The well plates then underwent two half-media changes (2 ml of EBCM is replaced with fresh media) per day. Care was taken during media changes with media added and removed slowly as to not disturb the EBs from their individual microwells. After 4 days of culture, a total of eight six-well microwell array plates could generate a sufficient volume of EBs to render ~1 ml of OBB slurry.

To generate cerebral organoids, PGP1 iPSCs were grown to ~80% confluency and, on day 0, were lifted from Matrigel-coated T225 flasks by incubation with gentle cell dissociation reagent for 12 min. The flasks were rinsed with DMEM/F12 containing Hepes and centrifuged for 5 min at 220g. Next, the cells were resuspended in EBCM and 10 M ROCKi and seeded into microwell arrays, prepared as described above, at a density of approximately 500 cells per microwell, and compacted by centrifugation at 100g for 3 min. The media was replaced on day 1 with fresh EBCM without ROCKi. On day 2, the EBs were measured to be approximately 200 m in diameter, and the media was replaced with neural induction medium (NIM) composed of DMEM/F12 with GlutaMAX, 1:1000 of heparin (10 mg/ml; Sigma-Aldrich), 1:100 N-2 supplement (Thermo Fisher Scientific), and 1 minimum essential medium nonessential amino acids (MEM-NEAA) (Thermo Fisher Scientific) with 1:2000 of 10 mM SB431542 in dimethyl sulfoxide (DMSO) (SB, EMD Millipore) and 1:2000 of 0.2 mM LDN 193189 in DMSO (LDN, EMD Millipore). Aliquots of SB and LDN were stored at 20C and added to the other media components immediately before adding NIM to cells. On day 3, the OBBs were removed from the microwells by pipetting up and down with DMEM/F12 containing 15 mM Hepes buffer, and the aggregates were transferred to a 50-ml centrifuge tube. The OBBs were allowed to settle under gravity, the supernatant was aspirated, and then, the OBBs were resuspended in 6 ml of NIM per two microwell arrays that were harvested. Next, they were seeded into untreated T25 flasks, and 6 ml of NIM was mixed with the OBBs per flask. After harvesting, an additional 5 ml of NIM was pipetted into each flask. The flasks were transferred to an orbital shaker (VWR Shaker, 1000 Standard) that was set to 60 rpm and placed in an incubator. The media was replaced on days 5 and 7 by placing the flasks on a custom-built sedimentation rack that held the flask at 45 such that cells settle into a single corner of the T25 flask. Once settled, the supernatant was aspirated and replaced with 6 ml of fresh NIM. On day 8, the medium was replaced with neural differentiation medium (NDM) composed of a 1:1 mix of DMEM/F12 with Neurobasal medium (Thermo Fisher Scientific), 1:200 N-2 supplement (Thermo Fisher Scientific), 1:100 B27 without vitamin A (Thermo Fisher Scientific), 1:4000 insulin solution (Sigma-Aldrich), 1:100 GlutaMAX, and 1:200 MEM-NEAA. A 1000 solution of -mercaptoethanol (1:1000; Thermo Fisher Scientific) was added to NDM immediately before adding the media to the OBBs. Thereafter, the media was replaced every other day, and the cerebral organoids were harvested for SWIFT vascularization on day 21.

To create the cardiac spheroids used for Fig. 1, D and E and fig. S2, we adapted a previous protocol by Lian et al. (28). Briefly, PGP1 iPSCs were grown to 80 to 95% confluency on T225 plates in mTeSR1 medium. Cardiac differentiation media was prepared by combining 500 ml of RPMI 1640 (Thermo Fisher Scientific) with either 10 ml of B27 without insulin (to form cardiomyocyte differentiation medium, CDM) or 10 ml of B27 (to form cardiomyocyte maturation medium, CMM). On day 0, differentiation was initiated by switching the mTeSR1 media to CDM and 6 M CHIR99021 for 24 hours. The following media changes were then made: day 1, CDM; day 3, CDM and 2 M iWR1 (BioGems; CDMI); days 5 and 7, CDM; day 9 (or by day of beating), CMM; days 11 and 13, RPMI without glucose supplemented with 2% B27 (to purify cardiomyocytes); and day 15, CMM. In a separate T225 flask, primary ventricular normal human cardiac fibroblasts (Lonza) were cultured in stromal cell growth medium (Lonza). Next, cardiomyocytes, at day 17, and fibroblasts were passaged by rinsing with PBS and then incubated with TrypLE for 9 min (for cardiomyocytes) or 5 min (for fibroblasts) at 37C. Next, the TrypLE was rinsed using 37C DMEM, and the cells were pelleted into separate 15-ml tubes via centrifugation at 220g for 5 min. Both cell types were resuspended in cardiac spheroid medium (CSM), comprising DMEM with high glucose (American Type Culture Collection) containing 10% fetal bovine serum (FBS) and penicillin-streptomycin (100 U/ml). The cells were counted and combined at 70% cardiomyocytes and 30% fibroblasts and then seeded into molded 800-m-diameter PDMS microwells at a cell density of 2000 cells per microwell. The microwell arrays were centrifuged at 100g for 3 min to form OBBs from the cells. CSM was half-changed daily, and the cardiac spheroids were observed to be beating by 2 to 4 days after aggregation, at which point they were harvested for SWIFT vascularization.

To demonstrate the versatility of this process, we also generated cardiac OBBs according to a protocol adapted from (27) and (28). BJFF iPSCs were formed into EBs on day 4, as described above. Two T225s (at ~80% confluency) were used to seed four six-well plates of microwells. On day 3, to remove ROCKi, the media was replaced with 4 ml of EBCM. On day 2, EBs from two wells were harvested, combined, and allowed to settle under gravity. The supernatant was aspirated, and the EBs from two wells were resuspended in 12 ml of EBCM and transferred to a single untreated T25 flask. Thus, EBs from four six-well plates (i.e., 24 wells in total) were harvested and placed into 12 T25 flasks. These flasks were incubated on an orbital shaker with a rotational speed of 60 rpm. To maintain and differentiate the cells in the T25 flasks, the following media changes were performed: day 1, 12 ml of EBCM; days 0 and 1, 10 ml of CDM and 5 M CHIR99021 (BioGems) (CDMC); day 2, 10 ml of CDM; days 3 and 4, 10 ml of CDMI; and day >5, 10 ml of CDM daily. On days 0, 2, and 5, a rinse was performed after aspiration, using 4 ml of the medium that the cells were being changed into, to increase the removal of growth factors or small molecules from the previous medium. After the cells resettled, the rinsing medium was removed before the new culture medium was added. Beating was typically initially observed by day 7, and almost all cardiac OBBs were visibly beating by day 9. SWIFT printing was performed using cardiac OBBs between days 9 and 11, after which the cells were maintained in CMM.

An ECM prepolymer gel was prepared as follows: A transglutaminase (TG) working solution was first prepared by adding a microbial tissue TG formulation (20 mg/ml; Moo Gloo TI) to cell culture medium (mTeSR1 for EBs, NIM for cerebral organoids) and sterile-filtering the solution via a 0.2-m syringe filter. For cardiac OBBs, TG was not added to the ECM. Next, the volume (x milliliters) of high-concentration rat tail collagen I (Corning) required to create a collagen concentration (4 mg/ml) in the final ECM solution (total volume, y milliliters) was calculated. The ECM prepolymer gel solution composed of (0.825y  1.135x) ml of mTeSR1, 0.125y ml of TG working solution (if applicable), x/10 ml of 10 PBS, y/100 ml of 250 mM CaCl2, ~x/40 ml of 1 N NaOH, x milliliters of high-concentration rat tail collagen type I, and 0.04y ml of Matrigel was added, in order, on ice. For example, to generate 8 ml of ECM prepolymer gel using a solution of collagen I (8.90 mg/ml), the following was added in order: 2.56 ml of mTeSR1, 1 ml of mTeSR1 containing TG (20 mg/ml), 0.36 ml of 10 PBS, 0.080 ml of CaCl2, 0.090 ml of 1 N NaOH, 3.596 ml of collagen I, and 0.320 ml of Matrigel. The volume of NaOH was adjusted, as necessary, to attain a final pH of ~7.4. To prevent premature gelation, the pipette was precooled using ice-cold PBS before pipetting the collagen and Matrigel. After addition and mixing of all components, the gel was centrifuged at 2000g for 3 min at 2C to remove air bubbles. The ECM prepolymer gel was stored on ice and used within ~1 hour.

The EBs were harvested from the microwell arrays once they grew 200 to 250 m in diameter (typically 2 to 4 days in microwell culture) by vigorous pipetting up and down in the wells with a P1000 pipette tip. The EB-containing media was transferred into a 50-ml tube, one tube per six wells. The wells were further rinsed with an excess of DMEM/F12 and Hepes using a 10-ml serological pipette to remove most (>90%) of the EBs. The media, with suspended EBs, was added to 50-ml centrifuge tubes, and the EBs were allowed to settle by incubating on ice for 10 min. All further cell handling was performed on ice to prevent gelation of the ECM. The supernatant was aspirated, and the EBs were collected into a single 50-ml tube. Next, 50 ml of fresh DMEM/F12 and Hepes was added to further wash away single cells. The EBs were then allowed to settle for another 10 min, and then, the supernatant was aspirated, and the EBs were resuspended up to 15 ml using DMEM/F12 and Hepes and transferred into a 15-ml tube and allowed to settle for 10 min. Next, to compact the EBs, the 15-ml tube was centrifuged at 100g for 3 min at 2C, the supernatant was aspirated, and the total volume of the EB matrix was estimated. To replace the DMEM/F12 and Hepes in the extracellular space with ECM, the cells were rinsed in the ECM solution at a 1:3 (EB matrix:ECM solution) volume ratio by pipetting with a P1000 pipette and then centrifuged at 100g, and the supernatant was aspirated. Last, the consolidated EB-ECM solution was resuspended in an equal volume of ECM solution, forming a dense EB-ECM slurry that could be handled and pipetted in subsequent steps. To estimate the cell density of the OBB matrix, the EB-ECM slurry was centrifuged at 150g for 3 min in an Eppendorf tube, and the OBB matrix volume was estimated visually by pipetting known volumes of water into a second Eppendorf tube to match the height of the pellet in the first tube. Next, the EBs were dissociated into single cells by incubating in Gentle Cell Dissociation Reagent for 15 min and triturating with a P1000 pipette, and the total cell number was counted. The cell density of the OBB matrix was then estimated by the ratio of cell number to OBB matrix volume.

Direct ink writing was used to manufacture customized silicone molds for printing and perfusing bulk tissues. First, a printable silicone ink was prepared by combining a 10:1 mass ratio of SE1700 base:curing agent (Dow Corning) and then adding 1% (w/w) black silicone pigment (Silc-Pig, Smooth-On). The components were then mixed using a speed mixer for 3 min at 2000 rpm (THINKY Inc.). Using a spatula, freshly mixed silicone ink was loaded into a 30-cm3 syringe (Nordson EFD) and centrifuged at >3000g to remove entrapped air. A tapered nozzle with a 0.41-mm inner diameter tip (Nordson EFD) was attached to the outlet of the syringe, and the syringe was affixed to a custom 3D printer (26). Ink was dispensed from the syringe by means of pressure generated by an Ultimus V pressure controller (Nordson EFD) onto an underlying glass substrate. The tissue molds are composed of silicone walls and small channels for tube-insertion sites, i.e., inlets and outlets. After printing, each mold was transferred to 80C for >2 hours to cure. After curing, fresh silicone ink was extruded onto the top layer of the mold (either manually or using the 3D printer), and a second glass slide was added to the top, such that the silicone gasket was sandwiched between two glass slides. For the mold used in Fig. 3, clear silicone was printed and cured onto the second glass slide before assembly, such that the assembled gasket had a slit in the middle of the top side for nozzle insertion. The mold was transferred to 80C to cure the silicone. Stainless steel tubes serving as inlets and outlets with an inner diameter of 0.83 mm were inserted into the inlet/outlet channels of the molds used for perfusion, and Versilon silicone tubing, with an inner diameter of 1/32 inches, was fitted to the external ends of the tubes and crimped with a hose clamp to seal the inlets and outlets. For the custom perfusion mold in Fig. 3, the stainless-steel tubes had a small ~1-mm notch cut into one end to facilitate a continuous connection from the printed vascular channels and the tube. Molds were autoclaved and stored at room temperature before use.

Alternately, polycarbonate chambers and front and back plates were machined to produce custom perfusion devices (fig. S7). The use of glass windows allowed for tissue visualization, and watertight seals were formed between the windows and chamber walls via o-rings. Metal tubes serving as fluidic inlets and outlets were press-fit on the side and bottom of the chamber. The perfusion chambers were autoclaved and stored at room temperature before use. Before printing vascular channels inside the OBB-based matrices, an anchoring structure was printed using a stereolithography 3D printer (Perfactory Aureus, Envisiontec. Inc), using E-Shell 300 clear resin, postcured under ultraviolet (UV) illumination, and positioned within the chamber. This structure had a pair of anchoring prongs for embedding inside the tissue, which served to hold the tissue in place during perfusion, while allowing for media perfusion around the outer surfaces of the tissue. Next, a tissue template (which ultimately forms the shape of the tissue) was stereolithographically printed, postcured under UV illumination, plasma-treated, and dipped into a 10% (w/w) aqueous solution of Pluronic F-127 (Sigma) and allowed to air-dry. This Pluronic F127coated template was inserted into the chamber and was fitted over the anchoring prongs. Next, 5% (w/v) molten gelatin, prepared by adding PBS to a stock solution of 15% (w/v) gelatin, was poured into the chamber and injected into the side metal tubes and allowed to gel around the template at 4C for 15 min. The template was then removed by pulling straight up, leaving a tissue mold with a custom tissue geometry and anchoring prongs protruding at its top. An o-ring was placed around the tissue inlet tube located at the bottom of the mold after which the printer was aligned to the chamber. To form a tissue-anchoring plug around the inlet, a solution of high-concentration ECM (hcECM) was prepared similar to the ECM described above, but with collagen (7.5 mg/ml) and 2% (w/v) TG [from a stock of TG (160 mg/ml) in medium]. This material was dispensed at the base of the tissue mold and around the O-ring before gelling at 22C for 15 min. Last, 5% (w/v) molten gelatin was injected into the bottom metal pin before filling the mold with an uncompacted mix of cell aggregate and ECM solution for SWIFT printing.

A sacrificial gelatin ink was prepared using the procedure described above to enable embedded 3D printing of vascular networks. The gelatin stock solution (solubilized at 85C, as described above) was melted at 80C for <30 min and then diluted in PBS to generate a 5% (w/v) solution. Next, red food coloring was added at 2% (v/v) to enable visualization of printed ink. The freshly prepared sacrificial gelatin ink was loaded into a 0.5-cm3 glass syringe and subsequently cooled to 4C for 15 min to induce gelation. The glass syringe was then loaded onto a custom-built syringe pump that was mounted onto the 3D printer. The syringe pump was operated by an Arduino microcontroller and a stepper-motor driver. The syringe barrel was housed in a water-cooled custom-built temperature controller to control the stiffness of the gelatin ink by maintaining a temperature of approximately 20C. A metal nozzle, with a length of 1.5 inches and an inner diameter of 0.25 mm (Nordson EFD), was fitted to the bottom of the syringe. Before printing, the syringe was incubated in its temperature-controlled housing for >15 min.

The OBB-ECM matrix was loaded into the silicone molds or perfusion chambers, and the following steps were performed rapidly to ensure that it remained cold and that the ECM did not undergo gelation before SWIFT was completed. First, for cardiac tissue fabrication, the volume of cardiac OBBs was estimated by compacting via centrifugation at 100g. Next, the ECM solution was mixed with 107 HNDF cells/ml of cardiac OBBs before resuspending the cardiac OBBs in the ECM:HNDF cell mixture to generate a slurry. The freshly prepared OBB slurry was pipetted into the molds, and the molds and their housing were centrifuged at 150g at 2C for 3 min to compact the aggregates, forming an OBB-based matrix, containing ~2.4 108 0.4 108 cells/ml and 6 105 2 105 HDNF cells/ml. Most of the supernatant ECM was removed leaving ~1 mm of the ECM solution above the top of the compacted bed of OBBs. To keep the matrix cool during the long-duration print in Fig. 3 (C and D), the mold was transferred onto the front of a custom-built water-cooled metal plate, comprising a metal plate, backed by an acrylic sheet that had milled channels through which ice-cold water was pumped. Both the metal plate and the acrylic sheet were cut to enable back illumination to facilitate visualization of the printing process. The mold was held on the water-cooled plate by a pair of clamps and was maintained in its housing in a vertical orientation by pressing it against two metal right-triangle prisms.

The desired vascular network was then printed within the OBB-based matrix at a print speed ranging from 0.5 to 4 mm/s. The extrusion rate of the syringe pump was set to define a certain diameter filament. The vascular geometry was designed in a custom-built MATLAB script that generated branched cubic splines, automatically assigned a printing order to prevent the nozzle from translating through previously printed features, and rendered each curvilinear line as piecewise linear segments in GCode. To vary the filament (and, hence, channel) diameter, the print speed was changed while maintaining a constant ink deposition. After each line segment, ink extrusion is stopped, and to avoid agitating the printed patterns, the nozzle was withdrawn vertically above the top of the tissue, translated to its new x and y coordinates, and reinserted to the correct new z coordinate by moving down vertically.

When templating a complex hierarchical vascular network (Fig. 3, C and D), the channel diameters were computed to ensure an equivalent shear stress within each segment of the network, independent of their diameter. Assuming Poiseuille flow, we determined the channel radius R and the flow Q, along with the pressure P at each node by solving the following set of equations (in MATLAB). On the basis of mass conservation, the sum of flow coming to a node is zeroiQi=0

The flow Q through each segment is proportional to the fourth power of its radius R4 and the difference of pressure across P, as given byQ=R48lPwhere is the viscosity of the culture medium and l is the length of a given segment. The shear stress exerted on a vessel wall, =4QR3, should be constant throughout the network, henceQiRi3=Q1R13

For the vascular network depicted in Fig. 3 (C and D), the ratio between the largest vessels (inlet and outlet) and the smallest one was found to be 1.3. To print these channels at a constant deposition rate, the nozzle velocity was modulated proportionally to the inverse of the square of the diameter of each channel (or to its cross-sectional area).

After SWIFT printing of cardiac tissues in machined perfusion chambers, the tissue was capped with hcECM to provide mechanical stability required to maintain patency for long-term culture. The mold was then transferred to a 37C, 5% CO2 incubator for 45 min to complete gelation of the OBB-ECM matrix, and to melt the printed sacrificial gelatin filaments and gelatin that surrounds the tissue. After incubation, silicone molds were sealed via an excess of autoclaved SE1700 silicone base (without cross-linker) into the gap at the top of the mold, while machined perfusion chambers were sealed by screwing on the polycarbonate lid. The sealed perfusion devices were then transferred to a custom-made tissue perfusion housing that held the tissue vertically and adjacent to a media reservoir, containing mTeSR1 media with PVA (4 mg/ml; for EB perfusion) or CMM (for cardiac tissue perfusion) with 1% antibiotic-antimycotic (Thermo Fisher Scientific). Next, the tissue was maintained at 37C, where the sacrificial gelatin ink is molten, and a peristaltic pump (Ismatec Inc.) was used to introduce fluid media that removes the gelatin and leaves behind open channels. For the cardiac tissue perfusion, after evacuating the printed gelatin at a vascular flow rate of 0.01 ml/min for 5 min, the inlet for the external perfusion was opened, and the molten gelatin that surrounds the cardiac tissue was evacuated at a high perfusion rate (0.8 ml/min) for 5 min. To ensure the removal of all gelatin, this process was repeated but with a vascular perfusion rate of 0.02 ml/min. Over the course of ~1 hour, the perfusion rate was then increased gradually to its final perfusion rate. Perfusion rates used are 250 l/min for Fig. 3 (C and D), 100 l/min for EB-tissue perfusion for viability measurements, and 500 l/min for cardiac tissue perfusion. For specific experiments, the media reservoir was hyperoxygenated by means of bubbling a mixture of sterile-filtered 95% O2/5% CO2 (Airgas Inc.) through the media at a volumetric flow rate of ~10 ml/min.

Tissues were cryosectioned by first fixing for 30 min in 4% formaldehyde and then rinsed three times in PBS containing 0.05% Tween 20 detergent (PBST). The fixed tissue was incubated overnight at 4C in PBS containing 30% sucrose (sucrose solution) and then transferred into a 1:1 solution of Optimal Cutting Temperature compound:sucrose (OCT:sucrose) for 1 hour and 30 min. Next, the tissue was gently placed into a cryostat tissue mold, which was subsequently filled with 100% OCT solution and frozen on a cryostat Peltier cooler. The tissue was sectioned using 40- to 60-m slices and transferred onto a Superfrost Plus slide (VWR Inc.). Sections were stored at 20C before immunostaining.

For immunostaining, the tissue cryosections were first permeabilized in a 0.125% solution of Triton X containing 0.5% bovine serum albumin (BSA) or 2% donkey serum (DS) for 10 min, rinsed three times in PBST, and blocked for at least 30 min in a solution of PBS containing 3% BSA or 2% DS. Next, primary antibodies [Oct 4 (1:200; sc-5279, Santa Cruz Biotechnology), Sox 2 (1:40; AF2018, R&D Systems), Tuj1 (1:1000; MAB1195, R&D Systems), collagen I (1:500; ab90395, Abcam), cTnT (1:200; ab45932, Abcam), and sarcomeric -actinin (1:200; ab9465, Abcam)] were added in 3% BSA or 2% DS in PBS and incubated overnight at 4C. On the following day, the primary antibody was removed by rinsing three times in PBST, and Alexa Fluor Plusconjugated secondary antibodies were added at 1:400 in 3% BSA or 2% DS in PBS for 45 min. The secondary antibodies were removed by rinsing three times in PBST. The nuclei were then labeled using a 1:5000 dilution of 42,6-diamidino-2-phenylindole in PBS for 3 min, followed by rinsing twice with PBST. The sections were imaged using a Zeiss confocal microscope.

To carry out flow cytometry, after measuring its volume (as described above for EB density estimation), the cardiac OBB matrix was rinsed thoroughly with cold DMEM to remove ECM. The remaining cardiac OBBs were dissociated by incubation for 30 min at 37C in a papain solution [papain (25 U/ml) in Earls balanced salt solution containing l-cystein and EDTA, per the manufacturers instructions (Worthington, Inc.)]. The cells were rinsed in fresh CDM, triturated via repeated pipetting with a P1000 pipette tip, and counted using a hemocytometer. The cells were then rinsed twice in PBS containing 3% BSA, fixed by incubating with CytoFix (BD Biosciences) for 15 min at 4C and then rinsed twice via centrifugation and addition of PBS containing 3% BSA. For flow cytometry, 1 106 cells were resuspended in perm-wash solution (BD Biosciences) and incubated for 15 min at room temperature before rinsing in fresh perm-wash solution. Next, samples were blocked for 15 min at room temperature by incubation in Pharmingen solution (BD Bioscience). The cells were then incubated for 1 hour in labeled primary antibodies (Vimentin PE, 562337, BD Pharmingen; cTnT Alexa Fluor 647, 565744, BD Pharmingen; 1:10 dilution in Pharmingen solution). After incubation, cells were washed twice in perm-wash and analyzed using an LSR II Analyzer (BD Biosciences). Gating of cell populations was performed as shown in fig. S10.

Acknowledgments: We thank J. Weaver for 3D printing cardiac molds and perfusion hardware, D. Mau for assistance in cell culture, S. Cotreau and the Harvard SEAS Machine Shop for the manufacturing of printing platforms, D. Foresti for constructing the temperature control system used for gelatin cooling, Srikanth Damaraju for consultations on coronary artery morphology, T. Ferrante for his assistance with calcium imaging, and L. K. Sanders for photography of printed tissue constructs. Funding: This research was primarily funded by the Office of Naval Research Vannevar Bush Faculty Fellowship program under award number N000141612823. Additional support was provided by the NSF CELL-MET (award number EEC-1647837), the National Human Genome Research Institute of the NIH under award number RM1HG008525 (cerebral organoids), GETTYLAB, and the 3D Organ Engineering Initiative at the Wyss Institute. Author contributions: M.A.S.-S., S.G.M.U., and J.A.L. designed research; M.A.S.-S., S.G.M.U., L.L.N., J.H.A., R.L.T., and S.D. performed research; M.A.S.-S. and S.G.M.U. analyzed data; and M.A.S.-S, S.G.M.U., and J.A.L wrote the manuscript. Competing interests: J.A.L., M.A.S.-S., and S.G.M.U. are inventors on a patent application related to this work filed by Harvard University (no. PCT/US18/51901, filed on 20 September 2018). The authors declare that they have no other competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.


Page 2

September 2019
Vol 5, Issue 9

    • By John Muthii Muriuki, Alexander J. Mentzer, Gavin Band, James J. Gilchrist, Tommy Carstensen, Swaib A. Lule, Morgan M. Goheen, Fatou Joof, Wandia Kimita, Reagan Mogire, Clare L. Cutland, Amidou Diarra, Anna Rautanen, Cristina Pomilla, Deepti Gurdasani, Kirk Rockett, Neema Mturi, Francis M. Ndungu, J. Anthony G. Scott, Sodiomon B. Sirima, Alireza Morovat, Andrew M. Prentice, Shabir A. Madhi, Emily L. Webb, Alison M. Elliott, Philip Bejon, Manjinder S. Sandhu, Adrian V. S. Hill, Dominic P. Kwiatkowski, Thomas N. Williams, Carla Cerami, Sarah H. Atkinson

    • By Daniel Schofield, Arsha Nagrani, Andrew Zisserman, Misato Hayashi, Tetsuro Matsuzawa, Dora Biro, Susana Carvalho

    • By Xuejian Wu, Zachary Pagel, Bola S. Malek, Timothy H. Nguyen, Fei Zi, Daniel S. Scheirer, Holger Müller

    • By Jiayang Li, Kejian Shi, Zeinab Farhadi Sabet, Wenjiao Fu, Huige Zhou, Shaoxin Xu, Tao Liu, Min You, Mingjing Cao, Mengzhen Xu, Xuejing Cui, Bin Hu, Ying Liu, Chunying Chen

    • By Suchan Chang, Dan Hyo Kim, Eun Young Jang, Seong Shoon Yoon, Young Seob Gwak, Yoo Jung Yi, Jun Yeon Lee, Song Hee Ahn, Jin Mook Kim, Yeon-Hee Ryu, Seung-Nam Kim, Hyo Sun Roh, Mi-Young Lee, Sang Chan Kim, Bong Hyo Lee, Hee Young Kim, Chae Ha Yang

    • By Folco Giomi, Alberto Barausse, Carlos M. Duarte, Jenny Booth, Susana Agusti, Vincent Saderne, Andrea Anton, Daniele Daffonchio, Marco Fusi

    • By Rebecca C. Ahrens-Nicklas, Christopher T. Pappas, Gerrie P. Farman, Rachel M. Mayfield, Tania M. Larrinaga, Livija Medne, Alyssa Ritter, Ian D. Krantz, Chaya Murali, Kimberly Y. Lin, Justin H. Berger, Sabrina W. Yum, Chrystalle Katte Carreon, Carol C. Gregorio

    • By Mark A. Skylar-Scott, Sebastien G. M. Uzel, Lucy L. Nam, John H. Ahrens, Ryan L. Truby, Sarita Damaraju, Jennifer A. Lewis

    • By Praneeth R. Kuninty, Ruchi Bansal, Susanna W. L. De Geus, Deby F. Mardhian, Jonas Schnittert, Joop van Baarlen, Gert Storm, Maarten F. Bijlsma, Hanneke W. van Laarhoven, Josbert M. Metselaar, Peter J. K. Kuppen, Alexander L. Vahrmeijer, Arne Östman, Cornelis F. M. Sier, Jai Prakash

    • By Manisha Pandey, Ainslie Calcutt, Victoria Ozberk, Zhenjun Chen, Matthew Croxen, Jessica Powell, Emma Langshaw, Jamie-Lee Mills, Freda E.-C. Jen, James McCluskey, Jenny Robson, Gregory J. Tyrrell, Michael F. Good

    • By Alexessander Couto Alves, N. Maneka G. De Silva, Ville Karhunen, Ulla Sovio, Shikta Das, H. Rob Taal, Nicole M. Warrington, Alexandra M. Lewin, Marika Kaakinen, Diana L. Cousminer, Elisabeth Thiering, Nicholas J. Timpson, Tom A. Bond, Estelle Lowry, Christopher D. Brown, Xavier Estivill, Virpi Lindi, Jonathan P. Bradfield, Frank Geller, Doug Speed, Lachlan J. M. Coin, Marie Loh, Sheila J. Barton, Lawrence J. Beilin, Hans Bisgaard, Klaus Bønnelykke, Rohia Alili, Ida J. Hatoum, Katharina Schramm, Rufus Cartwright, Marie-Aline Charles, Vincenzo Salerno, Karine Clément, Annique A. J. Claringbould, BIOS Consortium, Cornelia M. van Duijn, Elena Moltchanova, Johan G. Eriksson, Cathy Elks, Bjarke Feenstra, Claudia Flexeder, Stephen Franks, Timothy M. Frayling, Rachel M. Freathy, Paul Elliott, Elisabeth Widén, Hakon Hakonarson, Andrew T. Hattersley, Alina Rodriguez, Marco Banterle, Joachim Heinrich, Barbara Heude, John W. Holloway, Albert Hofman, Elina Hyppönen, Hazel Inskip, Lee M. Kaplan, Asa K. Hedman, Esa Läärä, Holger Prokisch, Harald Grallert, Timo A. Lakka, Debbie A. Lawlor, Mads Melbye, Tarunveer S. Ahluwalia, Marcella Marinelli, Iona Y. Millwood, Lyle J. Palmer, Craig E. Pennell, John R. Perry, Susan M. Ring, Markku J. Savolainen, Fernando Rivadeneira, Marie Standl, Jordi Sunyer, Carla M. T. Tiesler, Andre G. Uitterlinden, William Schierding, Justin M. O’Sullivan, Inga Prokopenko, Karl-Heinz Herzig, George Davey Smith, Paul O'Reilly, Janine F. Felix, Jessica L. Buxton, Alexandra I. F. Blakemore, Ken K. Ong, Vincent W. V. Jaddoe, Struan F. A. Grant, Sylvain Sebert, Mark I. McCarthy, Marjo-Riitta Järvelin, Early Growth Genetics (EGG) Consortium

    • By Yuqin Wang, Yu Wang, Xiaoyu Du, Shuanghong Yan, Panke Zhang, Hong-Yuan Chen, Shuo Huang

    • By A. Raveane, S. Aneli, F. Montinaro, G. Athanasiadis, S. Barlera, G. Birolo, G. Boncoraglio, A. M. Di Blasio, C. Di Gaetano, L. Pagani, S. Parolo, P. Paschou, A. Piazza, G. Stamatoyannopoulos, A. Angius, N. Brucato, F. Cucca, G. Hellenthal, A. Mulas, M. Peyret-Guzzon, M. Zoledziewska, A. Baali, C. Bycroft, M. Cherkaoui, J. Chiaroni, J. Di Cristofaro, C. Dina, J. M. Dugoujon, P. Galan, J. Giemza, T. Kivisild, S. Mazieres, M. Melhaoui, M. Metspalu, S. Myers, L. Pereira, F. X. Ricaut, F. Brisighelli, I. Cardinali, V. Grugni, H. Lancioni, V. L. Pascali, A. Torroni, O. Semino, G. Matullo, A. Achilli, A. Olivieri, C. Capelli

    • By Alexander T. Baker, Rosie M. Mundy, James A. Davies, Pierre J. Rizkallah, Alan L. Parker

    • By Jan Felix, Katharina Weinhäupl, Christophe Chipot, François Dehez, Audrey Hessel, Diego F. Gauto, Cecile Morlot, Olga Abian, Irina Gutsche, Adrian Velazquez-Campoy, Paul Schanda, Hugo Fraga

    • By E. Andrew Bennett, Isabelle Crevecoeur, Bence Viola, Anatoly P. Derevianko, Michael V. Shunkov, Thierry Grange, Bruno Maureille, Eva-Maria Geigl

    • By Chaojie Zhang, Chen-Kang Huang, Ken A. Marsh, Chris E. Clayton, Warren B. Mori, Chan Joshi

    • By J.A. Berenbeim, S. Boldissar, S. Owens, M.R. Haggmark, G. Gate, F.M. Siouri, T. Cohen, M. F. Rode, C. Schmidt Patterson, M.S. de Vries

    • By L. Daly, M. R. Lee, S. Piazolo, S. Griffin, M. Bazargan, F. Campanale, P. Chung, B. E. Cohen, A. E. Pickersgill, L. J. Hallis, P. W. Trimby, R. Baumgartner, L. V. Forman, G. K. Benedix

    • By Jennifer M. Ruddock, Haiwang Yong, Brian Stankus, Wenpeng Du, Nathan Goff, Yu Chang, Asami Odate, Andrés Moreno Carrascosa, Darren Bellshaw, Nikola Zotev, Mengning Liang, Sergio Carbajo, Jason Koglin, Joseph S. Robinson, Sébastien Boutet, Adam Kirrander, Michael P. Minitti, Peter M. Weber

    • By Roman Schuetz, Janille M. Maragh, James C. Weaver, Ira Rabin, Admir Masic

    • By Mandana T. Manzari, Grace R. Anderson, Kevin H. Lin, Ryan S. Soderquist, Merve Çakir, Mitchell Zhang, Chandler E. Moore, Rachel N. Skelton, Maréva Fèvre, Xinghai Li, Joseph J. Bellucci, Suzanne E. Wardell, Simone A. Costa, Kris C. Wood, Ashutosh Chilkoti


Page 3

A withdrawal-associated impairment in β-endorphin neurotransmission in the arcuate nucleus (ARC) of the hypothalamus is associated with alcohol dependence characterized by a chronic relapsing disorder. Although acupuncture activates β-endorphin neurons in the ARC projecting to the nucleus accumbens (NAc), a role for ARC β-endorphin neurons in alcohol dependence and acupuncture effects has not been examined. Here, we show that acupuncture at Shenmen (HT7) points attenuates behavioral manifestation of alcohol dependence by activating endorphinergic input to the NAc from the ARC. Acupuncture attenuated ethanol withdrawal tremor, anxiety-like behaviors, and ethanol self-administration in ethanol-dependent rats, which are mimicked by local injection of β-endorphin into the NAc. Acupuncture also reversed the decreased β-endorphin levels in the NAc and a reduction of neuronal activity in the ARC during ethanol withdrawal. These results suggest that acupuncture may provide a novel, potential treatment strategy for alcohol use disorder by direct activation of the brain pathway.

Alcohol use disorder (AUD) is characterized by chronic relapse after periods of abstinence and remains one of the world’s substantial public health problems. The most important method for successful treatment of AUD is to prevent the high rate of craving and relapse. Unfortunately, there is no satisfactory medical intervention to prevent alcohol relapse. Chronic ethanol consumption induces a reduction in β-endorphin neuron activity in the mediobasal hypothalamus, as indicated by the low levels of proopiomelanocortin messenger RNA in ethanol-dependent rats (1). The reduced β-endorphinergic activity in the hypothalamus may be responsible for the dysphoric and depressive states associated with ethanol withdrawal and may lead to the continued ethanol use (2). In addition, it has been proposed that β-endorphin in the nucleus accumbens (NAc) can modulate the behavioral response to stress (3) that can contribute to the development of alcohol dependence (4). Withdrawal from chronic ethanol administration causes the reduction of extracellular dopamine (DA) levels in the NAc, which is considered one of the critical underlying causes of a negative affect state and physical withdrawal signs associated with ethanol withdrawal (5, 6) and provokes relapse to alcohol drinking behavior (7). There is considerable evidence supporting the interaction between the endogenous opioid system and dopaminergic neurotransmission in the NAc (8). Thus, it is expected that β-endorphin may play an important role in attenuating ethanol dependence.

Acupuncture stimulates β-endorphinergic and enkephalinergic fibers in the arcuate nucleus (ARC) of the hypothalamus (9), and endorphinergic neurons originating from the ARC can, in turn, activate opioid receptors expressed on γ-aminobutyric acid (GABA) neurons in the NAc (10)—structures that are involved in ethanol reinforcement. It has been demonstrated that the endogenous opioid system is implicated in the acupuncture’s ability in suppressing the reinforcing effects of abused drugs (11). This is supported by studies showing that systemic administration of the nonselective opioid receptor antagonist naloxone inhibits ethanol intake in ethanol-dependent rats (12).

We have previously observed that acupuncture at Shenmen (HT7) point, located on the ulnar aspect of the wrist, may play an important role in modulating mesolimbic DA release and ethanol self-administration through activation of endorphinergic input to the NAc (13) and in reducing behavioral withdrawal signs by restoring DA levels in the NAc (14). As β-endorphin enhances DA release in the NAc (15), it was hypothesized that acupuncture might stimulate hypothalamic endorphinergic fibers and normalize the decreased β-endorphin levels in the NAc during withdrawal from chronic ethanol exposure, thus preventing physical withdrawal signs and negative emotional state induced by alcohol withdrawal by restoring DA levels in the NAc. To test this hypothesis, we evaluated (i) the effect of acupuncture on physical and psychological signs of ethanol withdrawal in physically ethanol-dependent rats and (ii) the roles of the endogenous opioid system in the NAc in acupuncture and ethanol withdrawal effects.

There were no significant differences in body weights between the control group and ethanol groups during the consumption period of total 16 days (fig. S1A). This suggests that their dietary intake of calories was not affected by the liquid ethanol diet. The average of ethanol consumption during the 9-day period was 10.3 ± 0.2 g kg−1 day−1 and 10.3 ± 0.4 g kg−1 day−1 in the Con-ethanol and HT7-ethanol group, respectively (fig. S1B). The mean blood ethanol concentrations (BECs) for the control group, ethanol group, and ethanol + HT7 group were 5.7 ± 1.3 mg/dl, 198.9 ± 9.0 mg/dl, and 180.7 ± 9.1 mg/dl, respectively. BEC was higher in the ethanol groups than in the control group (F2,18 = 122.976, P < 0.001; P < 0.001 versus control group; fig. S1C). These BECs were sufficiently high to produce physical and psychological withdrawal signs of dependence, as described by others (16). Moreover, rats receiving acupuncture at HT7 had BEC that did not differ significantly from the control rats, suggesting that acupuncture did not affect their ethanol metabolism. In the present study, it was shown that small amount of ethanol has been measured in the samples taken from the control rats. The commercial assay we used to measure the ethanol level actually measures the level of hydrogen peroxide, which is the product of ethanol oxidation. The assay includes alcohol oxidase, which oxidizes ethanol into hydrogen peroxide and measures the plasma level of hydrogen peroxide in terms of the ethanol level. Because small amount of hydrogen peroxide is present in every organism as a result of metabolism, the measurement is due to hydrogen peroxide naturally present in the rats and it is not the product of ethanol oxidation. Therefore, because hydrogen peroxide is the one that reacts with the assay probe, samples of the control rats show fluorescence.

Tremor was quantified in a real-time manner using automatic tremor activity monitoring system (ATAMS; fig. S2A). Power spectra for tremor were estimated by testing the effects of harmaline (10 mg/kg), the tremorgenic compound (17), on the distribution of frequency power generated from ATAMS force transducer in normal rats. Harmaline induced a motion power peak in the 10- to 22-Hz band (P < 0.05 versus pre-harmaline; fig. S2, B and C). Similar to harmaline data, the motion power between 10 and 22 Hz was greater in the ethanol diet rats than in the control diet rats during the withdrawal period (treatment F1,126 = 145.037, P < 0.001; time F20,126 = 8.826, P < 0.001; interaction F20,126 = 5.339, P < 0.001; fig. S2, D and E). On the basis of the results from validation of ATAMS using harmaline, ethanol withdrawal tremor was measured in motion power between 10 and 22 Hz for 5-min segment at each time point.

Tremor behaviors were monitored through ATAMS up to 24 hours after ethanol withdrawal. The ethanol-withdrawn rats exhibited significant increases in tremor ranging from 172 to 270% above baseline tremor (0 hour) between 2 and 5 hours after ethanol withdrawal (group F1,36 = 6.697, P = 0.061; time F9,36 = 2.489, P = 0.025; interaction F9,36 = 4.459, P < 0.001; P < 0.05 versus control group; Fig. 1, A and B). On the basis of these data, we chose the 2-hour time point after withdrawal to evaluate the effects of acupuncture on ethanol withdrawal tremor. The ethanol diet rats receiving HT7 acupuncture, but not those receiving LI5 acupuncture, showed a significant decrease in tremor compared with their pair-fed controls (F4,30 = 6.418, P < 0.001; P < 0.05, Con-ethanol group versus HT7-ethanol group; Fig. 1D). However, animals in the Con-HT7 group had tremor that did not differ from their control rats (P = 0.997, Con-control group versus HT7-control group).

Fig. 1 Effects of acupuncture on ethanol withdrawal tremor.

Rats on the chronic ethanol diet (ethanol group) received an ethanol-containing diet, and their pair-fed rats (control group) received an isocaloric diet with dextran substituted for ethanol. Rats received acupuncture at 2 hours after withdrawal from chronic ethanol consumption. Tremor was measured after acupuncture. (A and B) Representative force signals (filter, 10 to 22 Hz) at 0- and 2-hour time points in the control (n = 5) or ethanol (n = 5) group (A) and tremor activity expressed as a percentage of baseline tremor power at different time points (B). (C) Mechanical acupuncture instrument (MAI). (D) Effect of acupuncture on ethanol withdrawal tremor. Tremor was measured in rats given acupuncture at 2 hours after withdrawal from chronic consumption of an ethanol liquid diet. Another group of control diet rats (Con-control group) and ethanol diet rats (Con-ethanol group) did not receive any needle insertion. Con-control, n = 8; HT7-control, n = 7; Con-ethanol, n = 7; HT7-ethanol, n = 7; LI5-ethanol, n = 6. Graphs represent mean ± SEM. #P < 0.05 versus Con-control group, *P < 0.05 versus Con-ethanol group. Photo credit: Suchan Chang, College of Korean Medicine, Daegu Haany University.

Rats given yohimbine showed significantly lower percentage of time spent in the open arms than control rats given saline (P < 0.01 versus saline; Fig. 4A). Thus, the elevated plus maze apparatus can be used as an index of anxiety. While animals in the ethanol group had significantly higher percentage of time spent in open arm compared to animals in the control group (P < 0.001 versus Con-control or HT7-control group; Fig. 4B), rats (HT7-ethanol group) exhibited a marked elevation of percentage of time spent in the open arms following acupuncture (diet F1,31 = 21.362, P < 0.001; treatment F1,31 = 4.555, P = 0.041; interaction F1,31 = 3.554, P = 0.069; P < 0.05 versus Con-ethanol group; Fig. 4B). On the other hand, acupuncture at HT7 did not alter the time spent in open arm in the control diet rats, suggesting that acupuncture at HT7 may reduce anxiety-like behaviors without affecting the motor function. Similar to HT7 acupuncture, intra-NAc infusions of β-endorphin markedly increased the percentage of time spent in the open arms in ethanol diet–fed rats (F2,15 = 5.434, P < 0.05; Fig. 4C) or yohimbine-treated rats (fig. S4A). Intra-NAc infusions of β-endorphin itself did not affect anxiety-like behavior in control diet–fed rats (fig. S6). In addition, neither acupuncture (F3,31 = 2.049, P = 0.127; fig.S5, A and B) nor β-endorphin (F2,15 = 1.919, P = 0.181; fig.S5, C and D) affected the total distance of movement in the maze (fig. S5).

Fig. 4 Effects of acupuncture or intra-NAc infusions of β-endorphin on anxiety-like behavior in the elevated plus maze in ethanol-dependent rats.

(A) Effect of yohimbine on the time spent in the open arms. The percentage of time spent in the open arms was measured in ethanol-naïve rats 30 min after injection of yohimbine (5 mg/kg, intraperitoneally) (##P < 0.01 versus saline, n = 5 per group). (B) Effect of acupuncture on the time spent in the open arms. The percentage of time spent in the open arms was measured in rats given acupuncture at 2 hours after withdrawal from chronic consumption of an ethanol liquid diet. ###P < 0.001 versus Con-control group (n = 9) or HT7-control group (n = 9), *P < 0.05 versus Con-ethanol group (n = 9); HT7-ethanol group (n = 8). (C) Effect of intra-NAc infusions of β-endorphin on the time spent in the open arms. The percentage of time spent in the open arms was measured in rats given intra-NAc infusions of β-endorphin at 2 hours after withdrawal from chronic consumption of an ethanol liquid diet. #P < 0.05 Con-control group, *P < 0.05 versus Con-ethanol group, n = 6 per group. Graphs represent mean ± SEM of the percentage of time spent in the open arms.

In another set of animals, to see whether acupuncture effects are associated with hypothalamic-pituitary-adrenal axis, we performed an additional experiment concerning the effect of acupuncture on plasma corticosterone levels in yohimbine (an anxiogenic drug)–treated rats. The rats were randomized into control (saline, n = 8), yohimbine (yohimbine only, n = 7), and acupuncture (yohimbine + acupuncture, n = 8) groups. Acupuncture was given at HT7 acupuncture points 30 min after an intraperitoneal injection of yohimbine. Plasma levels of corticosterone were determined before (baseline) an intraperitoneal injection of yohimbine 10 and 40 min after acupuncture treatment. Yohimbine increased plasma corticosterone levels to about 240 ng/ml 40 min after injection, which was significantly suppressed by acupuncture at HT7 (group factor F2,12 = 48.73, P < 0.001; time factor F2,12 = 49.77, P < 0.001; interaction F4,24 = 25.55, P < 0.001; fig. S4B).

The numbers of active lever presses were almost two times higher in the ethanol diet rats than in the control diet rats. There was a significant difference in ethanol-reinforced responding between the ethanol diet rats and the control diet rats (Fig. 5B). Acupuncture significantly reduced ethanol-paired lever responding in the ethanol diet rats (diet F1,25 = 2.740, P < 0.110; treatment F1,25 = 14.124, P < 0.001; interaction F1,25 = 5.928, P = 0.022; P < 0.05 versus Con-ethanol group; Fig. 5B). In addition, a decreased number of inactive lever responses were observed following acupuncture at HT7 in the control or ethanol diet rats (diet F1,25 = 1.005, P = 0.326; treatment F1,25 = 11.172, P = 0.003; interaction F1,25 = 0.208, P = 0.653; Fig. 5C). Animals receiving intra-NAc infusions of β-endorphin showed a significant decrease in ethanol self-administration rate compared to animals given saline (F2,18 = 10.316, P < 0.001; P < 0.05 versus Con-ethanol group; Fig. 5D). The inactive lever responding was not altered in animals receiving intra-NAc infusions of β-endorphin (Fig. 5E).

Fig. 5 Effects of acupuncture or intra-NAc infusions of β-endorphin on ethanol self-administration in ethanol-dependent rats.

(A) Experimental procedures. (B and C) Effect of acupuncture on ethanol self-administration in dependent rats. Ethanol self-administration was measured in rats given acupuncture at 2 hours after withdrawal from chronic consumption of an ethanol liquid diet. #P < 0.05 versus Con-control group; *P < 0.05 versus Con-ethanol group, n = 6 per group. (D and E) Effect of intra-NAc infusions of β-endorphin on ethanol self-administration in ethanol-dependent rats. Ethanol self-administration was measured in rats given intra-NAc infusions of β-endorphin at 2 hours after withdrawal from chronic consumption of an ethanol liquid diet. #P < 0.05 versus Con-control group; *P < 0.05 versus Con-ethanol group; Con-control group, n = 6; Con-ethanol group, n = 8; β-ED-ethanol group, n = 7. Graphs represent mean ± SEM of lever presses.

To see whether acupuncture or β-endorphin itself influences general consummatory behaviors, we have performed further experiments with water self-administration using the same operant ethanol self-administration. Results showed that neither acupuncture nor intra-NAc infusions of β-endorphin affect the amount of water consumption in water-deprived rats. There was no significant difference (P = 0.671) in the baseline number of water infusions between the control group (n = 6) and the HT7 group (n = 6). Acupuncture at HT7 did not alter the number of water infusions compared to the control group (P = 0.501). In addition, there was no difference (P = 0.492) between control (n = 6) and β-endorphin (n = 6) in the baseline number of water infusions. Compared to saline, intra-NAc infusions of β-endorphin slightly, but not significantly, decreased the number of water infusions (P = 0.579) (fig.S7).

Our present data showed that acupuncture at HT7, but not at LI5, significantly decreased ethanol withdrawal tremor. This agrees with previous reports that acupuncture can attenuate AUD-associated withdrawal syndrome in the ethanol-withdrawn animals and alcoholics (14, 18). As insertion of needles into acupuncture points could cause motor impairment and produce generalized effects on tremor activity, the present study attempted to control for this possibility by comparing HT7 acupuncture with LI5 acupuncture. HT7 acupuncture reduced tremor activity in ethanol-dependent rats, whereas it had no effect on tremor activity in the control diet–fed rats. Similar trends were demonstrated for anxiety-like behavior, ethanol self-administration, and a decrease in extracellular β-endorphin levels in the NAc. This suggests that there is an interaction between HT7 acupuncture and ethanol withdrawal responses. In the present study, systemic administration of naloxone, a nonselective opioid receptor antagonist, blocked HT7 suppression of ethanol withdrawal tremor. Furthermore, HT7 acupuncture prevented a decrease in extracellular β-endorphin levels in the NAc in the ethanol-dependent rats. These findings imply that HT7 stimulation produces accumbal β-endorphin–mediated inhibition of ethanol withdrawal tremor. Our results showed that HT7 stimulation inhibited a reduction of neuron firing rate in the ARC of the hypothalamus in the ethanol-withdrawn rats. Similarly, HT7 acupuncture caused neuronal activation in the ARC projecting to the NAc as demonstrated by a marked increase in c-Fos expression in the ARC. These results provide strong evidence that acupuncture enhances opioid-mediated transmission in the ARC, which projects to the NAc. As β-endorphin–containing fibers highly exist in the ARC (19), innervate the NAc (10), and mediate acupuncture inhibition of morphine withdrawal syndromes (11), we speculated that β-endorphinergic neurons projecting to the NAc from the ARC might contribute to the ability of acupuncture to reduce ethanol dependence. Accordingly, a local injection of β-endorphin into the NAc mimicked acupuncture-mediated suppression of ethanol withdrawal symptoms including tremor, anxiety-like behavior, and ethanol self-administration. These results support the acupuncture’s ability in suppressing ethanol withdrawal symptoms by activating β-endorphin–mediated transmission in the ARC projecting to the NAc. We have previously proposed a possible mechanism that might explain how acupuncture at HT7 might activate β-endorphin neurons in the ARC of the hypothalamus (13). HT7 acupuncture may activate peripheral mechanoreceptors, and the signals are conveyed through large A-fibers within the ulnar nerve, which lies adjacent to HT7 points (20) and the dorsal column somatosensory pathway (21). Sensory stimulation induced by acupuncture triggers β-endorphin neurons in the ARC (9) projecting to the NAc, thereby activating μ-opioid receptors on accumbal GABA neurons (10).

The present study indicated that acupuncture attenuated the decrease in open-arm exploration of the elevated plus maze in ethanol-dependent rats, suggesting an anxiolytic-like function for acupuncture. In addition, acupuncture reduced ethanol self-administration during ethanol withdrawal. Results showed that neither acupuncture nor β-endorphin altered water self-administration and locomotor activity in the elevated plus maze, implying that suppression of ethanol-reinforced responding and anxiety-like behavior by acupuncture and β-endorphin may not be due to reduction in general consummatory behavior and motor impairment. Together with the previous observation that anxiety-like state causes relapse to alcohol drinking in dependent rats (16), the present results suggest that the decreased alcohol drinking behavior induced by acupuncture might result from a diminished anxiety-like state. In addition, although the precise relationship between physical withdrawal syndrome and alcohol drinking behavior remains to be elucidated, ethanol withdrawal tremor has been suggested to be a powerful stressor (22). This supports the idea that tremor may play an important role as a driving factor for ethanol relapse (22). Together, results suggest that acupuncture may attenuate relapse to ethanol seeking in dependent rats by reducing tremor and anxiety-like state. In addition, a local injection of β-endorphin into the NAc decreased anxiety-like behavior and acupuncture reduced plasma corticosterone levels in yohimbine-treated rats, respectively. Yohimbine has been suggested to increase stress through noradrenergic control of the hypothalamo-pituitary adrenal system (23). Thus, this implies the possibility that the reduced stress response produced by acupuncture is mediated via an inhibition of the hypothalamo-pituitary adrenal system in response to β-endorphin in the NAc. Several studies have suggested that β-endorphin in the NAc would play an important role in mediating an adaptive response to stress with DA-independent mechanism. For example, it has been shown that extracellular levels of β-endorphin in the NAc increased in response to the stressful component of extinction of heroin-reinforced behavior and an aversive stimulus (3). A similar result was obtained in another study in which stress facilitates aversion to novelty feeding in mice with lowered β-endorphin (24). The exact neurochemical mechanism mediating a role for β-endorphin is unknown, but some evidence suggests that β-endorphin attenuates stress response by suppressing the secretion of corticotrophin-releasing factor (CRF) (25). Thus, it can be suggested that β-endorphin in the NAc may reduce relapse to ethanol seeking via an inhibition of stress response. Similarly, reduction in β-endorphinergic activity in the brain after chronic exposure to ethanol was suggested to be responsible for the dysphoric symptoms and anxiety state of ethanol withdrawal in alcoholics that might contribute to negative reinforcement (2, 8). In addition, the reduction in hypothalamic β-endorphinergic neuronal activity during ethanol withdrawal might mediate physical withdrawal syndrome and anxiety-like behavior (26). Together, these findings support the hypothesis that β-endorphin neurons in the ARC may mediate the acupuncture’s ability in suppressing ethanol withdrawal symptoms and relapse to alcohol drinking.

Our previous study has shown that HT7 acupuncture prevented a decrease of DA release in the NAc and physical signs of withdrawal during ethanol withdrawal (14). This study raised the possibility that acupuncture may play a functional role in reducing ethanol dependence by normalizing the release of DA in the NAc. Withdrawal from chronic ethanol administration causes a reduction of DA transmission in the mesolimbic pathway, suggesting a likely neurochemical mechanism of the negative affective state and somatic withdrawal signs during ethanol abstinence (5, 6), the intense ethanol craving experienced by addicts (27), and ethanol-seeking behavior in ethanol-dependent rats (7). As the negative affective state could motivate continued ethanol-seeking behavior in ethanol-dependent rats (7), the reduced ethanol self-administration produced by acupuncture or infusions of β-endorphin into the NAc is likely mediated via an inhibition of DA depletion in the NAc, as revealed by the parallel decreases of anxiety-like behavior. Our results showed that intra-NAc infusions of β-endorphin restored deficient TH expression in the VTA and increased TH phosphorylation in the VTA during ethanol withdrawal. This suggests that infusions of β-endorphin into the NAc induce the reversal of the hypodopaminergic state during ethanol withdrawal. Conditions that induce synaptic activation of the cell appear to facilitate the expression and phosphorylation of TH in peripheral sympathetic neurons and central catecholaminergic neurons (28). As accumbal GABA input to VTA DA neurons is modulated by opioid receptors (29), we speculate that β-endorphin may enhance DA transmission in the mesolimbic pathway. Thus, the increased VTA TH expression and phosphorylation in response to β-endorphin is likely due to an enhancement of VTA DA neuron activity. Together with the present results that local injection of β-endorphin into the NAc mimics acupuncture effects on ethanol withdrawal symptoms, it is likely that HT7 acupuncture reduces the negative affective state and alcohol drinking behavior in ethanol-dependent rats via activation of DA neurons in the VTA, which is caused by activation of endorphinergic input into the NAc from the ARC.

The present results demonstrated a significant increase in tremor, anxiety-like behavior, and ethanol self-administration during withdrawal from chronic ethanol. In addition, significant β-endorphin reduction in the NAc was observed in the ethanol-dependent rats. Previous studies have shown that the hypothalamic β-endorphinergic activity was decreased in the ethanol-dependent animals (1, 26). These findings raise the possibility that the diminished β-endorphin release in the NAc may be attributed to suppression of endorphinergic input to the NAc during withdrawal from chronic ethanol. Similarly, it has been shown that NAc β-endorphin levels in the ethanol-intoxicated state were significantly decreased compared with water-intoxicated state (30). As hypothalamic endorphinergic neurons innervate the NAc (10), the decreased β-endorphinergic neuron activity during ethanol withdrawal would result in the excitation of NAc GABA neurons in the NAc. In support of this hypothesis, previous studies using whole-cell recording and in vivo microdialysis revealed that NAc GABA neurons become hyperexcitable during withdrawal from repeated cocaine exposure (31). Although how β-endorphin modulates DA levels in the NAc still remains elusive, it has been shown that morphine inhibits GABA neurons in the NAc and consequently disinhibits DA neurons in the VTA (29). In addition, NAc GABA neurons synapse into DA neurons in the VTA (32, 33). Thus, we can speculate that β-endorphin exerts its inhibitory action on GABA neurons in the NAc projecting to VTA DA neurons via opioid receptors. Accordingly, it can be expected that the lower level of NAc β-endorphin, endogenous opioid peptides, during ethanol withdrawal than normal states may disinhibit GABA neurons in the NAc and subsequently cause decrement of DA neuronal activity in the VTA. A marked deficiency of DA release in the NAc and reduced firing rates of DA neurons in the VTA was shown in the ethanol-dependent rats, and reversals of the hypodopaminergic state in the NAc during ethanol withdrawal attenuated affective and somatic withdrawal signs (5, 6). On the basis of interactions of the endogenous opioid system and DAergic neurotransmission, inhibition of β-endorphin neurons by chronic exposure to ethanol might lead to reduction of DAergic activity in the mesolimbic DA system. Thus, the present results suggest that a significant decrease of accumbal β-endorphinergic activity may represent the underlying mechanism of, at least in part, the dysphoric and neurological symptoms of ethanol abstinence. It has been proposed that the dysphoric and depressive state during ethanol withdrawal may result from the reduced β-endorphinergic activity in the hypothalamus (2). While our and other animal studies revealed consistent results that acupuncture suppresses alcohol-related behaviors such as alcohol craving or consumption (13, 34), human studies have shown mixed results in the effects of acupuncture on reducing clinical symptoms related to alcohol addiction. For examples, controlled studies showed that severe recidivist alcoholics treated with acupuncture at points specific for the treatment of substance abuse expressed less need for alcohol and less drinking episodes, compared to the control group given acupuncture at nonspecific acupoints (3537). However, Worner et al. (38) reported negative results of acupuncture treatment on alcohol-related symptoms. This discrepancy may be due to many vagaries including acupoint location, treatment duration, and stimulation frequency and intensity, leading to low reproducibility and high individual variations among acupuncturists in human studies. In support of this, our previous study demonstrated that considerable variations in intensity and frequency of acupuncture needle stimulation are observed within each acupuncturist or across acupuncturists (20).

In conclusion, we found that acupuncture inhibited tremor, anxiety-like behavior, and alcohol drinking behavior by restoring decreased β-endorphin levels in the NAc in physically ethanol-dependent rats. However, determination of the specific mechanisms involved in the neuronal links between acupuncture-induced sensory stimulation and hypothalamic β-endorphin neurons will require future study. Last, acupuncture may be an effective therapeutic intervention to treat AUD by direct activation of brain pathways.

The following drugs were used: β-endorphin [a μ-opioid receptor agonist; 0.25 and 0.5 μg per side, dissolved in sterile saline (site), Bachem, Switzerland]; yohimbine (an anxiogenic compound; 5 mg/kg, dissolved in distilled water, intraperitoneal injection; Sigma-Aldrich, USA); naloxone (a nonselective opioid receptor antagonist; 1.0 mg/kg, dissolved in sterile saline, intraperitoneal injection; Sigma-Aldrich, USA); harmaline (a tremogenic compound; 10 mg/kg, dissolved in sterile saline, subcutaneous injection; Sigma-Aldrich, USA); and a retrograde tracer, 1,1′-dioctadecyl-3,3,3′,3′-tetramethylindocarbocyanine (DiI) (2%, dissolved in methanol, Life Technologies, USA).

Rats were housed individually and given ad libitum access to a commercially available liquid diet (Dyets, Bethlehem, PA, USA), known as the Lieber-DeCarli diet. Rats in the chronic ethanol group (ethanol group) were exposed to 1 to 7% ethanol by gradually increasing ethanol concentration over the 7-day period and then maintained at a concentration of 7.2% for 9 days. Rats in the control group (control group) received ethanol-free isocaloric diets (or control diets) in which maltose dextrin replaced ethanol and the same amount of the standard liquid diet as their ethanol-paired rats consumed during the previous day.

Blood ethanol levels were measured 2 hours after removal of the ethanol diet on the last day of the ethanol liquid diet (16 days). To test the possibility that acupuncture might lower the blood ethanol level directly and attenuate alcohol dependence, we investigated whether acupuncture affected serum ethanol concentration in ethanol-dependent rats. Under isoflurane anesthesia, blood samples (0.5 ml) were taken by cardiac puncture 5 min after acupuncture stimulation. Serum ethanol levels were assayed using a commercially available colorimetric/fluorometric assay kit (BioVision, Milpitas, CA, USA).

This experiment examined the effects of acupuncture on ethanol withdrawal tremor in ethanol-dependent rats. Acupuncture was given for 20 s after withdrawals from chronic ethanol consumption. Ethanol withdrawal was precipitated by replacement of the ethanol-containing diet with the control liquid diet for 2 hours. To pharmacologically characterize the effects of acupuncture on ethanol withdrawal tremor, we gave rats an intraperitoneal injection of the nonselective opioid receptor antagonist naloxone (1.0 mg/kg) 30 min before acupuncture treatment. Rats also received intra-NAc infusions of β-endorphin (0.25 and 0.5 μg per side) to determine whether β-endorphin in the NAc can mimic acupuncture effect. As GABA neurons to the VTA projecting from the NAc were more abundant and more sensitive to μ-opioid receptor agonist in the shell than core (39), intra-infusions of β-endorphin were administered into the shell. Tremor was quantified in a real-time manner after acupuncture treatment by using a custom-made ATAMS (Fig. 2A), as previously described (40). In brief, the system was designed to monitor tremor activity through a force transducer (Grass Instruments, Braintree, MA, USA) mounted under a clear plastic cylinder (20 cm × 6 cm). The degree of tremor was measured for 15 min after placing each rat in the plastic cylinder. The second 5-min recording was performed for data. The signals from the force transducer were fed into bridged amplifiers (ETH-200, CB Sciences Inc., Dover, NH, USA), filtered between 10 and 22 Hz, and quantified using a LabChart and Scope program (ADInstruments). Power spectra for tremor were estimated by testing a tremogenic compound, harmaline. Harmaline, a derivative of β-carboline, is one of the most frequently used tremorgenic compound, and harmaline-induced tremor is regarded as a model for essential tremor in animals. It has been shown that subcutaneous or intraperitoneal injection of harmaline induces tremors within a few minutes in forelimb, hindlimb, head, or whole body, and the tremors last up to several hours (41). To find out specific oscillation frequencies to differentiate tremor-like activity from other movements, we assessed the motion powers during a 5-min pretreatment baseline and a 5-min posttreatment 10 min after subcutaneous injection of harmaline (10 mg/kg) in ethanol-naïve rats. On the basis of the results from validation of ATAMS using harmaline, ethanol withdrawal tremor was measured between 10 and 22 Hz. Tremor was expressed using the mean of motion power taken from a 5-min-long graphic recording.

Microdialysis was conducted as previously described (14). One week before the start of the ethanol diet, a single guide cannula was stereotaxically implanted into the NAc shell [anteroposterior (AP), 1.7; mediolateral (ML), 0.8; and dorsoventral (DV), 6.0 from bregma] under anesthesia (42). On the day of testing, a microdialysis probe (2 mm length, 20 kDa cutoff value, polyarylethersulfone membrane, CMA 12 Elite, Carnegie Medicine, Stockholm, Sweden) was inserted into the guide cannula in the brain of unanesthetized rats before starting the microdialysis session. Artificial cerebrospinal fluid (150 mM NaCl, 3.0 mM KCl, 1.4 mM CaCl2, and 0.8 mM MgCl2 in 10 mM phosphate buffer at pH 7.4) was perfused at a rate of 2.0 μl/min using a microinfusion pump (CMA 100, Carnegie Medicine, Stockholm, Sweden) 2 hours before ethanol withdrawal. The microdialysis samples were collected at 30-min intervals from 2 hours before and up to 2 hours after the start of ethanol abstinence and stored frozen at −80°C until assayed. Acupuncture was given 2 hours after the start of ethanol withdrawal. Levels of β-endorphin were measured using a commercially available enzyme-linked immunosorbent assay (ELISA) kit (Peninsula Laboratories Inc., Belmont, CA, USA) as described by others (3). Data in the microdialysates were reported from β-endorphin levels within the linear segment of the standard curve (0.05 to 5 ng/ml). The intra-assay coefficient was 4%, and the inter-assay coefficient was 13%. Baseline levels of β-endorphin were defined as the average value in four samples collected before acupuncture treatment.

To identify whether ARC neurons projecting to the NAc were activated by HT7 acupuncture, we assessed the expression of c-Fos in ARC neurons. DiI (2%) (a fluorescent retrograde neuronal tracer; Life Technologies, MA, USA) in 0.1 μl per side over 5 min was infused bilaterally into the NAc (AP, 1.7; ML, ±0.75; and DV, 5.7 from bregma) using a 10-μl Hamilton syringe under intraperitoneal sodium pentobarbital anesthesia (50 mg/kg). The tracer was allowed to spread for at least 7 days. After acupuncture treatment, rats were sacrificed and brains were prepared for c-Fos immunohistochemistry. Preparation of the brain tissue began with transcardial perfusion of 4% paraformaldehyde in 0.1 M phosphate buffer under pentobarbital anesthesia (80 mg/kg) 1 hour after acupuncture stimulation. Subsequently, brains were removed, immersed in 4% buffered paraformaldehyde for 2 hours, allowed to stand in 30% sucrose overnight, and cryosectioned at 30 μm (between 3.3 and 3.6 mm posterior to bregma). The brain slices were then mounted on gelatin-coated glass slides. The ARC labeled by DiI was identified with an Olympus AX70 fluorescence microscope (Olympus, Tokyo, Japan). For c-Fos immunohistochemistry study, the brain slices were incubated in blocking solution containing 3% bovine serum albumin for 1 hour at room temperature and incubated with anti–c-Fos rabbit polyclonal antibodies (1:1000, Santa Cruz Biotechnology, Santa Cruz, USA) for 24 hours at 4°C, followed by incubation with a biotinylated donkey anti-rabbit Alexa Fluor 488 (1:500, Invitrogen, Grand Island, NY, USA). Last, the sections were mounted onto gelatin-coated slides. The numbers of c-Fos immunoreactive neurons were counted under 200× with a confocal laser scanning microscope (LSM700, Carl Zeiss, Oberkochen, Germany) and averaged from three to four brain slices per rat.

This experiment was carried out to investigate the effects of acupuncture on neuronal firing rates in the ARC in ethanol-dependent rats. Two hours after withdrawals from the chronic ethanol diet, rats were placed in a stereotaxic apparatus under pentobarbital anesthesia (50 mg/kg, intraperitoneally) for in vivo extracellular recordings of ARC neurons. Anesthesia level was maintained by intravenous infusion of sodium pentobarbital (15 mg kg−1 hour−1) via an indwelling catheter into the jugular vein during the recording. The responsiveness of ARC neurons (AP, −3.3 to 3.6; ML, 2.0 to 2.2; and DV, 10 to 11.1 from bregma) was recorded by in vivo extracellular recording technique as previously described (43). ARC neurons were triggered by acupuncture stimulation. A single microcarbon filament-filled glass electrode (carbon fiber electrode; impedence, 0.4 to 1.2 megohms; Kation Scientific, Minneapolis, MN, USA) was aimed at the ARC (AP, −3.3 to −3.6; ML, 2.0 to 2.2; and DV, 10 to 11.1 from bregma) at an angle of 10° with a one-axis water hydraulic micromanipulator (Narishige, Amityville, NY, USA) following needle insertion into HT7 point. After identification of a single ARC neuron, single-cell activity was amplified (voltage gain, 104), filtered at 0.3 to 10 kHz with an ISO-80 amplifier (World Precision Instruments, Sarasota, FL, USA), fed directly into the data acquisition unit (CED Micro1401; Cambridge Electronic Design, UK), and stored on computer to construct the waveforms or plot peri-stimulus time histograms (spikes per second bin width, Spike2 software). To ensure that a single ARC neuron was held for the duration of recording, we used the Spike2 program to confirm the same action potential shape and amplitude. Single-cell activity was expressed using the average value in three samples obtained from a 10-s-long graphic recording during acupuncture at an interval of 20 s between samples. Baseline single-cell activity was recorded for a 20-s-long graphic recording at rest following stabilization of single-unit firing rate for at least 10 min.

The experimental procedure for operant ethanol self-administration in ethanol-dependent rats was based on a slight modification of methods used in earlier studies (7). Ethanol self-administration was carried out in daily 30-min sessions for 5 days a week during the dark cycle. After rats met criteria for stable responding with less than 20% variation in three consecutive sessions, rats were given 16 days of daily access to the ethanol liquid diet. Two hours after ethanol withdrawal on day 17, rats were given the opportunity to self-administer ethanol for 30 min after acupuncture treatment or intra-NAc infusions of β-endorphin. After completion of this test session, each rat was exposed again to the ethanol liquid diet for 4 to 5 days for the subsequent sessions identical to the first test. Each experiment was conducted using a within-subjects Latin square design, as described previously (7) (Fig. 5A). Operant test chambers (Med Associates Inc., Georgia, VT, USA) were equipped with two response levers and the house light that was illuminated during each self-administration session. Each press on the active lever delivered 0.1 ml of 10% (v/v) ethanol solution for 5 s to one of two drinking cups in the center panel of the operant chamber. Responses on the inactive lever had no programmed consequence but were recorded. During the lever training with sucrose solution, rats were subjected to a 22-hour water restriction schedule, with freely available food, in their home cage for successful acquisition of lever responding. The rats were trained in daily 30-min session to press a lever for a 20% (w/v) sucrose solution on a fixed ratio (FR1) reinforcement schedule until an acquisition criteria of 100 lever presses for three consecutive sessions was achieved. The rats were then trained to orally self-administer ethanol using a modified sucrose-fading method, as described previously (13). Following lever press training, rats self-administered 20% sucrose solution for the next 3 days. After establishment of stable lever responses, concentration of sucrose was gradually lowered to 0% and that of ethanol was raised to 10%. The detailed procedure was as follows: three to four sessions with 10% sucrose, three sessions with 2% ethanol in 10% sucrose, three sessions with 5% ethanol in 10% sucrose, and four sessions with 10% ethanol in 5% sucrose. Last, 10% ethanol without sucrose was introduced as the reinforcer. The rats met an established criterion for stable active lever responding with less than 20% variation in three consecutive responses. Typically, this required approximately 21 days following initiation of 10% ethanol self-administration. Different groups of animals were also subjected to water self-administration in an attempt to control for the effects of acupuncture or β-endorphin on general consummatory behaviors. Water self-administration took place using the same operant ethanol self-administration. Responses on the active lever produce a 0.1-ml drop of water, whereas responses on the inactive lever were recorded but had no consequence. Rats were trained to self-administer water for 30 min to facilitate lever pressing on an FR1 schedule at least three consecutive days. During the course of the entire sessions, rats were water-deprived for 16 hours. After rats exhibited stable water infusions in three consecutive responses with less than 20% variation, they were exposed to acupuncture or intra-NAc infusions of β-endorphin.

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/9/eaax1342/DC1

Fig. S1. Effects of chronic ethanol consumption on growth rate and BEC.

Fig. S2. Measurement of ethanol withdrawal tremor.

Fig. S3. Schematic localization of microdialysis probes and infusion sites.

Fig. S4. Effects of intra-NAc infusions of β-endorphin on anxiety-like behavior in the elevated plus maze and acupuncture at HT7 on plasma corticosterone levels in yohimbine-treated rats.

Fig. S5. Effects of HT7 acupuncture and intra-NAc infusions of β-endorphin on locomotor activity in the elevated plus maze.

Fig. S6. Effects of intra-NAc infusions of β-endorphin on anxiety-like behavior in the elevated plus maze in control diet–fed rats.

Fig. S7. Effects of HT7 acupuncture and intra-NAc infusions of β-endorphin on water self-administration.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank J. B. Kang for his helpful advice and B. Park for her English grammar correction. Funding: This work was supported by National Research Foundation of Korea (NRF) grant funded by the Korea government MSIT (no. 2018R1A5A2025272). Author contributions: H.Y.K. and C.H.Y. designed the experiments. S.C., D.H.K., E.Y.J., S.S.Y., Y.S.G., Y.J.Y., J.Y.L., S.H.A., J.M.K., B.H.L., Y.-H.R., S.-N.K., H.S.R., M.-Y.L., and S.C.K. carried out the experiments and the data analysis. S.C., H.Y.K., and C.H.Y. drafted the manuscript. H.Y.K. and C.H.Y. were responsible for the overall direction of the project and for edits to the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data need to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.


Page 4

September 2019
Vol 5, Issue 9

    • By John Muthii Muriuki, Alexander J. Mentzer, Gavin Band, James J. Gilchrist, Tommy Carstensen, Swaib A. Lule, Morgan M. Goheen, Fatou Joof, Wandia Kimita, Reagan Mogire, Clare L. Cutland, Amidou Diarra, Anna Rautanen, Cristina Pomilla, Deepti Gurdasani, Kirk Rockett, Neema Mturi, Francis M. Ndungu, J. Anthony G. Scott, Sodiomon B. Sirima, Alireza Morovat, Andrew M. Prentice, Shabir A. Madhi, Emily L. Webb, Alison M. Elliott, Philip Bejon, Manjinder S. Sandhu, Adrian V. S. Hill, Dominic P. Kwiatkowski, Thomas N. Williams, Carla Cerami, Sarah H. Atkinson

    • By Daniel Schofield, Arsha Nagrani, Andrew Zisserman, Misato Hayashi, Tetsuro Matsuzawa, Dora Biro, Susana Carvalho

    • By Xuejian Wu, Zachary Pagel, Bola S. Malek, Timothy H. Nguyen, Fei Zi, Daniel S. Scheirer, Holger Müller

    • By Jiayang Li, Kejian Shi, Zeinab Farhadi Sabet, Wenjiao Fu, Huige Zhou, Shaoxin Xu, Tao Liu, Min You, Mingjing Cao, Mengzhen Xu, Xuejing Cui, Bin Hu, Ying Liu, Chunying Chen

    • By Suchan Chang, Dan Hyo Kim, Eun Young Jang, Seong Shoon Yoon, Young Seob Gwak, Yoo Jung Yi, Jun Yeon Lee, Song Hee Ahn, Jin Mook Kim, Yeon-Hee Ryu, Seung-Nam Kim, Hyo Sun Roh, Mi-Young Lee, Sang Chan Kim, Bong Hyo Lee, Hee Young Kim, Chae Ha Yang

    • By Folco Giomi, Alberto Barausse, Carlos M. Duarte, Jenny Booth, Susana Agusti, Vincent Saderne, Andrea Anton, Daniele Daffonchio, Marco Fusi

    • By Rebecca C. Ahrens-Nicklas, Christopher T. Pappas, Gerrie P. Farman, Rachel M. Mayfield, Tania M. Larrinaga, Livija Medne, Alyssa Ritter, Ian D. Krantz, Chaya Murali, Kimberly Y. Lin, Justin H. Berger, Sabrina W. Yum, Chrystalle Katte Carreon, Carol C. Gregorio

    • By Mark A. Skylar-Scott, Sebastien G. M. Uzel, Lucy L. Nam, John H. Ahrens, Ryan L. Truby, Sarita Damaraju, Jennifer A. Lewis

    • By Praneeth R. Kuninty, Ruchi Bansal, Susanna W. L. De Geus, Deby F. Mardhian, Jonas Schnittert, Joop van Baarlen, Gert Storm, Maarten F. Bijlsma, Hanneke W. van Laarhoven, Josbert M. Metselaar, Peter J. K. Kuppen, Alexander L. Vahrmeijer, Arne Östman, Cornelis F. M. Sier, Jai Prakash

    • By Manisha Pandey, Ainslie Calcutt, Victoria Ozberk, Zhenjun Chen, Matthew Croxen, Jessica Powell, Emma Langshaw, Jamie-Lee Mills, Freda E.-C. Jen, James McCluskey, Jenny Robson, Gregory J. Tyrrell, Michael F. Good

    • By Alexessander Couto Alves, N. Maneka G. De Silva, Ville Karhunen, Ulla Sovio, Shikta Das, H. Rob Taal, Nicole M. Warrington, Alexandra M. Lewin, Marika Kaakinen, Diana L. Cousminer, Elisabeth Thiering, Nicholas J. Timpson, Tom A. Bond, Estelle Lowry, Christopher D. Brown, Xavier Estivill, Virpi Lindi, Jonathan P. Bradfield, Frank Geller, Doug Speed, Lachlan J. M. Coin, Marie Loh, Sheila J. Barton, Lawrence J. Beilin, Hans Bisgaard, Klaus Bønnelykke, Rohia Alili, Ida J. Hatoum, Katharina Schramm, Rufus Cartwright, Marie-Aline Charles, Vincenzo Salerno, Karine Clément, Annique A. J. Claringbould, BIOS Consortium, Cornelia M. van Duijn, Elena Moltchanova, Johan G. Eriksson, Cathy Elks, Bjarke Feenstra, Claudia Flexeder, Stephen Franks, Timothy M. Frayling, Rachel M. Freathy, Paul Elliott, Elisabeth Widén, Hakon Hakonarson, Andrew T. Hattersley, Alina Rodriguez, Marco Banterle, Joachim Heinrich, Barbara Heude, John W. Holloway, Albert Hofman, Elina Hyppönen, Hazel Inskip, Lee M. Kaplan, Asa K. Hedman, Esa Läärä, Holger Prokisch, Harald Grallert, Timo A. Lakka, Debbie A. Lawlor, Mads Melbye, Tarunveer S. Ahluwalia, Marcella Marinelli, Iona Y. Millwood, Lyle J. Palmer, Craig E. Pennell, John R. Perry, Susan M. Ring, Markku J. Savolainen, Fernando Rivadeneira, Marie Standl, Jordi Sunyer, Carla M. T. Tiesler, Andre G. Uitterlinden, William Schierding, Justin M. O’Sullivan, Inga Prokopenko, Karl-Heinz Herzig, George Davey Smith, Paul O'Reilly, Janine F. Felix, Jessica L. Buxton, Alexandra I. F. Blakemore, Ken K. Ong, Vincent W. V. Jaddoe, Struan F. A. Grant, Sylvain Sebert, Mark I. McCarthy, Marjo-Riitta Järvelin, Early Growth Genetics (EGG) Consortium

    • By Yuqin Wang, Yu Wang, Xiaoyu Du, Shuanghong Yan, Panke Zhang, Hong-Yuan Chen, Shuo Huang

    • By A. Raveane, S. Aneli, F. Montinaro, G. Athanasiadis, S. Barlera, G. Birolo, G. Boncoraglio, A. M. Di Blasio, C. Di Gaetano, L. Pagani, S. Parolo, P. Paschou, A. Piazza, G. Stamatoyannopoulos, A. Angius, N. Brucato, F. Cucca, G. Hellenthal, A. Mulas, M. Peyret-Guzzon, M. Zoledziewska, A. Baali, C. Bycroft, M. Cherkaoui, J. Chiaroni, J. Di Cristofaro, C. Dina, J. M. Dugoujon, P. Galan, J. Giemza, T. Kivisild, S. Mazieres, M. Melhaoui, M. Metspalu, S. Myers, L. Pereira, F. X. Ricaut, F. Brisighelli, I. Cardinali, V. Grugni, H. Lancioni, V. L. Pascali, A. Torroni, O. Semino, G. Matullo, A. Achilli, A. Olivieri, C. Capelli

    • By Alexander T. Baker, Rosie M. Mundy, James A. Davies, Pierre J. Rizkallah, Alan L. Parker

    • By Jan Felix, Katharina Weinhäupl, Christophe Chipot, François Dehez, Audrey Hessel, Diego F. Gauto, Cecile Morlot, Olga Abian, Irina Gutsche, Adrian Velazquez-Campoy, Paul Schanda, Hugo Fraga

    • By E. Andrew Bennett, Isabelle Crevecoeur, Bence Viola, Anatoly P. Derevianko, Michael V. Shunkov, Thierry Grange, Bruno Maureille, Eva-Maria Geigl

    • By Chaojie Zhang, Chen-Kang Huang, Ken A. Marsh, Chris E. Clayton, Warren B. Mori, Chan Joshi

    • By J.A. Berenbeim, S. Boldissar, S. Owens, M.R. Haggmark, G. Gate, F.M. Siouri, T. Cohen, M. F. Rode, C. Schmidt Patterson, M.S. de Vries

    • By L. Daly, M. R. Lee, S. Piazolo, S. Griffin, M. Bazargan, F. Campanale, P. Chung, B. E. Cohen, A. E. Pickersgill, L. J. Hallis, P. W. Trimby, R. Baumgartner, L. V. Forman, G. K. Benedix

    • By Jennifer M. Ruddock, Haiwang Yong, Brian Stankus, Wenpeng Du, Nathan Goff, Yu Chang, Asami Odate, Andrés Moreno Carrascosa, Darren Bellshaw, Nikola Zotev, Mengning Liang, Sergio Carbajo, Jason Koglin, Joseph S. Robinson, Sébastien Boutet, Adam Kirrander, Michael P. Minitti, Peter M. Weber

    • By Roman Schuetz, Janille M. Maragh, James C. Weaver, Ira Rabin, Admir Masic

    • By Mandana T. Manzari, Grace R. Anderson, Kevin H. Lin, Ryan S. Soderquist, Merve Çakir, Mitchell Zhang, Chandler E. Moore, Rachel N. Skelton, Maréva Fèvre, Xinghai Li, Joseph J. Bellucci, Suzanne E. Wardell, Simone A. Costa, Kris C. Wood, Ashutosh Chilkoti


Page 5

Since the early 1990s, the Greenland ice sheet (GrIS) has been losing mass at an accelerating rate, primarily due to enhanced meltwater runoff following atmospheric warming. Here, we show that a pronounced latitudinal contrast exists in the GrIS response to recent warming. The ablation area in north Greenland expanded by 46%, almost twice as much as in the south (+25%), significantly increasing the relative contribution of the north to total GrIS mass loss. This latitudinal contrast originates from a different response to the recent change in large-scale Arctic summertime atmospheric circulation, promoting southwesterly advection of warm air toward the GrIS. In the southwest, persistent high atmospheric pressure reduced cloudiness, increasing runoff through enhanced absorption of solar radiation; in contrast, increased early-summer cloudiness in north Greenland enhanced atmospheric warming through decreased longwave heat loss. This triggered a rapid snowline retreat, causing early bare ice exposure, amplifying northern runoff.

Since the 1990s, mass loss from the Greenland ice sheet (GrIS) has significantly accelerated (1) due to increased glacial discharge (2) and decreased surface mass balance (SMB), the latter representing the difference between mass gains from snowfall and mass losses from meltwater runoff and sublimation. While glacial discharge strongly reacts to variations in ocean temperatures (3, 4) and can regionally dominate mass loss in marine-terminating sectors of the ice sheet, i.e., particularly in northwest (NW) and southeast (SE) Greenland (5, 6), the recent (1991–2015) GrIS-wide mass loss is primarily (~61%) ascribed to a decrease in SMB resulting from increased meltwater runoff (7), concurrent with the atmospheric warming that followed a recent summertime circulation change (8). Historically, the wide ablation zone in southwest (SW) Greenland contributed most (~32%) to the ice sheet runoff total (table S1), mainly driven by absorbed solar radiation at the dark bare ice surface in summer (911). However, especially over highly reflective snow surfaces, clouds can also enhance melt and runoff through reduced surface heat loss by longwave radiation (12, 13). While the recent summertime circulation shift is responsible for reduced cloud cover over SW Greenland resulting in enhanced melt (14), to date, the physical processes responsible for the runoff changes in the north remain unclear.

To address this, here, we combine SMB of a dedicated, high-resolution (5.5 km) run from the Regional Atmospheric Climate Model (RACMO2.3p2; fig. S1A), statistically downscaled to 1 km (Fig. 1A), with bare ice extent (BIE) derived from remotely sensed broadband shortwave clear sky albedo from the moderate resolution imaging spectroradiometer (MODIS) MCD43A3 500-m daily albedo product. This high resolution is essential to resolve the narrow ablation zone and marginal glacier tongues in sufficient detail (15), e.g., applying the statistical downscaling technique increases GrIS runoff by 25% compared to the original model resolution of 5.5 km (from 238 to 297 Gt year−1 for 1958–2017). Simulated GrIS climate and SMB components are evaluated using in situ SMB measurements collected in the accumulation (SMB > 0; 182 sites in fig. S1B) and ablation zones (SMB < 0; 213 sites in figs. S1C and S4A); meteorological observations recorded at 37 automatic weather stations (AWSs) from the Programme for Monitoring of the Greenland Ice Sheet (PROMICE), Institute for Marine and Atmospheric Research Utrecht (IMAU), and Greenland Climate Network (GC-Net) networks (figs. S2 and S3); and meltwater discharge measured at the Watson River in west Greenland (fig. S4B). The results show that RACMO2.3p2 realistically represents near-surface temperature, specific humidity, wind speed and air pressure (0.73 < R2 < 0.98), cloud conditions through shortwave and longwave radiation components (0.85 < R2 < 0.96), SMB (R2 = 0.78), and Watson River meltwater discharge (R2 = 0.91). We distinguish seven Greenland sectors (Fig. 1A)—NW, NO (north), NE (northeast), CE (center east), SE, SW, and CW (center west)—based on their distinct climatologies while covering similar areas (16). Not considering the very rugged SE and CE sectors, where bare ice area variability is well captured, although the absolute extent is underestimated, there is very good agreement between modeled and MODIS remotely sensed bare ice area (R2 = 0.82; fig. S5D). The model evaluation is discussed in more detail in Materials and Methods.

Fig. 1 High-resolution runoff patterns and post-1990 changes.

(A) Annual mean runoff (RU) over the GrIS and peripheral glaciers and ice caps (1958–2017) as modeled by RACMO2.3p2 and statistically downscaled to 1-km horizontal resolution. Yellow stars locate 213 stake sites where SMB is measured. The seven selected sectors of the GrIS derived from (16) are overlaid. (B) Time series of annual cumulative runoff integrated over the GrIS (red). Panels (C), (D), and (E) show the same for the SW, NW, and NO sectors, respectively, of the GrIS [see (A)]. Dashed red lines show linear runoff trends over 1991–2017, and dashed black lines show the averaged runoff over the periods 1958–1990 and 1991–2017, respectively. Dashed yellow lines mark exceptional runoff years, i.e., 3 SDs (σ) above the 1958–1990 mean. Dashed gray lines highlight recent exceptional runoff years over the GrIS, i.e., 2010, 2012, and 2016 (see fig. S7). Individual trends (1991–2017), relative increase in runoff post-1990 (%), and SD (σ) for the period 1958–1990 are listed at the bottom of each subpanel.

Figure 1A shows annual mean runoff for the period 1958–2017. Before 1990, the GrIS mass balance was near zero (7, 17) or slightly negative (16), followed by accelerating mass loss. This mass loss is primarily driven by a 42% increase in runoff (1991–2017 versus 1958–1990), showing a positive trend of 6.6 ± 1.9 Gt year−2 (Fig. 1B), as compared to an estimated 4.4 ± 0.4 Gt year−2 trend in solid ice discharge (16). On the basis of a breakpoint analysis, we selected 1991 as the year after which the GrIS mass balance became predominantly negative (7). The frequency of extreme runoff years, i.e., at least 3 SDs (σ) above the 1958–1990 mean (dashed yellow lines in Fig. 1, B to E), increased from the 2000s onward when runoff reached a higher plateau (dashed black line). In GrIS sectors SW, NW, and NO, runoff during 1991–2017 increased by 33, 60, and 70%, respectively, relative to 1958–1990 (Fig. 1, C to E, and table S1). This shows that post-1990 runoff acceleration was relatively more pronounced in north Greenland than in south Greenland.

Figure 2A shows that before 1991, the SW sector contributed 32% to the total runoff, mostly originating from the relatively wide marginal ablation zone exposing dark bare ice in summer (Fig. 1A), whereas the relatively cold and dry NW and NO sectors only contributed 10 and 6%, respectively (Fig. 2A and table S1). After 1991, runoff is significantly enhanced, but mostly so in NW (+60%) and NO (+70%) compared to +33% in SW (Fig. 1, B to E). As a result, the relative runoff contributions from NW and NO significantly increased by 12 and 18% (stippled in Fig. 2B; table S1). Because NW-NO and CW-SW Greenland show similar responses (Fig. 2B), these sectors are merged into two regions named “North” and “South.” Figure 2C highlights how the post-1990 trend in runoff contribution is negative, although insignificant (P = 0.193), in south Greenland (red) and significantly positive (P = 0.001) in north Greenland (blue; table S1).

Fig. 2 Regional changes in runoff contribution.

(A) Average contribution of individual GrIS sectors to runoff totals for the period 1958–1990. (B) Post-1990 change in runoff contribution per sector (1991–2017 minus 1958–1990). Sectors experiencing significant change in runoff contribution, i.e., based on Student’s t test (t ≤ 0.10), are stippled with white dots. The post-1990 relative change in runoff contribution is also listed for each sector in white. (C) Time series of runoff contribution for North (blue; i.e., NO + NW) and South (red; i.e., CW + SW) Greenland. Post-1990 trends in runoff contribution are displayed as dashed lines.

We find that in concert with the runoff increase, the northern Greenland ablation zone expanded by 46% post-1990, almost twice as much as in the south (+25%; Fig. 3, A and B). At the same time, the bare ice zone, where the seasonal winter snow has completely melted at the end of summer, has grown more than twice as fast in the North (+33%) compared to the South (+14%), in good agreement with remotely sensed maximum BIE derived from MODIS (Fig. 3, C and D). Figure S6 shows how the SMB components changed in North and South Greenland for different ice sheet elevations. In the North, the snowline retreat moved the equilibrium line altitude (ELA; SMB = 0) upward by ~200 m from ~900 meters above sea level (masl) (fig. S6A) to ~1100 masl (fig. S6B), increasing runoff by 12 Gt year−1 (fig. S6C). In the South, the ELA migrated by ~100 m from ~1350 masl (fig. S6D) to ~1450 masl (fig. S6E), increasing runoff by ~17 Gt year−1 (red line in fig. S6F). We conclude that the faster snowline retreat in the North is responsible for its enhanced contribution to Greenland runoff totals (Fig. 2B). Note also how a saturated firn percolation zone with net ablation (SMB < 0; difference between the ablation zone and bare ice area in Fig. 3C) progressively forms in the North after the mid-1990s, resembling more and more the southern GrIS ablation zone (Fig. 3D).

Fig. 3 Rapid ablation zone expansion enhances runoff contribution from north Greenland.

(A) Map of SMB averaged for the period 1958–1990. Numbers refer to the ablation zone area for individual sectors (103 km2) and for the whole GrIS (bottom right). (B) Same as (A) for the period 1991–2017. Numbers refer to the relative ablation zone expansion (%) post-1990 for individual sectors and for the whole GrIS (bottom right). Time series of annual mean modeled ablation zone and summer bare ice area for (C) North Greenland (blue and cyan; i.e., NW + NO) and (D) South Greenland (red and orange; i.e., CW + SW sectors) over the period 1958–2017 compared to MODIS (black) bare ice area estimates (2000–2018). Dashed lines show the 1991–2017 trends. Numbers include trends, relative ablation/bare ice zone expansion, and change in ELA (i.e., SMB = 0) post-1990. In North and South Greenland, the modeled bare ice area is averaged over 10 and 5 days, respectively (see Materials and Methods). The cyan and yellow belt in (C) and (D) represents 1 SD of the 10 and 5 days used to estimate the modeled maximum bare ice area in North and South Greenland, respectively. The gray belt in (C) and (D) shows the uncertainty in measured MODIS bare ice area (see Eq. 1 in Materials and Methods).

We relate the latitudinal contrast in runoff increase to the recent change in the large-scale summertime atmospheric circulation. The latter change, observed in the Arctic from the late 1990s onward (8), promoted southwesterly advection of warm air from Baffin Bay toward the GrIS. Figure 4A illustrates the post-1990 change in June-July-August (JJA; 1991–2017 minus 1958–1990) large-scale circulation over Greenland (contours and arrows). A persistent high-pressure ridge results in an anticyclonic circulation anomaly over the GrIS (vectors in Fig. 4A). As a result, the NO and NW sectors experience anomalous westerly advection of warm and humid air from Baffin Bay, increasing summer (JJA) cloud content after 1991 by, respectively, 5 and 10% on average (3 and 8 ± 4 g m−2) and up to 20% locally (29 ± 4 g m−2). In contrast, northerly advection of relatively dry inland air reduces cloud content over the southern GrIS by 2 to 6% on average (−3 to −6 ± 4 g m−2), in line with (14). Clouds play a pivotal role in the GrIS melt climate by modulating downward shortwave (SWd) and longwave (LWd) radiation fluxes (18). Figure 4B shows post-1990 changes in summer SWd (JJA; 1991–2017 minus 1958–1990), mirroring the large-scale changes in clouds: SWd decreases in the NW (−2.3 ± 1.9 W m−2 on average) and NO (−3.2 ± 1.9 W m−2) sectors, while it mostly increases elsewhere (Fig. 4B). In contrast, LWd has increased everywhere by 3.5 W m−2 on average (Fig. 4C) due to overall higher free atmosphere temperatures (+0.9°C at 500 hPa; see inset in Fig. 4C). However, the LWd increase peaks in NW (4.7 ± 1.7 W m−2) and in NO (5.3 ± 1.7 W m−2) sectors, where increased clouds further enhance LWd compared to, e.g., SW Greenland (2.6 ± 1.7 W m−2). Therefore, the near-surface temperature increase is amplified in northern Greenland (e.g., +0.9 ± 0.1°C for NO) relative to other sectors (e.g., +0.5 ± 0.1°C for SW). Note the marked longitudinal contrast in cloud content change between NO-NW and NE sectors, where cloud content is reduced by 2% on average and down to 14% locally (−1 g m−2 down to −10 ± 4 g m−2). This is due to a large-scale Foehn effect, resulting in relatively warm and dry air flowing down the lee side of the main northern ice divide, after having generated enhanced orographic precipitation at the upwind slopes of the NW and NO sectors. The average circulation and radiative flux anomalies shown in Fig. 4 represent large interannual variability (fig. S7) that reflects the exact position of the high-pressure ridge in summer (19). As a result, runoff also shows increased interannual variability after 1991 (Fig. 1, B to E).

Fig. 4 Recent shift in summer atmospheric circulation and impact on cloudiness.

(A) Post-1990 change in summer cloud content (JJA; 1991–2017 minus 1958–1990) as modeled by RACMO2.3p2 at 5.5 km. Changes in large-scale circulation (black vectors; see inset for wind speed estimation) and in height of the 500-hPa geopotential (dashed black lines) are overlaid. (B) Change in modeled SWd and (C) LWd radiation. In (C), sector-averaged near-surface temperature (2 m) increase is displayed in black and white for southern and northern Greenland, respectively. Average GrIS-wide temperature increase at 500 hPa and at 2-m altitude is listed in the inset. The seven GrIS sectors are outlined in gray in (A) to (C). Stipples highlight regions showing a significant change, i.e., > 1 SD of the 1958-1990 period (σ) in (A) cloud content (σ = 4 g m−2), (B) SWd (σ = 1.9 W m−2), and (C) LWd (σ = 1.7 W m−2).

Figure 5 shows how anomalies in climate and runoff are interacting during the course of the melt season (April-September). Figure 5A shows the daily runoff contribution of North and South regions to GrIS total (% per day for the period April-September) averaged for 1958–1990 (cyan and orange lines) and 1991–2017 (blue and red lines). Figure 5B displays the cumulative anomalies (1991–2017 minus 1958–1990) in melt (orange; Gt) and runoff (red) as well as in cloud content (dashed gray; kg m−2) and refreezing capacity (dashed cyan; unitless), i.e., the fraction of meltwater and rainwater that is retained and/or refrozen within the firn pack. Figure S8 shows similar results for the South region. Figure 5 (C and D) shows the anomalies in June net longwave (LWn) and July-August net shortwave radiation (SWn). Figure 5E shows the daily fraction of the ablation zone area exposing bare ice in summer for North and South Greenland for the period before (cyan and orange) and after 1991 (blue and red). Note that, in a warming climate in which the equilibrium line migrates upward, the ablation zone comprises both the marginal bare ice area and the lower percolation zone, where saturated firn contributes to meltwater runoff in late summer.

Fig. 5 Increased early-summer cloudiness triggers runoff amplification in north Greenland.

(A) Time series of daily (April-September) mean runoff contribution to GrIS totals for periods 1958–1990 and 1991–2017 in South (orange and red; i.e., CW + SW sectors) and North (cyan and blue; i.e., NO + NW sectors) regions. (B) Time series of daily (April-September) cumulative anomalies (1991–2017 minus 1958–1990) in surface melt (orange) runoff (red) for the North region (NO + NW sectors; right y axis). Dashed gray and cyan lines (left y axis) show cumulative anomalies in cloud content and refreezing capacity, i.e., the fraction of meltwater and rainwater retained and/or refrozen in the firn. Anomalies in (C) June LWn and (D) July-August SWn. (E) Daily exposed bare ice area post-1990 and pre-1991 for the South (red and orange) and North (blue and cyan) regions expressed as a fraction (%) of the ablation zone area. Numbers at the bottom left corner of (E) express the relative increase (%) in maximum bare ice area post-1990 for North (blue) and South (red) Greenland. In (A), (B), and (E), the gray and yellow shades outline the period during which runoff contribution of North Greenland significantly increases (A) under high and low cloudiness successively (B).

In the South, the decrease in runoff contribution, i.e., the difference between the red and orange solid lines in Fig. 5A, occurs mainly during the peak melt season between July and early August (yellow shade in fig. S8). In that period, the cumulative cloud content anomaly gradually decreases and becomes negative (i.e., less clouds compared to pre-1991) from mid-May and throughout the melt season (fig. S8), enhancing SWd (Fig. 4B). The post-1990 melt and runoff increase, estimated at 46 and 37 Gt, is driven by increasing SWn (Fig. 5D): Earlier removal of the seasonal snow cover exposes dark bare ice over the extensive South ablation zone, enhancing the absorption of incoming solar radiation (10). In contrast, increased runoff contribution in North Greenland after 1991, i.e., the difference between the blue and cyan solid lines in Fig. 5A, starts early in the melt season in June and ends in late July. This increase is concurrent with high cloud content in June (gray shade in Fig. 5B) followed by low cloud content until late July (yellow shade in Fig. 5B). In early summer, increased cloudiness and atmospheric temperatures in North Greenland (Fig. 4, A and C) act in concert to warm the snow-covered surface through reduced longwave heat loss (Fig. 5C). This warming triggers an early melting of the shallow snow cover, rapidly migrating the snowline further inland. As a result, early exposure of dark bare ice causes a rapid decline in surface albedo (fig. S9A) and reduces meltwater retention (fig. S9B). At higher elevations, the presence of summer clouds maintains firn temperatures close to the melting point during nighttime (13), preventing refreezing of percolating meltwater; note the sharp decline in firn refreezing capacity that mirrors the cloud content increase (gray shade; Fig. 5B and fig. S9B). Figure 5E highlights the rapid snowline retreat in the North, resulting in an expansion of the (maximum) daily bare ice area by 34% compared to 20% in the South. Early-summer exposure of bare ice further leads to anomalous high absorption of incoming shortwave radiation for reduced cloud content in July (yellow shade; Fig. 5, B and D). This chain of events significantly amplifies the relative runoff increase in North Greenland (+63%; table S1) compared to South Greenland (+34%) and other sectors (fig. S9C).

The whole of Greenland has warmed since the early 1990s, but mostly in the north, amplifying northern mass loss through enhanced meltwater runoff. Using a combination of remote sensing and output of a regional climate model, we show that this latitudinal gradient can be coupled to a circulation shift that brings more clouds to northern Greenland. This triggered early snowline retreat in the dry north, causing a rapid expansion of the bare ice (+33%) and ablation zones (+46%), twice as fast as in the south (+14 and +25%). As a result, North Greenland has significantly increased its relative contribution to total GrIS mass loss through enhanced runoff. Superimposed on increased melt, rising temperatures and increased cloudiness in the North also led to a summer rainfall increase of ~42%, contributing 5 Gt year−1 or 8% to North Greenland runoff totals (fig. S6, A to C). In late summer (September), rainfall in North Greenland occasionally even equals runoff over bare ice, in line with (20). If the current trends in ablation area expansion were to continue, the northern ablation zone would equal the size of the southern counterpart in another ~45 years. This would require the current circulation anomaly to persist, but predicting this is highly uncertain, as Earth System Models (ESMs) from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) fail to reproduce the contemporary large-scale Arctic circulation change (21) that led to the latitudinal gradient in runoff response reported here. This highlights the importance of better resolving Arctic circulation variability and cloud microphysics in ESMs (12, 13), e.g., cloud phase, water content, and optical thickness, to obtain reliable GrIS mass change projections.

We used a new run at 5.5-km horizontal resolution of the polar (p) version of RACMO2.3p2 for the period 1958–2017. For detailed description of the model and recent updates, we refer to (22). In brief, RACMO2.3p2 incorporates the dynamical core of the High Resolution Limited Area Model (23) and the physics from the European Centre for Medium-Range Weather Forecasts–Integrated Forecast System (ECMWF-IFS cycle CY33r1) (24). RACMO2.3p2 includes a multilayer snow module that simulates melt, water percolation, and retention in snow, refreezing, and runoff (25). The model also accounts for dry snow densification (26), and drifting snow erosion and sublimation (27). Snow albedo is calculated on the basis of snow grain size, cloud optical thickness, solar zenith angle, and impurity concentration in snow (28). Compared to (22), no model physics have been changed. However, increased horizontal resolution of the host model, i.e., 5.5 km instead of 11 km, better resolves gradients in SMB components over the topographically complex ice sheet margins and neighboring peripheral glaciers and ice caps.

We evaluated the modeled GrIS present-day climate by comparing modeled meteorological variables (fig. S2), i.e., 2-m air temperature and specific humidity, 10-m wind speed and surface pressure, and radiative fluxes (fig. S3), i.e., short/longwave down/upward radiation (SWd, LWd, SWu, and LWu, respectively), to daily measurements collected at 37 AWS (green dots in fig. S1A). These AWSs are operated by the PROMICE (18 sites for 2007–2016) (33), by the IMAU (5 sites for 2004–2016) (34), and by the GC-Net (14 sites for 1995–2017) (35). Erroneous radiation measurements, caused by, e.g., sensor riming in winter and sensor heating in summer, were discarded from the analysis following the method described in (22). Sites showing an elevation difference of >100 m compared to the GIMP DEM down-sampled to 5.5-km resolution were not used (nine sites). Because of frequent data gaps and sensor issues, the GC-Net results for air temperature, specific humidity, wind speed, and surface pressure are included in fig. S2 for completeness but were not used for evaluation statistics.

Figures S2 and S3 show that RACMO2.3p2 at 5.5 km agrees well when compared to daily measurements of meteorological variables and radiative fluxes. Notably, 2-m temperature and specific humidity only show small biases of 0.14°C and −0.11 g kg−1, respectively. Modeled and observed radiative fluxes also agree well with small biases, e.g., 4.3 and −7.8 W m−2 for SWd and LWd, respectively. Note that these biases can be partly ascribed to sensor uncertainty, estimated at 2.7 and 1.2% for daily mean SWd (4.8 W m−2) and LWd (2.7 W m−2) radiation (34). Additional uncertainty originates from window heating in summer, i.e., shortwave-induced sensor heating, potentially overestimating LWd records by 10 to 15 W m−2 on a 6-min basis (34) or less than 5 W m−2 on daily average. Similar evaluations in four individual GrIS sectors (22) demonstrated that modeled SWd and LWd were as well represented in northern Greenland as in other regions. The strong correlations with remotely sensed (13) and in situ SWd (fig. S3A) and LWd (fig. S3C) demonstrate that RACMO2.3p2 realistically represents Greenland cloud characteristics. Furthermore, errors in modeled SWd/LWd [root mean square error (RMSE) of ~20 W m−2] were much smaller than the difference between clear sky and overcast conditions, reaching ~100 W m−2 (12). Note that surface pressure shows systematic biases at several stations, resulting from an up to 100-m difference between the model and station elevation.

For SMB evaluation, we used 182 accumulation measurements from stakes, firn pits, and cores (36) derived from airborne radar campaign (37) in the GrIS accumulation zone (white dots in fig. S1A) as well as 1073 ablation measurements from 213 stake sites (yellow dots in fig. S1A) (38). Figure S1B compares modeled (5.5 km in blue and 11 km in red) and observed SMB in the accumulation zone. Compared to the previous simulation at 11-km resolution, the 5.5-km run slightly improves the SMB representation in the accumulation zone, with a smaller bias and RMSE reduced by 4.5 and 2.7 mm we (millimeters water equivalent) or 21 and 4%, respectively. Relative to (22), significant improvements were found in the ablation zone, where RACMO2.3p2 at 5.5 km now explains 48% of the observed SMB variance (R2) instead of 41% previously, with a smaller bias and RMSE reduced by 230 and 140 mm we or 38 and 11%, respectively. This is due to the better resolved SMB patterns over narrow glaciers and marginal ablation zones at a resolution of 5.5 km (fig. S1A) compared to 11 km [figure 1 of (22)].

Although improving on previous model versions, significant SMB biases persist locally in the low ablation zone, where surface ablation exceeds 4 m we year−1 (fig. S1C). For instance, at stake QAS_L (291 masl), i.e., the lowest site of the 30-km-long Qagssimiut transect at the southern tip of the GrIS, strong winter precipitation and summer ablation gradients are not well resolved, resulting in SMB biases of up to 8 m we for summer 2010 (fig. S1C). In these rugged marginal regions, a spatial resolution of 5.5 km remains insufficient to resolve narrow glaciers and confined ablation zones that are prime contributors to GrIS runoff totals (15). To address this, we applied statistical downscaling to the 5.5-km product as described in (15, 39), which corrects runoff for biases in elevation and bare ice albedo. This allows us to reproduce more accurately the high runoff rates observed at the GrIS margins, significantly improving the agreement with SMB measurements (fig. S4A). The downscaled product explains 78% of the SMB variance (R2), with a positive bias of 70 mm we, suggesting a small runoff underestimation. In addition, a comparison between modeled and observed meltwater discharge from the Watson River catchment in SW Greenland (gray contour in sector SW of Fig. 1A) also shows good agreement despite a small negative bias of, on average, ~0.4 Gt year−1 (or km3) (1976–2016) and up to ~2 Gt year−1 for exceptional melt episodes that occurred in summers of 2010 and 2012 (fig. S4B). The average runoff bias reaches 5% for the period 1976–2016, peaking at ~20% for extreme runoff years (2010 and 2012), in line with (40). On the basis of this, we used a conservative 20% runoff uncertainty for annual basin integrated values. Good agreement between modeled and in situ surface ablation (fig. S4A), i.e., primarily driven by runoff, as well as with measured meltwater discharge from the Watson River in west Greenland (fig. S4B), confirms that meltwater runoff is well reproduced by RACMO2.3p2 (0.78 < R2 < 0.91). Last, a recent comparison between the period 2003–2017 mass changes derived from the Gravity Recovery and Climate Experiment (GRACE) and from a combination of modeled RACMO2.3p2 SMB at 1-km resolution and estimated solid ice discharge showed good agreement on both GrIS-wide and individual basin scales (6).

Remotely sensed annual bare ice area was derived from broadband shortwave clear sky albedo data from the MCD43A3 MODIS 500-m daily albedo product (http://dx.doi.org/10.5067/MODIS/MCD43A3.006). To discard invalid albedo records, we masked erratic albedo grid cells from the GrIS-wide, daily 16-day MCD43A3 MODIS product (2000–2018), resorting to full Bidirectional Reflectance Distribution Function (BRDF) inversions. Valid daily MODIS data were classified as ice- or snow-covered grid cells using an upper threshold for shortwave albedo of 0.60 (SW0.60), slightly higher than the assumed maximum albedo of bright bare ice (0.55). The justification for this is to ensure that the maximum bare ice area is fully captured by MODIS. Subsequently, the daily ice/snow cells were converted to annual bare ice area assuming that the bare ice is exposed in summer if (i) the current pixel is classified as ice at least 5 days in that year (5th percentile) and (ii) that pixel is located within the long-term modeled ablation zone of the GrIS (SMB < 0; 1991–2017), i.e., including both the bare ice and lower percolation zone below the long-term ELA of ~1450 masl. These criteria allow the elimination of spurious bare ice conditions, e.g., meltwater lakes, meltwater runoff streams, and superimposed ice, that are detected at higher elevations in the percolation zone of the GrIS. For instance, in summer 2012, MODIS detected bare ice conditions at AWS S9 (1520 masl) located close to the equilibrium line along the K-transect in western Greenland (67°N), whereas fieldwork reported superimposed ice conditions resulting from firn saturation. Masked pixels were then filled on the basis of the recurrence of ice/snow observations over 2000–2018 for that cell. In brief, masked pixels exposing ice more than 50% of the time for the period 2000–2018 are assumed bare ice. This method may result in bare ice area overestimation in early years of MODIS operation, which includes more invalid albedo estimates. Alternatively, daily bare ice conditions could have been extrapolated onto masked pixels based on surrounding grid cells. However, we discarded this approach because (i) local MODIS pixels may not be representative of surrounding areas and (ii) a large number of MODIS records were masked due to, e.g., a large solar zenith angle and local cloud cover, deteriorating the quality of local MODIS estimates. Extrapolation may thus strongly overestimate remotely sensed (observed) bare ice area for the entire 2000–2018 period. That is why we used the ice recurrence method to obtain a spatially continuous annual BIE product covering the whole of the GrIS for the period 2000–2018 (fig. S5A). The MODIS product is sensitive to the shortwave albedo threshold used to discriminate between ice and snow conditions. To reflect this, we estimated an uncertainty in maximum MODIS BIE by repeating the analysis described above using the assumed maximum albedo of bright bare ice of 0.55 (SW0.55) as an upper threshold. The BIE uncertainty was estimated GrIS wide and for individual sectors as followsUncertainty = BIESW0.60BIESW0.55(1)

The latter uncertainty was used to derive a symmetric error band around the MODIS bare ice area estimates in Fig. 3 (C and D) and fig. S5 (B to D).

For model evaluation that is consistent with the MODIS product, we estimated a 10-day bare ice area in northern sectors (NW, NO, and NE) and a 5-day one for southern regions. Daily modeled bare ice area was estimated as the integrated area of pixels showing a surface albedo of ≤0.55 on the original 5.5-km grid, bilinearly interpolated to the 1-km GIMP ice mask. The 10 and 5 largest daily estimates of bare ice area are averaged for North and South Greenland, respectively, to obtain a maximum annual bare ice area comparable to MODIS data. This is deemed reasonable as MODIS records 16-day albedo and because less valid MODIS albedo estimates are available for the North (55% of valid records, on average, for the period 2000–2018) compared to the South (79%). This difference stems from large solar zenith angles at high latitudes, negatively affecting the quality of the satellite measurements. Figure S5A shows minimum and maximum annual BIE derived from MODIS for the period 2000–2018, i.e., summer of 2000 (black) and 2012 (red). Modeled GrIS bare ice area (orange) agrees well with MODIS (black; fig. S5B), although with a negative bias of ~8100 km2 (fig. S5C) peaking for exceptional melt years, e.g., 2012, 2014, and 2016. This negative bias, i.e., up to 16,100 km2 in 2014, mainly originates from the CE (~25%) and SE (~25%) sectors of the GrIS (fig. S5D). One possible explanation is that, in CE and SE sectors, RACMO2.3p2 simulates too much snowfall in winter, delaying or preventing the exposure of bare ice in summer. At the same time, these two sectors also experience frequent overcast conditions, negatively affecting the quality of the MODIS bare ice product. In other sectors, RACMO2.3p2 agrees well with observed bare ice area (R2 = 0.82; fig. S5D) with a small negative bias (140 km2).

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/9/eaaw0123/DC1

Table S1. Changes in runoff production and contribution per sector.

Fig. S1. RACMO2.3p2 integration domain and SMB evaluation.

Fig. S2. Evaluation of modeled meteorological variables at 5.5 km.

Fig. S3. Evaluation of radiative fluxes at 5.5 km.

Fig. S4. Evaluation of the downscaled product using in situ and catchment measurements.

Fig. S5. Evaluation of the modeled bare ice area using remote sensing.

Fig. S6. Post-1990 upward migration of the equilibrium line.

Fig. S7. High interannual variability in summer atmospheric circulation and impacts on the cloudiness.

Fig. S8. Reduced summer cloud cover enhances melt in southern Greenland.

Fig. S9. Post-1990 changes in surface conditions.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: Funding: B.N., W.J.v.d.B., and M.R.v.d.B. acknowledge funding from the Polar Program of the Netherlands Organization for Scientific Research (NWO; grant no. 866.15.202) and the Netherlands Earth System Science Centre (NESSC). Author contributions: B.N. prepared the manuscript, carried out the RACMO2.3p2 simulation, and produced the downscaled dataset at 1 km. B.N., W.J.v.d.B., and M.R.v.d.B. designed the study and analyzed the data. S.L. processed the MODIS albedo data and produced MODIS BIE estimates. All authors commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The annual data sets required to reproduce the figures shown in the manuscript and Supplementary Material are available at https://doi.pangaea.de/10.1594/PANGAEA.904428. The daily downscaled SMB dataset presented in this study is freely available from the authors without conditions upon request. Besides SMB, the dataset includes daily total precipitation (snow and rain), snowfall, total melt (snow and ice), meltwater runoff, retention and refreezing, total sublimation (surface and drifting snow), and snow drift erosion at 1-km horizontal resolution for the period 1958–2017. The dataset also includes daily 2-m air temperature at 1-km resolution. Data from RACMO2.3p2 at 5.5-km spatial resolution and MODIS BIE are also available upon request.
« Go back