To explore the physics of large-scale flows in solar-like stars, we perform 3D anelastic simulations of rotating convection for global models with stratification resembling the solar interior. The numerical method is based on an implicit large-eddy simulation approach designed to capture effects from non-resolved small scales. We obtain two regimes of differential rotation, with equatorial zonal flows accelerated either in the direction of rotation (solar-like) or in the opposite direction (anti-solar). While the models with the solar-like differential rotation tend to produce multiple cells of meridional circulation, the models with anti-solar differential rotation result in only one or two meridional cells. Our simulations indicate that the rotation and large-scale flow patterns critically depend on the ratio between buoyancy and Coriolis forces. By including a subadiabatic layer at the bottom of the domain, corresponding to the stratification of a radiative zone, we reproduce a layer of strong radial shear similar to the solar tachocline. Similarly, enhanced superadiabaticity at the top results in a near-surface shear layer located mainly at lower latitudes. The models reveal a latitudinal entropy gradient localized at the base of the convection zone and in the stable region, which however does not propagate across the convection zone. In consequence, baroclinicity effects remain small and the rotation iso-contours align in cylinders along the rotation axis. Our results confirm the alignment of large convective cells along the rotation axis in the deep convection zone, and suggest that such banana-cell pattern can be hidden beneath the supergranulation layer.