We study the abundance of subhaloes in the hydrodynamical cosmological simulation Illustris, which includes both baryons and dark matter in a $\Lambda$CDM volume 106.5 Mpc a side. We compare Illustris to its dark matter-only (DMO) analogue, Illustris-Dark, and quantify the effects of baryonic processes on the demographics of subhaloes in the host mass range $10^{11}$ to $3 \times 10^{14} \msun$. We focus on both the evolved ($z=0$) subhalo cumulative mass functions (SHMF) and the statistics of subhaloes ever accreted, i.e. infall subhalo mass function. We quantify the variance in subhalo abundance at fixed host mass and investigate the physical reasons responsible for such scatter. We find that in Illustris, baryonic physics impacts both the infall and $z=0$ subhalo abundance by tilting the DMO function and suppressing the abundance of low-mass subhaloes. The breaking of self-similarity in the subhalo abundance at $z=0$ is enhanced by the inclusion of baryonic physics. The non-monotonic alteration of the evolved subhalo abundances can be explained by the modification of the concentration--mass relation of Illustris hosts compared to Illustris-Dark.Interestingly, the baryonic implementation in Illustris does not lead to an increase in the halo-to-halo variation compared to Illustris-Dark. In both cases, the normalized intrinsic scatter today is larger for Milky Way-like haloes than for cluster-sized objects. For Milky Way-like haloes, it increases from about eight per cent at infall to about 25 per cent at the current epoch. In both runs, haloes of fixed mass formed later host more subhaloes than early formers.

We study the effect of baryonic processes on the shapes of dark matter (DM) haloes from Illustris, a suite of hydrodynamical (Illustris) and DM-only (Illustris-Dark) cosmological simulations performed with the moving-mesh code {\sc arepo}. DM halo shapes are determined using an iterative method based on the inertia tensor for a wide range of $z=0$ masses ($M_{200} = 1 \times 10^{11} - 3 \times 10^{14} M_\odot$). Convergence tests shows that the local DM shape profiles are converged only for $r > 9\epsilon$, $\epsilon$ being the Plummer-equivalent softening length, larger than expected. Haloes from non-radiative simulations (i.e. neglecting radiative processes, star formation, and feedback) exhibit no alteration in shapes from their DM-only counterparts: thus moving-mesh hydrodynamics alone is insufficient to cause differences in DM shapes. With the full galaxy-physics implementation, condensation of baryons results in significantly rounder and more oblate haloes, with the median minor-to-major axis ratio $\left< s\equiv c/a \right> \approx 0.7$, almost invariant throughout the halo and across halo masses. This somewhat improves the agreement between simulation predictions and observational estimates of the Milky Way halo shape. Consistently, the velocity anisotropy of DM is also reduced in Illustris, across halo masses and radii. Within the inner halo ($r=0.15 R_{200}$), both $s$ and $q$ (intermediate-to-major axis ratio) exhibit non-monotonicity with galaxy mass, peaking at $m_* \approx 10^{10.5-11} M_\odot$, which we find is due to the strong dependence of inner halo shape with galaxy formation efficiency. Baryons in Illustris affect the correlation of halo shape with halo properties, leading to a positive correlation of sphericity of MW-mass haloes with halo formation time and concentration, the latter being mildly more pronounced than in Illustris-Dark.