We employ three-dimensional, fully resolved numerical simulations using the volume-of-fluid method to study the motion and interaction of two in-line bubbles ascending in both Newtonian and shear-thinning fluids. Additionally, we explore passive scalar transfer between the fluid phases across a variety of fluidic conditions, modeling shear-thinning behavior in non-Newtonian fluids through the Carreau model. The impact of the Galilei (Ga) and Bond (Bo) numbers, the bubble pair radius ratio, the inelastic time constant (λ), and the flow index (n) on the bubbles dynamics and the transient Sherwood number (Sht) and the surface-averaged Sherwood number (⟨Sh⟩) are reported. Using the well-known Ga–Bo regime phase diagram for a single rising bubble in a Newtonian ambient fluid, the present numerical experiments are used to study the departure from this reference case due to the presence and characteristics of a second bubble and the non-Newtonian nature of the ambient fluid. When categorized based on the single bubble phase diagram, we found that in regimes I (axisymmetric) and III (oscillatory), a pair of bubbles does not breakup or merge during our simulations. However, their behaviors vary due to the second bubble and change in non-Newtonian fluid parameters like the inelastic time constant and flow index. Likewise, we explored this parameter space for regime II (skirted), where the two bubbles eventually merge, and regimes IV (peripheral breakup) and V (central breakup), known for multiple bubble breakups. Additionally, we present results on differently sized bubbles, showing that their merging tendency depends on their arrangement as leading or trailing positions in the pair.